new Laplacian for exchange BC

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H

@@ -223,6 +223,156 @@ 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 i, int const j, int const k, int const ncomp=0 ) {
+
+        amrex::Real const inv_dx = coefs_x[0];
+        if ( F(i+1,j,k) == 0.) //non-magnetic on right
+        {
+           if (F(i,j+1,k) == 0. || F(i,j,k+1) == 0.) //Top corner is non-magnetic -- left sided 2nd order derivative
+           {
+              return inv_dx*inv_dx*(2.*F(i,j,k) - 5.*F(i-1,j,k) + 4.*F(i-2,j,k) - F(i-3,j,k));
+           }
+           else if (F(i,j-1,k) == 0. || F(i,j,k-1) == 0.) //Bottom corner is non-magnetic -- right sided 2nd order derivative
+           {
+              return inv_dx*inv_dx*(2.*F(i,j,k) - 5.*F(i+1,j,k) + 4.*F(i+2,j,k) - F(i+3,j,k));
+           }
+           else
+           {
+           return inv_dx*(0. - DownwardDx(F, coefs_x, i, j, k, ncomp));
+           }
+        }
+        else if (F(i-1,j,k) == 0.) //non-magnetic on left
+        {
+           if (F(i,j+1,k) == 0. || F(i,j,k+1) == 0.) //Top corner is non-magnetic -- left sided 2nd order derivative
+           {
+              return inv_dx*inv_dx*(2.*F(i,j,k) - 5.*F(i-1,j,k) + 4.*F(i-2,j,k) - F(i-3,j,k));
+           }
+           else if (F(i,j-1,k) == 0. || F(i,j,k-1) == 0.) //Bottom corner is non-magnetic -- right sided 2nd order derivative
+           {
+              return inv_dx*inv_dx*(2.*F(i,j,k) - 5.*F(i+1,j,k) + 4.*F(i+2,j,k) - F(i+3,j,k));
+           }
+           else
+           {
+              return inv_dx*(UpwardDx(F, coefs_x, i, j, k, ncomp) - 0.);
+           }
+        }
+        else
+        {
+           return inv_dx*(UpwardDx(F, coefs_x, i, j, k, ncomp) - DownwardDx(F, 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 i, int const j, int const k, int const ncomp=0 ) {
+
+        amrex::Real const inv_dy = coefs_y[0];
+        if ( F(i,j+1,k) == 0.) //non-magnetic on right
+        {
+           if (F(i+1,j,k) == 0. || F(i,j,k+1) == 0.) //Top corner is non-magnetic -- left sided 2nd order derivative
+           {
+              return inv_dy*inv_dy*(2.*F(i,j,k) - 5.*F(i,j-1,k) + 4.*F(i,j-2,k) - F(i,j-3,k));
+           }
+           else if (F(i-1,j,k) == 0. || F(i,j,k-1) == 0.) //Bottom corner is non-magnetic -- right sided 2nd order derivative
+           {
+              return inv_dy*inv_dy*(2.*F(i,j,k) - 5.*F(i,j+1,k) + 4.*F(i,j+2,k) - F(i,j+3,k));
+           }
+           else
+           {
+           return inv_dy*(0. - DownwardDy(F, coefs_y, i, j, k, ncomp));
+           }
+        }
+        else if (F(i,j-1,k) == 0.) //non-magnetic on left
+        {
+           if (F(i+1,j,k) == 0. || F(i,j,k+1) == 0.) //Top corner is non-magnetic -- left sided 2nd order derivative
+           {
+              return inv_dy*inv_dy*(2.*F(i,j,k) - 5.*F(i,j-1,k) + 4.*F(i,j-2,k) - F(i,j-3,k));
+           }
+           else if (F(i-1,j,k) == 0. || F(i,j,k-1) == 0.) //Bottom corner is non-magnetic -- right sided 2nd order derivative
+           {
+              return inv_dy*inv_dy*(2.*F(i,j,k) - 5.*F(i,j+1,k) + 4.*F(i,j+2,k) - F(i,j+3,k));
+           }
+           else
+           {
+              return inv_dy*(UpwardDy(F, coefs_y, i, j, k, ncomp) - 0.);
+           }
+        }
+        else
+        {
+           return inv_dy*(UpwardDy(F, coefs_y, i, j, k, ncomp) - DownwardDy(F, 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 i, int const j, int const k, int const ncomp=0 ) {
+
+        amrex::Real const inv_dz = coefs_z[0];
+        if ( F(i,j,k+1) == 0.) //non-magnetic on right
+        {
+           if (F(i+1,j,k) == 0. || F(i,j+1,k) == 0.) //Top corner is non-magnetic -- left sided 2nd order derivative
+           {
+              return inv_dz*inv_dz*(2.*F(i,j,k) - 5.*F(i,j,k-1) + 4.*F(i,j,k-2) - F(i,j,k-3));
+           }
+           else if (F(i-1,j,k) == 0. || F(i,j-1,k) == 0.) //Bottom corner is non-magnetic -- right sided 2nd order derivative
+           {
+              return inv_dz*inv_dz*(2.*F(i,j,k) - 5.*F(i,j,k+1) + 4.*F(i,j,k+2) - F(i,j,k+3));
+           }
+           else
+           {
+           return inv_dz*(0. - DownwardDz(F, coefs_z, i, j, k, ncomp));
+           }
+        }
+        else if (F(i,j,k-1) == 0.) //non-magnetic on left
+        {
+           if (F(i+1,j,k) == 0. || F(i,j+1,k) == 0.) //Top corner is non-magnetic -- left sided 2nd order derivative
+           {
+              return inv_dz*inv_dz*(2.*F(i,j,k) - 5.*F(i,j,k-1) + 4.*F(i,j,k-2) - F(i,j,k-3));
+           }
+           else if (F(i-1,j,k) == 0. || F(i,j-1,k) == 0.) //Bottom corner is non-magnetic -- right sided 2nd order derivative
+           {
+              return inv_dz*inv_dz*(2.*F(i,j,k) - 5.*F(i,j,k+1) + 4.*F(i,j,k+2) - F(i,j,k+3));
+           }
+           else
+           {
+              return inv_dz*(UpwardDz(F, coefs_z, i, j, k, ncomp) - 0.);
+           }
+        }
+        else
+        {
+           return inv_dz*(UpwardDz(F, coefs_z, i, j, k, ncomp) - DownwardDz(F, 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 i, int const j, int const k, int const ncomp=0 ) {
+
+        return LaplacianDx_Mag(F, coefs_x, i, j, k, ncomp) + LaplacianDy_Mag(F, coefs_y, i, j, k, ncomp) + LaplacianDz_Mag(F, coefs_z, i, j, k, ncomp);
+    }
+
+#endif
+
 };
 
 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM.cpp

@@ -185,9 +185,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian(
                         // H_exchange
                         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_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                     }
 
                     if (mag_anisotropy_coupling == 1){
@@ -298,9 +298,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian(
                         // H_exchange
                         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_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                     }
 
                     if (mag_anisotropy_coupling == 1){
@@ -411,9 +411,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian(
                         // H_exchange
                         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_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                     }
 
                     if (mag_anisotropy_coupling == 1){

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHM_2nd.cpp

@@ -194,9 +194,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd(
                         // H_exchange
                         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_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                     }
 
                     if (mag_anisotropy_coupling == 1){
@@ -275,9 +275,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd(
                         // H_exchange
                         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_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                     }
 
                     if (mag_anisotropy_coupling == 1){
@@ -356,9 +356,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd(
                         // H_exchange
                         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_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                        Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                        Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                        Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                     }
 
                     if (mag_anisotropy_coupling == 1){
@@ -510,9 +510,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd(
                             // H_exchange
                             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_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                            Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                            Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                            Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                            Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                            Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_xface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                         }
 
                         if (mag_anisotropy_coupling == 1){
@@ -620,9 +620,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd(
                             // H_exchange
                             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_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                            Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                            Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                            Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                            Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                            Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_yface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                         }
 
                         if (mag_anisotropy_coupling == 1){
@@ -731,9 +731,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveHMCartesian_2nd(
                             // H_exchange
                             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_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
-                            Hy_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
-                            Hz_eff += H_exchange_coeff * T_Algo::Laplacian(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
+                            Hx_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 0);
+                            Hy_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 1);
+                            Hz_eff += H_exchange_coeff * T_Algo::Laplacian_Mag(M_zface, coefs_x, coefs_y, coefs_z, i, j, k, 2);
                         }
 
                         if (mag_anisotropy_coupling == 1){