Macroscopic pml with mur in H PML

Docs/source/running_cpp/parameters.rst

@@ -1216,7 +1216,9 @@ Numerics and algorithms
     The medium for evaluating the Maxwell solver. Available options are :
 
     - ``vacuum``: vacuum properties are used in the Maxwell solver.
-    - ``macroscopic``: macroscopic Maxwell equation is evaluated. If this option is selected, then the corresponding properties of the medium must be provided using ``macroscopic.sigma``, ``macroscopic.epsilon``, and ``macroscopic.mu`` for each case where the initialization style is ``constant``.  Otherwise if the initialization style uses the parser, ``macroscopic.sigma_function(x,y,z)``, ``macroscopic.epsilon_function(x,y,z)`` and/or ``macroscopic.mu_function(x,y,z)`` must be provided using the parser initialization style for spatially varying macroscopic properties.
+    - ``macroscopic``: macroscopic Maxwell equation is evaluated. If this option is selected, then the corresponding properties of the medium must be provided using ``macroscopic.sigma``, 
+    ``macroscopic.epsilon``, and ``macroscopic.mu`` for each case where the initialization style is ``constant``.  Otherwise if the initialization style uses the parser, ``macroscopic.sigma_function(x,y,z)``, ``macroscopic.epsilon_function(x,y,z)`` 
+    and/or ``macroscopic.mu_function(x,y,z)`` must be provided using the parser initialization style for spatially varying macroscopic properties.
 
     If ``algo.em_solver_medium`` is not specified, ``vacuum`` is the default.
 
@@ -1236,7 +1238,10 @@ Numerics and algorithms
 * ``macroscopic.sigma``, ``macroscopic.epsilon``, ``macroscopic.mu`` (`double`)
     To initialize a constant conductivity, permittivity, and permeability of the
     computational medium, respectively. The default values are the corresponding values
-    in vacuum.
+    in vacuum. 
+    In particular, if `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms` is non-zero, ``macroscopic.mu`` must be set to the vacuum value (``mu0``, i.e. 1.25663706212e-06). 
+    If `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms == 0`, LLG evolution is turned off and the macroscopic magnetic permeability, i.e. ``mur * mu0`` should be parsed in ``macroscopic.mu``, with 
+    the dimensionless relative permeability ``mur >=1``.
 
 * ``macroscopic.mag_normalized_error`` (`double`; default: `0.1`)
     The maximum relative amount we let M deviate from Ms before aborting for the LLG equation for saturated cases, i.e., `mag_M_normalization>0`.
@@ -2219,6 +2224,9 @@ Solving magnetization using LLG equation
     ``macroscopic.mag_Ms_function(x,y,z)`` to initialize the saturation magnetization.
     If ``algo.em_solver_medium`` is set to macroscopic, and ``USE_LLG = TRUE``,
     then this input property must be provided.
+    In addition, if `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms` is non-zero, ``macroscopic.mu`` must be set to the vacuum value (``mu0``, i.e. 1.25663706212e-06). 
+    If `USE_LLG=TRUE` in the GNUMakefile and `mag_Ms == 0`, LLG evolution is turned off and the macroscopic magnetic permeability, i.e. ``mur * mu0`` should be parsed in ``macroscopic.mu``, with 
+    the dimensionless relative permeability ``mur >=1``.
 
 * ``macroscopic.mag_alpha_init_style`` (string) optional (default is "default")
     This parameter determines the type of initialization for the Gilbert damping factor

Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_noMs

@@ -1,10 +1,9 @@
 # This is a E+H (use USE_LLG=TRUE) simulation of a
 # plane wave in a periodic domain
-
 ################################
 ####### GENERAL PARAMETERS ######
 #################################
-max_step = 500
+max_step = 1000
 amr.n_cell = 16 16 128
 amr.max_grid_size = 32
 amr.blocking_factor = 16
