new Laplacian for exchange BC

Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_exchangeBC

@@ -0,0 +1,88 @@
+################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 60000
+amr.n_cell = 64 64 64 # number of cells spanning the domain in each coordinate direction at level 0
+amr.max_grid_size = 32 # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 16
+geometry.coord_sys = 0
+
+geometry.prob_lo = -10e-9 -10e-9 -10.0e-9
+geometry.prob_hi =  10e-9  10e-9  10.0e-9
+boundary.field_lo = periodic periodic periodic
+boundary.field_hi = periodic periodic periodic
+amr.max_level = 0
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.verbose = 1
+warpx.use_filter = 0
+warpx.cfl = 1000
+warpx.mag_time_scheme_order = 2 # default 1
+warpx.mag_M_normalization = 1 # 1 is saturated
+warpx.mag_LLG_coupling = 0
+warpx.mag_LLG_exchange_coupling = 1
+
+algo.em_solver_medium = macroscopic # vacuum/macroscopic
+algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler
+
+macroscopic.sigma_function(x,y,z) = "0.0"
+macroscopic.epsilon_function(x,y,z) = "8.8541878128e-12"
+macroscopic.mu_function(x,y,z) = "1.25663706212e-06"
+
+my_constants.mag_lo_x = -5.e-9
+my_constants.mag_hi_x = 5.e-9
+my_constants.mag_lo_y = -5.e-9
+my_constants.mag_hi_y = 5.e-9
+my_constants.mag_lo_z = -5.e-9
+my_constants.mag_hi_z = 5.e-9
+
+#unit conversion: 1 Gauss = (1000/4pi) A/m
+macroscopic.mag_Ms_init_style = "parse_mag_Ms_function" # parse or "constant"
+macroscopic.mag_Ms_function(x,y,z) = "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"  # in unit A/m, equal to 1750 Gauss; Ms must be nonzero for LLG
+
+macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant"
+macroscopic.mag_alpha_function(x,y,z) = "0.058 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"  # alpha is unitless, calculated from linewidth Delta_H = 40 Oersted
+
+macroscopic.mag_gamma_init_style = "parse_mag_gamma_function" # parse or "constant"
+macroscopic.mag_gamma_function(x,y,z) = "-1.759e11 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # gyromagnetic ratio is constant for electrons in all materials
+
+macroscopic.mag_exchange_init_style = "parse_mag_exchange_function" # parse or "constant"
+macroscopic.mag_exchange_function(x,y,z) = "3.76e-12 * (x >= mag_lo_x) * (x <= mag_hi_x) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # exchange coupling constanit; Must be non-zero when exchange coupling is ON
+
+macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step
+macroscopic.mag_tol = 1.e-6 # M magnitude relative error tolerance compared to previous iteration
+macroscopic.mag_normalized_error = 0.1 # if M magnitude relatively changes more than this value, raise a red flag
+
+#################################
+############ FIELDS #############
+#################################
+warpx.E_ext_grid_init_style = parse_E_ext_grid_function
+warpx.Ex_external_grid_function(x,y,z) = 0.
+warpx.Ey_external_grid_function(x,y,z) = 0.
+warpx.Ez_external_grid_function(x,y,z) = 0.
+
+warpx.H_ext_grid_init_style = parse_H_ext_grid_function
+warpx.Hx_external_grid_function(x,y,z)= 0.
+warpx.Hy_external_grid_function(x,y,z) = 0.
+warpx.Hz_external_grid_function(x,y,z) = 0.
+
+#unit conversion: 1 Gauss = 1 Oersted = (1000/4pi) A/m
+#calculation of H_bias: H_bias (oe) = frequency / 2.8e6
+warpx.H_bias_ext_grid_init_style = parse_H_bias_ext_grid_function
+warpx.Hx_bias_external_grid_function(x,y,z)= 0.
+warpx.Hy_bias_external_grid_function(x,y,z)= "3.7e4" # in A/m, equal to 464 Oersted
+warpx.Hz_bias_external_grid_function(x,y,z)= 0.
+
+warpx.M_ext_grid_init_style = parse_M_ext_grid_function
+warpx.Mx_external_grid_function(x,y,z)= "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y < 0.) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"
+warpx.My_external_grid_function(x,y,z) = 0.
+warpx.Mz_external_grid_function(x,y,z) = "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= 0.) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"
+
+#Diagnostics
+diagnostics.diags_names = plt
+plt.intervals = 10
+plt.diag_type = Full
+plt.fields_to_plot = Ex Ey Ez Hx Hy Hz Bx By Bz Mx_xface My_xface Mz_xface Mx_yface My_yface Mz_yface Mx_zface My_zface Mz_zface
+plt.plot_raw_fields = 0

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H

