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/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
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
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
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
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
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
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.);