Spontaneous generation of elastic velocity waves leading to zeroing parameters during the course of simulation

[kbrzusz]

Hello, I spotted following issue during preparation of simulations on 1.1.0. Later I replicated this issue on 1.1.2.

Here's MWE (this code follows the example on magnetoelastics):

cx, cy, cz = 5e-9, 5e-9, 20e-9
sx, sy, sz = 240e-9, 240e-9, 20e-9
nx, ny, nz = int(round(sx / cx)) + 1, int(round(sy / cy)) + 1, int(round(sz / cz))

time: float = 2.5e-9
steps: int = 250
dt = time / steps

circle_shape = mumaxplus.util.Cylinder(sx, sz)
circle_shape = circle_shape.translate(sx / 2, sy / 2, 0)

cellsize = (cx, cy, cz)
grid = Grid((nx, ny, nz))
world = World(cellsize)
magnet = Ferromagnet(world, grid, geometry=circle_shape)

magnet.msat = 4.9e5
magnet.aex = 8e-12
magnet.alpha = 0
magnet.ku1 = 0
magnet.anisU = (0, 0, 1)
magnet.magnetization = (1, 0, 0)

magnet.enable_elastodynamics = True
magnet.rho = 8.9e3
magnet.B1 = -3.4e6
magnet.B2 = -1e6
magnet.C11 = 246e9
magnet.C44 = 124e9
magnet.C12 = 147e9
magnet.eta = 1e10
magnet.elastic_displacement = (0, 0, 0)

world.minimize()

parameter_array = np.zeros((steps + 1, 3, 3, *magnet.geometry.shape))  # step, parameter, direction, Z, Y, X

from tqdm import tqdm
print("running...")
for i in tqdm(range(1, steps + 1)):
    parameter_array[i, 0] = magnet.full_magnetization.eval()
    parameter_array[i, 1] = magnet.elastic_displacement.eval()
    parameter_array[i, 2] = magnet.elastic_velocity.eval()

    world.timesolver.run(dt)

It seems there occurs a generation and amplification of elastics waves -- it's especially destructive in case of curved shapes -- in the case of cx (cy)=2nm, simulation silently becomes malformed (eval()s saved into parameter_array consist only of 0s). Below are snapshots of magnetization and elastic velocity for t=0.26ns (colors depict direction and magnitude):

On magnetization snapshot ripples can be seen, meanwhile elastic velocity has giant values in diagonal regions (up to ~1.8e3) while center remain at <10.

Same effect happens for bigger cx (cy), for 5nm it occurs at ~0.49ns, 10nm -> ~1ns, 20nm -> ~0.89ns.

I suspect that the presence of gaps of non-magnetic medium plays crucial role in this effect -- Cuboid remains stable in much longer timescale (for 5nm it's up to ~3.8ns), but even in such case there are elastic waves generated. Moreover speed of the simulation is elevated -- it takes larger number of steps to drop below 1it/s.

Am I missing something? How can I make magnetoelastics work correctly?

[kbrzusz]

Hello,

apologies, part of reported issue is due to my own code. The zero part of eval() is due to:

1. initialization of parameter_array as np.zeros and
2. termination (which I didn't include to keep MWE truly Minimal) of the loop due to low dt per

from collections import deque
timestep_deque = deque(maxlen=8)

outside the loop and

if np.mean(timestep_deque) < 2.5e-14:         
    break

inside the loop.

But the issue of generation of elastic waves leading to slowdown of calculations and abrupt changes of direction of vectors persists.
Given the fact the system is minimized and no external factors are applied to it, I suppose it remain at equilibrium.

I forgot to mention in my initial post, but I checked both single and double precision -- I didn't notice any difference.

Apologies once more!

[ilateur]

Hi, sorry for the delayed response. I'd like to debug this issue myself, but I'm having some technical difficulties at the moment (broken GPU), so I'll try to give some suggestions instead.

• I'd recommend using fixed time steps for (magneto)elastic simulations, as in the example. Adaptive time stepping is only really supported for magnetic equations of motion (unless using the code from ancient pull request Add adaptive time stepping to elastic sims via enum in Variable #60).

world.timesolver.adaptive_timestep = False
world.timesolver.timestep = 1e-12  # play around with this value, smaller should be more stable (hopefully)

• You could try to increase the damping parameters.
  • eta could be higher, but probably won't help much here. 1e10 is very low, 1e14 is very high.
  • you could increase the hidden stiffness_damping, this could have the largest effect as it's good at damping high frequency oscillations.
  • non-zero magnetic damping alpha might also help? Purely magnetic equations also tend to be unstable without damping.
• I also do not know to what extent minimize minimizes a coupled system, as it was written for magnetization specifically (same for relax). Otherwise, (manually) try to get as close to the relaxed state as possible. This might not be the expected relaxed magnetization with 0 displacement.

If these suggestions do not help, I fear that the boundary conditions might be inherently unstable for (combinations of) corners, of which there are a lot here, even though they are technically written to be general. In that case, my apologies...

P.S. The middle of the magnet is confusingly located at c*(n - 1)/2 (not s/2), because the origin is at the middle of the bottom left cell, not the corner. If you change this in your circle_shape translation, it should get rid of the asymmetric single cells at the right and bottom of the circle.

circle_shape.translate(cx * (nx - 1) / 2, cy * (ny - 1) / 2, 0)