@@ -15,6 +15,9 @@
 #include <AMReX_Array4.H>
 #include <AMReX_Gpu.H>
 #include <AMReX_REAL.H>
+#include "Utils/WarpXUtil.H"
+#include "WarpX.H"
+#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
 
 #include <array>
 #include <cmath>
@@ -223,6 +226,131 @@ struct CartesianYeeAlgorithm {
 
 #endif
 
+#ifdef WARPX_MAG_LLG
+
+    /**
+     * Perform divergence of gradient along x on M field when exchange coupling is on */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real LaplacianDx_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_x, int const n_coefs_x, amrex::Real const Ms_lo_x, amrex::Real const Ms_hi_x,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        amrex::Real const inv_dx = coefs_x[0];
+    if(ncomp != nodality){
+          if ( Ms_hi_x == 0.) //non-magnetic on right
+          {
+                return inv_dx*inv_dx*(2.*F(i,j,k,ncomp) - 5.*F(i-1,j,k,ncomp) + 4.*F(i-2,j,k,ncomp) - F(i-3,j,k,ncomp));
+          }
+          else if (Ms_lo_x == 0.) //non-magnetic on left
+          {
+                return inv_dx*inv_dx*(2.*F(i,j,k,ncomp) - 5.*F(i+1,j,k,ncomp) + 4.*F(i+2,j,k,ncomp) - F(i+3,j,k,ncomp));
+      } else
+          {
+             return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp));
+          }
+    } else {
+          if ( Ms_hi_x == 0.) //non-magnetic on right
+          {
+             return inv_dx*(0. - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp));
+          }
+          else if (Ms_lo_x == 0.) //non-magnetic on left
+          {
+             return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - 0.);
+      } else
+          {
+             return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp));
+          }
+    }
+    }
+
+    /**
+     * Perform divergence of gradient along y on M field when exchange coupling is on*/
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real LaplacianDy_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_y, int const n_coefs_y, amrex::Real const Ms_lo_y, amrex::Real const Ms_hi_y,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        amrex::Real const inv_dy = coefs_y[0];
+    if(ncomp != nodality){
+          if ( Ms_hi_y == 0.) //non-magnetic on right
+          {
+                return inv_dy*inv_dy*(2.*F(i,j,k,ncomp) - 5.*F(i,j-1,k,ncomp) + 4.*F(i,j-2,k,ncomp) - F(i,j-3,k,ncomp));
+          }
+          else if (Ms_lo_y == 0.) //non-magnetic on left
+          {
+                return inv_dy*inv_dy*(2.*F(i,j,k,ncomp) - 5.*F(i,j+1,k,ncomp) + 4.*F(i,j+2,k,ncomp) - F(i,j+3,k,ncomp));
+      } else
+          {
+             return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp));
+          }
+    } else {
+          if ( Ms_hi_y == 0.) //non-magnetic on right
+          {
+             return inv_dy*(0. - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp));
+          }
+          else if (Ms_lo_y == 0.) //non-magnetic on left
+          {
+             return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - 0.);
+      } else
+          {
+             return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp));
+          }
+    }
+    }
+
+     /**
+     * Perform divergence of gradient along z on M field when exchange coupling is on*/
+   AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real LaplacianDz_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_z, int const n_coefs_z, amrex::Real const Ms_lo_z, amrex::Real const Ms_hi_z,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        amrex::Real const inv_dz = coefs_z[0];
+    if(ncomp != nodality){
+          if ( Ms_hi_z == 0.) //non-magnetic on right
+          {
+                return inv_dz*inv_dz*(2.*F(i,j,k,ncomp) - 5.*F(i,j,k-1,ncomp) + 4.*F(i,j,k-2,ncomp) - F(i,j,k-3,ncomp));
+          }
+          else if (Ms_lo_z == 0.) //non-magnetic on left
+          {
+                return inv_dz*inv_dz*(2.*F(i,j,k,ncomp) - 5.*F(i,j,k+1,ncomp) + 4.*F(i,j,k+2,ncomp) - F(i,j,k+3,ncomp));
+      } else
+          {
+             return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp));
+          }
+    } else {
+          if ( Ms_hi_z == 0.) //non-magnetic on right
+          {
+             return inv_dz*(0. - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp));
+          }
+          else if (Ms_lo_z == 0.) //non-magnetic on left
+          {
+             return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - 0.);
+      } else
+          {
+             return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp));
+          }
+    }
+    }
+
+     /**
+     * Compute the sum to get Laplacian of M field when exchange coupling is on*/
+   AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real Laplacian_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_x, amrex::Real const * const coefs_y, amrex::Real const * const coefs_z,
+        int const n_coefs_x, int const n_coefs_y, int const n_coefs_z,
+        amrex::Real const Ms_lo_x, amrex::Real const Ms_hi_x, amrex::Real const Ms_lo_y, amrex::Real const Ms_hi_y, amrex::Real const Ms_lo_z, amrex::Real const Ms_hi_z,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        return LaplacianDx_Mag(F, coefs_x, n_coefs_x, Ms_lo_x, Ms_hi_x, i, j, k, ncomp, nodality) + LaplacianDy_Mag(F, coefs_y, n_coefs_y, Ms_lo_y, Ms_hi_y, i, j, k, ncomp, nodality) + LaplacianDz_Mag(F, coefs_z, n_coefs_z, Ms_lo_z, Ms_hi_z, i, j, k, ncomp, nodality);
+    }
+
+#endif
+
 };
 
 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM.cpp