@@ -21,22 +20,23 @@ warpx.verbose = 0
 warpx.use_filter = 0
 warpx.cfl = 0.9
 warpx.do_pml = 1
-
 algo.em_solver_medium = macroscopic # vacuum/macroscopic
 algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler
 macroscopic.sigma_init_style = "parse_sigma_function" # parse or "constant"
 macroscopic.sigma_function(x,y,z) = "0.0"
 macroscopic.epsilon_init_style = "parse_epsilon_function" # parse or "constant"
 macroscopic.epsilon_function(x,y,z) = "8.8541878128e-12"
+macroscopic.epsilon_function(x,y,z) = "epsilon_r*8.8541878128e-12"
 macroscopic.mu_init_style = "parse_mu_function" # parse or "constant"
 macroscopic.mu_function(x,y,z) = "1.25663706212e-06"
-
+macroscopic.mu_function(x,y,z) = "mu_r*1.25663706212e-06"
 particles.nspecies = 0
 
 #################################
 ############ FIELDS #############
 #################################
-
+my_constants.mu_r = 4
+my_constants.epsilon_r = 1
 my_constants.pi = 3.14159265359
 my_constants.L = 141.4213562373095e-6
 my_constants.c = 299792458.
@@ -45,10 +45,10 @@ my_constants.wavelength = 64.e-6
 warpx.E_ext_grid_init_style = parse_E_ext_grid_function
 warpx.Ez_external_grid_function(x,y,z) = 0.
 warpx.Ex_external_grid_function(x,y,z) = 0.
-warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)"
-
+#warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)"
+warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))"
 warpx.H_ext_grid_init_style = parse_H_ext_grid_function
-warpx.Hx_external_grid_function(x,y,z)= "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)/(c*1.25663706212e-06)"
+warpx.Hx_external_grid_function(x,y,z)= "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r)) / sqrt( mu_r*1.25663706212e-06 / (epsilon_r*8.8541878128e-12) )"
 warpx.Hy_external_grid_function(x,y,z)= 0.
 warpx.Hz_external_grid_function(x,y,z) = 0.
 
@@ -65,4 +65,4 @@ macroscopic.mag_gamma = 0.
 diagnostics.diags_names = plt
 plt.period = 10
 plt.fields_to_plot = Ex Ey Ez Hx Hy Hz Bx By Bz
-plt.diag_type = Full
+plt.diag_type = Full

Examples/Tests/Macroscopic_Maxwell/inputs_3d_noMs

@@ -1,10 +1,9 @@
-# This is a E+B (use USE_LLG=FALSE, ie., no M field) simulation of a
+# This is a E+B (use USE_LLG=FALSE) simulation of a
 # plane wave in a periodic domain
-
 ################################
 ####### GENERAL PARAMETERS ######
 #################################
-max_step = 500
+max_step = 1000
 amr.n_cell = 16 16 128
 amr.max_grid_size = 32
 amr.blocking_factor = 16
@@ -21,22 +20,23 @@ warpx.verbose = 0
 warpx.use_filter = 0
 warpx.cfl = 0.9
 warpx.do_pml = 1
-
 algo.em_solver_medium = macroscopic # vacuum/macroscopic
 algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler
 macroscopic.sigma_init_style = "parse_sigma_function" # parse or "constant"
 macroscopic.sigma_function(x,y,z) = "0.0"
 macroscopic.epsilon_init_style = "parse_epsilon_function" # parse or "constant"
 macroscopic.epsilon_function(x,y,z) = "8.8541878128e-12"
+macroscopic.epsilon_function(x,y,z) = "epsilon_r*8.8541878128e-12"
 macroscopic.mu_init_style = "parse_mu_function" # parse or "constant"
 macroscopic.mu_function(x,y,z) = "1.25663706212e-06"
-
+macroscopic.mu_function(x,y,z) = "mu_r*1.25663706212e-06"
 particles.nspecies = 0
 
 #################################
 ############ FIELDS #############
 #################################
