From 126c2a296c1389fcb1b3a516573056e9edd51089 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 09:34:39 -0700
Subject: [PATCH 01/30] docs

---
 README.md | 36 +++++++++++++++++++++++++++++++++++-
 1 file changed, 35 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index ffe91c0..9f986db 100644
--- a/README.md
+++ b/README.md
@@ -1 +1,35 @@
-# Micromagnetics
\ No newline at end of file
+# MagneX
+MagneX is a massively parallel, 3D micromagnetics solver for modeling magnetic materials.
+MagneX solves the Landau-Lifshitz-Gilbert (LLG) equations, including exchange, anisotropy, demagnetization, and Dzyaloshinskii-Moriya interaction (DMI) coupling.
+The algorithm is implemented using Exascale Computing Project software framework, AMReX, which provides effective scalability on manycore and GPU-based supercomputing architectures.
+
+# Installation
+## Download AMReX Repository
+``` git clone git@github.com:AMReX-Codes/amrex.git ```
+## Download MagneX Repository
+``` git@github.com:AMReX-Microelectronics/MagneX.git ```
+## Build
+Make sure that the AMReX and MagneX are cloned in the same location in their filesystem. Navogate to the Exec folder of MagneX and execute
+```make -j 4```
+
+# Running MagneX
+Example input scripts are located in `Exec/standard_problem_inputs/` directory. 
+## Simple Testcase
+You can run the following to simulate muMAG Standard Problem 4 dynamics:
+## For pure MPI build (but with a single MPI rank)
+```./main3d.gnu.MPI.ex standard_problem_inputs/inputs_std4```
+# Visualization and Data Analysis
+Refer to the following link for several visualization tools that can be used for AMReX plotfiles. 
+
+[Visualization](https://amrex-codes.github.io/amrex/docs_html/Visualization_Chapter.html)
+
+### Data Analysis in Python using yt 
+You can extract the data in numpy array format using yt (you can refer to this for installation and usage of [yt](https://yt-project.org/). After you have installed yt, you can do something as follows, for example, to get variable 'Pz' (z-component of polarization)
+```
+import yt
+ds = yt.load('./plt00001000/') # for data at time step 1000
+ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+P_array = ad0['Pz'].to_ndarray()
+```
+# Publications
+1. Z. Yao, P. Kumar, J. C. LePelch, and A. Nonaka, MagneX: An Exascale-Enabled Micromagnetics Solver for Spintronic Systems, in preparation.

From bfe6903e8fb19b5935a8e17cbaa38cec2a3536b1 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 09:42:30 -0700
Subject: [PATCH 02/30] docs update

---
 README.md | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index 9f986db..7c40f02 100644
--- a/README.md
+++ b/README.md
@@ -4,12 +4,15 @@ MagneX solves the Landau-Lifshitz-Gilbert (LLG) equations, including exchange, a
 The algorithm is implemented using Exascale Computing Project software framework, AMReX, which provides effective scalability on manycore and GPU-based supercomputing architectures.
 
 # Installation
-## Download AMReX Repository
+## Download AMReX and MagneX Repositories
+Make sure that AMReX and MagneX are cloned at the same root location.
 ``` git clone git@github.com:AMReX-Codes/amrex.git ```
-## Download MagneX Repository
-``` git@github.com:AMReX-Microelectronics/MagneX.git ```
+``` git clone git@github.com:AMReX-Microelectronics/MagneX.git ```
+## Dependencies
+Beyond a standard Ubuntu22 installation, the Ubuntu packages libfftw3-dev and libfftw3-mpi-dev are required.
+Also heFFTe is required.
 ## Build
-Make sure that the AMReX and MagneX are cloned in the same location in their filesystem. Navogate to the Exec folder of MagneX and execute
+ Navogate to the Exec folder of MagneX and execute
 ```make -j 4```
 
 # Running MagneX

From 84f19dcc4660c5e337d33786e7548c38485a16bc Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 09:43:32 -0700
Subject: [PATCH 03/30] docs update

---
 README.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 7c40f02..f4f7401 100644
--- a/README.md
+++ b/README.md
@@ -5,8 +5,8 @@ The algorithm is implemented using Exascale Computing Project software framework
 
 # Installation
 ## Download AMReX and MagneX Repositories
-Make sure that AMReX and MagneX are cloned at the same root location.
-``` git clone git@github.com:AMReX-Codes/amrex.git ```
+Make sure that AMReX and MagneX are cloned at the same root location. \
+``` git clone git@github.com:AMReX-Codes/amrex.git ``` \
 ``` git clone git@github.com:AMReX-Microelectronics/MagneX.git ```
 ## Dependencies
 Beyond a standard Ubuntu22 installation, the Ubuntu packages libfftw3-dev and libfftw3-mpi-dev are required.

From d71eb798575603413bf53335b5c4919cfd76ea77 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 09:49:12 -0700
Subject: [PATCH 04/30] docs update

---
 README.md | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/README.md b/README.md
index f4f7401..08ac0b7 100644
--- a/README.md
+++ b/README.md
@@ -4,15 +4,23 @@ MagneX solves the Landau-Lifshitz-Gilbert (LLG) equations, including exchange, a
 The algorithm is implemented using Exascale Computing Project software framework, AMReX, which provides effective scalability on manycore and GPU-based supercomputing architectures.
 
 # Installation
+Here are instructions for a basic, pure-MPI (no GPU) installation.  More detailed instructions for GPU systems are in the full documentation.
 ## Download AMReX and MagneX Repositories
 Make sure that AMReX and MagneX are cloned at the same root location. \
-``` git clone git@github.com:AMReX-Codes/amrex.git ``` \
-``` git clone git@github.com:AMReX-Microelectronics/MagneX.git ```
+``` git clone https://github.com/AMReX-Codes/amrex.git ``` \
+``` git clone https://AMReX-Microelectronics/MagneX.git ```
 ## Dependencies
-Beyond a standard Ubuntu22 installation, the Ubuntu packages libfftw3-dev and libfftw3-mpi-dev are required.
-Also heFFTe is required.
+Beyond a standard Ubuntu22 installation, the Ubuntu packages libfftw3-dev, libfftw3-mpi-dev, and cmake are required.
+Also heFFTe is required.  At the same level that AMReX and MagneX are cloned, run: \
+``` git clone https://github.com/icl-utk-edu/heffte.git ```\
+``` cd heffte ```\
+``` mkdir build ```\
+``` cd build ```\
+``` cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_STANDARD=17 -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=. -DHeffte_ENABLE_FFTW=ON -DHeffte_ENABLE_CUDA=OFF .. ```\
+``` make -j4 ```\
+``` make install```
 ## Build
- Navogate to the Exec folder of MagneX and execute
+ Navigate to the Exec folder of MagneX and execute
 ```make -j 4```
 
 # Running MagneX

From 7fe43b2c412a3281fa9492d56959b2863e787510 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 09:56:54 -0700
Subject: [PATCH 05/30] use forward euler to simple build works

---
 Exec/standard_problem_inputs/inputs_std4 | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/Exec/standard_problem_inputs/inputs_std4 b/Exec/standard_problem_inputs/inputs_std4
index 7bc126c..47f4ea5 100644
--- a/Exec/standard_problem_inputs/inputs_std4
+++ b/Exec/standard_problem_inputs/inputs_std4
@@ -6,10 +6,10 @@ max_grid_size_x = 160
 max_grid_size_y = 40
 max_grid_size_z = 1
 dt = 5.e-15  # forward Euler (685s), 1.e-14 unstable
-dt = 2.5e-13 # Predictor-corrector (70s), 5.e-13 unstable
-dt = 1.e-13 # trapezoidal (63s), 2.5e-13 unstable
-dt = 2.5e-13 # SSPRK3 (36s), 5.e-13 unstable
-dt = 5.e-13 # RK4 (24s), 1.e-12 unstable
+#dt = 2.5e-13 # Predictor-corrector (70s), 5.e-13 unstable
+#dt = 1.e-13 # trapezoidal (63s), 2.5e-13 unstable
+#dt = 2.5e-13 # SSPRK3 (36s), 5.e-13 unstable
+#dt = 5.e-13 # RK4 (24s), 1.e-12 unstable
 
 ##############
 # 1.5625nm case
@@ -73,10 +73,10 @@ DMI_coupling = 0
 
 # INTEGRATION
 
-#TimeIntegratorOption = 1 #Forward Euler
+TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-TimeIntegratorOption = 4 #amrex/sundials backend integrators
+#TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
@@ -98,7 +98,7 @@ integration.type = RungeKutta
 ### 2 = Trapezoid Method
 ### 3 = SSPRK3 Method
 ### 4 = RK4 Method
-integration.rk.type = 4
+integration.rk.type = 1
 
 ## If using a user-specified Butcher Tableau, then
 ## set nodes, weights, and table entries here:

From d80cf56ff3025f24c2c46306d01631c07b1bd83c Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 10:00:05 -0700
Subject: [PATCH 06/30] Update README.md

---
 README.md | 31 +++++++++++++++----------------
 1 file changed, 15 insertions(+), 16 deletions(-)

diff --git a/README.md b/README.md
index 08ac0b7..bef5634 100644
--- a/README.md
+++ b/README.md
@@ -7,28 +7,27 @@ The algorithm is implemented using Exascale Computing Project software framework
 Here are instructions for a basic, pure-MPI (no GPU) installation.  More detailed instructions for GPU systems are in the full documentation.
 ## Download AMReX and MagneX Repositories
 Make sure that AMReX and MagneX are cloned at the same root location. \
-``` git clone https://github.com/AMReX-Codes/amrex.git ``` \
-``` git clone https://AMReX-Microelectronics/MagneX.git ```
+``` >> git clone https://github.com/AMReX-Codes/amrex.git ``` \
+``` >> git clone https://AMReX-Microelectronics/MagneX.git ```
 ## Dependencies
 Beyond a standard Ubuntu22 installation, the Ubuntu packages libfftw3-dev, libfftw3-mpi-dev, and cmake are required.
 Also heFFTe is required.  At the same level that AMReX and MagneX are cloned, run: \
-``` git clone https://github.com/icl-utk-edu/heffte.git ```\
-``` cd heffte ```\
-``` mkdir build ```\
-``` cd build ```\
-``` cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_STANDARD=17 -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=. -DHeffte_ENABLE_FFTW=ON -DHeffte_ENABLE_CUDA=OFF .. ```\
-``` make -j4 ```\
-``` make install```
+``` >> git clone https://github.com/icl-utk-edu/heffte.git ```\
+``` >> cd heffte ```\
+``` >> mkdir build ```\
+``` >> cd build ```\
+``` >> cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_STANDARD=17 -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=. -DHeffte_ENABLE_FFTW=ON -DHeffte_ENABLE_CUDA=OFF .. ```\
+``` >> make -j4 ```\
+``` >> make install```
 ## Build
- Navigate to the Exec folder of MagneX and execute
-```make -j 4```
+ Navigate to MagneX/Exec/ and run:\
+```>> make -j4```
 
 # Running MagneX
 Example input scripts are located in `Exec/standard_problem_inputs/` directory. 
 ## Simple Testcase
-You can run the following to simulate muMAG Standard Problem 4 dynamics:
-## For pure MPI build (but with a single MPI rank)
-```./main3d.gnu.MPI.ex standard_problem_inputs/inputs_std4```
+You can run the following to simulate muMAG Standard Problem 4 dynamics:\
+```>> ./main3d.gnu.MPI.ex standard_problem_inputs/inputs_std4```
 # Visualization and Data Analysis
 Refer to the following link for several visualization tools that can be used for AMReX plotfiles. 
 
@@ -38,9 +37,9 @@ Refer to the following link for several visualization tools that can be used for
 You can extract the data in numpy array format using yt (you can refer to this for installation and usage of [yt](https://yt-project.org/). After you have installed yt, you can do something as follows, for example, to get variable 'Pz' (z-component of polarization)
 ```
 import yt
-ds = yt.load('./plt00001000/') # for data at time step 1000
+ds = yt.load('./plt00010000/') # for data at time step 10000
 ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
-P_array = ad0['Pz'].to_ndarray()
+Px_array = ad0['Mx'].to_ndarray()
 ```
 # Publications
 1. Z. Yao, P. Kumar, J. C. LePelch, and A. Nonaka, MagneX: An Exascale-Enabled Micromagnetics Solver for Spintronic Systems, in preparation.

From 09256d6b8929583e72ce027b659d854d451a7e7b Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 10:04:25 -0700
Subject: [PATCH 07/30] Update README.md

---
 README.md | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index bef5634..32e2939 100644
--- a/README.md
+++ b/README.md
@@ -10,8 +10,9 @@ Make sure that AMReX and MagneX are cloned at the same root location. \
 ``` >> git clone https://github.com/AMReX-Codes/amrex.git ``` \
 ``` >> git clone https://AMReX-Microelectronics/MagneX.git ```
 ## Dependencies
-Beyond a standard Ubuntu22 installation, the Ubuntu packages libfftw3-dev, libfftw3-mpi-dev, and cmake are required.
-Also heFFTe is required.  At the same level that AMReX and MagneX are cloned, run: \
+Beyond a standard Ubuntu22 installation, the Ubuntu packages libfftw3-dev, libfftw3-mpi-dev, and cmake are required.\
+SUNDIALS is optional an enabled Runge-Kutta, implicit, and multirate integrators (more detailed instructions in the full documentation).\
+heFFTe is required.  At the same level that AMReX and MagneX are cloned, run: \
 ``` >> git clone https://github.com/icl-utk-edu/heffte.git ```\
 ``` >> cd heffte ```\
 ``` >> mkdir build ```\

From fe8ecfdc3faa6b3886b49c0130ca1769c3d86594 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 22 May 2024 10:08:43 -0700
Subject: [PATCH 08/30] Update README.md

---
 README.md | 2 --
 1 file changed, 2 deletions(-)

diff --git a/README.md b/README.md
index 32e2939..0a89797 100644
--- a/README.md
+++ b/README.md
@@ -25,8 +25,6 @@ heFFTe is required.  At the same level that AMReX and MagneX are cloned, run: \
 ```>> make -j4```
 
 # Running MagneX
-Example input scripts are located in `Exec/standard_problem_inputs/` directory. 
-## Simple Testcase
 You can run the following to simulate muMAG Standard Problem 4 dynamics:\
 ```>> ./main3d.gnu.MPI.ex standard_problem_inputs/inputs_std4```
 # Visualization and Data Analysis

From 1df8bdbab9d6b21a4da39600df354e57eebca386 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Thu, 23 May 2024 11:06:49 -0700
Subject: [PATCH 09/30] whitespace

---
 README.md | 2 --
 1 file changed, 2 deletions(-)

diff --git a/README.md b/README.md
index 220f6e4..ea52767 100644
--- a/README.md
+++ b/README.md
@@ -2,11 +2,9 @@
 MagneX is a massively parallel, 3D micromagnetics solver for modeling magnetic materials.
 MagneX solves the Landau-Lifshitz-Gilbert (LLG) equations, including exchange, anisotropy, demagnetization, and Dzyaloshinskii-Moriya interaction (DMI) coupling.
 The algorithm is implemented using Exascale Computing Project software framework, AMReX, which provides effective scalability on manycore and GPU-based supercomputing architectures.
-
 # Documentation and Getting Help
 More extensive documentation is available [HERE](https://amrex-microelectronics.github.io).
 Our community is here to help. Please report installation problems or general questions about the code in the github [Issues](https://github.com/AMReX-Microelectronics/MagneX/issues) tab above.
-
 # Installation
 Here are instructions for a basic, pure-MPI (no GPU) installation.  More detailed instructions for GPU systems are in the full documentation.
 ## Download AMReX and MagneX Repositories

From 15a993ce8f8f9dd2b6be8e231f86a4e4b978daed Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 31 May 2024 11:57:44 -0700
Subject: [PATCH 10/30] enable DMI and anis in scaling test

---
 Exec/inputs_performance | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Exec/inputs_performance b/Exec/inputs_performance
index 3b291c2..3512d27 100644
--- a/Exec/inputs_performance
+++ b/Exec/inputs_performance
@@ -31,16 +31,16 @@ Ms_parser(x,y,z) =         "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.
 gamma_parser(x,y,z) =      "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -1.759e11"
 exchange_parser(x,y,z) =   "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 1.3e-11"
 anisotropy_parser(x,y,z) = "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -139.26"
-DMI_parser(x,y,z) = "0."
+DMI_parser(x,y,z) =        "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -4.5e-3"
 
 precession = 1
 demag_coupling = 1
 FFT_solver = 1
 M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
 exchange_coupling = 1
-anisotropy_coupling = 0
-anisotropy_axis = 0.0 1.0 0.0
-DMI_coupling = 0
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
 
 # INTEGRATION
 

From 535bee5a0a2bfe1cceff0a5ca5c6227be2c5a73c Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 12 Jun 2024 18:10:17 -0700
Subject: [PATCH 11/30] skyrmion diagnostic post processing utility

---
 Exec/DMI_diagnostic/DMI_diagnostic.cpp | 144 +++++++++++++++++++++++++
 Exec/DMI_diagnostic/GNUmakefile        |  33 ++++++
 Exec/DMI_diagnostic/Make.package       |  20 ++++
 3 files changed, 197 insertions(+)
 create mode 100644 Exec/DMI_diagnostic/DMI_diagnostic.cpp
 create mode 100644 Exec/DMI_diagnostic/GNUmakefile
 create mode 100644 Exec/DMI_diagnostic/Make.package

diff --git a/Exec/DMI_diagnostic/DMI_diagnostic.cpp b/Exec/DMI_diagnostic/DMI_diagnostic.cpp
new file mode 100644
index 0000000..70ffc58
--- /dev/null
+++ b/Exec/DMI_diagnostic/DMI_diagnostic.cpp
@@ -0,0 +1,144 @@
+
+#include <fstream>
+#include <iostream>
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_PlotFileUtil.H>
+
+using namespace std;
+using namespace amrex;
+
+static
+void
+PrintUsage (const char* progName)
+{
+    Print() << std::endl
+            << "Use this utility to extract a 1D profile along the x-axis of Mz and theta for a DMI skyrmion test." << std::endl << std::endl;
+    Print() << "Usage:" << '\n';
+    Print() << progName << " infile=inputFileName" << '\n' << '\n';
+
+    exit(1);
+}
+
+int
+main (int   argc,
+      char* argv[])
+{
+    amrex::Initialize(argc,argv);
+
+    if (argc == 1) {
+        PrintUsage(argv[0]);
+    }
+
+    // plotfile name
+    std::string iFile;
+
+    // read in parameters from inputs file
+    ParmParse pp;
+
+    // read in plotfile name
+    pp.query("infile", iFile);
+    if (iFile.empty())
+        amrex::Abort("You must specify `infile'");
+
+    // for the Header
+    std::string Header = iFile;
+    Header += "/Header";
+
+    // open header
+    ifstream x;
+    x.open(Header.c_str(), ios::in);
+
+    // read in first line of header
+    string str;
+    x >> str;
+
+    // read in number of components from header
+    int ncomp;
+    x >> ncomp;
+
+    // read in variable names from header
+    for (int n=0; n<ncomp; ++n) {
+        x >> str;
+    }
+
+    // read in dimensionality from header
+    int dim;
+    x >> dim;
+
+    if (dim != AMREX_SPACEDIM) {
+        Print() << "\nError: you are using a " << AMREX_SPACEDIM << "D build to open a "
+                << dim << "D plotfile\n\n";
+        Abort();
+    }
+
+    int lev = 0;
+
+    do {
+
+        if (lev > 9) {
+            Abort("Utility only works for 10 levels of refinement or less");
+        }
+
+        // storage for the MultiFab
+        MultiFab mf;
+
+        std::string iFile_lev = iFile;
+
+        std::string levX  = "/Level_"+to_string(lev)+"/Cell";
+        std::string levXX = "/Level_0"+to_string(lev)+"/Cell";
+
+        // now read in the plotfile data
+        // check to see whether the user pointed to the plotfile base directory
+        // or the data itself
+        if (amrex::FileExists(iFile+levX+"_H")) {
+            iFile_lev += levX;
+        } else if (amrex::FileExists(iFile+levXX+"_H")) {
+            iFile_lev += levXX;
+        } else {
+            break; // terminate while loop
+        }
+
+        // read in plotfile to MultiFab
+        VisMF::Read(mf, iFile_lev);
+
+        if (lev == 0) {
+            ncomp = mf.nComp();
+            Print() << "Number of components in the plotfile = " << ncomp << std::endl;
+            Print() << "Nodality of plotfile = " << mf.ixType().toIntVect() << std::endl;
+        }
+
+        // get boxArray to compute number of grid points at the level
+        BoxArray ba = mf.boxArray();
+        Print() << "Number of grid points at level " << lev << " = " << ba.numPts() << std::endl;
+
+        for ( MFIter mfi(mf,false); mfi.isValid(); ++mfi ) {
+
+            const Box& bx = mfi.validbox();
+            const auto lo = amrex::lbound(bx);
+            const auto hi = amrex::ubound(bx);
+
+            const Array4<Real>& mfdata = mf.array(mfi);
+
+            Real offset = 0.;
+            
+            int k = (hi.z+1)/2;
+            int j = (hi.y+1)/2;
+            for (auto i = (hi.x+1)/2; i <= hi.x; ++i) {
+                std::cout << i << " " << mfdata(i,j,k,3)/1.1e6 << " " << mfdata(i,j,k,16)+offset << "\n";
+                // 2pi cyclic fix
+                if (i<hi.x) {
+                    if (offset==0. && mfdata(i,j,k,16)>5. && mfdata(i+1,j,k,16)<1.) {
+                        offset = 2.*M_PI;
+                    }
+                }
+            }
+           
+        } // end MFIter
+
+        // proceed to next level of refinement
+        ++lev;
+
+    } while(true);
+
+}
diff --git a/Exec/DMI_diagnostic/GNUmakefile b/Exec/DMI_diagnostic/GNUmakefile
new file mode 100644
index 0000000..7ead0dd
--- /dev/null
+++ b/Exec/DMI_diagnostic/GNUmakefile
@@ -0,0 +1,33 @@
+AMREX_HOME ?= ../../../amrex
+
+DEBUG	= TRUE
+DIM	= 3
+
+COMP    = gcc
+
+PRECISION = DOUBLE
+
+USE_MPI   = FALSE
+USE_OMP   = FALSE
+
+###################################################
+
+EBASE = DMI_diagnostic
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+include ./Make.package
+include $(AMREX_HOME)/Src/Base/Make.package
+
+vpath %.c   : . $(vpathdir)
+vpath %.h   : . $(vpathdir)
+vpath %.cpp : . $(vpathdir)
+vpath %.H   : . $(vpathdir)
+vpath %.F   : . $(vpathdir)
+vpath %.f   : . $(vpathdir)
+vpath %.f90 : . $(vpathdir)
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules
+
+clean::
+	$(SILENT) $(RM) particle_compare.exe
diff --git a/Exec/DMI_diagnostic/Make.package b/Exec/DMI_diagnostic/Make.package
new file mode 100644
index 0000000..159417d
--- /dev/null
+++ b/Exec/DMI_diagnostic/Make.package
@@ -0,0 +1,20 @@
+CEXE_sources += ${EBASE}.cpp
+
+
+ifneq ($(EBASE), particle_compare)
+ifneq ($(EBASE), WritePlotfileToASCII)
+
+  INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Base
+  include $(AMREX_HOME)/Src/Base/Make.package
+  vpathdir += $(AMREX_HOME)/Src/Base
+
+  INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/amrdata
+  include $(AMREX_HOME)/Src/Extern/amrdata/Make.package
+  vpathdir += $(AMREX_HOME)/Src/Extern/amrdata
+
+  ifeq ($(NEEDS_f90_SRC),TRUE)
+    f90EXE_sources += ${EBASE}_nd.f90
+  endif
+
+endif
+endif

From a32bd3184978107a6a250cb763df301117277388 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Thu, 13 Jun 2024 13:47:04 -0700
Subject: [PATCH 12/30] renname file

---
 Exec/{input_2D_DMI => standard_problem_inputs/inputs_dmi1} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename Exec/{input_2D_DMI => standard_problem_inputs/inputs_dmi1} (100%)

diff --git a/Exec/input_2D_DMI b/Exec/standard_problem_inputs/inputs_dmi1
similarity index 100%
rename from Exec/input_2D_DMI
rename to Exec/standard_problem_inputs/inputs_dmi1

From c5f33493b3f0efcf8a788b74ec8de8ab8be040be Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Thu, 13 Jun 2024 17:37:05 -0700
Subject: [PATCH 13/30] restart bugfix

---
 Source/main.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Source/main.cpp b/Source/main.cpp
index d10a091..519c072 100644
--- a/Source/main.cpp
+++ b/Source/main.cpp
@@ -90,7 +90,7 @@ void main_main ()
     int start_step = 1;
     Real time = 0.0;
 
-    if (restart > 0) {
+    if (restart >= 0) {
 
         start_step = restart+1;
 

From 23f20f61243bbaef53da10ceeebd37beaf822579 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 15 Jul 2024 20:12:55 -0700
Subject: [PATCH 14/30] sundials 7.1.1

---
 Exec/README_sundials | 50 +++++++++++++++++++++++++++-----------------
 1 file changed, 31 insertions(+), 19 deletions(-)

diff --git a/Exec/README_sundials b/Exec/README_sundials
index 7d7356d..3fb088e 100644
--- a/Exec/README_sundials
+++ b/Exec/README_sundials
@@ -3,13 +3,19 @@ https://computing.llnl.gov/projects/sundials/faq#inst
 
 Installation
 
-# Download sundials-x.y.z.tar.gz file for SUNDIALS and extract it at the same level as amrex.
-# Update 6/11/24 - v7.0.0. tested
-https://computing.llnl.gov/projects/sundials/sundials-software
-
+# You need SUNDIALS v7.1.1 or later.
+# Check https://computing.llnl.gov/projects/sundials/sundials-software to see if it's available for download.
+# If so, download sundials-x.y.z.tar.gz and extract it at the same level as amrex using
 >> tar -xzvf sundials-x.y.z.tar.gz # where x.y.z is the version of sundials
+>> mv sundials-x.y.z sundials-src
+
+# If v7.1.1. is not available on the website, clone the git repo directly and use the latest version
+# At the same level that amrex is cloned, do:
+
+>> git clone https://github.com/LLNL/sundials.git
+>> mv sundials sundials-src
 
-# at the same level that amrex is cloned, do:
+# Next
 
 >> mkdir sundials
 >> cd sundials
@@ -22,7 +28,7 @@ HOST BUILD
 >> mkdir builddir
 >> cd builddir
 
->> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir/examples -DENABLE_MPI=ON ../../sundials-x.y.z
+>> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir/examples -DENABLE_MPI=ON ../../sundials-src
 
 ######################
 NVIDIA/CUDA BUILD
@@ -34,7 +40,7 @@ NVIDIA/CUDA BUILD
 >> mkdir builddir_cuda
 >> cd builddir_cuda
 
->> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir_cuda -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir_cuda/examples -DENABLE_CUDA=ON -DENABLE_MPI=ON ../../sundials-x.y.z
+>> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir_cuda -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir_cuda/examples -DENABLE_CUDA=ON -DENABLE_MPI=ON ../../sundials-src
 
 ######################
 
@@ -59,16 +65,14 @@ NVIDIA/CUDA BUILD
 
 TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
-## *** Selecting the integrator backend ***
-## integration.type can take on the following string or int values:
-## (without the quotation marks)
-## "ForwardEuler" or "0" = Native Forward Euler Integrator
-## "RungeKutta" or "1"   = Native Explicit Runge Kutta
-## "SUNDIALS" or "2"     = SUNDIALS ARKODE Integrator
-## for example:
-integration.type = <pick an option>
+# INTEGRATION
+## integration.type can take on the following values:
+## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator
+## 1 or "RungeKutta"   => Native AMReX Explicit Runge Kutta controlled by integration.rk.type
+## 2 or "SUNDIALS"     => SUNDIALS backend controlled by integration.sundials.strategy
+integration.type = SUNDIALS
 
-## *** Parameters Needed For Native Explicit Runge-Kutta ***
+## Native AMReX Explicit Runge-Kutta parameters
 #
 ## integration.rk.type can take the following values:
 ### 0 = User-specified Butcher Tableau
@@ -76,7 +80,15 @@ integration.type = <pick an option>
 ### 2 = Trapezoid Method
 ### 3 = SSPRK3 Method
 ### 4 = RK4 Method
-integration.rk.type = <pick an option>
+integration.rk.type = 1
+
+# Set the SUNDIALS method type:
+# ERK      = Explicit Runge-Kutta method
+# DIRK     = Diagonally Implicit Runge-Kutta method
+#
+# Optionally select a specific SUNDIALS method by name, see the SUNDIALS
+# documentation for the supported method names
 
-integration.sundials.strategy = ERK
-#integration.sundials.erk.method = SSPRK3
+# Use forward Euler (fixed step sizes only)
+integration.sundials.type = ERK
+integration.sundials.method = ARKODE_FORWARD_EULER_1_1

From 437c39b0e8d508f88e1c25268008b6267b035c54 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 2 Aug 2024 16:25:22 -0700
Subject: [PATCH 15/30] inputs file for paper demag plot

---
 Exec/inputs_demagviz | 115 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 115 insertions(+)
 create mode 100644 Exec/inputs_demagviz

diff --git a/Exec/inputs_demagviz b/Exec/inputs_demagviz
new file mode 100644
index 0000000..28b6b22
--- /dev/null
+++ b/Exec/inputs_demagviz
@@ -0,0 +1,115 @@
+n_cell = 64 64 64
+max_grid_size_x = 64
+max_grid_size_y = 64
+max_grid_size_z = 64
+
+dt = 1.e-15
+nsteps = 1
+plot_int = 1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 256.e-9 256.e-9 256.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>64.e-9)*(x<192.e-9)*(y>64.e-9)*(y<192.e-9)*(z>64.e-9)*(z<192.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>64.e-9)*(x<192.e-9)*(y>64.e-9)*(y<192.e-9)*(z>64.e-9)*(z<192.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>64.e-9)*(x<192.e-9)*(y>64.e-9)*(y<192.e-9)*(z>64.e-9)*(z<192.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>64.e-9)*(x<192.e-9)*(y>64.e-9)*(y<192.e-9)*(z>64.e-9)*(z<192.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>64.e-9)*(x<192.e-9)*(y>64.e-9)*(y<192.e-9)*(z>64.e-9)*(z<192.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>64.e-9)*(x<192.e-9)*(y>64.e-9)*(y<192.e-9)*(z>64.e-9)*(z<192.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>64.e-9)*(x<192.e-9)*(y>64.e-9)*(y<192.e-9)*(z>64.e-9)*(z<192.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 2 #Predictor-corrector
+#TimeIntegratorOption = 3 #2nd order artemis way
+#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+
+# tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
+iterative_tolerance = 1.e-9
+
+## amrex/sundials backend integrators
+## *** Selecting the integrator backend ***
+## integration.type can take on the following string or int values:
+## (without the quotation marks)
+## "ForwardEuler" or "0" = Native Forward Euler Integrator
+## "RungeKutta" or "1"   = Native Explicit Runge Kutta controlled by integration.rk.type
+## "SUNDIALS" or "2"     = SUNDIALS Integrators controlled by integration.sundials.strategy
+integration.type = RungeKutta
+
+## *** Parameters Needed For Native Explicit Runge-Kutta ***
+#
+## integration.rk.type can take the following values:
+### 0 = User-specified Butcher Tableau
+### 1 = Forward Euler
+### 2 = Trapezoid Method
+### 3 = SSPRK3 Method
+### 4 = RK4 Method
+integration.rk.type = 3
+
+## If using a user-specified Butcher Tableau, then
+## set nodes, weights, and table entries here:
+#
+## The Butcher Tableau is read as a flattened,
+## lower triangular matrix (but including the diagonal)
+## in row major format.
+integration.rk.weights = 1
+integration.rk.nodes = 0
+integration.rk.tableau = 0.0
+
+## *** Parameters Needed For SUNDIALS Integrators ***
+## integration.sundials.strategy specifies which ARKODE strategy to use.
+## The available options are (without the quotations):
+## "ERK" = Explicit Runge Kutta
+## "MRI" = Multirate Integrator
+## "MRITEST" = Tests the Multirate Integrator by setting a zero-valued fast RHS function
+## for example:
+integration.sundials.strategy = ERK
+
+## *** Parameters Specific to SUNDIALS ERK Strategy ***
+## (Requires integration.type=SUNDIALS and integration.sundials.strategy=ERK)
+## integration.sundials.erk.method specifies which explicit Runge Kutta method
+## for SUNDIALS to use. The following options are supported:
+## "SSPRK3" = 3rd order strong stability preserving RK (default)
+## "Trapezoid" = 2nd order trapezoidal rule
+## "ForwardEuler" = 1st order forward euler
+## for example:
+integration.sundials.erk.method = SSPRK3
+
+## *** Parameters Specific to SUNDIALS MRI Strategy ***
+## (Requires integration.type=SUNDIALS and integration.sundials.strategy=MRI)
+## integration.sundials.mri.implicit_inner specifies whether or not to use an implicit inner solve
+## integration.sundials.mri.outer_method specifies which outer (slow) method to use
+## integration.sundials.mri.inner_method specifies which inner (fast) method to use
+## The following options are supported for both the inner and outer methods:
+## "KnothWolke3" = 3rd order Knoth-Wolke method (default for outer method)
+## "Trapezoid" = 2nd order trapezoidal rule
+## "ForwardEuler" = 1st order forward euler (default for inner method)
+## for example:
+integration.sundials.mri.implicit_inner = false
+integration.sundials.mri.outer_method = KnothWolke3
+integration.sundials.mri.inner_method = Trapezoid

From 2fe1a87e413328bb9c1f66ec1080a40ffe566c2a Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 19 Aug 2024 16:44:21 -0700
Subject: [PATCH 16/30] license and copyright

---
 Legal.txt   | 17 +++++++++++++++++
 license.txt | 43 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 60 insertions(+)
 create mode 100644 Legal.txt
 create mode 100644 license.txt

diff --git a/Legal.txt b/Legal.txt
new file mode 100644
index 0000000..3fe4c70
--- /dev/null
+++ b/Legal.txt
@@ -0,0 +1,17 @@
+*** Copyright Notice ***
+
+Exascale-Enabled Models and Algorithms for Microelectronics Applications
+(MicroEleX) Copyright (c) 2024, The Regents of the University of California,
+through Lawrence Berkeley National Laboratory (subject to receipt of any
+required approvals from the U.S. Dept. of Energy). All rights reserved.
+
+If you have questions about your rights to use or distribute this software,
+please contact Berkeley Lab's Intellectual Property Office at
+IPO@lbl.gov.
+
+NOTICE.  This Software was developed under funding from the U.S. Department
+of Energy and the U.S. Government consequently retains certain rights.  As
+such, the U.S. Government has been granted for itself and others acting on
+its behalf a paid-up, nonexclusive, irrevocable, worldwide license in the
+Software to reproduce, distribute copies to the public, prepare derivative 
+works, and perform publicly and display publicly, and to permit others to do so.
diff --git a/license.txt b/license.txt
new file mode 100644
index 0000000..e3eed5c
--- /dev/null
+++ b/license.txt
@@ -0,0 +1,43 @@
+*** License Agreement ***
+
+Exascale-Enabled Models and Algorithms for Microelectronics Applications
+(MicroEleX) Copyright (c) 2024, The Regents of the University of California,
+through Lawrence Berkeley National Laboratory (subject to receipt of any
+required approvals from the U.S. Dept. of Energy). All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+(1) Redistributions of source code must retain the above copyright notice,
+this list of conditions and the following disclaimer.
+
+(2) Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+(3) Neither the name of the University of California, Lawrence Berkeley
+National Laboratory, U.S. Dept. of Energy nor the names of its contributors
+may be used to endorse or promote products derived from this software
+without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+You are under no obligation whatsoever to provide any bug fixes, patches,
+or upgrades to the features, functionality or performance of the source
+code ("Enhancements") to anyone; however, if you choose to make your
+Enhancements available either publicly, or directly to Lawrence Berkeley
+National Laboratory, without imposing a separate written license agreement
+for such Enhancements, then you hereby grant the following license: a
+non-exclusive, royalty-free perpetual license to install, use, modify,
+prepare derivative works, incorporate into other computer software,
+distribute, and sublicense such enhancements or derivative works thereof,
+in binary and source code form.

From 24a4e0988b0caa7fbdcbf6ab3cb608b3076d636e Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 30 Sep 2024 08:24:33 -0700
Subject: [PATCH 17/30] patches to get regression tests running again some
 benchmarks will be changing

---
 Exec/regression_inputs/inputs_PSSW              | 2 +-
 Exec/regression_inputs/inputs_compare_subdomain | 2 +-
 Source/main.cpp                                 | 4 ++--
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Exec/regression_inputs/inputs_PSSW b/Exec/regression_inputs/inputs_PSSW
index 8141fff..0f1bcd7 100644
--- a/Exec/regression_inputs/inputs_PSSW
+++ b/Exec/regression_inputs/inputs_PSSW
@@ -32,7 +32,7 @@ anisotropy_parser(x,y,z) = "(x>-2587.e-9)*(x<2587.e-9)*(y>-2587.e-9)*(y<2587.e-9
 DMI_parser(x,y,z) = "0.0"
 
 precession = 1
-demag_coupling = 1
+demag_coupling = 0
 FFT_solver = 1
 M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
 exchange_coupling = 1
diff --git a/Exec/regression_inputs/inputs_compare_subdomain b/Exec/regression_inputs/inputs_compare_subdomain
index 0e826d0..527b929 100644
--- a/Exec/regression_inputs/inputs_compare_subdomain
+++ b/Exec/regression_inputs/inputs_compare_subdomain
@@ -32,7 +32,7 @@ anisotropy_parser(x,y,z) = "(x>-1.5625e-8)*(x<1.5625e-8)*(y>-15.625e-9)*(y<15.62
 DMI_parser(x,y,z) = "0.0"
 
 precession = 1
-demag_coupling = 1
+demag_coupling = 0
 FFT_solver = 1
 M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
 exchange_coupling = 1
diff --git a/Source/main.cpp b/Source/main.cpp
index 0e9162b..c9c429a 100644
--- a/Source/main.cpp
+++ b/Source/main.cpp
@@ -211,11 +211,11 @@ void main_main ()
     if (demag_coupling == 1) {
 
         const Real* dx = geom.CellSize();
-        if (dx[0] != dx[1]) {
+        if (!almostEqual(dx[0], dx[1], 10)) {
             Abort("Demag requires dx=dy");
         }
 #if (AMREX_SPACEDIM==3)
-        if (dx[0] != dx[2]) {
+        if (!almostEqual(dx[0], dx[2], 10)) {
             Abort("Demag requires dx=dy=dz");
         }
 #endif

From afe92a244f3c91d64da29bf16f34b5f3c8ba7dd0 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 2 Oct 2024 11:59:32 -0700
Subject: [PATCH 18/30] fix sundials settings

---
 Exec/inputs_performance | 79 ++++++++++++++++++-----------------------
 1 file changed, 35 insertions(+), 44 deletions(-)

diff --git a/Exec/inputs_performance b/Exec/inputs_performance
index 3512d27..9381c61 100644
--- a/Exec/inputs_performance
+++ b/Exec/inputs_performance
@@ -43,7 +43,6 @@ anisotropy_axis = 0.0 0.0 1.0
 DMI_coupling = 1
 
 # INTEGRATION
-
 TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
@@ -52,24 +51,19 @@ TimeIntegratorOption = 1 #Forward Euler
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
 
-## amrex/sundials backend integrators
-## *** Selecting the integrator backend ***
-## integration.type can take on the following string or int values:
-## (without the quotation marks)
-## "ForwardEuler" or "0" = Native Forward Euler Integrator
-## "RungeKutta" or "1"   = Native Explicit Runge Kutta controlled by integration.rk.type
-## "SUNDIALS" or "2"     = SUNDIALS Integrators controlled by integration.sundials.strategy
+# for TimeIntegratorOption = 4, integration.type can take on the following values:
+## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator
+## 1 or "RungeKutta"   => Native AMReX Explicit Runge Kutta controlled by integration.rk.type
+## 2 or "SUNDIALS"     => SUNDIALS backend controlled by integration.sundials.type
 integration.type = RungeKutta
 
-## *** Parameters Needed For Native Explicit Runge-Kutta ***
-#
-## integration.rk.type can take the following values:
+## for integration.type = RungeKutta, integration.rk.type can take the following values:
 ### 0 = User-specified Butcher Tableau
 ### 1 = Forward Euler
 ### 2 = Trapezoid Method
 ### 3 = SSPRK3 Method
 ### 4 = RK4 Method
-integration.rk.type = 3
+integration.rk.type = 4
 
 ## If using a user-specified Butcher Tableau, then
 ## set nodes, weights, and table entries here:
@@ -81,35 +75,32 @@ integration.rk.weights = 1
 integration.rk.nodes = 0
 integration.rk.tableau = 0.0
 
-## *** Parameters Needed For SUNDIALS Integrators ***
-## integration.sundials.strategy specifies which ARKODE strategy to use.
-## The available options are (without the quotations):
-## "ERK" = Explicit Runge Kutta
-## "MRI" = Multirate Integrator
-## "MRITEST" = Tests the Multirate Integrator by setting a zero-valued fast RHS function
-## for example:
-integration.sundials.strategy = ERK
-
-## *** Parameters Specific to SUNDIALS ERK Strategy ***
-## (Requires integration.type=SUNDIALS and integration.sundials.strategy=ERK)
-## integration.sundials.erk.method specifies which explicit Runge Kutta method
-## for SUNDIALS to use. The following options are supported:
-## "SSPRK3" = 3rd order strong stability preserving RK (default)
-## "Trapezoid" = 2nd order trapezoidal rule
-## "ForwardEuler" = 1st order forward euler
-## for example:
-integration.sundials.erk.method = SSPRK3
-
-## *** Parameters Specific to SUNDIALS MRI Strategy ***
-## (Requires integration.type=SUNDIALS and integration.sundials.strategy=MRI)
-## integration.sundials.mri.implicit_inner specifies whether or not to use an implicit inner solve
-## integration.sundials.mri.outer_method specifies which outer (slow) method to use
-## integration.sundials.mri.inner_method specifies which inner (fast) method to use
-## The following options are supported for both the inner and outer methods:
-## "KnothWolke3" = 3rd order Knoth-Wolke method (default for outer method)
-## "Trapezoid" = 2nd order trapezoidal rule
-## "ForwardEuler" = 1st order forward euler (default for inner method)
-## for example:
-integration.sundials.mri.implicit_inner = false
-integration.sundials.mri.outer_method = KnothWolke3
-integration.sundials.mri.inner_method = Trapezoid
+## for integration.type = SUNDIALS, set the SUNDIALS method type:
+# ERK      = Explicit Runge-Kutta method
+# DIRK     = Diagonally Implicit Runge-Kutta method
+# IMEX-RK  = Implicit-Explicit Additive Runge-Kutta method
+# EX-MRI   = Explicit Multirate Infatesimal method
+# IM-MRI   = Implicit Multirate Infatesimal method
+# IMEX-MRI = Implicit-Explicit Multirate Infatesimal method
+integration.sundials.type = EX-MRI
+
+# *** Select a specific SUNDIALS ERK or DIRK method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#explicit-butcher-tables
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#implicit-butcher-tables
+#integration.sundials.method = ARKODE_FORWARD_EULER_1_1
+#integration.sundials.method = ARKODE_BACKWARD_EULER_1_1
+
+# *** Select a specific SUNDIALS IMEX-RK method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#additive-butcher-tables
+#integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2
+#integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2
+
+# *** Select a specific SUNDIALS EX-MRI, IM-MRI, or IMEX-MRI method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Usage/MRIStep/MRIStepCoupling.html#mri-coupling-tables
+#  *** Select a specific SUNDIALS MRI slow method coupling table ***
+integration.sundials.method = ARKODE_MIS_KW3
+#  *** Select a specific SUNDIALS ERK or DIRK fast method ***
+integration.sundials.fast_type = ERK  # ERK or DIRK
+integration.sundials.fast_method = ARKODE_KNOTH_WOLKE_3_3
+
+fast_exchange = 1

From 52cf85cc2874758e19f61b0aab604bf5bedf087f Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 1 Nov 2024 10:09:28 -0700
Subject: [PATCH 19/30] cleanup per Weiqun

---
 Exec/GNUmakefile           | 12 +-----------
 Source/Demagnetization.cpp |  2 +-
 2 files changed, 2 insertions(+), 12 deletions(-)

diff --git a/Exec/GNUmakefile b/Exec/GNUmakefile
index 8cbbc1e..6c9ceb2 100644
--- a/Exec/GNUmakefile
+++ b/Exec/GNUmakefile
@@ -8,6 +8,7 @@ USE_CUDA     = FALSE
 USE_HIP      = FALSE
 COMP         = gnu
 DIM          = 3
+USE_FFT      = TRUE
 
 USE_SUNDIALS = FALSE
 
@@ -47,17 +48,6 @@ ifeq ($(USE_SUNDIALS),TRUE)
 
 endif
 
-ifeq ($(USE_CUDA),TRUE)
-  libraries += -lcufft
-else ifeq ($(USE_HIP),TRUE)
-  # Use rocFFT.  ROC_PATH is defined in amrex
-  INCLUDE_LOCATIONS += $(ROC_PATH)/rocfft/include
-  LIBRARY_LOCATIONS += $(ROC_PATH)/rocfft/lib
-  LIBRARIES += -L$(ROC_PATH)/rocfft/lib -lrocfft
-else
-  libraries += -lfftw3_mpi -lfftw3f -lfftw3
-endif
-
 include $(AMREX_HOME)/Src/Base/Make.package
 include $(AMREX_HOME)/Src/FFT/Make.package
 include $(AMREX_HOME)/Src/Boundary/Make.package 
diff --git a/Source/Demagnetization.cpp b/Source/Demagnetization.cpp
index c495a34..d24a6ff 100644
--- a/Source/Demagnetization.cpp
+++ b/Source/Demagnetization.cpp
@@ -60,7 +60,7 @@ void Demagnetization::define()
     Geometry cgeom_large(cdomain_large, real_box_large, CoordSys::cartesian, is_periodic);
     auto cba_large = amrex::decompose(cdomain_large, ParallelContext::NProcsSub(),
                                       {AMREX_D_DECL(true,true,false)});
-    DistributionMapping cdm_large = amrex::FFT::detail::make_iota_distromap(cba_large.size());
+    DistributionMapping cdm_large(cba_large);
     
     Kxx_fft.define(cba_large, cdm_large, 1, 0);
     Kxy_fft.define(cba_large, cdm_large, 1, 0);

From 97d152a61fe686babbe6f0c4bf45bf4324d0332d Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 10 Mar 2025 21:25:49 -0700
Subject: [PATCH 20/30] logic for different physics wasn't right for
 single-rate case

---
 Source/main.cpp | 48 ++++++++++++++++++++++++++++--------------------
 1 file changed, 28 insertions(+), 20 deletions(-)

diff --git a/Source/main.cpp b/Source/main.cpp
index 7ad879b..4893c33 100644
--- a/Source/main.cpp
+++ b/Source/main.cpp
@@ -558,39 +558,47 @@ void main_main ()
                 }
 
                 // exchange
-                if ( (using_MRI==1 && fast_exchange) || (using_IMEX==1 && implicit_exchange) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_exchangefield[idim].setVal(0.);
+                if (exchange_coupling == 1) {
+                    if ( (using_MRI==1 && fast_exchange) || (using_IMEX==1 && implicit_exchange) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_exchangefield[idim].setVal(0.);
+                        }
+                    } else {
+                        CalculateH_exchange(ar_state, H_exchangefield, Ms, exchange, DMI, geom);
                     }
-                } else {
-                    CalculateH_exchange(ar_state, H_exchangefield, Ms, exchange, DMI, geom);
                 }
 
                 // DMI
-                if ( (using_MRI==1 && fast_DMI) || (using_IMEX==1 && implicit_DMI) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_DMIfield[idim].setVal(0.);
+                if (DMI_coupling == 1) {
+                    if ( (using_MRI==1 && fast_DMI) || (using_IMEX==1 && implicit_DMI) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_DMIfield[idim].setVal(0.);
+                        }
+                    } else {
+                        CalculateH_DMI(ar_state, H_DMIfield, Ms, exchange, DMI, geom);
                     }
-                } else {
-                    CalculateH_DMI(ar_state, H_DMIfield, Ms, exchange, DMI, geom);
                 }
 
                 // anisotropy
-                if ( (using_MRI==1 && fast_anisotropy) || (using_IMEX==1 && implicit_anisotropy) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_anisotropyfield[idim].setVal(.0);
+                if (anisotropy_coupling == 1) {
+                    if ( (using_MRI==1 && fast_anisotropy) || (using_IMEX==1 && implicit_anisotropy) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_anisotropyfield[idim].setVal(.0);
+                        }
+                    } else {
+                        CalculateH_anisotropy(ar_state, H_anisotropyfield, Ms, anisotropy);
                     }
-                } else {
-                    CalculateH_anisotropy(ar_state, H_anisotropyfield, Ms, anisotropy);
                 }
 
                 // H_demag
-                if ( (using_MRI==1 && fast_demag) || (using_IMEX==1 && implicit_demag) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_demagfield[idim].setVal(0.);
+                if (demag_coupling == 1) {
+                    if ( (using_MRI==1 && fast_demag) || (using_IMEX==1 && implicit_demag) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_demagfield[idim].setVal(0.);
+                        }
+                    } else {
+                        demag_solver.CalculateH_demag(ar_state, H_demagfield);
                     }
-                } else {
-                    demag_solver.CalculateH_demag(ar_state, H_demagfield);
                 }
 
                 // Compute f^n = f(M^n, H^n)

From 7ec1149367e2b0526d49c4bccd6bf11a975cca32 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Tue, 29 Jul 2025 13:56:49 -0700
Subject: [PATCH 21/30] add SUNDIALS_HOME

---
 Exec/GNUmakefile | 31 +------------------------------
 1 file changed, 1 insertion(+), 30 deletions(-)

diff --git a/Exec/GNUmakefile b/Exec/GNUmakefile
index 6c9ceb2..0379b6f 100644
--- a/Exec/GNUmakefile
+++ b/Exec/GNUmakefile
@@ -11,6 +11,7 @@ DIM          = 3
 USE_FFT      = TRUE
 
 USE_SUNDIALS = FALSE
+SUNDIALS_HOME ?= ../../sundials/instdir
 
 include $(AMREX_HOME)/Tools/GNUMake/Make.defs
 
@@ -18,36 +19,6 @@ include ../Source/Make.package
 VPATH_LOCATIONS  += ../Source
 INCLUDE_LOCATIONS += ../Source
 
-ifeq ($(USE_SUNDIALS),TRUE)
-  ifeq ($(USE_CUDA),TRUE)
-    SUNDIALS_ROOT ?= $(TOP)../../sundials/instdir_cuda
-  else
-    SUNDIALS_ROOT ?= $(TOP)../../sundials/instdir
-  endif
-  ifeq ($(NERSC_HOST),perlmutter)
-    SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib64
-  else
-    SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib
-  endif
-
-  USE_CVODE_LIBS ?= TRUE
-  USE_ARKODE_LIBS ?= TRUE
-
-  DEFINES += -DAMREX_USE_SUNDIALS
-  INCLUDE_LOCATIONS += $(SUNDIALS_ROOT)/include
-  LIBRARY_LOCATIONS += $(SUNDIALS_LIB_DIR)
-
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_cvode
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_arkode
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_nvecmanyvector
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_core
-
-  ifeq ($(USE_CUDA),TRUE)
-    LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_nveccuda
-  endif
-
-endif
-
 include $(AMREX_HOME)/Src/Base/Make.package
 include $(AMREX_HOME)/Src/FFT/Make.package
 include $(AMREX_HOME)/Src/Boundary/Make.package 

From b95ceba2359aa2d9f47bb741cf3b7b698738964b Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 30 Jul 2025 09:56:08 -0700
Subject: [PATCH 22/30] sundials update

---
 Exec/README_sundials | 96 ++------------------------------------------
 1 file changed, 3 insertions(+), 93 deletions(-)

diff --git a/Exec/README_sundials b/Exec/README_sundials
index 3fb088e..3b6830a 100644
--- a/Exec/README_sundials
+++ b/Exec/README_sundials
@@ -1,94 +1,4 @@
-SUNDIALS installation guide:
-https://computing.llnl.gov/projects/sundials/faq#inst
+Refer to https://amrex-codes.github.io/amrex/docs_html/TimeIntegration_Chapter.html
+for SNUDIALS installation, build, and usage instructions.
 
-Installation
-
-# You need SUNDIALS v7.1.1 or later.
-# Check https://computing.llnl.gov/projects/sundials/sundials-software to see if it's available for download.
-# If so, download sundials-x.y.z.tar.gz and extract it at the same level as amrex using
->> tar -xzvf sundials-x.y.z.tar.gz # where x.y.z is the version of sundials
->> mv sundials-x.y.z sundials-src
-
-# If v7.1.1. is not available on the website, clone the git repo directly and use the latest version
-# At the same level that amrex is cloned, do:
-
->> git clone https://github.com/LLNL/sundials.git
->> mv sundials sundials-src
-
-# Next
-
->> mkdir sundials
->> cd sundials
-
-######################
-HOST BUILD 
-######################
-
->> mkdir instdir
->> mkdir builddir
->> cd builddir
-
->> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir/examples -DENABLE_MPI=ON ../../sundials-src
-
-######################
-NVIDIA/CUDA BUILD
-######################
-
-# Navigate back to the 'sundials' directory and do:
-
->> mkdir instdir_cuda
->> mkdir builddir_cuda
->> cd builddir_cuda
-
->> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir_cuda -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir_cuda/examples -DENABLE_CUDA=ON -DENABLE_MPI=ON ../../sundials-src
-
-######################
-
->> make -j4
->> make install
-
-# in your .bashrc or preferred configuration file, add the following (and then "source ~/.bashrc")
-
-# If you have a CPU build:
-
->> export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/pathto/sundials/instdir/lib/
-
-# If you have a NVIDIA/CUDA build:
-
->> export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/pathto/sundials/instdir_cuda/lib/
-
-# now you are ready to compile MagneX with:
-
->> make -j4 USE_SUNDIALS=TRUE # optional to have 'USE_CUDA=TRUE' as well
-
-# in your inputs file, you will need to have:
-
-TimeIntegratorOption = 4 #amrex/sundials backend integrators
-
-# INTEGRATION
-## integration.type can take on the following values:
-## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator
-## 1 or "RungeKutta"   => Native AMReX Explicit Runge Kutta controlled by integration.rk.type
-## 2 or "SUNDIALS"     => SUNDIALS backend controlled by integration.sundials.strategy
-integration.type = SUNDIALS
-
-## Native AMReX Explicit Runge-Kutta parameters
-#
-## integration.rk.type can take the following values:
-### 0 = User-specified Butcher Tableau
-### 1 = Forward Euler
-### 2 = Trapezoid Method
-### 3 = SSPRK3 Method
-### 4 = RK4 Method
-integration.rk.type = 1
-
-# Set the SUNDIALS method type:
-# ERK      = Explicit Runge-Kutta method
-# DIRK     = Diagonally Implicit Runge-Kutta method
-#
-# Optionally select a specific SUNDIALS method by name, see the SUNDIALS
-# documentation for the supported method names
-
-# Use forward Euler (fixed step sizes only)
-integration.sundials.type = ERK
-integration.sundials.method = ARKODE_FORWARD_EULER_1_1
+Make sure that SUNDIALS_HOME in the GNUmakefile points to the installation directory.

From a689f78c3ee836dfc38f1ee85b31a076019de56b Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 2 Jan 2026 08:30:31 -0800
Subject: [PATCH 23/30] comments in inputs file

---
 Exec/standard_problem_inputs/inputs_std4_IMEX | 120 ++++++++++++++++++
 Exec/standard_problem_inputs/inputs_std4_MRI  |   2 +-
 Exec/standard_problem_inputs/inputs_std4_RK4  |   4 +-
 3 files changed, 123 insertions(+), 3 deletions(-)
 create mode 100644 Exec/standard_problem_inputs/inputs_std4_IMEX

diff --git a/Exec/standard_problem_inputs/inputs_std4_IMEX b/Exec/standard_problem_inputs/inputs_std4_IMEX
new file mode 100644
index 0000000..c0cbe9f
--- /dev/null
+++ b/Exec/standard_problem_inputs/inputs_std4_IMEX
@@ -0,0 +1,120 @@
+amrex.use_gpu_aware_mpi=1
+
+##############
+# 0.78125nm case
+##############
+n_cell = 640 160 4
+max_grid_size_x = 640
+max_grid_size_y = 160
+max_grid_size_z = 4
+dt = 5.e-15 # IMEX, 1.e-14 unstable
+
+stop_time = 4.e-9
+
+plot_int = 25000
+chk_int = -1
+restart = -1
+diag_type = 4
+
+prob_lo = 0. 0. 0.
+prob_hi = 500.e-9 125.e-9 3.125e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 1
+Hx_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 - (t>1.e-9)*19576"
+Hy_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 + (t>1.e-9)*3422"
+Hz_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16"
+
+# Field 2:  mu_0 Hx=-35.5 mT, mu_0 Hy=-6.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 36 mT, directed 190 degrees counterclockwise from the positive x axis
+#timedependent_Hbias = 1
+#Hx_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 - (t>1.e-9)*28259"
+#Hy_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 - (t>1.e-9)*5013"
+#Hz_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16"
+
+timedependent_alpha = 1
+alpha_parser(x,y,z,t) =  " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9) * (0.5*(t<=1.e-9) + 0.02*(t>1.e-9))"
+Ms_parser(x,y,z) =       " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9) * 8.e5"
+gamma_parser(x,y,z) =    " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9) * -1.759e11"
+exchange_parser(x,y,z) = " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9)* 1.3e-11"
+anisotropy_parser(x,y,z) = "0."
+DMI_parser(x,y,z) = "0."
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 0
+anisotropy_axis = 0.0 1.0 0.0
+DMI_coupling = 0
+
+# INTEGRATION
+#TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 2 #Predictor-corrector
+#TimeIntegratorOption = 3 #2nd order artemis way
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
+
+# tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
+iterative_tolerance = 1.e-9
+
+# for TimeIntegratorOption = 4, integration.type can take on the following values:
+## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator
+## 1 or "RungeKutta"   => Native AMReX Explicit Runge Kutta controlled by integration.rk.type
+## 2 or "SUNDIALS"     => SUNDIALS backend controlled by integration.sundials.type
+integration.type = SUNDIALS
+
+## for integration.type = RungeKutta, integration.rk.type can take the following values:
+### 0 = User-specified Butcher Tableau
+### 1 = Forward Euler
+### 2 = Trapezoid Method
+### 3 = SSPRK3 Method
+### 4 = RK4 Method
+integration.rk.type = 4
+
+## If using a user-specified Butcher Tableau, then
+## set nodes, weights, and table entries here:
+#
+## The Butcher Tableau is read as a flattened,
+## lower triangular matrix (but including the diagonal)
+## in row major format.
+integration.rk.weights = 1
+integration.rk.nodes = 0
+integration.rk.tableau = 0.0
+
+## for integration.type = SUNDIALS, set the SUNDIALS method type:
+# ERK      = Explicit Runge-Kutta method
+# DIRK     = Diagonally Implicit Runge-Kutta method
+# IMEX-RK  = Implicit-Explicit Additive Runge-Kutta method
+# EX-MRI   = Explicit Multirate Infatesimal method
+# IM-MRI   = Implicit Multirate Infatesimal method
+# IMEX-MRI = Implicit-Explicit Multirate Infatesimal method
+integration.sundials.type = IMEX-RK
+
+# *** Select a specific SUNDIALS ERK or DIRK method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#explicit-butcher-tables
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#implicit-butcher-tables
+#integration.sundials.method = ARKODE_FORWARD_EULER_1_1
+#integration.sundials.method = ARKODE_BACKWARD_EULER_1_1
+
+# *** Select a specific SUNDIALS IMEX-RK method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#additive-butcher-tables
+integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2
+integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2
+
+# *** Select a specific SUNDIALS EX-MRI, IM-MRI, or IMEX-MRI method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Usage/MRIStep/MRIStepCoupling.html#mri-coupling-tables
+#  *** Select a specific SUNDIALS MRI slow method coupling table ***
+#integration.sundials.method = ARKODE_MIS_KW3
+#  *** Select a specific SUNDIALS ERK or DIRK fast method ***
+#integration.sundials.fast_type = ERK  # ERK or DIRK
+#integration.sundials.fast_method = ARKODE_KNOTH_WOLKE_3_3
+
+implicit_exchange = 1
diff --git a/Exec/standard_problem_inputs/inputs_std4_MRI b/Exec/standard_problem_inputs/inputs_std4_MRI
index 9b15ce3..b9770da 100644
--- a/Exec/standard_problem_inputs/inputs_std4_MRI
+++ b/Exec/standard_problem_inputs/inputs_std4_MRI
@@ -7,7 +7,7 @@ n_cell = 640 160 4
 max_grid_size_x = 640
 max_grid_size_y = 160
 max_grid_size_z = 4
-dt = 1.25e-13 # SUNDIALS MRI
+dt = 1.25e-13 # SUNDIALS MRI, 2.5e-13 unstable
 
 stop_time = 4.e-9
 
diff --git a/Exec/standard_problem_inputs/inputs_std4_RK4 b/Exec/standard_problem_inputs/inputs_std4_RK4
index becd5c2..13b08db 100644
--- a/Exec/standard_problem_inputs/inputs_std4_RK4
+++ b/Exec/standard_problem_inputs/inputs_std4_RK4
@@ -7,11 +7,11 @@ n_cell = 640 160 4
 max_grid_size_x = 640
 max_grid_size_y = 160
 max_grid_size_z = 4
-dt = 2.5e-14 # RK4 (438s) 5.e-14 unstable
+dt = 2.5e-14 # RK4, 5.e-14 unstable
 
 stop_time = 4.e-9
 
-plot_int = 10000
+plot_int = 5000
 chk_int = -1
 restart = -1
 diag_type = 4

From c8db8c5fabd5dc723259fd38d1ada5e282b4cc57 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Thu, 12 Feb 2026 18:09:20 -0800
Subject: [PATCH 24/30] scaling tests

---
 Exec/scaling_inputs/inputs_scaling_cpu_0.25   | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_cpu_0.5    | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_cpu_1      | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_cpu_128    | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_cpu_16     | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_cpu_2      | 52 +++++++++++++++++++
 .../inputs_scaling_gpu_0.25}                  |  0
 Exec/scaling_inputs/inputs_scaling_gpu_0.5    | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_gpu_1      | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_gpu_128    | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_gpu_16     | 52 +++++++++++++++++++
 Exec/scaling_inputs/inputs_scaling_gpu_2      | 52 +++++++++++++++++++
 12 files changed, 572 insertions(+)
 create mode 100644 Exec/scaling_inputs/inputs_scaling_cpu_0.25
 create mode 100644 Exec/scaling_inputs/inputs_scaling_cpu_0.5
 create mode 100644 Exec/scaling_inputs/inputs_scaling_cpu_1
 create mode 100644 Exec/scaling_inputs/inputs_scaling_cpu_128
 create mode 100644 Exec/scaling_inputs/inputs_scaling_cpu_16
 create mode 100644 Exec/scaling_inputs/inputs_scaling_cpu_2
 rename Exec/{inputs_scaling => scaling_inputs/inputs_scaling_gpu_0.25} (100%)
 create mode 100644 Exec/scaling_inputs/inputs_scaling_gpu_0.5
 create mode 100644 Exec/scaling_inputs/inputs_scaling_gpu_1
 create mode 100644 Exec/scaling_inputs/inputs_scaling_gpu_128
 create mode 100644 Exec/scaling_inputs/inputs_scaling_gpu_16
 create mode 100644 Exec/scaling_inputs/inputs_scaling_gpu_2

diff --git a/Exec/scaling_inputs/inputs_scaling_cpu_0.25 b/Exec/scaling_inputs/inputs_scaling_cpu_0.25
new file mode 100644
index 0000000..6c27b5e
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_cpu_0.25
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 128 128 128
+max_grid_size_x = 32
+max_grid_size_y = 32
+max_grid_size_z = 64
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 512.e-9 512.e-9 512.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<512.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_cpu_0.5 b/Exec/scaling_inputs/inputs_scaling_cpu_0.5
new file mode 100644
index 0000000..cc4359a
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_cpu_0.5
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 256 128 128
+max_grid_size_x = 32
+max_grid_size_y = 32
+max_grid_size_z = 64
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 1024.e-9 512.e-9 512.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_cpu_1 b/Exec/scaling_inputs/inputs_scaling_cpu_1
new file mode 100644
index 0000000..93041fd
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_cpu_1
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 256 256 128
+max_grid_size_x = 32
+max_grid_size_y = 32
+max_grid_size_z = 64
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 1024.e-9 1024.e-9 512.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_cpu_128 b/Exec/scaling_inputs/inputs_scaling_cpu_128
new file mode 100644
index 0000000..4a0ae57
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_cpu_128
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 1024 1024 1024
+max_grid_size_x = 32
+max_grid_size_y = 32
+max_grid_size_z = 64
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 4096.e-9 4096.e-9 4096.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_cpu_16 b/Exec/scaling_inputs/inputs_scaling_cpu_16
new file mode 100644
index 0000000..4c02f6f
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_cpu_16
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 512 512 512
+max_grid_size_x = 32
+max_grid_size_y = 32
+max_grid_size_z = 64
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 2048.e-9 2048.e-9 2048.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_cpu_2 b/Exec/scaling_inputs/inputs_scaling_cpu_2
new file mode 100644
index 0000000..f04eaa1
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_cpu_2
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 256 256 256
+max_grid_size_x = 32
+max_grid_size_y = 32
+max_grid_size_z = 64
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 1024.e-9 1024.e-9 1024.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/inputs_scaling b/Exec/scaling_inputs/inputs_scaling_gpu_0.25
similarity index 100%
rename from Exec/inputs_scaling
rename to Exec/scaling_inputs/inputs_scaling_gpu_0.25
diff --git a/Exec/scaling_inputs/inputs_scaling_gpu_0.5 b/Exec/scaling_inputs/inputs_scaling_gpu_0.5
new file mode 100644
index 0000000..aced596
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_gpu_0.5
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 256 128 128
+max_grid_size_x = 128
+max_grid_size_y = 128
+max_grid_size_z = 128
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 1024.e-9 512.e-9 512.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<1024.e-9)*(y>0.)*(y<512.e-9)*(z>0.)*(z<512.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_gpu_1 b/Exec/scaling_inputs/inputs_scaling_gpu_1
new file mode 100644
index 0000000..a36e82a
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_gpu_1
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 256 256 128
+max_grid_size_x = 128
+max_grid_size_y = 128
+max_grid_size_z = 128
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 1024.e-9 1024.e-9 512.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<512.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_gpu_128 b/Exec/scaling_inputs/inputs_scaling_gpu_128
new file mode 100644
index 0000000..a977856
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_gpu_128
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 1024 1024 1024
+max_grid_size_x = 128
+max_grid_size_y = 128
+max_grid_size_z = 128
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 4096.e-9 4096.e-9 4096.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<4096.e-9)*(y>0.)*(y<4096.e-9)*(z>0.)*(z<4096.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_gpu_16 b/Exec/scaling_inputs/inputs_scaling_gpu_16
new file mode 100644
index 0000000..d490dcd
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_gpu_16
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 512 512 512
+max_grid_size_x = 128
+max_grid_size_y = 128
+max_grid_size_z = 128
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 2048.e-9 2048.e-9 2048.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<2048.e-9)*(y>0.)*(y<2048.e-9)*(z>0.)*(z<2048.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1
diff --git a/Exec/scaling_inputs/inputs_scaling_gpu_2 b/Exec/scaling_inputs/inputs_scaling_gpu_2
new file mode 100644
index 0000000..488cadb
--- /dev/null
+++ b/Exec/scaling_inputs/inputs_scaling_gpu_2
@@ -0,0 +1,52 @@
+amrex.use_gpu_aware_mpi=1
+
+n_cell = 256 256 256
+max_grid_size_x = 128
+max_grid_size_y = 128
+max_grid_size_z = 128
+
+dt = 1.e-15
+nsteps = 100
+plot_int = -1
+chk_int = -1
+restart = -1
+
+prob_lo = 0. 0. 0.
+prob_hi = 1024.e-9 1024.e-9 1024.e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 0
+Hx_bias_parser(x,y,z,t) = "1.e5"
+Hy_bias_parser(x,y,z,t) = "1.e5"
+Hz_bias_parser(x,y,z,t) = "1.e5"
+
+timedependent_alpha = 0
+alpha_parser(x,y,z,t) =    "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * 0.5"
+Ms_parser(x,y,z) =         "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * 8.e5"
+gamma_parser(x,y,z) =      "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * -1.759e11"
+exchange_parser(x,y,z) =   "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * 1.3e-11"
+anisotropy_parser(x,y,z) = "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * -139.26"
+DMI_parser(x,y,z) =        "(x>0.)*(x<1024.e-9)*(y>0.)*(y<1024.e-9)*(z>0.)*(z<1024.e-9) * -4.5e-3"
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 1
+anisotropy_axis = 0.0 0.0 1.0
+DMI_coupling = 1
+
+# INTEGRATION
+
+TimeIntegratorOption = 1 #Forward Euler
+
+integration.type = RungeKutta
+integration.rk.type = 1

From d6d49ce292740c12b51000ee5dc87e970d3ac984 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Fri, 13 Feb 2026 13:43:44 -0800
Subject: [PATCH 25/30] turn on RK4 for STD3

---
 Exec/standard_problem_inputs/inputs_std3_10cell | 4 ++--
 Exec/standard_problem_inputs/inputs_std3_20cell | 4 ++--
 Exec/standard_problem_inputs/inputs_std3_40cell | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/Exec/standard_problem_inputs/inputs_std3_10cell b/Exec/standard_problem_inputs/inputs_std3_10cell
index 45821b9..57072e0 100644
--- a/Exec/standard_problem_inputs/inputs_std3_10cell
+++ b/Exec/standard_problem_inputs/inputs_std3_10cell
@@ -52,10 +52,10 @@ anisotropy_axis = 0.0 0.0 1.0
 DMI_coupling = 0
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std3_20cell b/Exec/standard_problem_inputs/inputs_std3_20cell
index eec0c4f..999c15f 100644
--- a/Exec/standard_problem_inputs/inputs_std3_20cell
+++ b/Exec/standard_problem_inputs/inputs_std3_20cell
@@ -52,10 +52,10 @@ anisotropy_axis = 0.0 0.0 1.0
 DMI_coupling = 0
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std3_40cell b/Exec/standard_problem_inputs/inputs_std3_40cell
index 5da47bb..9755969 100644
--- a/Exec/standard_problem_inputs/inputs_std3_40cell
+++ b/Exec/standard_problem_inputs/inputs_std3_40cell
@@ -52,10 +52,10 @@ anisotropy_axis = 0.0 0.0 1.0
 DMI_coupling = 0
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9

From e9fa1ba67f2bf48153f9bd27d70b49bb932f7e2f Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Sat, 14 Feb 2026 07:35:01 -0800
Subject: [PATCH 26/30] tweaks for final paper

---
 Exec/standard_problem_inputs/inputs_std4_IMEX | 4 ++--
 Exec/standard_problem_inputs/inputs_std4_MRI  | 4 ++--
 Exec/standard_problem_inputs/inputs_std4_RK4  | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/Exec/standard_problem_inputs/inputs_std4_IMEX b/Exec/standard_problem_inputs/inputs_std4_IMEX
index c0cbe9f..f493b29 100644
--- a/Exec/standard_problem_inputs/inputs_std4_IMEX
+++ b/Exec/standard_problem_inputs/inputs_std4_IMEX
@@ -9,9 +9,9 @@ max_grid_size_y = 160
 max_grid_size_z = 4
 dt = 5.e-15 # IMEX, 1.e-14 unstable
 
-stop_time = 4.e-9
+stop_time = 1.25e-12
 
-plot_int = 25000
+plot_int = -1
 chk_int = -1
 restart = -1
 diag_type = 4
diff --git a/Exec/standard_problem_inputs/inputs_std4_MRI b/Exec/standard_problem_inputs/inputs_std4_MRI
index b9770da..b728d7b 100644
--- a/Exec/standard_problem_inputs/inputs_std4_MRI
+++ b/Exec/standard_problem_inputs/inputs_std4_MRI
@@ -9,9 +9,9 @@ max_grid_size_y = 160
 max_grid_size_z = 4
 dt = 1.25e-13 # SUNDIALS MRI, 2.5e-13 unstable
 
-stop_time = 4.e-9
+stop_time = 1.25e-12
 
-plot_int = 1000
+plot_int = -1
 chk_int = -1
 restart = -1
 diag_type = 4
diff --git a/Exec/standard_problem_inputs/inputs_std4_RK4 b/Exec/standard_problem_inputs/inputs_std4_RK4
index 13b08db..4284390 100644
--- a/Exec/standard_problem_inputs/inputs_std4_RK4
+++ b/Exec/standard_problem_inputs/inputs_std4_RK4
@@ -9,9 +9,9 @@ max_grid_size_y = 160
 max_grid_size_z = 4
 dt = 2.5e-14 # RK4, 5.e-14 unstable
 
-stop_time = 4.e-9
+stop_time = 1.25e-12
 
-plot_int = 5000
+plot_int = -1
 chk_int = -1
 restart = -1
 diag_type = 4

From 55ba7df3433ba2187a627a40e073949728b5d0ad Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Sun, 15 Feb 2026 14:25:07 -0800
Subject: [PATCH 27/30] bugfix - for hysteresis do not compute Hbias each time
 step!!!

---
 Source/main.cpp | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/Source/main.cpp b/Source/main.cpp
index 8392b9f..3c91de2 100644
--- a/Source/main.cpp
+++ b/Source/main.cpp
@@ -554,7 +554,9 @@ void main_main ()
                         H_biasfield[idim].setVal(0.);
                     }
                 } else {
-                    ComputeHbias(H_biasfield, time_in, geom);
+                    if (timedependent_Hbias) {
+                        ComputeHbias(H_biasfield, time_in, geom);
+                    }
                 }
 
                 // exchange
@@ -627,7 +629,9 @@ void main_main ()
 
                 // H_bias
                 if (fast_H_bias==1) {
-                    ComputeHbias(H_biasfield, time_in, geom);
+                    if (timedependent_Hbias) {
+                        ComputeHbias(H_biasfield, time_in, geom);
+                    }
                 } else {
                     for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
                         H_biasfield[idim].setVal(0.);
@@ -705,7 +709,9 @@ void main_main ()
 
                 // H_bias
                 if (implicit_H_bias==1) {
-                    ComputeHbias(H_biasfield, time_in, geom);
+                    if (timedependent_Hbias) {
+                        ComputeHbias(H_biasfield, time_in, geom);
+                    }
                 } else {
                     for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
                         H_biasfield[idim].setVal(0.);

From 6475e14bc5f69575580706cb8fc311e4a9b9137a Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Sun, 15 Feb 2026 15:44:44 -0800
Subject: [PATCH 28/30] fix std2 inputs to match text in paper (and get proper
 results)

---
 Exec/standard_problem_inputs/inputs_std2_dlex12 |  8 ++++----
 Exec/standard_problem_inputs/inputs_std2_dlex15 |  8 ++++----
 Exec/standard_problem_inputs/inputs_std2_dlex18 |  8 ++++----
 Exec/standard_problem_inputs/inputs_std2_dlex21 |  8 ++++----
 Exec/standard_problem_inputs/inputs_std2_dlex24 |  8 ++++----
 Exec/standard_problem_inputs/inputs_std2_dlex3  | 14 +++++++-------
 Exec/standard_problem_inputs/inputs_std2_dlex6  | 14 +++++++-------
 Exec/standard_problem_inputs/inputs_std2_dlex9  | 14 +++++++-------
 8 files changed, 41 insertions(+), 41 deletions(-)

diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex12 b/Exec/standard_problem_inputs/inputs_std2_dlex12
index 4d0fae8..0d1fc5b 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex12
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex12
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex15 b/Exec/standard_problem_inputs/inputs_std2_dlex15
index 59e8fa4..09f9f2b 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex15
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex15
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex18 b/Exec/standard_problem_inputs/inputs_std2_dlex18
index fcec880..408868e 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex18
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex18
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex21 b/Exec/standard_problem_inputs/inputs_std2_dlex21
index 94bc5e3..5ca1d34 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex21
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex21
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex24 b/Exec/standard_problem_inputs/inputs_std2_dlex24
index 9590508..dcc9630 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex24
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex24
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex3 b/Exec/standard_problem_inputs/inputs_std2_dlex3
index 14c255c..a869eb4 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex3
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex3
@@ -29,9 +29,9 @@ Mz_parser(x,y,z) = "8.e5/sqrt(1.5) * (x>0.)*(x<75.e-9)*(y>0.)*(y<15.e-9)*(z>0.)*
 # Field Ensemble for hysteresis:  [Hx/Ms= +.08 -> -.08 -> +.08], [Hy/Ms= +.08 -> -.08 -> +.08], mu_0 [Hz/Ms= +.08 -> -.08 -> +.08]
 # which is a field directed about 57 degrees offset from each coordinate axis
 timedependent_Hbias = 0
-Hx_bias_parser(x,y,z,t) = "1.54e4"
-Hy_bias_parser(x,y,z,t) = "1.54e4"
-Hz_bias_parser(x,y,z,t) = "1.54e4"
+Hx_bias_parser(x,y,z,t) = "6.4e4"
+Hy_bias_parser(x,y,z,t) = "6.4e4"
+Hz_bias_parser(x,y,z,t) = "6.4e4"
 
 timedependent_alpha = 0
 alpha_parser(x,y,z,t) =  " (x>0.)*(x<75.e-9)*(y>0.)*(y<15.e-9)*(z>0.)*(z<1.5e-9) * 0.5"
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex6 b/Exec/standard_problem_inputs/inputs_std2_dlex6
index 2aef3d3..96b0a7c 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex6
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex6
@@ -29,9 +29,9 @@ Mz_parser(x,y,z) = "8.e5/sqrt(3.) * (x>0.)*(x<150.e-9)*(y>0.)*(y<30.e-9)*(z>0.)*
 # Field Ensemble for hysteresis:  [Hx/Ms= +.08 -> -.08 -> +.08], [Hy/Ms= +.08 -> -.08 -> +.08], mu_0 [Hz/Ms= +.08 -> -.08 -> +.08]
 # which is a field directed about 57 degrees offset from each coordinate axis
 timedependent_Hbias = 0
-Hx_bias_parser(x,y,z,t) = "3.4e4"
-Hy_bias_parser(x,y,z,t) = "3.4e4"
-Hz_bias_parser(x,y,z,t) = "3.4e4"
+Hx_bias_parser(x,y,z,t) = "6.4e4"
+Hy_bias_parser(x,y,z,t) = "6.4e4"
+Hz_bias_parser(x,y,z,t) = "6.4e4"
 
 timedependent_alpha = 0
 alpha_parser(x,y,z,t) =  " (x>0.)*(x<150.e-9)*(y>0.)*(y<30.e-9)*(z>0.)*(z<3.e-9) * 0.5"
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_std2_dlex9 b/Exec/standard_problem_inputs/inputs_std2_dlex9
index a680a0b..e9861ff 100644
--- a/Exec/standard_problem_inputs/inputs_std2_dlex9
+++ b/Exec/standard_problem_inputs/inputs_std2_dlex9
@@ -29,9 +29,9 @@ Mz_parser(x,y,z) = "8.e5/sqrt(3.) * (x>0.)*(x<225.e-9)*(y>0.)*(y<45.e-9)*(z>0.)*
 # Field Ensemble for hysteresis:  [Hx/Ms= +.08 -> -.08 -> +.08], [Hy/Ms= +.08 -> -.08 -> +.08], mu_0 [Hz/Ms= +.08 -> -.08 -> +.08]
 # which is a field directed about 57 degrees offset from each coordinate axis
 timedependent_Hbias = 0
-Hx_bias_parser(x,y,z,t) = "4.54e4"
-Hy_bias_parser(x,y,z,t) = "4.54e4"
-Hz_bias_parser(x,y,z,t) = "4.54e4"
+Hx_bias_parser(x,y,z,t) = "6.4e4"
+Hy_bias_parser(x,y,z,t) = "6.4e4"
+Hz_bias_parser(x,y,z,t) = "6.4e4"
 
 timedependent_alpha = 0
 alpha_parser(x,y,z,t) =  " (x>0.)*(x<225.e-9)*(y>0.)*(y<45.e-9)*(z>0.)*(z<4.5e-9) * 0.5"
@@ -54,16 +54,16 @@ DMI_coupling = 0
 equilibrium_tolerance = 1.e-9
 
 # hysteresis increment for H_bias
-increment_size = 1280
+increment_size = 128
 
 # number of increments ONE WAY: Must = 2*(Hx_bias/increment_size)
-nsteps_hysteresis = 100
+nsteps_hysteresis = 1000
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9

From e1e2708f40c4758ce2d00256670f30095924b16a Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Sun, 15 Feb 2026 19:10:15 -0800
Subject: [PATCH 29/30] fix up DMI1 problem

---
 Exec/DMI_diagnostic/DMI_diagnostic.cpp   |  4 ++--
 Exec/standard_problem_inputs/inputs_dmi1 | 11 ++++++-----
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/Exec/DMI_diagnostic/DMI_diagnostic.cpp b/Exec/DMI_diagnostic/DMI_diagnostic.cpp
index dab694b..5e86909 100644
--- a/Exec/DMI_diagnostic/DMI_diagnostic.cpp
+++ b/Exec/DMI_diagnostic/DMI_diagnostic.cpp
@@ -118,10 +118,10 @@ main (int   argc,
             int k = (hi.z+1)/2;
             int j = (hi.y+1)/2;
             for (auto i = (hi.x+1)/2; i <= hi.x; ++i) {
-                std::cout << i << " " << mfdata(i,j,k,3)/1.1e6 << " " << mfdata(i,j,k,16)+offset << "\n";
+                std::cout << i << " " << mfdata(i,j,k,2)/1.1e6 << " " << mfdata(i,j,k,3)+offset << "\n";
                 // 2pi cyclic fix
                 if (i<hi.x) {
-                    if (offset==0. && mfdata(i,j,k,16)>5. && mfdata(i+1,j,k,16)<1.) {
+                    if (offset==0. && mfdata(i,j,k,3)>5. && mfdata(i+1,j,k,3)<1.) {
                         offset = 2.*M_PI;
                     }
                 }
diff --git a/Exec/standard_problem_inputs/inputs_dmi1 b/Exec/standard_problem_inputs/inputs_dmi1
index 962538a..40e5370 100644
--- a/Exec/standard_problem_inputs/inputs_dmi1
+++ b/Exec/standard_problem_inputs/inputs_dmi1
@@ -15,7 +15,7 @@ plot_exchange = 0
 plot_DMI = 0
 plot_anisotropy = 0
 plot_demag = 0
-plot_theta = 0
+plot_theta = 1
 
 prob_lo = -50.0e-9 -50.0e-9 -1.0e-9
 prob_hi =  50.0e-9  50.0e-9  1.0e-9
@@ -40,10 +40,11 @@ DMI_parser(x,y,z) = "        -4.5e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9)
 anisotropy_parser(x,y,z) = " 510.0e3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^3
 anisotropy_axis = 0.0 0.0 1.0
 
-# quasi-uniform initialization: change the value before pi to 0; or Mx = My = 0; Mz(i,j,k) = (sqrt(x*x + y*y) <= 50.0e-9) ? Ms_arr(i,j,k):0.0;
-# skyrmion initialization: sin(pi * sqrt(x*x + y*y) / 50.0e-9)
-# 2pi state initialization: sin(2*pi * sqrt(x*x + y*y) / 50.0e-9)
-# 3pi state initialization: sin(3*pi * sqrt(x*x + y*y) / 50.0e-9)
+# change coefficient in front of pi as follows:
+# quasi-uniform initialization: sin and cos(0*pi ...
+# skyrmion initialization     : sin and cos(1*pi ...
+# 2pi state initialization    : sin and cos(2*pi ...
+# 3pi state initialization    : sin and cos(3*pi ...
 Mx_parser(x,y,z) = "1.1e6 * sin(2*3.141592654 * sqrt(x*x + y*y) / 50.0e-9) * x/sqrt(x*x + y*y) *(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)"
 My_parser(x,y,z) = "1.1e6 * sin(2*3.141592654 * sqrt(x*x + y*y) / 50.0e-9) * y/sqrt(x*x + y*y) *(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)"
 Mz_parser(x,y,z) = "1.1e6 * cos(2*3.141592654 * sqrt(x*x + y*y) / 50.0e-9) *(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # for the skyrmion state

From 87a43f93e581bb99b78224b7a8795af2e6633bee Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Sun, 15 Feb 2026 21:13:30 -0800
Subject: [PATCH 30/30] paper update

---
 Exec/standard_problem_inputs/inputs_dmi2_100 | 16 ++++++----------
 Exec/standard_problem_inputs/inputs_dmi2_25  | 16 ++++++----------
 Exec/standard_problem_inputs/inputs_dmi2_50  | 16 ++++++----------
 Exec/standard_problem_inputs/inputs_dmi2_75  | 16 ++++++----------
 4 files changed, 24 insertions(+), 40 deletions(-)

diff --git a/Exec/standard_problem_inputs/inputs_dmi2_100 b/Exec/standard_problem_inputs/inputs_dmi2_100
index 6d4b8e7..ab8e0be 100644
--- a/Exec/standard_problem_inputs/inputs_dmi2_100
+++ b/Exec/standard_problem_inputs/inputs_dmi2_100
@@ -7,9 +7,11 @@ max_grid_size_y = 256
 max_grid_size_z = 4
 dt = 2.5e-14
 nsteps = 100000
-plot_int = 1000
+plot_int = 100000
 chk_int = 100000
 
+restart = -1
+
 prob_lo = -100.0e-9 -100.0e-9 -1.0e-9
 prob_hi =  100.0e-9  100.0e-9  1.0e-9
 
@@ -32,9 +34,7 @@ anisotropy_parser(x,y,z) = " 510.0e3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9
 anisotropy_axis = 0.0 0.0 1.0
 
 # D_c = 4 sqrt(AK) / pi = 3.6371e-3
-#DMI_parser(x,y,z) = "       -1.25e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -1.5e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -2.0e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
+#DMI_parser(x,y,z) = "        -2.3e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -2.5e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -3.0e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "       -3.25e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
@@ -46,10 +46,6 @@ DMI_parser(x,y,z) = "      -3.625e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9
 #DMI_parser(x,y,z) = "       -4.25e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -4.5e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -5.0e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -5.5e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.0e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.5e-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "         -7.0-3*(sqrt(x*x + y*y) <= 100.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 
 # quasi-uniform initialization: change the value before pi to 0; or Mx = My = 0; Mz(i,j,k) = (sqrt(x*x + y*y) <= 100.0e-9) ? Ms_arr(i,j,k):0.0;
 # skyrmion initialization: sin(pi * sqrt(x*x + y*y) / 100.0e-9)
@@ -65,10 +61,10 @@ Hy_bias_parser(x,y,z,t) = "0.0"
 Hz_bias_parser(x,y,z,t) = "0.0"
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_dmi2_25 b/Exec/standard_problem_inputs/inputs_dmi2_25
index 063a169..aa088bf 100644
--- a/Exec/standard_problem_inputs/inputs_dmi2_25
+++ b/Exec/standard_problem_inputs/inputs_dmi2_25
@@ -7,9 +7,11 @@ max_grid_size_y = 64
 max_grid_size_z = 4
 dt = 2.5e-14
 nsteps = 100000
-plot_int = 1000
+plot_int = 100000
 chk_int = 100000
 
+restart = -1
+
 prob_lo = -25.0e-9 -25.0e-9 -1.0e-9
 prob_hi =  25.0e-9  25.0e-9  1.0e-9
 
@@ -32,9 +34,7 @@ anisotropy_parser(x,y,z) = " 510.0e3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9)
 anisotropy_axis = 0.0 0.0 1.0
 
 # D_c = 4 sqrt(AK) / pi = 3.6371e-3
-#DMI_parser(x,y,z) = "       -1.25e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -1.5e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -2.0e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
+#DMI_parser(x,y,z) = "        -2.3e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -2.5e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -3.0e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "       -3.25e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
@@ -46,10 +46,6 @@ DMI_parser(x,y,z) = "      -3.625e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9)
 #DMI_parser(x,y,z) = "       -4.25e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -4.5e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -5.0e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -5.5e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.0e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.5e-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "         -7.0-3*(sqrt(x*x + y*y) <= 25.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 
 # quasi-uniform initialization: change the value before pi to 0; or Mx = My = 0; Mz(i,j,k) = (sqrt(x*x + y*y) <= 25.0e-9) ? Ms_arr(i,j,k):0.0;
 # skyrmion initialization: sin(pi * sqrt(x*x + y*y) / 25.0e-9)
@@ -65,10 +61,10 @@ Hy_bias_parser(x,y,z,t) = "0.0"
 Hz_bias_parser(x,y,z,t) = "0.0"
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_dmi2_50 b/Exec/standard_problem_inputs/inputs_dmi2_50
index 1a8f769..f012ea7 100644
--- a/Exec/standard_problem_inputs/inputs_dmi2_50
+++ b/Exec/standard_problem_inputs/inputs_dmi2_50
@@ -7,9 +7,11 @@ max_grid_size_y = 128
 max_grid_size_z = 4
 dt = 2.5e-14
 nsteps = 100000
-plot_int = 1000
+plot_int = 100000
 chk_int = 100000
 
+restart = -1
+
 prob_lo = -50.0e-9 -50.0e-9 -1.0e-9
 prob_hi =  50.0e-9  50.0e-9  1.0e-9
 
@@ -32,9 +34,7 @@ anisotropy_parser(x,y,z) = " 510.0e3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9)
 anisotropy_axis = 0.0 0.0 1.0
 
 # D_c = 4 sqrt(AK) / pi = 3.6371e-3
-#DMI_parser(x,y,z) = "       -1.25e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -1.5e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -2.0e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
+#DMI_parser(x,y,z) = "        -2.3e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -2.5e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -3.0e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "       -3.25e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
@@ -46,10 +46,6 @@ DMI_parser(x,y,z) = "      -3.625e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9)
 #DMI_parser(x,y,z) = "       -4.25e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -4.5e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -5.0e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -5.5e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.0e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.5e-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "         -7.0-3*(sqrt(x*x + y*y) <= 50.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 
 # quasi-uniform initialization: change the value before pi to 0; or Mx = My = 0; Mz(i,j,k) = (sqrt(x*x + y*y) <= 50.0e-9) ? Ms_arr(i,j,k):0.0;
 # skyrmion initialization: sin(pi * sqrt(x*x + y*y) / 50.0e-9)
@@ -65,10 +61,10 @@ Hy_bias_parser(x,y,z,t) = "0.0"
 Hz_bias_parser(x,y,z,t) = "0.0"
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9
diff --git a/Exec/standard_problem_inputs/inputs_dmi2_75 b/Exec/standard_problem_inputs/inputs_dmi2_75
index 903dce4..84fd27b 100644
--- a/Exec/standard_problem_inputs/inputs_dmi2_75
+++ b/Exec/standard_problem_inputs/inputs_dmi2_75
@@ -7,9 +7,11 @@ max_grid_size_y = 192
 max_grid_size_z = 4
 dt = 2.5e-14
 nsteps = 100000
-plot_int = 1000
+plot_int = 100000
 chk_int = 100000
 
+restart = -1
+
 prob_lo = -75.0e-9 -75.0e-9 -1.0e-9
 prob_hi =  75.0e-9  75.0e-9  1.0e-9
 
@@ -32,9 +34,7 @@ anisotropy_parser(x,y,z) = " 510.0e3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9)
 anisotropy_axis = 0.0 0.0 1.0
 
 # D_c = 4 sqrt(AK) / pi = 3.6371e-3
-#DMI_parser(x,y,z) = "       -1.25e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -1.5e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -2.0e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
+#DMI_parser(x,y,z) = "        -2.3e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -2.5e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -3.0e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "       -3.25e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
@@ -46,10 +46,6 @@ DMI_parser(x,y,z) = "      -3.625e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9)
 #DMI_parser(x,y,z) = "       -4.25e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -4.5e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 #DMI_parser(x,y,z) = "        -5.0e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -5.5e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.0e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "        -6.5e-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
-#DMI_parser(x,y,z) = "         -7.0-3*(sqrt(x*x + y*y) <= 75.0e-9) * (z>=-1.0e-9) * (z<=1.0e-9)" # unit is J/m^2
 
 # quasi-uniform initialization: change the value before pi to 0; or Mx = My = 0; Mz(i,j,k) = (sqrt(x*x + y*y) <= 75.0e-9) ? Ms_arr(i,j,k):0.0;
 # skyrmion initialization: sin(pi * sqrt(x*x + y*y) / 75.0e-9)
@@ -65,10 +61,10 @@ Hy_bias_parser(x,y,z,t) = "0.0"
 Hz_bias_parser(x,y,z,t) = "0.0"
 
 # INTEGRATION
-TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 1 #Forward Euler
 #TimeIntegratorOption = 2 #Predictor-corrector
 #TimeIntegratorOption = 3 #2nd order artemis way
-#TimeIntegratorOption = 4 #amrex/sundials backend integrators
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
 
 # tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
 iterative_tolerance = 1.e-9