diff --git a/.github/workflows/build-wheels.yml b/.github/workflows/build-wheels.yml index ff5585a1..e24399e6 100644 --- a/.github/workflows/build-wheels.yml +++ b/.github/workflows/build-wheels.yml @@ -19,7 +19,7 @@ jobs: matrix: os: - {name: "ubuntu", version: "22.04"} - - {name: "windows", version: "latest"} + - {name: "windows", version: "2022"} cuda: ["12.4.1"] python: ["3.11", "3.12", "3.13", "3.14"] @@ -129,7 +129,7 @@ jobs: steps: - name: Download all artifacts - uses: actions/download-artifact@v4 + uses: actions/download-artifact@v7 with: path: artifacts @@ -158,6 +158,6 @@ jobs: with: tag_name: ${{ github.event.release.tag_name || github.event.inputs.tag }} files: dist-gh-release/* - overwrite: true + overwrite_files: true env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/deploy-docs.yml b/.github/workflows/deploy-docs.yml index c9527695..d996ea39 100644 --- a/.github/workflows/deploy-docs.yml +++ b/.github/workflows/deploy-docs.yml @@ -14,7 +14,7 @@ jobs: steps: - name: Checkout repository - uses: actions/checkout@v4 + uses: actions/checkout@v6 - name: Deploy docs/_build/html to GitHub Pages uses: peaceiris/actions-gh-pages@v4 diff --git a/docs/_static/custom.css b/docs/_static/custom.css index 45da8fc6..e85a365b 100644 --- a/docs/_static/custom.css +++ b/docs/_static/custom.css @@ -21,3 +21,7 @@ html[data-theme="dark"] .navbar-brand img { display: block; width: 100%; } + +.bd-sidebar-primary .navbar-icon-links { + justify-content: space-evenly !important; +} diff --git a/docs/api.rst b/docs/api.rst index c54b9285..fe4b3b21 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -28,6 +28,7 @@ API PoissonSystem DmiTensor BoundaryTraction + Window .. toctree:: :maxdepth: 1 diff --git a/docs/examples.rst b/docs/examples.rst index 2ef5a810..d1954739 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -13,3 +13,4 @@ Here we show some example Python scripts to run mumax⁺ simulations. examples/DW_SAW examples/voronoi examples/Bloch_wall_altermagnet + examples/moving_window diff --git a/docs/examples/moving_window.rst b/docs/examples/moving_window.rst new file mode 100644 index 00000000..8b6b12dc --- /dev/null +++ b/docs/examples/moving_window.rst @@ -0,0 +1,22 @@ +:nosearch: + +Moving simulation window +======================== + +In this example we move a domain wall in a ferromagnet using a Zhang-Li STT. We let the simulation +window move together with the wall, keeping the domain wall centered in the simulation space. +Using this, we can virtually simulate an infinitely long magnetic nanowire using a limited number +of simulation cells. + +Note: +The moving window functionality only works properly if +- The magnet parameters are uniform +- The magnet has no geometry +- The magnet has no regions + +.. literalinclude:: ../../examples/moving_window.py + :language: python + :lines: 13- + +.. video:: ../images/moving_window.mp4 + :align: center \ No newline at end of file diff --git a/docs/images/moving_window.mp4 b/docs/images/moving_window.mp4 new file mode 100644 index 00000000..fb09abed Binary files /dev/null and b/docs/images/moving_window.mp4 differ diff --git a/examples/Bloch_wall_altermagnet.py b/examples/Bloch_wall_altermagnet.py index e3aa0c90..3d626579 100644 --- a/examples/Bloch_wall_altermagnet.py +++ b/examples/Bloch_wall_altermagnet.py @@ -21,13 +21,10 @@ A12 = A0/2 length = 256e-9 -width = 64e-9 # ----------- Create altermagnet ----------- Nx = int(length / cs) -Ny = int(width / cs) - world = World((cs, cs, cs)) grid = Grid((Nx, 1, 1)) magnet = Altermagnet(world, grid) diff --git a/examples/moving_window.py b/examples/moving_window.py new file mode 100644 index 00000000..1eb535a9 --- /dev/null +++ b/examples/moving_window.py @@ -0,0 +1,146 @@ +""" +In this example we move a domain wall in a ferromagnet using a Zhang-Li STT. We let the simulation +window move together with the wall, keeping the domain wall centered in the simulation space. +Using this, we can virtually simulate an infinitely long magnetic nanowire using a limited number +of simulation cells. + +Note: +The moving window functionality only works properly if +- The magnet parameters are uniform +- The magnet has no geometry +- The magnet has no regions +""" +from mumaxplus import World, Grid, Ferromagnet +from mumaxplus.util import twodomain, plot_field + +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation, FFMpegWriter +import numpy as np + +# ----------- Material and simulation parameters ----------- +cs = 1e-9 +length = 256e-9 +width = 64e-9 +thickness = 1e-9 + +Ms = 600e3 +aex = 10e-12 +alpha = 0.02 +ku = 6e5 +anisU = (1, 0, 0) + +# ----------- Create magnet ----------- +world = World((cs, cs, cs)) +grid = Grid((int(length / cs), int(width / cs), int(thickness / cs))) + +magnet = Ferromagnet(world, grid) + +magnet.msat = Ms +magnet.aex = aex +magnet.alpha = alpha +magnet.ku1 = ku +magnet.anisU = anisU + +magnet.enable_demag = False +magnet.enable_openbc = True + + +# ----------- Create two domain state ----------- +magnet.magnetization = twodomain((1, 0, 0), (0, 1, 0), (-1, 0, 0), magnet.center[0], 5e-9) +magnet.minimize() + +# ----------- Add current and simulate ----------- +magnet.jcur = (-1e12, 0, 0) +magnet.xi = 0.2 +magnet.pol = 1 + +# Center the simulation window, keeping component 0 (x) close to zero. +# We expect motion along the x axis. +world.center_domain_wall(comp=0, axis=0) + +tmax = 0.5e-9 +timepoints = np.linspace(0, tmax, 100) + +def DW_pos(field): + mx = field[0, 0, int(width / cs/2), :] + sign = np.sign(mx) + # find cell where magnetization component crosses zero + i = np.where(np.diff(sign) != 0)[0][0] + # linearly interpolate + m1, m2 = mx[i], mx[i + 1] + return (i - m1 / (m2 - m1)) * cs + +outputquantities = {"mag": magnet.magnetization, + "window position": lambda: world.window.position[0], + "wall position": lambda: world.window.position[0] + DW_pos(magnet.magnetization())} + +output = world.timesolver.solve(timepoints, outputquantities, tqdm=True) + +p0 = output["wall position"][0] + +# ----------- Create movie ----------- +print("Creating animation...") +fig, axes = plt.subplots(2, 1, figsize=(10, 7)) + +# Time trace subplot +ax_trace = axes[0] +lines = {} +keys = ["window position", "wall position"] +for key in keys: + lines[key], = ax_trace.plot([], [], '-', label=key) + +ax_trace.set_xlim(0, tmax * 1e9) +ax_trace.set_ylim(min(0, min(output["window position"]) * 1e9), + max(0, max(output["wall position"] - p0) * 1e9)) +ax_trace.set_xlabel("Time $t$ (ns)") +ax_trace.set_ylabel("position (nm)") +ax_trace.legend() +ax_trace.grid() + +# Magnetization image subplot +ax_image = axes[1] +plot_field(output["mag"][0], ax=ax_image, arrow_size=8) + +ax_image.set_xlim(0, int(length/cs)) +ticks = ax_image.get_xticks() + +ax_image.set_title("$t$ = 0.000 ns") +ax_image.set_xlabel("$x$ (nm)") +ax_image.set_ylabel("$y$ (nm)") + +fig.tight_layout() + +# --- Animation Function --- +def update(frame): + # Update image + ax_image.clear() + plot_field(output["mag"][frame], ax=ax_image, arrow_size=8) + ax_image.set_xlabel("$x$ (nm)") + ax_image.set_ylabel("$y$ (nm)") + + # update x-ticks + shift = output["window position"][frame] * 1e9 + tick_positions = (ticks - shift) + mask = (0 <= tick_positions) & (tick_positions <= int(length/cs)) + + ax_image.set_xticks(tick_positions[mask]) + ax_image.set_xticklabels([f"{t:.0f}" for t in ticks[mask]]) + ax_image.set_title(f"$t$ = {output['time'][frame] * 1e9:.3f} ns") + + # Update time trace + lines["window position"].set_data(np.array(output["time"][:frame+1]) * 1e9, + np.array(output["window position"][:frame+1]) * 1e9) + lines["wall position"].set_data(np.array(output["time"][:frame+1]) * 1e9, + np.array(output["wall position"][:frame+1] - p0) * 1e9) + + return [ax_image] + list(lines.values()) + +# Animation parameters +fps = 15 +anim = FuncAnimation(fig, update, frames=len(output["time"]), + interval=1000 / fps, repeat_delay=5000 / fps) + +# --- Save the animation --- +save_filename = "moving_window.mp4" +writer = FFMpegWriter(fps=fps) +anim.save(save_filename, writer=writer) \ No newline at end of file diff --git a/examples/standardproblem2.py b/examples/standardproblem2.py index f7fd75b1..dd8b46ab 100644 --- a/examples/standardproblem2.py +++ b/examples/standardproblem2.py @@ -30,7 +30,7 @@ def get_next_power_of_2(x): y *= 2 return y -def get_gridsize(L, d, t, l_ex=l_ex): +def get_gridsize(L, d, t): """Cell length should at least be < l_ex/2. The number of cells is best a power of 2 for FFT. This results in cell sizes between 0.25*l_ex and 0.5*l_ex.""" @@ -42,7 +42,7 @@ def get_gridsize(L, d, t, l_ex=l_ex): L = L_p_d * d # dimensionless length L = length/l_ex t = t_p_d * d # dimensionless thickness t = thickness/l_ex - nx, ny, nz = get_gridsize(L, d, t, l_ex=l_ex) + nx, ny, nz = get_gridsize(L, d, t) world = World(cellsize=(L*l_ex/nx, d*l_ex/ny, t*l_ex/nz)) magnet = Ferromagnet(world, Grid((nx, ny, nz))) magnet.msat = msat diff --git a/mumaxplus/__init__.py b/mumaxplus/__init__.py index 8ebba134..ef21d0e5 100644 --- a/mumaxplus/__init__.py +++ b/mumaxplus/__init__.py @@ -55,6 +55,7 @@ from .traction import BoundaryTraction from .variable import Variable from .world import World +from .window import Window from . import util FP_PRECISION = {1: "SINGLE", 2: "DOUBLE"}.get(_cpp.FP_PRECISION, "UNKNOWN") # Use _cpp value, as that is certainly the correct one @@ -77,6 +78,7 @@ "TimeSolver", "Variable", "World", + "Window", "PoissonSystem", "util", "FP_PRECISION" diff --git a/mumaxplus/ferromagnet.py b/mumaxplus/ferromagnet.py index 45a79405..dcda2c3b 100644 --- a/mumaxplus/ferromagnet.py +++ b/mumaxplus/ferromagnet.py @@ -653,7 +653,7 @@ def B1(self) -> Parameter: See Also -------- - B2 + B2, B_chiral """ return Parameter(self._impl.B1) @@ -678,7 +678,7 @@ def B2(self) -> Parameter: See Also -------- - B1 + B1, B_chiral """ return Parameter(self._impl.B2) @@ -697,6 +697,36 @@ def B2(self, value): + " is set to a positive value, instead of negative (or zero)." + " Make sure this is intentional!", UserWarning) + @property + def B_chiral(self) -> Parameter: + r"""Chiral magnetoelastic coupling constant (J/m³). + + Notes + ----- + Materials of the cubic point group 23 (or B20 compounds) can have a + chiral magnetoelastic coupling, with the following energy density. + + .. math:: \mathcal{E} = B_\text{chiral} \sum_{i, j, k} \epsilon_{ijk} \varepsilon_{ii} m_j^2 + + Here :math:`\epsilon_{ijk}` is the Levi-Civita symbol and + :math:`\varepsilon_{ii}` denotes the normal strain components. + This energy density comes from equations (8.12) and (8.16) in Ref. [1], where + B_chiral corresponds to :math:`\lambda_{12}`. Magnetoelastic coupling constants + B1 and B2 correspond to :math:`\lambda_{11}` and :math:`2 \lambda_{44}` respectively. + These lambdas are not the usual magnetostrictive coefficients. + + .. [1] L\ . Franke, “Elastic Coupling at Quantum Phase Transitions and in Chiral Magnets,” Das Karlsruher Institut für Technologie, Karlsruhe, 2025. doi: 10.5445/IR/1000184834. + + See Also + -------- + B1, B2 + """ + return Parameter(self._impl.B_chiral) + + @B_chiral.setter + def B_chiral(self, value): + self.B_chiral.set(value) + # ----- POISSON SYSTEM ---------------------- @property @@ -1128,7 +1158,7 @@ def magnetoelastic_field(self) -> FieldQuantity: See Also -------- - B1, B2 + B1, B2, B_chiral Magnet.strain_tensor, Magnet.rigid_norm_strain, Magnet.rigid_shear_strain magnetoelastic_force """ @@ -1160,7 +1190,7 @@ def magnetoelastic_force(self) -> FieldQuantity: See Also -------- - B1, B2 + B1, B2, B_chiral Magnet.effective_body_force, magnetoelastic_field """ return FieldQuantity(_cpp.magnetoelastic_force(self._impl)) \ No newline at end of file diff --git a/mumaxplus/util/__init__.py b/mumaxplus/util/__init__.py index f18e1552..eedd1aea 100644 --- a/mumaxplus/util/__init__.py +++ b/mumaxplus/util/__init__.py @@ -27,7 +27,7 @@ "show_magnet_geometry", "show_field_3D", "show_regions", - # voronoi + # misc "VoronoiTessellator", "MFM" ] diff --git a/mumaxplus/util/voronoi.py b/mumaxplus/util/voronoi.py index ac8c0655..b5fb0219 100644 --- a/mumaxplus/util/voronoi.py +++ b/mumaxplus/util/voronoi.py @@ -6,6 +6,7 @@ from .. import _cpp from mumaxplus.world import World from mumaxplus.grid import Grid + class VoronoiTessellator: diff --git a/mumaxplus/window.py b/mumaxplus/window.py new file mode 100644 index 00000000..dda1e4c7 --- /dev/null +++ b/mumaxplus/window.py @@ -0,0 +1,60 @@ +class Window: + """Simulation window of the world. + + Each world already has its own Window. This Window can be accessed through + the world.window property. + + Windows should not be initialized by the end user. + """ + + def __init__(self, impl): + self._impl = impl + + def _check_boundary(self, boundary): + if not isinstance(boundary, int) or boundary not in (0, 1): + raise ValueError(f"Invalid boundary: {boundary}. Must be one of (0, 1).") + + def insert_magnetization(self, boundary, value=(0, 0, 0)): + """ + Set magnetization value at a given boundary. + If set to zero (default), the current edge value is used. + + For multi-sublattice magnets, the magnetization value of `sub1` should + be provided. + + Parameters + ---------- + boundary : int + Boundary index: 0 = Left/Bottom, 1 = Right/Top + value : tuple of 3 floats + Magnetization vector (x, y, z) + """ + self._check_boundary(boundary) + self._impl.insert_magnetization(boundary, value) + + def disable_motion(self): + """Disable the motion of the simulation window. + + See Also + -------- + World.center_domain_wall + """ + self._impl.disable_motion() + + @property + def position(self): + """Returns the current position of the simulation window (m). + The origin of the window coincides (when unmoved) with the origin of a `Grid` instance, + i.e. it is determined by the coordinate of the lower left cell. + """ + return self._impl.position + + @property + def velocity(self): + """Returns the current velocity of the simulation window (m/s).""" + return self._impl.velocity + + @property + def total_shift(self): + """Returns the total amount shifted by the simulation window (m).""" + return self._impl.total_shift \ No newline at end of file diff --git a/mumaxplus/world.py b/mumaxplus/world.py index 0373dd85..d6b59bfb 100644 --- a/mumaxplus/world.py +++ b/mumaxplus/world.py @@ -4,11 +4,15 @@ from .timesolver import TimeSolver from .grid import Grid +from .magnet import Magnet from .ferromagnet import Ferromagnet from .antiferromagnet import Antiferromagnet +from .altermagnet import Altermagnet from .ncafm import NcAfm +from .window import Window import warnings +import numpy as np class World: """Construct a world with a given cell size.""" @@ -96,21 +100,31 @@ def get_ncafm(self, name): raise KeyError(f"No magnet named {name}") return NcAfm._from_impl(magnet_impl) + @property + def magnets(self) -> dict[str,Magnet]: + """Get a dictionary of all magnet names.""" + return {**self.ferromagnets, **self.antiferromagnets, **self.altermagnets, **self.ncafms} + @property def ferromagnets(self) -> dict[str,Ferromagnet]: - """Get a dictionairy of :class:`Ferromagnet` names.""" + """Get a dictionary of :class:`Ferromagnet` names.""" return {key: Ferromagnet._from_impl(impl) for key, impl in self._impl.ferromagnets.items()} @property def antiferromagnets(self) -> dict[str,Antiferromagnet]: - """Get a dictionairy of :class:`Antiferromagnet` names.""" + """Get a dictionary of :class:`Antiferromagnet` names.""" return {key: Antiferromagnet._from_impl(impl) for key, impl in self._impl.antiferromagnets.items()} + @property + def altermagnets(self) -> dict[str,Altermagnet]: + """Get a dictionary of :class:`Altermagnet` names.""" + return {key: Altermagnet._from_impl(impl) for key, impl in + self._impl.altermagnets.items()} @property def ncafms(self): - """Get a dictionairy of non-collinear antiferromagnets by name.""" + """Get a dictionary of non-collinear antiferromagnets by name.""" return {key: NcAfm._from_impl(impl) for key, impl in self._impl.ncafms.items()} @@ -334,3 +348,62 @@ def unset_pbc(self): set_pbc """ self._impl.unset_pbc() + + @property + def window(self) -> Window: + """Simulation window for this world.""" + return Window(self._impl.window) + + def center_domain_wall(self, comp, axis=None): + """Move the simulation window along with the domain wall. This function + should be called before the `TimeSolver`. + + Note + ---- + The domain wall should be centered already for this function to work properly. + + Parameters + ---------- + comp : int + The magnetization direction of the domains. The average of this + magnetization component will be kept close to zero. Possible values + are 0 (x-component), 1 (y-component) and 2 (z-component). + axis : int + The axis along which the domain wall moves (i.e. the wall normal). + Possible values are 0 (x-direction), 1 (y-direction) and 2 (z-direction). + If axis is `None` (default), then the first axis with the most number of + grid cells is chosen. + + warning + ------- + `center_domain_wall` will not work properly with non-uniform parameters. + + See Also + -------- + Window.disable_motion + """ + if comp not in (0, 1, 2): + raise ValueError("The component `comp` should be 0 (x), 1 (y) or 2 (z).") + + if len(self.magnets) != 1: + raise RuntimeError("The moving window functionality only works if exactly 1 magnet exists.") + + magnet = list(self.magnets.values())[0] + if not np.all(magnet.geometry) or np.any(magnet.regions): + raise RuntimeError("The moving window functionality doesn't work well with geometry" + " or regions as of yet.") + + # If no axis is given, return first axis with most number of cells + if axis is None: + axis = np.argmax(self.bounding_grid.size) # bounding grid is safe if only 1 magnet + warnings.warn("There is no axis provided in the moving simulation window." + + f" The {('x', 'y', 'z')[axis]}-direction is used as normal to the" + + " domain wall", UserWarning) + + av = magnet.magnetization.average() if isinstance(magnet, Ferromagnet) else magnet.sub1.magnetization.average() + if np.abs(av[comp]) > 4 / magnet.grid.size[0]: + raise RuntimeError(f"The domain wall does not seem centered (average {('x', 'y', 'z')[comp]}-" + + f"component is {av[comp]:.2e}). `center_domain_wall` only works properly " + + "if the wall is initialized near the center of the magnet.") + + self._impl.center_domain_wall(comp, axis) \ No newline at end of file diff --git a/src/bindings/CMakeLists.txt b/src/bindings/CMakeLists.txt index f0ca66ab..955c4fdb 100644 --- a/src/bindings/CMakeLists.txt +++ b/src/bindings/CMakeLists.txt @@ -30,6 +30,7 @@ pybind11_add_module(${MUMAX_MODULE_NAME} wrap_poissonsystem.cpp wrap_linsolver.cpp wrap_voronoi.cpp + wrap_window.cpp wrap_traction.cpp ) diff --git a/src/bindings/main.cpp b/src/bindings/main.cpp index 142f9c5d..7532cd4a 100644 --- a/src/bindings/main.cpp +++ b/src/bindings/main.cpp @@ -32,5 +32,6 @@ PYBIND11_MODULE(MUMAX_MODULE_NAME, m) { wrap_system(m); wrap_dmitensor(m); wrap_voronoi(m); + wrap_window(m); wrap_traction(m); } diff --git a/src/bindings/wrap_ferromagnet.cpp b/src/bindings/wrap_ferromagnet.cpp index 91ec4a94..927a9914 100644 --- a/src/bindings/wrap_ferromagnet.cpp +++ b/src/bindings/wrap_ferromagnet.cpp @@ -72,6 +72,7 @@ void wrap_ferromagnet(py::module& m) { .def_readonly("poisson_system", &Ferromagnet::poissonSystem) .def_readonly("B1", &Ferromagnet::B1) .def_readonly("B2", &Ferromagnet::B2) + .def_readonly("B_chiral", &Ferromagnet::BChiral) .def("reset_noise_generator", &Ferromagnet::resetNoiseGenerator) .def("minimize", &Ferromagnet::minimize, py::arg("tol"), py::arg("nsamples")) diff --git a/src/bindings/wrap_window.cpp b/src/bindings/wrap_window.cpp new file mode 100644 index 00000000..433c30ea --- /dev/null +++ b/src/bindings/wrap_window.cpp @@ -0,0 +1,15 @@ +#include "window.hpp" +#include "wrappers.hpp" + +void wrap_window(py::module& m) { + + py::class_(m, "Window") + .def("insert_magnetization", [](Window& self, int side, real3 value) { + self.setMagValue(static_cast(side), value); + }, py::arg("side"), py::arg("value")) + + .def_property_readonly("position", &Window::position) + .def_property_readonly("velocity", &Window::velocity) + .def_property_readonly("total_shift", &Window::totalShift) + .def("disable_motion", &Window::disableMotion); + } \ No newline at end of file diff --git a/src/bindings/wrap_world.cpp b/src/bindings/wrap_world.cpp index 00119623..f02c260f 100644 --- a/src/bindings/wrap_world.cpp +++ b/src/bindings/wrap_world.cpp @@ -11,6 +11,7 @@ #include "system.hpp" #include "timesolver.hpp" #include "wrappers.hpp" +#include "window.hpp" /* Helper function to add any magnet instance to the world*/ template @@ -120,7 +121,8 @@ void wrap_world(py::module& m) { .def_property_readonly("antiferromagnets", &MumaxWorld::antiferromagnets, "get a map of all antiferromagnets in this world") - + .def_property_readonly("altermagnets", &MumaxWorld::altermagnets, + "get a map of all altermagnets in this world") .def_property_readonly("ncafms", &MumaxWorld::ncafms, "get a map of all non-collinear antiferromagnets in this world") @@ -144,5 +146,10 @@ void wrap_world(py::module& m) { &MumaxWorld::setMastergrid, "mastergrid of the world") .def_property("pbc_repetitions", &MumaxWorld::pbcRepetitions, &MumaxWorld::setPbcRepetitions, "PBC repetitions of the world") + + // Moving simulation window + .def_property_readonly("window", &MumaxWorld::window, + py::return_value_policy::reference) + .def("center_domain_wall", &MumaxWorld::centerDomainWall, py::arg("comp"), py::arg("axis")) ; } diff --git a/src/bindings/wrappers.hpp b/src/bindings/wrappers.hpp index 9558a2ce..19df7483 100644 --- a/src/bindings/wrappers.hpp +++ b/src/bindings/wrappers.hpp @@ -43,4 +43,5 @@ void wrap_linsolver(py::module& m); void wrap_system(py::module& m); void wrap_dmitensor(py::module& m); void wrap_voronoi(py::module& m); -void wrap_traction(py::module& m); \ No newline at end of file +void wrap_window(py::module& m); +void wrap_traction(py::module& m); diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 92a30b62..a6f46282 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -23,7 +23,9 @@ add_library(core STATIC rungekutta.cpp rungekutta.hpp scalarquantity.cpp - scalarquantity.hpp + scalarquantity.hpp + shift.cu + shift.hpp stepper.cpp stepper.hpp system.cpp @@ -34,6 +36,8 @@ add_library(core STATIC variable.cpp voronoi.cpp voronoi.hpp + window.cpp + window.hpp world.cpp world.hpp ) @@ -41,4 +45,4 @@ add_library(core STATIC target_include_directories(core PRIVATE ${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES}) target_include_directories(core PUBLIC .) -target_link_libraries(core PUBLIC cudautil) \ No newline at end of file +target_link_libraries(core PUBLIC cudautil physics) \ No newline at end of file diff --git a/src/core/shift.cu b/src/core/shift.cu new file mode 100644 index 00000000..8f8da4ab --- /dev/null +++ b/src/core/shift.cu @@ -0,0 +1,75 @@ +#include "cudaerror.hpp" +#include "cudalaunch.hpp" +#include "cudastream.hpp" +#include "field.hpp" +#include "shift.hpp" + +__global__ void k_shift_field(CuField result, + CuField field, + int dir, + int axis, + real3 leftValue, + real3 rightValue) { + const int idx = blockIdx.x * blockDim.x + threadIdx.x; + + if (!field.cellInGeometry(idx)) { return; } + + Grid grid = field.system.grid; + int3 gridsize = grid.size(); + + real3 value; + + const int3 dst = grid.index2coord(idx); + + int3 direction{0, 0, 0}; + (&direction.x)[axis] = dir; + + int3 src = dst - direction; + + // TODO: this only works because non-trivial geometries are not allowed at this point + if (field.cellInGeometry(src)) + value = field.vectorAt(src); + else { + if (dir == 1) + value = (leftValue != real3{0, 0, 0}) ? leftValue : field.vectorAt(dst); + else + value = (rightValue != real3{0, 0, 0}) ? rightValue : field.vectorAt(dst); + } + result.setVectorInCell(idx, value); +} + +real sgn(real value) { + return (value > 0.) ? 1. : -1.; +} + +int calculateShiftDirection(const Field& field, int comp, int axis, real3 leftValue, real3 rightValue) { + real av = field.average()[comp]; + int3 gridsize = field.grid().size(); + real tolerance = 4.0 / (&gridsize.x)[axis]; + + if (abs(av) > tolerance) { + // If left insertion value is absent, deduce sign of 'left' domain + if (leftValue == real3{0,0,0}) { + auto grid = field.grid(); + int3 coo = {grid.size().x / 2, grid.size().y / 2, grid.size().z / 2}; + (&coo.x)[axis] = 0; + + // TODO: What if coo not in geometry? + int idx = grid.coord2index(coo); + real value; + checkCudaError(cudaMemcpy(&value, field.device_ptr(comp) + idx, + sizeof(real), cudaMemcpyDeviceToHost)); + (&leftValue.x)[comp] = value; + } + return - sgn((&leftValue.x)[comp]) * sgn(av); + } + else { + return 0; + } +} + +Field shift(const Field& field, int dir, int comp, int axis, real3 leftValue, real3 rightValue) { + Field result(field.system(), 3); + cudaLaunch(field.grid().ncells(), k_shift_field, result.cu(), field.cu(), dir, axis, leftValue, rightValue); + return result; +} \ No newline at end of file diff --git a/src/core/shift.hpp b/src/core/shift.hpp new file mode 100644 index 00000000..ac34c847 --- /dev/null +++ b/src/core/shift.hpp @@ -0,0 +1,8 @@ +#pragma once + +#include "datatypes.hpp" + +class Field; + +int calculateShiftDirection(const Field& field, int comp, int axis, real3 leftValue, real3 rightValue); +Field shift(const Field& field, int dir, int comp, int axis, real3 leftValue, real3 RightValue); \ No newline at end of file diff --git a/src/core/timesolver.cpp b/src/core/timesolver.cpp index c2494565..81ff6b81 100644 --- a/src/core/timesolver.cpp +++ b/src/core/timesolver.cpp @@ -108,6 +108,7 @@ void TimeSolver::step() { "Timesolver can not make a step because the timestep is smaller than " "or equal to zero."); stepper_->step(); + postStep(); } void TimeSolver::steps(unsigned int nSteps) { @@ -135,5 +136,6 @@ void TimeSolver::run(real duration) { real oldTimestep = timestep(); setTimeStep(stoptime - time_); step(); + postStep(); if (fixedTimeStep_) setTimeStep(oldTimestep); } diff --git a/src/core/timesolver.hpp b/src/core/timesolver.hpp index 58069cbc..775fb348 100644 --- a/src/core/timesolver.hpp +++ b/src/core/timesolver.hpp @@ -55,6 +55,8 @@ class TimeSolver { void setUpperBound(real upperBound) { upperBound_ = upperBound; } void enableAdaptiveTimeStep() { fixedTimeStep_ = false; } void disableAdaptiveTimeStep() { fixedTimeStep_ = true; } + void setPostStepFunction(std::function func) { postStep_ = func; } + void clearPostStepFunction() { postStep_ = nullptr; } //------------- EXECUTING THE SOLVER ----------------------------------------- @@ -62,6 +64,7 @@ class TimeSolver { void steps(unsigned int nsteps); void runwhile(std::function); void run(real duration); + void postStep() { if (postStep_) return postStep_(); } //------------- HELPER FUNCTIONS FOR ADAPTIVE TIMESTEPPING ------------------- @@ -80,6 +83,7 @@ class TimeSolver { real timestep_ = 0.0; real upperBound_ = 2.0; bool fixedTimeStep_ = false; + std::function postStep_ = nullptr; std::vector eqs_; //------------- THE INTERNAL STEPPER ----------------------------------------- diff --git a/src/core/window.cpp b/src/core/window.cpp new file mode 100644 index 00000000..a6fa6b98 --- /dev/null +++ b/src/core/window.cpp @@ -0,0 +1,28 @@ +#include "datatypes.hpp" +#include "mumaxworld.hpp" +#include "shift.hpp" +#include "timesolver.hpp" +#include "window.hpp" + +Window::Window(MumaxWorld& world) + : world_(world), + position_(real3{0, 0, 0}), + velocity_(real3{0, 0, 0}), + total_dist_(real3{0, 0, 0}) { + magValues_.fill({0, 0, 0}); + } + +void Window::move(int dir, int axis, int comp) { + real3 cs = world_.cellsize(); + // Multiply dir by -1 because the field is moved to dir iff the window is moved to -dir + (&position_.x)[axis] += -1. * dir * (&cs.x)[axis]; + (&velocity_.x)[axis] = -1. * dir * (&cs.x)[axis] / world_.timesolver().timestep(); + (&total_dist_.x)[axis] += (&cs.x)[axis]; +} +Field Window::centerOnExcitation(const Field& field, int dir, int axis, int comp) { + return shift(field, dir, comp, axis, magValues_[0], magValues_[1]); +} + +void Window::disableMotion() { + world_.timesolver().clearPostStepFunction(); +} diff --git a/src/core/window.hpp b/src/core/window.hpp new file mode 100644 index 00000000..ce6c84dc --- /dev/null +++ b/src/core/window.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "datatypes.hpp" +#include "field.hpp" + +#include +#include + +enum class Boundary { + Left, + Right +}; + +class MumaxWorld; +class Window { + public: + explicit Window(MumaxWorld& world); + ~Window() = default; + + // Set origin of the simulation window + void setOrigin(real3 origin) { origin_ = origin; } + // Set values to be inserted at the boundaries + void setMagValue(Boundary side, real3 value) { magValues_[idx(side)] = value; } + // Get values to be inserted at the boundaries + std::array getMagValues() { return magValues_; } + + void move(int dir, int axis, int comp); + Field centerOnExcitation(const Field& field, int dir, int axis, int comp); + void disableMotion(); + + // Get total amount shifted + real3 position() const { return origin_ + position_; } + real3 velocity() const { return velocity_; } + real3 totalShift() const { return total_dist_; } + + private: + // Keep reference of the world to which this window belongs + MumaxWorld& world_; + real3 origin_; + // Values to insert at the boundaries + std::array magValues_; + // Current window position and velocity + real3 position_; + real3 velocity_; + // Total distance travelled by window + real3 total_dist_; + + static constexpr size_t idx(Boundary b) { return static_cast(b); } +}; \ No newline at end of file diff --git a/src/physics/ferromagnet.cpp b/src/physics/ferromagnet.cpp index 397d788d..474bb6d1 100644 --- a/src/physics/ferromagnet.cpp +++ b/src/physics/ferromagnet.cpp @@ -59,7 +59,8 @@ Ferromagnet::Ferromagnet(std::shared_ptr system_ptr, poissonSystem(this), // magnetoelasticity B1(system(), 0.0, name + ":B1", "J/m3"), - B2(system(), 0.0, name + ":B1", "J/m3") { + B2(system(), 0.0, name + ":B2", "J/m3"), + BChiral(system(), 0.0, name + ":BChiral", "J/m3") { {// Initialize random magnetization // TODO: this can be done much more efficient somewhere else int nvalues = 3 * this->grid().ncells(); diff --git a/src/physics/ferromagnet.hpp b/src/physics/ferromagnet.hpp index 1768ac0c..1e325506 100644 --- a/src/physics/ferromagnet.hpp +++ b/src/physics/ferromagnet.hpp @@ -103,4 +103,5 @@ class Ferromagnet : public Magnet { // Magnetoelasticity Parameter B1; // First magnetoelastic coupling constant Parameter B2; // Second magnetoelastic coupling constant + Parameter BChiral; // Chiral magnetoelastic coupling constant }; \ No newline at end of file diff --git a/src/physics/magnetoelasticfield.cu b/src/physics/magnetoelasticfield.cu index cbdc4e71..ebdde257 100644 --- a/src/physics/magnetoelasticfield.cu +++ b/src/physics/magnetoelasticfield.cu @@ -23,7 +23,8 @@ bool dynamicMagnetoelasticAssuredZero(const Ferromagnet* magnet) { } return (!enableElastodynamics || magnet->msat.assuredZero() || - (magnet->B1.assuredZero() && magnet->B2.assuredZero())); + (magnet->B1.assuredZero() && magnet->B2.assuredZero() && + magnet->BChiral.assuredZero())); } bool rigidMagnetoelasticAssuredZero(const Ferromagnet* magnet) { @@ -38,7 +39,8 @@ bool rigidMagnetoelasticAssuredZero(const Ferromagnet* magnet) { } return (!appliedStrain || magnet->msat.assuredZero() || - (magnet->B1.assuredZero() && magnet->B2.assuredZero())); + (magnet->B1.assuredZero() && magnet->B2.assuredZero() && + magnet->BChiral.assuredZero())); } __global__ void k_dynamicMagnetoelasticField(CuField hField, @@ -46,6 +48,7 @@ __global__ void k_dynamicMagnetoelasticField(CuField hField, const CuField strain, const CuParameter B1, const CuParameter B2, + const CuParameter BChiral, const CuParameter msat) { const int idx = blockIdx.x * blockDim.x + threadIdx.x; const CuSystem system = hField.system; @@ -66,7 +69,9 @@ __global__ void k_dynamicMagnetoelasticField(CuField hField, if (ip2 >= 3) ip2 -= 3; hField.setValueInCell(idx, i, - 2 / msat.valueAt(idx) * - (B1.valueAt(idx) * strain.valueAt(idx, i) * mField.valueAt(idx, i) + + ((B1.valueAt(idx) * strain.valueAt(idx, i) + + BChiral.valueAt(idx) * (- strain.valueAt(idx, ip1) + strain.valueAt(idx, ip2)) + ) * mField.valueAt(idx, i) + B2.valueAt(idx) * (strain.valueAt(idx, i+ip1+2) * mField.valueAt(idx, ip1) + strain.valueAt(idx, i+ip2+2) * mField.valueAt(idx, ip2)))); } @@ -78,6 +83,7 @@ __global__ void k_rigidMagnetoelasticField(CuField hField, const CuVectorParameter shearStrain, const CuParameter B1, const CuParameter B2, + const CuParameter BChiral, const CuParameter msat) { const int idx = blockIdx.x * blockDim.x + threadIdx.x; const CuSystem system = hField.system; @@ -98,7 +104,9 @@ __global__ void k_rigidMagnetoelasticField(CuField hField, if (ip2 >= 3) ip2 -= 3; hField.setValueInCell(idx, i, - 2 / msat.valueAt(idx) * - (B1.valueAt(idx) * normStrain.valueAt(idx, i) * mField.valueAt(idx, i) + + ((B1.valueAt(idx) * normStrain.valueAt(idx, i) + + BChiral.valueAt(idx) * (- normStrain.valueAt(idx, ip1) + normStrain.valueAt(idx, ip2)) + ) * mField.valueAt(idx, i) + B2.valueAt(idx) * (shearStrain.valueAt(idx, i+ip1-1) * mField.valueAt(idx, ip1) + shearStrain.valueAt(idx, i+ip2-1) * mField.valueAt(idx, ip2)))); } @@ -115,6 +123,7 @@ Field evalMagnetoelasticField(const Ferromagnet* magnet) { CuField mField = magnet->magnetization()->field().cu(); CuParameter B1 = magnet->B1.cu(); CuParameter B2 = magnet->B2.cu(); + CuParameter BChiral = magnet->BChiral.cu(); CuParameter msat = magnet->msat.cu(); if (!rigidMagnetoelasticAssuredZero(magnet)) { // maybe use rigid strain @@ -123,13 +132,13 @@ Field evalMagnetoelasticField(const Ferromagnet* magnet) { CuVectorParameter shearStrain = magnet->hostMagnet()->rigidShearStrain.cu(); cudaLaunch(ncells, k_rigidMagnetoelasticField, hField.cu(), mField, - normStrain, shearStrain, B1, B2, msat); + normStrain, shearStrain, B1, B2, BChiral, msat); } else { // independent magnet CuVectorParameter normStrain = magnet->rigidNormStrain.cu(); CuVectorParameter shearStrain = magnet->rigidShearStrain.cu(); cudaLaunch(ncells, k_rigidMagnetoelasticField, hField.cu(), mField, - normStrain, shearStrain, B1, B2, msat); + normStrain, shearStrain, B1, B2, BChiral, msat); } return hField; @@ -144,7 +153,7 @@ Field evalMagnetoelasticField(const Ferromagnet* magnet) { } cudaLaunch(ncells, k_dynamicMagnetoelasticField, hField.cu(), mField, - strain.cu(), B1, B2, msat); + strain.cu(), B1, B2, BChiral, msat); return hField; } diff --git a/src/physics/magnetoelasticforce.cu b/src/physics/magnetoelasticforce.cu index 903fe4b7..e1522e54 100644 --- a/src/physics/magnetoelasticforce.cu +++ b/src/physics/magnetoelasticforce.cu @@ -11,6 +11,7 @@ __global__ void k_magnetoelasticForce(CuField fField, const CuField m, const CuParameter B1, const CuParameter B2, + const CuParameter BChiral, const real3 w, // w = 1/cellsize const Grid mastergrid) { const int idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -89,6 +90,7 @@ __global__ void k_magnetoelasticForce(CuField fField, real f_i = 2 * B1.valueAt(idx) * m_here[i] * der[i][i]; f_i += B2.valueAt(idx) * m_here[i] * (der[ip1][ip1] + der[ip2][ip2]); f_i += B2.valueAt(idx) * (m_here[ip1] * der[ip1][i] + m_here[ip2] * der[ip2][i]); + f_i += BChiral.valueAt(idx) * (m_here[ip1] * der[i][ip1] - m_here[ip2] * der[i][ip2]); fField.setValueInCell(idx, i, f_i); } } @@ -104,9 +106,10 @@ Field evalMagnetoelasticForce(const Ferromagnet* magnet) { CuField m = magnet->magnetization()->field().cu(); CuParameter B1 = magnet->B1.cu(); CuParameter B2 = magnet->B2.cu(); + CuParameter BChiral = magnet->BChiral.cu(); real3 w = 1 / magnet->cellsize(); Grid mastergrid = magnet->world()->mastergrid(); - cudaLaunch(ncells, k_magnetoelasticForce, fField.cu(), m, B1, B2, w, mastergrid); + cudaLaunch(ncells, k_magnetoelasticForce, fField.cu(), m, B1, B2, BChiral, w, mastergrid); return fField; } diff --git a/src/physics/mumaxworld.cpp b/src/physics/mumaxworld.cpp index 60a672c4..85080425 100644 --- a/src/physics/mumaxworld.cpp +++ b/src/physics/mumaxworld.cpp @@ -16,20 +16,24 @@ #include "minimizer.hpp" #include "ncafm.hpp" #include "relaxer.hpp" +#include "shift.hpp" #include "system.hpp" #include "thermalnoise.hpp" #include "timesolver.hpp" #include "torque.hpp" +#include "window.hpp" MumaxWorld::MumaxWorld(real3 cellsize) : World(cellsize), biasMagneticField({0, 0, 0}), - RelaxTorqueThreshold(-1.0) {} + RelaxTorqueThreshold(-1.0), + window_(std::make_unique(*this)) {} MumaxWorld::MumaxWorld(real3 cellsize, Grid mastergrid, int3 pbcRepetitions) : World(cellsize, mastergrid, pbcRepetitions), biasMagneticField({0, 0, 0}), - RelaxTorqueThreshold(-1.0) {} + RelaxTorqueThreshold(-1.0), + window_(std::make_unique(*this)) {} MumaxWorld::~MumaxWorld() {} @@ -327,3 +331,42 @@ void MumaxWorld::unsetPBC() { } // -------------------------------------------------- +// Moving simulation window + +void MumaxWorld::centerDomainWall(int comp, int axis) { + if (magnets_.size() > 1) + throw std::runtime_error("Moving the simulation window is only possible when only one " + "magnet lives in the world."); + if (magnets_.size() < 1) + throw std::runtime_error("Moving the simulation window is not possible when there is no " + "magnet in the world."); + + Magnet* magnet = magnets_.begin()->second; + timesolver_->setPostStepFunction([this, magnet, comp, axis]() { + const Field& mag = magnet->asHost() ? magnet->asHost()->sublattices()[0]->magnetization()->field() + : magnet->asFM()->magnetization()->field(); + int dir = calculateShiftDirection(mag, + comp, axis, + window_->getMagValues()[0], + window_->getMagValues()[1]); + if (dir != 0) { + window_->move(dir, axis, comp); + // Shift magnetization + auto shifted = window_->centerOnExcitation(mag, dir, axis, comp); + + // Multi-sublattice systems + if (auto host = magnet->asHost()) { + auto sub0 = host->sublattices()[0]; + sub0->magnetization()->set(shifted); + for (auto sub : host->getOtherSublattices(sub0)) { + auto shifted = window_->centerOnExcitation(sub->magnetization()->field(), dir, axis, comp); + sub->magnetization()->set(shifted); + } + } + + // Ferromagnet + else + magnet->asFM()->magnetization()->set(shifted); + } + }); +} \ No newline at end of file diff --git a/src/physics/mumaxworld.hpp b/src/physics/mumaxworld.hpp index ffe90386..50604e7a 100644 --- a/src/physics/mumaxworld.hpp +++ b/src/physics/mumaxworld.hpp @@ -11,6 +11,7 @@ #include "gpubuffer.hpp" #include "grid.hpp" #include "torque.hpp" +#include "window.hpp" #include "world.hpp" class Altermagnet; @@ -128,6 +129,7 @@ class MumaxWorld : public World { magnets_[name] = raw; handleNewStrayfield(raw); + window_->setOrigin(int3_to_real3(grid.origin()) * this->cellsize()); return raw; } @@ -259,6 +261,10 @@ class MumaxWorld : public World { // -------------------------------------------------- + // Moving simulation window + Window& window() const { return *window_; } + void centerDomainWall(int comp, int axis); + private: std::map magnets_; @@ -267,4 +273,6 @@ class MumaxWorld : public World { std::map> antiferromagnets_; std::map> altermagnets_; std::map> ncafms_; -}; + + std::unique_ptr window_; +}; \ No newline at end of file diff --git a/test/test_magnetoelasticfield.py b/test/test_magnetoelasticfield.py index 2d2652e2..c300a59b 100644 --- a/test/test_magnetoelasticfield.py +++ b/test/test_magnetoelasticfield.py @@ -45,14 +45,14 @@ def create_magnet(d_comp, m_comp): return magnet -def sine_displacement(magnet, i_comp, j_comp, B1, B2): - """Creates the displacement following a sine in the d_comp direction - and cosine in another direction it then calculates the magnetoelasticforce. +def sine_displacement(magnet, i_comp, j_comp, B1=0, B2=0, Bc=0): + """Creates a displacement in the j_comp direction following a sine along the + i_comp direction. Then calculates and compares analytical and numerical + magnetoelastic field. """ - magnet.enable_elastodynamics = True # just in case - magnet.B1 = B1 magnet.B2 = B2 + magnet.B_chiral = Bc L = N*cellsize[i_comp] k = P*2*math.pi/L @@ -90,64 +90,93 @@ def displacement_func(x, y, z): B_anal[i,...] = - 2 / msat * ( B1 * strain_anal[i,...] * m[i,...] + + Bc * (-strain_anal[ip1,...] + strain_anal[ip2,...]) * m[i,...] + B2 * (strain_anal[i+ip1+2,...] * m[ip1,...] + strain_anal[i+ip2+2,...] * m[ip2,...])) assert max_semirelative_error(B_num, B_anal) < RTOL -def test_x_Exx_B1(): +# --- B1 --- + +def test_mx_Exx_B1(): m_comp, i, j = 0, 0, 0 - B1, B2 = B, 0 magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B1=B) -def test_y_Eyy_B1(): +def test_my_Eyy_B1(): m_comp, i, j = 1, 1, 1 - B1, B2 = B, 0 magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B1=B) -def test_z_Ezz_B1(): +def test_mz_Ezz_B1(): m_comp, i, j = 2, 2, 2 - B1, B2 = B, 0 magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B1=B) + +# --- B2 --- -def test_y_Exy_B2(): +def test_my_Exy_B2(): m_comp, i, j = 1, 0, 1 - B1, B2 = 0, B magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B2=B) -def test_z_Exz_B2(): +def test_mz_Exz_B2(): m_comp, i, j = 2, 0, 2 - B1, B2 = 0, B magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B2=B) -def test_x_Exy_B2(): +def test_mx_Exy_B2(): m_comp, i, j = 0, 0, 1 - B1, B2 = 0, B magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B2=B) -def test_z_Eyz_B2(): +def test_mz_Eyz_B2(): m_comp, i, j = 2, 1, 2 - B1, B2 = 0, B magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B2=B) -def test_x_Exz_B2(): +def test_mx_Exz_B2(): m_comp, i, j = 0, 0, 2 - B1, B2 = 0, B magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B2=B) -def test_y_Eyz_B2(): +def test_my_Eyz_B2(): m_comp, i, j = 1, 1, 2 - B1, B2 = 0, B magnet = create_magnet(i, m_comp) - sine_displacement(magnet, i, j, B1, B2) + sine_displacement(magnet, i, j, B2=B) + +# --- Bc --- + +def test_mx_Eyy_Bc(): + m_comp, i, j = 0, 1, 1 + magnet = create_magnet(i, m_comp) + sine_displacement(magnet, i, j, Bc=B) + +def test_mx_Ezz_Bc(): + m_comp, i, j = 0, 2, 2 + magnet = create_magnet(i, m_comp) + sine_displacement(magnet, i, j, Bc=B) + +def test_my_Exx_Bc(): + m_comp, i, j = 1, 0, 0 + magnet = create_magnet(i, m_comp) + sine_displacement(magnet, i, j, Bc=B) + +def test_my_Ezz_Bc(): + m_comp, i, j = 1, 2, 2 + magnet = create_magnet(i, m_comp) + sine_displacement(magnet, i, j, Bc=B) + +def test_mz_Exx_Bc(): + m_comp, i, j = 2, 0, 0 + magnet = create_magnet(i, m_comp) + sine_displacement(magnet, i, j, Bc=B) + +def test_mz_Eyy_Bc(): + m_comp, i, j = 2, 1, 1 + magnet = create_magnet(i, m_comp) + sine_displacement(magnet, i, j, Bc=B) + def test_random(): """Test with a random magnetization and displacement. @@ -163,12 +192,10 @@ def test_random(): magnet.enable_elastodynamics = True magnet.msat = msat - B1, B2 = B, 0.5*B + B1, B2, Bc = B, 0.5*B, 0.3*B magnet.B1 = B1 magnet.B2 = B2 - - L = N*cellsize[0] - k = P*2*math.pi/L + magnet.B_chiral = Bc def displacement_func(x, y, z): return tuple(np.random.rand(3)) @@ -186,7 +213,44 @@ def displacement_func(x, y, z): B_anal[i,...] = - 2 / msat * ( B1 * strain[i,...] * m[i,...] + + Bc * (-strain[ip1,...] + strain[ip2,...]) * m[i,...] + B2 * (strain[i+ip1+2,...] * m[ip1,...] + strain[i+ip2+2,...] * m[ip2,...])) - assert max_semirelative_error(B_num, B_anal) < RTOL \ No newline at end of file + assert max_semirelative_error(B_num, B_anal) < 1e-6 + + +def test_rigid_magnetoelastic_field(): + """Test with a random magnetization and rigid norm and shear strain. + """ + nx, ny, nz = N, 4, 2 # any shape would do + + world = World(cellsize) + magnet = Ferromagnet(world, Grid((nx, ny, nz))) + + magnet.msat = msat + B1, B2, Bc = B, 0.5*B, 0.3*B + magnet.B1 = B1 + magnet.B2 = B2 + magnet.B_chiral = Bc + + norm_strain = np.random.rand(3, nz, ny, nx) + shear_strain = np.random.rand(3, nz, ny, nx) + magnet.rigid_norm_strain = norm_strain + magnet.rigid_shear_strain = shear_strain + + B_num = magnet.magnetoelastic_field.eval() + B_anal = np.zeros(shape=B_num.shape) + + m = magnet.magnetization.eval() + for i in range(3): + ip1 = (i+1)%3 + ip2 = (i+2)%3 + + B_anal[i,...] = - 2 / msat * ( + B1 * norm_strain[i,...] * m[i,...] + + Bc * (-norm_strain[ip1,...] + norm_strain[ip2,...]) * m[i,...] + + B2 * (shear_strain[i+ip1-1,...] * m[ip1,...] + + shear_strain[i+ip2-1,...] * m[ip2,...])) + + assert max_semirelative_error(B_num, B_anal) < 1e-6 \ No newline at end of file diff --git a/test/test_magnetoelasticforce.py b/test/test_magnetoelasticforce.py index 2490f1ad..3dabbb0d 100644 --- a/test/test_magnetoelasticforce.py +++ b/test/test_magnetoelasticforce.py @@ -22,7 +22,7 @@ def max_semirelative_error(result, wanted): return max_absolute_error(result, wanted) / np.max(abs(wanted)) -def create_magnet(d_comp, B1, B2): +def create_magnet(d_comp): """Makes a world with a 1D magnet in the d_comp direction. """ gridsize, gridsize_magnet, pbc_repetitions = [0, 0, 0], [1, 1, 1], [0, 0, 0] @@ -34,17 +34,16 @@ def create_magnet(d_comp, B1, B2): magnet = Ferromagnet(world, Grid(gridsize_magnet)) magnet.enable_elastodynamics = True - magnet.B1 = B1 - magnet.B2 = B2 - - return world, magnet + return magnet -def sine_force(magnet, d_comp, comp_cos, B1, B2): +def sine_force(magnet, d_comp, comp_cos, B1=0, B2=0, Bc=0): """Creates the magnetization following a sine in the d_comp direction and cosine in another direction it then calculates the magnetoelasticforce. """ - magnet.enable_elastodynamics = True # just in case + magnet.B1 = B1 + magnet.B2 = B2 + magnet.B_chiral = Bc L = N*cellsize[d_comp] k = P*2*math.pi/L @@ -57,66 +56,95 @@ def magnetization_func(x, y, z): return tuple(m) magnet.magnetization = magnetization_func + mag = magnet.magnetization.eval() force_num = magnet.magnetoelastic_force.eval() force_anal = np.zeros(shape=force_num.shape) - force_anal[d_comp,...] = 2 * B1 * k * magnet.magnetization.eval()[d_comp,...] * magnet.magnetization.eval()[comp_cos,...] - force_anal[comp_cos,...] = B2 * k * (magnet.magnetization.eval()[comp_cos,...]**2 - magnet.magnetization.eval()[d_comp,...]**2) + chiral_sign = -1 if comp_cos == (d_comp + 1)%3 else 1 + force_anal[d_comp,...] = 2 * B1 * k * mag[d_comp,...] * mag[comp_cos,...] + \ + chiral_sign * Bc * k * mag[comp_cos,...] * mag[d_comp,...] + force_anal[comp_cos,...] = B2 * k * (mag[comp_cos,...]**2 - mag[d_comp,...]**2) assert max_semirelative_error(force_num, force_anal) < RTOL +# --- B1 --- def test_sinx_cosx_0_B1(): d_comp = 0 - B1, B2 = B, 0 - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+1)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, B1=B) def test_0_siny_cosy_B1(): d_comp = 1 - B1, B2 = B, 0 - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+1)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, B1=B) def test_cosz_0_sinz_B1(): d_comp = 2 - B1, B2 = B, 0 - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+1)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, B1=B) + +# --- B2 --- def test_sinx_cosx_0_B2(): d_comp = 0 - B1, B2 = 0, B - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+1)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, B2=B) def test_0_siny_cosy_B2(): d_comp = 1 - B1, B2 = 0, B - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+1)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, B2=B) def test_cosz_0_sinz_B2(): d_comp = 2 - B1, B2 = 0, B - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+1)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, B2=B) def test_sinx_0_cosx_B2(): d_comp = 0 - B1, B2 = 0, B - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+2)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+2)%3, B2=B) def test_cosy_siny_0_B2(): d_comp = 1 - B1, B2 = 0, B - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+2)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+2)%3, B2=B) def test_0_cosz_sinz_B2(): d_comp = 2 - B1, B2 = 0, B - world, magnet = create_magnet(d_comp, B1, B2) - sine_force(magnet, d_comp, (d_comp+2)%3, B1, B2) + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+2)%3, B2=B) + +# --- Bc --- + +def test_sinx_cosx_0_Bc(): + d_comp = 0 + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, Bc=B) + +def test_0_siny_cosy_Bc(): + d_comp = 1 + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, Bc=B) + +def test_cosz_0_sinz_Bc(): + d_comp = 2 + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+1)%3, Bc=B) + +def test_sinx_0_cosx_Bc(): + d_comp = 0 + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+2)%3, Bc=B) + +def test_cosy_siny_0_Bc(): + d_comp = 1 + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+2)%3, Bc=B) + +def test_0_cosz_sinz_Bc(): + d_comp = 2 + magnet = create_magnet(d_comp) + sine_force(magnet, d_comp, (d_comp+2)%3, Bc=B) diff --git a/test/test_moving_window.py b/test/test_moving_window.py new file mode 100644 index 00000000..544a273c --- /dev/null +++ b/test/test_moving_window.py @@ -0,0 +1,262 @@ +import numpy as np +import pytest +from mumaxplus import Antiferromagnet, Ferromagnet, Grid, World + +cs = 1e-9 +runtime = 1e-12 + +def make_ferromagnet(nx=64, ny=1, nz=1): + world = World((cs, cs, cs)) + grid = Grid((nx, ny, nz)) + magnet = Ferromagnet(world, grid) + return world, magnet + +def set_parameters(magnet, comp): + magnet.enable_demag = False + magnet.enable_openbc = True + magnet.msat = 800e3 + magnet.aex = 13e-12 + magnet.alpha = 0.1 + magnet.ku1 = 5e5 + anisU = [0.0, 0.0, 0.0] + anisU[comp] = 1.0 + magnet.anisU = tuple(anisU) + magnet.xi = 0.2 + magnet.pol = 1 + + +def dw_profile(magnet, comp, axis, width=5, center=None, minimize=True): + nz, ny, nx = magnet.grid.shape + + n = (nx, ny, nz)[axis] + if center is None: + center = int(0.5 * n) + + def slc(start, stop): + s = [slice(None), slice(None), slice(None)] + s[2 - axis] = slice(start, stop) + return tuple(s) + + m = np.zeros((3, nz, ny, nx)) + + m[comp][slc(None, center - width)] = 1 # 'left' domain + m[comp][slc(center + width, None)] = -1 # 'right' domain + m[(comp + 1) % 3][slc(center - width, center + width)] = 1 # domain wall + + magnet.magnetization = m + if minimize: + magnet.minimize() + +class TestValidArguments: + def test_valid_boundary(self): + world, _ = make_ferromagnet() + for b in (0, 1): + world.window.insert_magnetization(b, (1, 0, 0)) # no raise + + def test_invalid_boundary(self): + world, _ = make_ferromagnet() + for b in (-1, 4, 1.5): + with pytest.raises((ValueError)): + world.window.insert_magnetization(b, (1, 0, 0)) + + def test_invalid_comp(self): + world, _ = make_ferromagnet() + for c in (-1, 4, 1.5): + with pytest.raises((ValueError)): + world.center_domain_wall(c) +class TestInitialConditions: + def test_geometry(self): + w, g = World((1, 1, 1)), Grid((10, 10, 10)) + geo = np.ones(g.shape) + geo[0, 0, 0] = 0 + magnet = Ferromagnet(w, g, geometry=geo) + with pytest.raises((RuntimeError)): + w.center_domain_wall(0, 0) + + def test_regions(self): + w, g = World((1, 1, 1)), Grid((10, 10, 10)) + reg = np.zeros(g.shape) + reg[0, 0, 0] = 1 + magnet = Ferromagnet(w, g, regions=reg) + with pytest.raises((RuntimeError)): + w.center_domain_wall(0, 0) + + def test_multiple_magnets(self): + world = World((1, 1, 1)) + magnet_1 = Ferromagnet(world, Grid((10, 10, 1))) + magnet_2 = Antiferromagnet(world, Grid((10, 10, 1), origin=(0, 0, 1))) + with pytest.raises((RuntimeError)): + world.center_domain_wall(0, 0) + + def test_origin(self): + origin = (5, 6, 7) + cellsize = (2e-9, 3e-9, 4e-9) + origin_phys = np.array(origin) * np.array(cellsize) + world = World(cellsize) + magnet = Ferromagnet(world, Grid((11, 12, 13), origin=origin)) + assert np.all(np.isclose(world.window.position, origin_phys, atol=cs/10)) + + def test_offset_domain_wall(self): + world, magnet = make_ferromagnet(nx=100) + dw_profile(magnet, 0, 0, center=20) + with pytest.raises((RuntimeError)): + world.center_domain_wall(0, 0) + + def test_zero_average(self): + Nx = 128 + world, magnet = make_ferromagnet(Nx, 1, 1) + set_parameters(magnet, 0) + dw_profile(magnet, 0, 0) + + world.center_domain_wall(0, 0) + magnet.jcur = (1e12, 0, 0) + world.timesolver.run(runtime) + av = magnet.magnetization.average()[0] + assert np.abs(av) < 4 / Nx + +class TestShift: + def setup_and_shift(self, comp=0, axis=0, nx=64, ny=1, nz=1, current=1e14): + world, magnet = make_ferromagnet(nx, ny, nz) + set_parameters(magnet, comp) + dw_profile(magnet, comp, axis) + + jcur = [0, 0, 0] + jcur[axis] = current + magnet.jcur = jcur + + world.center_domain_wall(comp, axis) + world.timesolver.run(runtime) + return world.window.position + + def test_axis0(self): + shift = self.setup_and_shift(comp=0, axis=0, nx=64, ny=1, nz=1) + assert abs(shift[0]) >= cs + assert shift[1] == 0.0 + assert shift[2] == 0.0 + + def test_axis1(self): + shift = self.setup_and_shift(comp=0, axis=1, nx=1, ny=64, nz=1) + assert abs(shift[1]) >= cs + assert shift[0] == 0.0 + assert shift[2] == 0.0 + + def test_axis2(self): + shift = self.setup_and_shift(comp=0, axis=2, nx=1, ny=1, nz=64) + assert abs(shift[2]) >= cs + assert shift[0] == 0.0 + assert shift[1] == 0.0 + + def test_comp0(self): + shift = self.setup_and_shift(comp=0, axis=0) + assert abs(shift[0]) >= cs + + def test_comp1(self): + shift = self.setup_and_shift(comp=1, axis=0) + assert abs(shift[0]) >= cs + + def test_comp2(self): + shift = self.setup_and_shift(comp=2, axis=0) + assert abs(shift[0]) >= cs + + def test_no_current(self): + shift = self.setup_and_shift(comp=0, axis=0, current=0) + assert np.all(np.isclose(shift, 0.0, atol=cs/10)) + + def test_opposite_current(self): + shift_right = self.setup_and_shift(axis=0, current=1e14) + shift_left = self.setup_and_shift(axis=0, current=-1e14) + assert np.sign(shift_right[0]) != np.sign(shift_left[0]) + + def test_disable_motion(self): + Nx = 64 + world, magnet = make_ferromagnet(Nx, 1, 1) + set_parameters(magnet, 0) + dw_profile(magnet, 0, 0) + magnet.jcur = (1e14, 0, 0) + + world.center_domain_wall(0, 0) + world.timesolver.run(runtime) + assert magnet.magnetization.average()[0] < 4 / Nx + + world.window.disable_motion() + + world.timesolver.run(runtime) + av = magnet.magnetization.average()[0] + assert np.abs(av) > 4 / Nx + +class TestInsertion: + def setup(self, comp, axis=0): + nx, ny, nz = 64, 1, 1 + world, magnet = make_ferromagnet(nx, ny, nz) + set_parameters(magnet, comp) + dw_profile(magnet, comp, axis) + return world, magnet + + def test_carry_magnetization(self): + comp, axis = 2, 0 + world, magnet = self.setup(comp, axis) + left_init = magnet.magnetization()[comp, :, :, 0].squeeze() + right_init = magnet.magnetization()[comp, :, :, -1].squeeze() + + world.center_domain_wall(comp, axis) + + magnet.jcur = (1e12, 0, 0) + world.timesolver.run(runtime) + + left = magnet.magnetization()[comp, :, :, 0].squeeze() + right = magnet.magnetization()[comp, :, :, -1].squeeze() + + assert np.allclose((left_init, right_init), (left, right), atol=cs/10) + + def test_insertion_values(self): + comp, axis = 2, 0 + world, magnet = self.setup(comp, axis) + left_init = magnet.magnetization()[comp, :, :, 0].squeeze() + right_init = magnet.magnetization()[comp, :, :, -1].squeeze() + + world.window.insert_magnetization(0, np.array([0, 0, left_init])) + world.window.insert_magnetization(1, np.array([0, 0, right_init])) + + world.center_domain_wall(comp, axis) + magnet.jcur = (10e12, 0, 0) + world.timesolver.run(runtime) + + left = magnet.magnetization()[comp, :, :, 0].squeeze() + right = magnet.magnetization()[comp, :, :, -1].squeeze() + assert np.allclose((left_init, right_init), (left, right), atol=cs/10) + +class TestAntiferromagnet: + + def setup(self, comp=0, axis=0): + nx, ny, nz = 64, 1, 1 + world = World((cs, cs, cs)) + magnet = Antiferromagnet(world, Grid((nx, ny, nz))) + set_parameters(magnet, comp) + + magnet.afmex_cell = -10e-12 + magnet.afmex_nn = -5e-12 + + dw_profile(magnet.sub1, comp, axis, minimize=False) + magnet.sub2.magnetization = - magnet.sub1.magnetization() + magnet.minimize() + + jcur = [0.0, 0.0, 0.0] + jcur[axis] = 5e12 + magnet.sub1.jcur = tuple(jcur) + magnet.sub2.jcur = tuple(jcur) + + return world, magnet + + def test_same_shift(self): + comp, axis = 0, 0 + world, magnet = self.setup(comp, axis) + m1_before = magnet.sub1.magnetization() + m2_before = magnet.sub2.magnetization() + + world.center_domain_wall(comp, axis) + world.timesolver.run(runtime) + + diff1 = np.abs(magnet.sub1.magnetization() - m1_before) + diff2 = np.abs(magnet.sub2.magnetization() - m2_before) + + assert np.allclose(diff1, diff2, atol=cs/10) \ No newline at end of file