-
+my_constants.mu_r = 1
+my_constants.epsilon_r = 1
 my_constants.pi = 3.14159265359
 my_constants.L = 141.4213562373095e-6
 my_constants.c = 299792458.
@@ -45,10 +45,11 @@ my_constants.wavelength = 64.e-6
 warpx.E_ext_grid_init_style = parse_E_ext_grid_function
 warpx.Ez_external_grid_function(x,y,z) = 0.
 warpx.Ex_external_grid_function(x,y,z) = 0.
-warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)"
+#warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)"
+warpx.Ey_external_grid_function(x,y,z) = "1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))"
 
 warpx.B_ext_grid_init_style = parse_B_ext_grid_function
-warpx.Bx_external_grid_function(x,y,z)= "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)/c"
+warpx.Bx_external_grid_function(x,y,z)= "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))/c * sqrt(epsilon_r*mu_r)"
 warpx.By_external_grid_function(x,y,z)= 0.
 warpx.Bz_external_grid_function(x,y,z) = 0.
 

Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt

@@ -19,7 +19,7 @@ if(WarpX_MAG_LLG)
         MacroscopicEvolveHM.cpp
         MacroscopicEvolveM_2nd.cpp
         MacroscopicEvolveM.cpp
-        EvolveHPML.cpp
+        MacroscopicEvolveHPML.cpp
     )
 endif()
 add_subdirectory(MacroscopicProperties)

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H

@@ -145,9 +145,11 @@ class FiniteDifferenceSolver
 
 
 #ifdef WARPX_MAG_LLG
-        void EvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield,
+        void MacroscopicEvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield,
                       std::array< amrex::MultiFab*, 3 > const Efield,
-                      amrex::Real const dt );
+                      amrex::Real const dt,
+                      std::unique_ptr<MacroscopicProperties> const& macroscopic_properties,
+                      amrex::MultiFab* const mu_mf );
 #endif
 
     private:
@@ -300,10 +302,12 @@ class FiniteDifferenceSolver
 
 #ifdef WARPX_MAG_LLG
         template< typename T_Algo >
-        void EvolveHPMLCartesian (
-            std::array< amrex::MultiFab*, 3 > Bfield,
+        void MacroscopicEvolveHPMLCartesian (
+            std::array< amrex::MultiFab*, 3 > Hfield,
             std::array< amrex::MultiFab*, 3 > const Efield,
-            amrex::Real const dt );
+            amrex::Real const dt,
+            std::unique_ptr<MacroscopicProperties> const& macroscopic_properties,
+            amrex::MultiFab* const mu_mf );
 #endif
 
 #endif

Source/FieldSolver/FiniteDifferenceSolver/EvolveHPML.cpp → Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp

@@ -1,9 +1,3 @@
-/* Copyright 2020 Remi Lehe
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
 
 #include "Utils/WarpXAlgorithmSelection.H"
 #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
@@ -17,36 +11,39 @@
 #include "BoundaryConditions/PMLComponent.H"
 #include <AMReX_Gpu.H>
 #include <AMReX.H>
+#include "Utils/CoarsenIO.H"
 
 using namespace amrex;
 
 #ifdef WARPX_MAG_LLG
 
 /**
- * \brief Update the H field, over one timestep
+ * \brief Update Hfield in PML region
  */
