Magnetoelastic Sims

[jplauzie]

Hi,

I am working on extending the minimizer to allow for magnetoelastics. The initial results are promising and it works as expected, I just need to do some final cleanup and I hope to be able to share it soon. To verify the results, I am comparing it to some long runs using the dynamics solver. I ran into a few issues with the runs():

1. In Add adaptive time stepping to elastic sims via enum in Variable #60 , is the adaptive timestep known to collapse to very small values (~1e-15 to ~1e-16, it basically stops) for certain problems? The PR mentions looking at sensibletimestep, so you may already be aware of this. I can send a test script, I'm not sure if this is expected or not.

2. How are boundary conditions handled when you have a geometry that doesn't fill the simulation box? I was peeking through, and Magnetoelasticity #4 mentions only free BC's are currently implemented. However, the displacements outside the geometry are being pinned to 0,0,0 , so it seems to be doing something.

3. Uniform displacements and rotations seem to build up over long run()'s, is there a way to subtract them off? Similar to how you would normalize the magnetization. I suppose this isn't a huge deal since it shouldn't affect the dynamics, it just messes with the plots unless some post-processing is done?

edit:
It looks like the minimizer's lastDu is sensitive to those uniform rigid modes

[ilateur]

Woah, before holiday I was also looking into extending the minimizer to include magnetoelastics, but it seems you are further along than I am.

1. No, I was not aware of this. Although depending on the problem, a time step of 1e-15 s to 1e-16 s can still be normal. For example we have used time steps in that range with (magnetoelastic) spin waves in antiferromagnets before. You should be able to tweak the "adaptivity" by playing with displacement_max_error, magnetization_max_error and velocity_max_error. The displacement and/or velocity max errors are perhaps too strict.
   sensibleTimestep is only used at the beginning of a run to get the time step somewhere close to the correct order of magnitude, and the adaptive time stepping should be able to quickly tune the time step further. This can go wrong it it is multiple orders of magnitude off. The current sensible timestep formula is nonsense for elastic equations, but I don't know how much it matters. It can be bypassed by calling

world.timesolver.timestep = choose_starting_timestep
world.timesolver._impl.run(duration)  # notice _impl

