Error in extracting dispersion curves in post processing

[jpla...@rams.colostate.edu]

Hello,

I am trying to reproduce the dispersion curves of a CoFeB waveguide, shown in Figure 3 of "Confined magnetoelastic waves in thin waveguides" by Vanderkeven et al.

They provided their code in another paper, "Finite difference magnetoelastic simulator", so i am using it essentially verbatim, except tweaking the excitation to be a sinc function- the sample code is for a specific mode. To do the postprocessing, I am following the procedure given in the mumax workshop for the nanowire. However, when i plot the curve, all i see is a large spike at k=0 and f=0.

The excitation is defined as:

B_ext.setregion(2, vector(5e-3,(1e-3)*sinc(2*pi*50e9*t), 0))

where region 2 is a 100 nm strip of the waveguide. So it applies a 5mT field along the waveguide, and a 50 GHz sinc in the y direction

Can anyone give advice? I've tried a few tweaks, like making it a sinc pulse at 50 GHz to ensure it's exciting all of the modes. However, i suspect my issue is with the post processing in python, perhaps how i'm defining the x component of the magnetization as 'mx = m[:,0,0,0,:]' . The relevant bit to generate the 2-D FFT, following the workshop is:

    table, fields = run_mumax3(script,"spinwaves")

    m = np.stack([fields[key] for key in sorted(fields.keys())])

    mx = m[:,0,0,0,:]

    mx_fft = np.fft.fft2(mx)

    mx_fft = np.fft.fftshift(mx_fft)

    plt.figure(figsize=(10,6))

    extent = [ -(2*np.pi)/(2*dx), (2*np.pi)/(2*dx), -1/(2*dt), 1/(2*dt)] 

  plt.imshow(np.abs(mx_fft)**2,extent=extent,aspect='auto',origin='lower',cmap="inferno")

    plt.xlim([-2e8,2e8])

    plt.ylim([1e9,fmax])

    plt.ylabel("$f$ (Hz)")

    plt.xlabel("$k$ (1/m)")

    plt.show()

If i understand correctly, this should take the 2-D FFT along the x-direction and the time domain. It seems very similar to the nanowire example of the workshop, except the y-dimension is bigger than 1 cell,  the bias field is in the x direction, and excitation in y. 

I can post the full code below if it is helpful to solve the issue. I can also post the normal mumax3 equivalent, without the magneto-elastic effects, the error is the same. The message board was not allowing attachments/links.

Thanks,

Josh L.

[Felipe Garcia]

Hi Josh L.,

With your configuration, one has to analyze my.

This is the part that is wrong

mx = m[:,0,0,0,:]

mx_fft = np.fft.fft2(mx)

In your case one has to analyze my. With the field along x, it is normal that the biggest peak is for k=0 because this is the continuous part because the magnetization is oriented along x, so one will have second order spin waves which are much smaller than the unity.

I think you have to change to (please check)

my = m[:,1,0,0,:]

    

Well actually, the name of the variable is not important but one has to select my. In the example the sinc pulse is applied along x.

Best regards,

Felipe

--
You received this message because you are subscribed to the Google Groups "mumax2" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mumax2+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mumax2/4dcb1e4a-a089-4a4e-ba12-f41a333d5f7an%40googlegroups.com.

[jpla...@rams.colostate.edu]

Hi Felipe,

Thanks! It makes sense to take my, since that the direction of the excitation. And i think your form is correct. my = m[:,1,0,0,:]. The first index is time, the second gets the y component of magnetization, and the last 3 are z,y,x coordinates respectively (and I want along the x, so : in the x index). (I guess my = m[:,1,0,19,:] would be better, to select along the middle of the waveguide (Since it is 40 cells long in the y direction), but that should be mostly equivalent in this case)

When i do that, i got something that didn't follow the expected dispersion of k^2 . It looks like this

The response is still being dominated by the static magnetization (this is also basically what i get if i FFT with no excitation field running). I suppose not surprising, since the waveguide has regions of high alpha to prevent reflections at the boundaries, so significant portions of the waveguide aren't excited.2 follow up questions, if you don't mind. 

I guess i probably have to subtract this off to see the detail of the second order spin waves. Is their a standard procedure for this? I would guess something like substract off the static magnetization (either pre or post FFT), or a semi-log plot?

How come in the wire example, the static mode is not pronounced? I would've expected a similar issue for that problem as well, but they don't do anything to subtract this off and it is not an issue. The k^2 relation dominates quite strongly.

Best,

Josh Lauzier

[Felipe Garcia]

Hi Josh Lauzier,

Regarding your questions