-void FiniteDifferenceSolver::EvolveHPML (
+void FiniteDifferenceSolver::MacroscopicEvolveHPML (
     std::array< amrex::MultiFab*, 3 > Hfield,
     std::array< amrex::MultiFab*, 3 > const Efield,
-    amrex::Real const dt ) {
+    amrex::Real const dt,
+    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties,
+    amrex::MultiFab* const mu_mf ) {
 
    // Select algorithm (The choice of algorithm is a runtime option,
    // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
-    amrex::ignore_unused(Hfield, Efield, dt);
+    amrex::ignore_unused(Hfield, Efield, dt, macroscopic_properties, mu_mf);
     amrex::Abort("PML are not implemented in cylindrical geometry.");
 #else
     if (m_do_nodal) {
 
-        EvolveHPMLCartesian <CartesianNodalAlgorithm> ( Hfield, Efield, dt );
+        MacroscopicEvolveHPMLCartesian <CartesianNodalAlgorithm> ( Hfield, Efield, dt, macroscopic_properties, mu_mf );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
 
-        EvolveHPMLCartesian <CartesianYeeAlgorithm> ( Hfield, Efield, dt );
+        MacroscopicEvolveHPMLCartesian <CartesianYeeAlgorithm> ( Hfield, Efield, dt, macroscopic_properties, mu_mf );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
-        EvolveHPMLCartesian <CartesianCKCAlgorithm> ( Hfield, Efield, dt );
+        MacroscopicEvolveHPMLCartesian <CartesianCKCAlgorithm> ( Hfield, Efield, dt, macroscopic_properties, mu_mf );
 
     } else {
         amrex::Abort("Unknown algorithm");
@@ -58,10 +55,18 @@ void FiniteDifferenceSolver::EvolveHPML (
 #ifndef WARPX_DIM_RZ
 
 template<typename T_Algo>
-void FiniteDifferenceSolver::EvolveHPMLCartesian (
+void FiniteDifferenceSolver::MacroscopicEvolveHPMLCartesian (
     std::array< amrex::MultiFab*, 3 > Hfield,
     std::array< amrex::MultiFab*, 3 > const Efield,
-    amrex::Real const dt ) {
+    amrex::Real const dt,
+    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties,
+    amrex::MultiFab* const mu_mf ) {
+
+    amrex::GpuArray<int, 3> const& mu_stag  = macroscopic_properties->mu_IndexType;
+    amrex::GpuArray<int, 3> const& Hx_stag  = macroscopic_properties->Hx_IndexType;
+    amrex::GpuArray<int, 3> const& Hy_stag  = macroscopic_properties->Hy_IndexType;
+    amrex::GpuArray<int, 3> const& Hz_stag  = macroscopic_properties->Hz_IndexType;
+    amrex::GpuArray<int, 3> const& macro_cr = macroscopic_properties->macro_cr_ratio;
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef _OPENMP
@@ -90,39 +95,51 @@ void FiniteDifferenceSolver::EvolveHPMLCartesian (
         Box const& tby  = mfi.tilebox(Hfield[1]->ixType().ixType());
         Box const& tbz  = mfi.tilebox(Hfield[2]->ixType().ixType());
 
-        amrex::Real mu0_inv = 1._rt/PhysConst::mu0;
+        // starting component to interpolate macro properties to Hx, Hy, Hz locations
+        const int scomp = 0;
+        // mu_mf will be imported but will only be called at grids where Ms == 0
+        Array4<Real> const& mu_arr = mu_mf->array(mfi);
 
         // Loop over the cells and update the fields
         amrex::ParallelFor(tbx, tby, tbz,
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                Hx(i, j, k, PMLComp::xz) += mu0_inv * dt * (
+
+                Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, scomp);
+
+                Hx(i, j, k, PMLComp::xz) += mu_inv * dt * (
                     T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yx)
                   + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yy)
                   + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) );
-                Hx(i, j, k, PMLComp::xy) -= mu0_inv * dt * (
+                Hx(i, j, k, PMLComp::xy) -= mu_inv * dt * (
                     T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zx)
                   + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zy)
                   + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz) );
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                Hy(i, j, k, PMLComp::yx) += mu0_inv * dt * (
+
+                Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, scomp);
+
+                Hy(i, j, k, PMLComp::yx) += mu_inv * dt * (
                     T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zx)
                   + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zy)
                   + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz) );
