On the problem of simulating the dispersion of magnetic moments in the z-direction with magnetoelastic coupling.

[LQ]

On the problem of simulating the dispersion of magnetic moments in the z-direction with magnetoelastic coupling.

I used the mumax extension package of the magneto-elastic coupling (https://github.com/Fredericvdv/Magnetoelasticity_MuMax3)  to simulate the dispersion of the magnetic moment in the z direction with the parameters of the prb (https://open-research-europe.ec.europa.eu/articles/1-35/v1)  article, but the magnetic moment will become very messy when running only. 

I am not sure if it is a bug in this software package or if I set it wrong? Does anyone know why, thanks a lot.

 this is may mumax code:

OutputFormat = OVF2_TEXT
//mesh
dx := 5e-9
dy := 5e-9
dz := 20e-9
Nx := 4000
Ny := 40
Nz := 1
SetMesh(Nx, Ny, Nz, dx, dy, dz, 0, 0, 0)
SetPBC(1, 1, 0);
//Parameters
f := 50e9
Hac := 1e-3
Hdc := 5e-3

//Elastics
C11 = 283e9
C12 = 166e9
C44 = 58e9
rho = 8e3
eta = 0
//B1=-8.8e6
B2=-8.8e6


//Magnetisation
Msat = 1.2e6
Aex = 18e-12
alpha = 4e-3

//##############################
//Left
defregion(24,xrange(-5.2e-6,-5e-6))
eta.setregion(24,1e12)
alpha.setregion(24,1e-2)
defregion(25,xrange(-5.4e-6,-5.2e-6))
eta.setregion(25,3e12)
alpha.setregion(25,5e-2)
defregion(26,xrange(-5.6e-6,-5.4e-6))
eta.setregion(26,6e12)
alpha.setregion(26,1e-1)
defregion(27,xrange(-Inf,-5.6e-6))
eta.setregion(27,5e13)
alpha.setregion(27,0.5)
//right
defregion(14,xrange(5e-6,5.2e-6))
eta.setregion(14,1e12)
alpha.setregion(14,5e-2)
defregion(15,xrange(5.2e-6,5.4e-6))
eta.setregion(15,3e12)
alpha.setregion(15,5e-2)
defregion(16,xrange(5.4e-6,5.6e-6))
eta.setregion(16,6e12)
alpha.setregion(16,1e-1)
defregion(17,xrange(5.6e-6,Inf))
eta.setregion(17,5e13)
alpha.setregion(17,0.5)
//##############################

//Static field
B_ext = vector(0, 0, Hdc)

//Excitation
defregion(2,xrange(-50e-9,50e-9))
B_ext.setregion(2, vector(0, Hac*sinc(2*pi*f*(t-2e-12)), Hdc))

//Initial state
u=uniform(0,0,0)
m=uniform(0,0,1)
relax();
//Solver
SetSolver(9)
dtt := 1e-13
fixdt= dtt

//savings
//autosave(B_mel, 5e-9)
//autosave(F_mel, 0.5e-9)
autosave(m, 1e-11)
autosave(u, 1e-11)
run(10e-9)

[LQ]

Strangely, I found that in Magnetoelasticity_MuMax3/cuda/shearstrainkernel.cu, the z-direction counting all comments are not implemented. Is this a bug in the module?

[Josh Lauzier]

Hello,

Can you explain what you mean by "very messy"? I have run the code before (both the ME extension, and this particular example script), so I can confirm that it works. However, as you probably noticed, the default code is for the 9.8 Ghz mode (If you have read their PRB paper (here), also Ref 71 of the open-EU paper, it corresponds to Fig 7b, not the overall dispersion curve). So a couple small tweaks need to be made to the provided script:

One, you should run() for a period of time before recording, to allow the system to relax. There is some slight relaxation at the edges/corners due to demag, which generate some transient ME waves. I found between ~10-50ns to be good. 10 is probably enough, but 50 is very very clean. What I did was run the system once for 50ns, save m and u, and then load them into the script in place of "u=uniform(0,0,0)

m=uniform(0,0,1)"). I don't recall, but relax() wasn't sufficient (I think because relax only does the magneto portion, not the full coupled equations, but I didn't check), I do remember trying it and still getting some transient behavior.  Probably this is what you mean by "messy", it does interfere a lot with the FFT signal if you do not do this.

Second, you need to tweak the excitation. You already changed it to a sinc, which is good. Ideally, for a dispersion curve you want this centered. To center a pulse that occurs in 10ns, you want it shifted by 5 ns. So "sinc(2*pi*f*(t-2e-12))" should be "sinc(2*pi*f*(t-5e-9))". Otherwise the sync is centered at 2e-12, essentially the start of the simulation. This will cut off a good portion of the sinc. In their PRB paper (here) they mention only pulsing for 20ps. So i sharply cut off the pulse instead of using a sinc. The difference shouldn't matter I think, i have seen both work. What matters is exciting all the modes.

Third, there is a minor typo- Bac and Bdc should be Hac and Hdc, respectively (but it looks like you caught that).

My version of the script is attached (I was studying a slightly different material, so please double check the material parameters, i did make some changes, including to the bumper areas). But I was able to reproduce their figures before tweaking it for my own use.

As for the kernel- I am not sure offhand if that is correct or not. My initial gut reaction is that is strange. That might be worth opening an issue up on the github to confirm. Looking at existing issues, it looks like a bug that was introduced when Frederic corrected a previous error in the kernel: https://github.com/Fredericvdv/Magnetoelasticity_MuMax3/issues/1 . It may not pop up as an issue in simulations where the shear strains happens to be zero

Best,

Josh L.

[LQ]

thanks, I read the author's documentation carefully, ‘u’ only considers displacement in the x and y direction (No z).   

But there is still a very obvious problem when the cell size setting is relatively small. For example 1nm. Regardless of the magnetic order, it will become disordered.  You can have a try. It has nothing to do with long-term relaxation.

Generally, the smaller the mesh, the more accurate it is, but this magnetoelastic module has obvious problems. This is very strange maybe a bug of the module.

[Josh Lauzier]

Hello,

Unfortunately I'm not currently at a workstation with the ME extension to test, but I can take a look tomorrow. Can you upload 2 example scripts, one with it occurring and not occurring? Preferably also with outputs (either OVF or pictures).

Also, I would try a couple sanity tests. For example, a too large dt can cause this sort of 'noise' issue, even with normal micromagnetics (see for instance this discussion). You could also compare by running it with the ME effects disabled, and/or to the base mumax, to verify that it is in fact the ME module that is the issue. 

For example, just running your above script briefly on normal mumax (with only the ME lines removed) with normal mesh sizes posted above, i get something like the attached picture. Which is fairly disordered- but that probably has more to do with the shape anisotropy. I hadn't previously noticed, but it looks like you changed the configuration to the z direction, which is very different than the one shown in the paper. If you disable demag ( like in this example in the mumax summer tutorial example 5, which makes a nice direct comparison if you tweak a few lines, when looking at just the magnetic response.), you get a nice clean response, since the external field is enough to orient the magnetization

I don't have it linked handy, but I do recall looking at smaller mesh sizes, and not noticing any issues. However, I was studying with the dc field along the x-direction and the excitation in the y, as in Vanderveken et al's paper, where the demag field acts to stabilize the magnetization.

Also, your discretization in the z direction is maybe too large if you're interested in z-direction dynamics. 20nm is much larger than the exchange length. You can often get away with a large z discretization it if everything interesting is occurring in plane, but that is not true in your case, so you can't just change x and y discretizations. It may not matter since it's one cell thick, but I would be wary, and discretize as the colab example does to be safe. (I did get a much cleaner response when switched 1 20nm layer in the z to 4 5nm layers, so I suspect this may be at least one of your culprits).

Cheers,

Josh L.

[LQ]

Hello,

I now confirm that the cell size is fine, but the issue of magnetic order becoming disordered still persists, not sure where the problem is. The attachments were run with magnetoelastic (B1, B2 = 0) and the original mumax program respectively, and the results were very different. You can try it. 

thank you.

cheers,

Lei

[Josh Lauzier]

Hi,

I see immediately what you mean, when running it :)