I guess i probably have to subtract this off to see the detail of the second order spin waves. Is their a standard procedure for this? I would guess something like substract off the static magnetization (either pre or post FFT), or a semi-log plot?

I recommend you try to get the linear spin waves first because they will have larger intensity. To avoid that one can remove the equilibrium magnetization before the FFT. The equilibrium magnetization has to be transversal to the rf field to have linear spin waves (or more precisely the reversed sentence). The best is to be sure that this happens and remove the equilibrium value. It is possible to have an intermediate situation when one applies a tilted field but there in almost all the cases there will be one component which is transversal to the equilibrium magnetization. One can also do after the FFT is done removing the k=0 value but this is not nice if one wants to plot for example the full parabola.. I mentioned the second order spin waves because they will appear analyzing in the same direction like the applied dc field but they will be much smaller than the linear spin waves.

How come in the wire example, the static mode is not pronounced? I would've expected a similar issue for that problem as well, but they don't do anything to subtract this off and it is not an issue. The k^2 relation dominates quite strongly.

In that case, the magnetization in equilibrium is not oriented along X. The field is strong enough to saturate the magnetization in the z direction because there is no dipolar field. The DC in the FFT of mx is basically the average mx. There may be a DC component during the dynamics but it is small compared to the case of magnetization saturated in that direction, which is a huge value. Therefore, there is no need to do that because there is no such a problem.  Then when one does this the linear spin waves are directly visible.

Best regards,

Felipe

To view this discussion on the web visit https://groups.google.com/d/msgid/mumax2/ccd030db-c47d-4b0f-8a44-e1ee4becad7fn%40googlegroups.com.

[jpla...@rams.colostate.edu]

Hi Felipe,

Thanks, that makes sense. For the above picture, that was with the dc field and equilibrium magnetization in the x direction, along the long direction of the wave guide (10microns in x, 200nm in y, 20nm in z). The excitation field is in the y direction, and for that picture i took the FFT of my. Shouldn't that satisfy the transverse condition for linear spin waves to dominate the response? 