-                Hy(i, j, k, PMLComp::yz) -= mu0_inv * dt * (
+                Hy(i, j, k, PMLComp::yz) -= mu_inv * dt * (
                     T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xx)
                   + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xy)
                   + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xz) );
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                Hz(i, j, k, PMLComp::zy) += mu0_inv * dt * (
+
+                Real mu_inv = 1._rt/CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, scomp);
+
+                Hz(i, j, k, PMLComp::zy) += mu_inv * dt * (
                     T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xx)
                   + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xy)
                   + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) );
-                Hz(i, j, k, PMLComp::zx) -= mu0_inv * dt * (
+                Hz(i, j, k, PMLComp::zx) -= mu_inv * dt * (
                     T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yx)
                   + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yy)
                   + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) );

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H

@@ -58,6 +58,12 @@ public:
      amrex::GpuArray<int, 3> My_IndexType;
      /** Gpu Vector with index type of the Mz multifab */
      amrex::GpuArray<int, 3> Mz_IndexType;
+     /** Gpu Vector with index type of the Hx multifab */
+     amrex::GpuArray<int, 3> Hx_IndexType;
+     /** Gpu Vector with index type of the Hy multifab */
+     amrex::GpuArray<int, 3> Hy_IndexType;
+     /** Gpu Vector with index type of the Hz multifab */
+     amrex::GpuArray<int, 3> Hz_IndexType;
 #endif
      /** Gpu Vector with index type of coarsening ratio with default value (1,1,1) */
      amrex::GpuArray<int, 3> macro_cr_ratio;

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp

@@ -235,6 +235,9 @@ MacroscopicProperties::InitData ()
     IntVect Mx_stag = warpx.getMfield_fp(0,0).ixType().toIntVect(); // face-centered
     IntVect My_stag = warpx.getMfield_fp(0,1).ixType().toIntVect();
     IntVect Mz_stag = warpx.getMfield_fp(0,2).ixType().toIntVect();
+    IntVect Hx_stag = warpx.getHfield_fp(0,0).ixType().toIntVect(); // face-centered
+    IntVect Hy_stag = warpx.getHfield_fp(0,1).ixType().toIntVect();
+    IntVect Hz_stag = warpx.getHfield_fp(0,2).ixType().toIntVect();
 #endif
     for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
         sigma_IndexType[idim]   = sigma_stag[idim];
@@ -250,6 +253,9 @@ MacroscopicProperties::InitData ()
         Mx_IndexType[idim]        = Mx_stag[idim];
         My_IndexType[idim]        = My_stag[idim];
         Mz_IndexType[idim]        = Mz_stag[idim];
+        Hx_IndexType[idim]        = Hx_stag[idim];
+        Hy_IndexType[idim]        = Hy_stag[idim];
+        Hz_IndexType[idim]        = Hz_stag[idim];
 #endif
         macro_cr_ratio[idim]    = 1;
     }
@@ -267,6 +273,9 @@ MacroscopicProperties::InitData ()
         Mx_IndexType[2]        = 0;
         My_IndexType[2]        = 0;
         Mz_IndexType[2]        = 0;
+        Hx_IndexType[2]        = 0;
+        Hy_IndexType[2]        = 0;
+        Hz_IndexType[2]        = 0;
 #endif
         macro_cr_ratio[2]    = 1;
 #endif

Source/FieldSolver/FiniteDifferenceSolver/Make.package

@@ -8,7 +8,7 @@ CEXE_sources += MacroscopicEvolveE.cpp
 #ifdef WARPX_MAG_LLG
 CEXE_sources += MacroscopicEvolveHM.cpp
 CEXE_sources += MacroscopicEvolveHM_2nd.cpp
-CEXE_sources += EvolveHPML.cpp
+CEXE_sources += MacroscopicEvolveHPML.cpp
 #endif
 
 CEXE_sources += EvolveBPML.cpp