It is the timestep. It still occurs for fixdt= 0.5e-13, but for 1e-14 or 1e-15, the problem resolves itself. So it looks like for your particular set up, you will need to use a smaller step. The necessary time step is related to a lot of factors, but one of them is mesh-size. Finer mesh-sizes in spatial coordinates will require finer timesteps. This will cause it to take a bit longer to run, but otherwise it is fine. Ideally, you want as large a time step as possible (so it finishes quicker) but still small enough to be accurate.

And this also occurs in normal mumax/micromagnetics as well, if you set the solver identically (either setting ME coefficients to 0, using the purely magnetic solver in the ME, or just a non-ME mumax installation), you get the same problem if you use fixdt=1e-13, and the same solver (i think the corresponding solver is RK4 not RK45, it's 4th order Runge-Kutta. So setSolver(4). Mumax by default use RK45 or setsolver(5)). You still see the problem with RK45 as well in this case, though. If you want to recreate the behavior, simply leave in fixdt and setsolver(4) (or setsolver(5) ) in your second script.

If you don't set fixdt and the solver explicitly, mumax will default to RK45 and an adaptive time step (in your latter script, you have both commented out), so it will pick appropriate steps for you. So  in your latter script, it is switching to smaller steps. The adaptive solver seems to go for ~3e-14, so you can probably get away with around that size step. I'd maybe go with 1e-14 to be safe, since RK4 is not quite as the same as RK45. 

The different solvers have different tolerances for the time step (there are trade offs, higher order methods usually allow larger time steps but might require more evaluations of the torque per time step, so the simulation still ends up taking longer. If you're curious, the original paper for mumax3 goes into some detail into the various solvers and their pros/cons in terms of order/convergence etc.). Further, since the default solver has built in adaptive stepping, so it essentially will pick an appropriate timestep for you if you don't manually override it. Hopefully in a future release someone can port over an adaptive timestepping scheme, it is very convenient instead of having to fiddle with it manually. But without that, you just need to judiciously pick an appropriate time step. 

Cheers,

Josh L.

[LQ]

hi，

thank you very much. The problem of magnetic moment disorder has been solved, and it is indeed caused by insufficient time step.   However, no matter how small the time step is, the displacement (u) can never be relaxed to a stable state. The Fourier transform of u is shown in the figure below, and the equally spaced frequencies are moving, but this is not in line with physics. It should be stationary without any disturbance. I don't know what you think, I attach the mumax program reference.

Cheers,

LQ

[Josh Lauzier]

Hello,

It should settle into a steady state, but a few things to consider- 

Your values for the damping (both alpha and eta, but mostly eta) seem relatively low. I don't know normal values for eta offhand (and it will of course vary a lot by material), but in the example paper, they use eta values ~1e12 and higher for the bumpers (the higher damping regions to dissipate the waves without reflecting). Your eta value is several orders of magnitude smaller than that. Even with the 1e12 bumpers, it still took ~10-50ns for the quasi-elastic transient excitations generated by the edges to fully damp out and get to a steady state. So I would expect yours to take much much longer, although it should technically get there eventually. For your value of eta, probably it will be too long to be reasonable to wait for it.

Second, even with realistic etas, the attenuation time for the phonon system can be much much longer than magnon systems. So you might have to go to very large times for it to fully damp out for realistic eta. If you're only interested in achieving the steady state, you can crank up the damping parameters to damp out transitory behavior quicker. Unfortunately, I don't think relax() or minimize() were updated, but having a high damping achieves a very similar effect to relax() (in normal mumax, relax disables the precession completely. But high damping can achieve a similar effect. There is an option to completely disable precession in run() , but I'm not sure how it interacts with the magnetoelastics).

Third, related to the bumpers, you will get reflections at the boundaries without them. This is why the example program has increasing eta/alpha as you get away from the center of the wire. I would recommend introducing bumpers around the edges to absorb these transients, if you don't just use a high eta overall. Since u=(0,0,0) is not a relaxed state (due to magnetostrictive effect), you'll get some excitations, particularly on the edges, as the system relaxes. I am not quite sure what you mean by the frequencies are moving', but I suspect it might be because these reflections. They will appear in the FFT, and they can give some strange looking behavior, especially as they start to interfere with each other and with the non-uniform steady state (with low damping, they will continually bounce back and forth between edges). They are essentially just artifacts. If you watch the animation as it runs, you can see them bounce/interfere. 

I didn't have time to test extensively, but if I crank eta up it settles quite fast. I attached a picture of the steady state. For uniform eta ~1e14,  u settled into steady state less than 1ns. With bumpers (but leaving your eta of 1e6), it took ~3-5ns or so to damp out the transient signal. Probably uniform eta ~1e12 should settle in some nanoseconds (maybe more than 20 though). You can actually see this same steady state form in your original program, it just gets overlaid with the transient reflections behavior that doesn't damp out fast enough.

Best regards,

Josh L.