(Granted, I have not subtracted off the equilibrium contribution yet, but it sounds like i shouldn't need to, just to see the linear waves in this configuration)

The Bfield is as

B_ext = vector(5e-3, 0, 0)

defregion(2,xrange(-50e-9,50e-9))

B_ext.setregion(2, vector( 5e-3,  (1e-3)*sinc(2*pi*50e9*t), 0))

with the FFT now as

my = m[:,1,0,19,:]

 my_fft = np.fft.fft2(my)

 my_fft = np.fft.fftshift(my_fft)

What you say makes perfect sense, but I must be making a simple mistake somewhere, since I'm not seeing the linear waves, even if i increase the dc or rf fields. I've also uploaded the full files (file extension changed to .txt to allow the filter to post them) in case it helps. run_mumax3 is a .py python script (includes the mx3, running it, and FFT post processing), and equivalentscript is just the .mx3 portion

Best,

Josh Lauzier

[Felipe Garcia]

Hi Josh Lauzier,

I send you a version of your script working and a figure of the results, which agree well with the expected results from your configuration.

I have made a number of changes.

-You have changed some things that should not be changed because they are automatically calculated. Examples of that are autosave time and run time. One has to keep this order:

autosave(m,{dt})

This part I don't think is necessary to change because the original script automatically calculates dt so that the maximum frequency of the FFT coincides with the maximum frequency of your sinc pulse. Same for other parameters.

More importantly, you modified the center of the pulse and moved it to t=0. This is a bad idea because you remove plenty of the AC components of the signal (you can check using mumax3-fft when the pulse is in the center and at the beginning). Ideally one has to have run and the pulse like:

B_ext.setregion(2, vector(Hdc,Hac*sinc(2*pi*f*(t-{T}/2)), 0))

run({T})

This way the pulse is automatically centered in the middle of the time window of the simulation.

-The system has IP anisotropy. The conditions are normal micromagnetic. This makes that mz at the boundary for the equilibrium configuration is different from zero. For that reason the best variable is mz (or alternatively consider removing the equilibrium part before doing the analysis). This is the main reason you observe mostly DC in your script. The others are to improve the results. In the case of the tutorial the mx values is exactly zero, which is the direction of the rf field, so there is no need to analyze other components. 

Then one has to add:

 mz = m[:,2,0,19,:]

-The DC field is not strong enough to orient the magnetization in the direction of the DC field but, moreover, it is close to the DC component of the excitation. Therefore, you observe reorientation due to the dc part of the sinc function. To avoid that I have increased the DC field and then this effect is removed. The tutorial case has a strong field which removes this effect. You can see this because my does not tend to zero in your original simulations.

Hdc := 50e-3

-The equilibrium state is different from your initial condition which makes then that you see the transient to the equilibrium state. To do that is better to relax and minimize before the sinc pulse.   Otherwise you increase the k=0 intensity. This is not necessary in the tutorial because the uniform state is already the equilibrium.

Then one has to include in the script

relax()
minimize()

-I also modified the value of alpha but I think that is not necessary and maybe it is to assign a smaller value to that parameter in order to obtain larger intensities.

Best regards,

Felipe

--
You received this message because you are subscribed to the Google Groups "mumax2" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mumax2+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mumax2/064db422-5f0b-4034-b0f1-9fcb887bd237n%40googlegroups.com.

[Felipe Garcia]

Actually, while editing I have removed one of the comments. The working script is this one.

[Eloi Junior]

Hello everyone, could someone please help me with a similar error? This is my code:

m = np.stack([fields[key] for key in sorted(fields.keys())])

# Select the y component and the (only) layer y=0
my = m[:,1,0,:,:]

# Apply the FFT for every cell
my_fft = np.fft.fft(my, axis=0)

# Select the the two modes

mode3 = my_fft[mode3_idx]


# Plot the result
plt.figure(figsize=(10,4))
plt.subplot(1,3,1)
plt.title("$m_z$")
plt.imshow(my[0])

plt.title("Mode 3")
plt.imshow(np.abs(mode3)**2)

plt.show()

is returning the following error, but I don't know exactly how to solve it:

ValueError Traceback (most recent call last) <ipython-input-22-423c5fd653a5> in <module> 1 # Stack all snapshots (4D arrays) of the magnetization on top of each other 2 # The results in a single 5D array (first index is the snapshot index) ----> 3 m = np.stack([fields[key] for key in sorted(fields.keys())]) 4 5 # Select the y component and the (only) layer y=0 <__array_function__ internals> in stack(*args, **kwargs) ~/.local/lib/python3.8/site-packages/numpy/core/shape_base.py in stack(arrays, axis, out) 424 shapes = {arr.shape for arr in arrays} 425 if len(shapes) != 1: --> 426 raise ValueError('all input arrays must have the same shape') 427 428 result_ndim = arrays[0].ndim + 1 ValueError: all input arrays must have the same shape

I appreciate any help

[Josh Lauzier]

Hello,

You should check the fields variable and see what size each array is (alternatively, just check the OVF files themselves generated by mumax). Somewhere in there, you likely have a bad (mis-sized) OVF file. Most likely, you have a left over OVF sitting in the folder from some other task that is the wrong size.

For example if you were generating 45 OVF files, but you cancelled it and shortened it to 25 to test faster. The new script will overwrite the first 25 OVFs automatically, but the second 20 old files will just sit in the directory. And if you cancelled midway through say the 45th run, that 45 OVF file will only be half-full or empty. 

The read_mumax3_ovffiles python script from the colab just grabs all OVFs in a folder based on the filetype

[Eloi Junior]

Hello, Josh

I ran the script again and saved it in another folder to avoid the possible errors you mentioned. At the end I checked the sizes of the OVF files and they all have the same size. However, the error continues to occur. Is there any other solution, or reason why this error is happening?

[David Cortés-Ortuño]

Hi, as mentioned before you should check your arrays are the same size before stacking them:

for key in sorted(fields.keys()):

    print(fields[key].shape)

or even better, just print which elements have different size:

ks = list(fields.keys())

ref = fields[ks[0]].shape

for k in ks:

    for i in range(4):   # check each size of every field (4 dim array)

        if fields[k].shape[i] != ref[i]:

            print(f"Different size for key {key} at shape index {i}")

[Josh Lauzier]

The reason it's giving you that error is that the arrays need to be the same shape in order to use np.stack . For example, if you're stacking N 2x2x2x1 arrays it creates a 5-D array 2x2x2x1xN. But if one of them is 3x3x1, it doesn't know what to do with it, because that doesn't easily map/stack onto the 2x2x2x1 arrays

The reason i suggested checking the OVF file is we know that the python script (assuming you didn't make any changes to it) should be pretty robust for loading values from the OVF files. So likely the issue is in the format of the OVFs themselves you're feeding into python, for some reason. Which then gets carried over to the fields variable in python, with some of the fields arrays not being the same size.

I would next recommend checking the shape of the fields variable next. Either just printing them out, or checking out their .shape property (or both). The code snippets David listed are a good ways to do so. For that error, there should be at least 1 array (maybe more) that is not the same shape as the others, or has a weird shape (like (n,) instead of (n,1)). You'll then have to track down in your mx3 script why it's not the same size as the others (and since one of the elements is in the time domain, it could be something as simple as you ran it for a shorter simulation time than the others).