2. The boundary conditions are now traction(-free) boundary conditions, $\bf{n} \cdot \sigma = \bf{0}$, with $\bf{n}$ the normal to the relevant surface. They should behave exactly the same for partially-filled simulation boxes. Although, I can imagine that strangely shaped geometries can become unstable much more quickly. Cells outside the geometry are completely ignored and thereby also pinned. The implementation sneaks its way into the calculation of the internal body force (see Improve accuracy of traction-free BC to 2nd order #68).
3. I have noticed those too. There is no built-in way to remove them, as I was a bit unsure as of how to cleanly do that, and perhaps one wants to study that behavior (but probably not). And like you said, it shouldn't affect the dynamics.
   To fix the translation, the center of mass should stay fixed. Multiplying displacement with the mass density (rho), summing and dividing by the "total mass density" should give the position of the center of mass, which you can then subtract from the displacement. Or if mass density is uniform, simply subtracting the average displacement from the displacement works too. There's probably a similar formula to undo the total rotation. You can do this in Python before plotting if you'd like.

By the way, how are you extending the minimizer? I was thinking of setting the velocity to 0, then using the force density directly as the step direction for the displacement, as that is the energy descent direction, similar to torque for magnetization. Magnetic and elastic equations would get their own Barzilai Borwein step sizes.

[jplauzie]

No, I was not aware of this. Although depending on the problem, a time step of 1e-15 s to 1e-16 s can still be normal. For example we have used time steps in that range with (magnetoelastic) spin waves in antiferromagnets before. You should be able to tweak the "adaptivity" by playing with displacement_max_error, magnetization_max_error and velocity_max_error. The displacement and/or velocity max errors are perhaps too strict.

I'm not as familiar with AFMs, but I think it makes sense to have a smaller timestep for those, because their characteristic frequency is higher, typically THz compared to GHz for FM. I've only been testing with FMs for now, the main case I ran into this was a donut shape with a vortex initial state (nothing particularly special, I just wanted a non-uniform test). I tried fixing a timestep of 1e-14 to force it, and that blows up the results around 10ns or so. So it seems to be something real. (Strangely, when I tested a similar configuration using Frederic's mumax3 ME module, it didn't blow up?? Maybe some slight difference in how the BC's are defined? I know that for instance he used u_z had to be uniform/zero). This test script was without absorbing bumpers, but bumpers also didn't seem to help.

minitestv3_donut.txt

And like you said, it shouldn't affect the dynamics.

Do you ever have to worry if they get large? I noticed for large runs, you start to see net displacements around ~1-3nm or so. So I was wondering if it matters specifically if they get close to the cellsize? Particularly, the displacements pre step tend to be pretty small, like ~1e-15, is there a noise floor concern? Even if in general it shouldn't matter

To fix the translation, the center of mass should stay fixed. Multiplying displacement with the mass density (rho), summing and dividing by the "total mass density" should give the position of the center of mass, which you can then subtract from the displacement. Or if mass density is uniform, simply subtracting the average displacement from the displacement works too. There's probably a similar formula to undo the total rotation. You can do this in Python before plotting if you'd like.

Cool, that makes sense. I found your Python example and copied that, but I noticed it was making the minimizer take many extra steps, so I just turned it into a Cuda function. Rotation was a bit more involved, I basically just did a least-squares fit for uniform rotation and subtracted it off, which should be the equivalent. I'm not sure if this is the best way, but it seems to work ok. I think solvers like PetsC seem to do something a bit different 1 2, but I think that's more because it's general purpose and they need to handle arbitrary equations/operators, not because it's better?

By the way, how are you extending the minimizer? I was thinking of setting the velocity to 0, then using the force density directly as the step direction for the displacement, as that is the energy descent direction, similar to torque for magnetization. Magnetic and elastic equations would get their own Barzilai Borwein step sizes.

That is pretty much exactly what I did, and then I just alternated magnetic/elastic steps separately until they both hit convergence. If you want to take a look, it is here. I'll also throw up a PR in a bit. It mostly works as you would expect, I only had to make one or two changes; a scaling trick because the elastic BB1 numerator had a habit of underflowing in float32, and the subtracting off of the rigid modes (I used a lastDu criteria to match the magnetic minimizer's lastDmsamples. Maybe it would be better to check if the force stagnates instead of Du? I'm not sure).

There's basically just 2 files with major changes, minimizer.cu and removerigidbodymodes.cu. Minimizer.cu I think is mostly done (but feel free to adjust/change whatever if something seems silly), I just need to do some minor things like remove the diagnostic messages, and figure out if rigid modes only need to be removed every ~100 steps or so, for performance. It converges to the same states that a long run does. The minimizer portion will work without removerigidbodymodes, it just inflates the number of steps it takes if it's not used.

So far, I've only tested it with FMs, with a uniform 2-D geometry, and a few starting magnetic configurations (uniform, vortex, etc). It also handles a user-applied body-force. The function hooks are there for it to work with AFMs/altermagnets, but I don't know as much about them for magnetoelastics so I haven't tested them yet.

removerigidbodymodes works mechanically (at least for simple test cases of uniform rho and a plain 2-D geometry), but I haven't had a chance to clean it up yet so it's a bit of a mess. There's definitely some nonsense in terms of readability, unnecessary use of doubles on the GPU, etc. In principle it has the code to handle non-uniform rho/other geometries, but I've only tested against uniform 2-D sheets after i ran into the donut issue. I have family visiting for a few days, but once I'm back it shouldn't take that long to get it cleaned up

[jplauzie]

In regards to 1):

I found your (very nice) plotting code, and it seems the stepsize problem is actually because the no-traction boundary conditions are causing the displacement to break for some simulations? This is running magnetoelastic.py from #60, but with 2 minor changes:

No PBC, by changing
world = World(cellsize, mastergrid=Grid((nx, ny, 0)), pbc_repetitions=(2, 2, 0))
to
world = World(cellsize)

and setting the driving field to 0, so there is no excitation (just to rule out that it was pumping energy into the system or something). This is also with the same 5e4 amplification in the original script for the graphic

magnetoelastic.mp4

magnetoelastictest.txt

I did try removing the rigidmodes from the dynamic solver as well, just to make sure translation/rotation weren't causing the weird graphics, but that doesn't really seems to change this output. So it seems to be something real, not just due to the sample translating/rotating.

The vortex/donut has a similar break, except a bit more violent:

magnetoelastic.mp4

magnetoelastictest2.txt