Macroscopic pml with mur in H PML

Docs/source/running_cpp/parameters.rst

@@ -1265,7 +1265,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.
 
@@ -1285,7 +1287,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`.
@@ -2300,6 +2305,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_MsBoundary

@@ -0,0 +1,73 @@
+# This is a E+H (use USE_LLG=TRUE) simulation of a
+# plane wave in a periodic domain
+################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1000
+amr.n_cell = 16 16 128
+amr.max_grid_size = 32
+amr.blocking_factor = 16
+geometry.coord_sys   = 0
+geometry.is_periodic = 1  1  0
+geometry.prob_lo     =  -16.e-6 -16.e-6 -128.e-6
+geometry.prob_hi     =   16.e-6  16.e-6  128.e-6
+amr.max_level = 0
+
+#################################
+############ NUMERICS ###########
+#################################
+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) = "epsilon_r*8.8541878128e-12"
+macroscopic.mu_init_style = "parse_mu_function" # parse or "constant"
+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.
+my_constants.wavelength = 64.e-6
+
+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) = "1.e7*exp(-z**2/L**2)*cos(2*pi*z/wavelength * sqrt(epsilon_r*mu_r))"
+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) = "-1.e7*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.
+
+# If you want to use a USE_LLG build with Ms=0 everywhere, uncomment these settings
+warpx.mag_M_normalization = 1
+macroscopic.mag_Ms_init_style = parse_mag_Ms_function
+macroscopic.mag_Ms_function(x,y,z) = "1.4e5 * (z>64e-6)"
+
+macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant"
+macroscopic.mag_alpha_function(x,y,z) = "0.0058 * (z>64e-6)" # 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 * (z>64e-6)" # gyromagnetic ratio is constant for electrons in all materials
+
+warpx.M_ext_grid_init_style = parse_M_ext_grid_function
+warpx.Mx_external_grid_function(x,y,z)= 0.
+warpx.My_external_grid_function(x,y,z)= "1.4e5 * (z>64e-6)"
+warpx.Mz_external_grid_function(x,y,z) = 0.
+
+# Diagnostics
+diagnostics.diags_names = plt
+plt.period = 10
+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.diag_type = Full

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,35 +20,34 @@ 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.
 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 * sqrt(epsilon_r*mu_r))"
+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)= "-1.e5*exp(-z**2/L**2)*cos(2*pi*z/wavelength)/(c*1.25663706212e-06)"
-warpx.Hy_external_grid_function(x,y,z)= 0.
+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.
 
 # If you want to use a USE_LLG build with Ms=0 everywhere, uncomment these settings

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,35 +20,34 @@ 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.
 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 * sqrt(epsilon_r*mu_r))"
+warpx.Ez_external_grid_function(x,y,z) = 0.
 
 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.By_external_grid_function(x,y,z)= 0.
+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.
 
 # Diagnostics

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

@@ -146,9 +146,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,
+                      std::unique_ptr<MacroscopicProperties> const& macroscopic_properties,
+                      amrex::MultiFab* const mu_mf,
                       const bool dive_cleaning);
 #endif
 
@@ -303,10 +305,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,
+            std::unique_ptr<MacroscopicProperties> const& macroscopic_properties,
+            amrex::MultiFab* const mu_mf,
             const bool dive_cleaning);
 #endif