@@ -141,9 +141,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian(
         Box const &tbz = mfi.tilebox(Hfield[2]->ixType().toIntVect());
 
         // Extract stencil coefficients for calculating the exchange field H_exchange and the anisotropy field H_anisotropy
-        amrex::Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
-        amrex::Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr();
-        amrex::Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
+        amrex::Real const *const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
+        int const n_coefs_x = m_stencil_coefs_x.size();
+        amrex::Real const *const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr();
+        int const n_coefs_y = m_stencil_coefs_y.size();
+        amrex::Real const *const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
+        int const n_coefs_z = m_stencil_coefs_z.size();
 
         const auto dx = warpx.Geom(lev).CellSizeArray();
         const auto problo = warpx.Geom(lev).ProbLoArray();
@@ -185,9 +188,23 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian(
                         // H_exchange - use M^(old_time)
                         if (mag_exchange_arrx == 0._rt) amrex::Abort("The mag_exchange_arrx is 0.0 while including the exchange coupling term H_exchange for H_eff");
                         amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_arrx / PhysConst::mu0 / mag_Ms_arrx / mag_Ms_arrx;
-                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+
+                        WarpXUtilAlgo::getCellCoordinates(i-1, j, k, Mx_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_x = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i+1, j, k, Mx_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_x = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j-1, k, Mx_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_y = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j+1, k, Mx_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_y = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j, k-1, Mx_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_z = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j, k+1, Mx_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_z = mag_parser(x,y,z);
+
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 0); //Last argument is nodality -- xface = 0
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 0); //Last argument is nodality -- xface = 0
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_xface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 0); //Last argument is nodality -- xface = 0
                     }
 
                     if (mag_anisotropy_coupling == 1){
@@ -298,9 +315,23 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian(
                         // H_exchange - use M^(old_time)
                         if (mag_exchange_arry == 0._rt) amrex::Abort("The mag_exchange_arry is 0.0 while including the exchange coupling term H_exchange for H_eff");
                         amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_arry / PhysConst::mu0 / mag_Ms_arry / mag_Ms_arry;
-                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+
+                        WarpXUtilAlgo::getCellCoordinates(i-1, j, k, My_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_x = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i+1, j, k, My_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_x = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j-1, k, My_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_y = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j+1, k, My_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_y = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j, k-1, My_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_z = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j, k+1, My_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_z = mag_parser(x,y,z);
+
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 1); //Last argument is nodality -- yface = 1
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 1); //Last argument is nodality -- yface = 1
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_yface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 1); //Last argument is nodality -- yface = 1
                     }
 
                     if (mag_anisotropy_coupling == 1){
@@ -412,9 +443,23 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian(
                         // H_exchange - use M^(old_time)
                         if (mag_exchange_arrz == 0._rt) amrex::Abort("The mag_exchange_arrz is 0.0 while including the exchange coupling term H_exchange for H_eff");
                         amrex::Real const H_exchange_coeff = 2.0 * mag_exchange_arrz / PhysConst::mu0 / mag_Ms_arrz / mag_Ms_arrz;
-                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_old_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+
+                        WarpXUtilAlgo::getCellCoordinates(i-1, j, k, Mz_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_x = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i+1, j, k, Mz_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_x = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j-1, k, Mz_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_y = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j+1, k, Mz_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_y = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j, k-1, Mz_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_lo_z = mag_parser(x,y,z);
+                        WarpXUtilAlgo::getCellCoordinates(i, j, k+1, Mz_stag, problo, dx, x, y, z);
+                        amrex::Real Ms_hi_z = mag_parser(x,y,z);
+
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 0, 2); //Last argument is nodality -- zface = 2
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 1, 2); //Last argument is nodality -- zface = 2
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_old_zface, coefs_x, coefs_y, coefs_z, n_coefs_x, n_coefs_y, n_coefs_z, Ms_lo_x, Ms_hi_x, Ms_lo_y, Ms_hi_y, Ms_lo_z, Ms_hi_z, i, j, k, 2, 2); //Last argument is nodality -- zface = 2
                     }
 
                     if (mag_anisotropy_coupling == 1){