Macroscopic pml with mur in H PML

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,15 +1,14 @@
-# This is a E+B (use USE_LLG=FALSE, ie., no M field) simulation of a
+# This is a E+H (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
 geometry.coord_sys   = 0
-geometry.is_periodic = 1  1  0
+geometry.is_periodic = 1  1  1
 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
@@ -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,15 +45,16 @@ 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.
 
 # Diagnostics
 diagnostics.diags_names = plt
 plt.period = 10
 plt.fields_to_plot = Ex Ey Ez Bx By Bz
-plt.diag_type = Full
+plt.diag_type = Full

Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt

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

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H

@@ -148,6 +148,12 @@ class FiniteDifferenceSolver
         void EvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield,
                       std::array< amrex::MultiFab*, 3 > const Efield,
                       amrex::Real const dt );
+
+        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 );
 #endif
 
     private:
@@ -301,9 +307,17 @@ class FiniteDifferenceSolver
 #ifdef WARPX_MAG_LLG
         template< typename T_Algo >
         void EvolveHPMLCartesian (
-            std::array< amrex::MultiFab*, 3 > Bfield,
+            std::array< amrex::MultiFab*, 3 > Hfield,
             std::array< amrex::MultiFab*, 3 > const Efield,
             amrex::Real const dt );
+
+        template< typename T_Algo >
+        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 );
 #endif
 
 #endif

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveHPML.cpp

@@ -0,0 +1,162 @@
+/* Copyright 2020 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
+#ifdef WARPX_DIM_RZ
+#   include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#else
+#   include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
+#   include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
+#   include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H"
+#endif
+#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
+ */
+void FiniteDifferenceSolver::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 ) {
+
+   // 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, macroscopic_properties, mu_mf);
+    amrex::Abort("PML are not implemented in cylindrical geometry.");
+#else
+    if (m_do_nodal) {
+
+        MacroscopicEvolveHPMLCartesian <CartesianNodalAlgorithm> ( Hfield, Efield, dt, macroscopic_properties, mu_mf );
+
+    } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+
+        MacroscopicEvolveHPMLCartesian <CartesianYeeAlgorithm> ( Hfield, Efield, dt, macroscopic_properties, mu_mf );
+
+    } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
+
+        MacroscopicEvolveHPMLCartesian <CartesianCKCAlgorithm> ( Hfield, Efield, dt, macroscopic_properties, mu_mf );
+
+    } else {
+        amrex::Abort("Unknown algorithm");
+    }
+#endif
+}
+
+
+#ifndef WARPX_DIM_RZ
+
+template<typename T_Algo>
+void FiniteDifferenceSolver::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 ) {
+
+    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
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+
+    for ( MFIter mfi(*Hfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+
+        // Extract field data for this grid/tile
+        Array4<Real> const& Hx = Hfield[0]->array(mfi);
+        Array4<Real> const& Hy = Hfield[1]->array(mfi);
+        Array4<Real> const& Hz = Hfield[2]->array(mfi);
+        Array4<Real> const& Ex = Efield[0]->array(mfi);
+        Array4<Real> const& Ey = Efield[1]->array(mfi);
+        Array4<Real> const& Ez = Efield[2]->array(mfi);
+
+        // Extract stencil coefficients
+        Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
+        int const n_coefs_x = m_stencil_coefs_x.size();
+        Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr();
+        int const n_coefs_y = m_stencil_coefs_y.size();
+        Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
+        int const n_coefs_z = m_stencil_coefs_z.size();
+
+        // Extract tileboxes for which to loop
+        Box const& tbx  = mfi.tilebox(Hfield[0]->ixType().ixType());
+        Box const& tby  = mfi.tilebox(Hfield[1]->ixType().ixType());
+        Box const& tbz  = mfi.tilebox(Hfield[2]->ixType().ixType());
+        // 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){
+
+                Real mu_arrx = CoarsenIO::Interp( mu_arr, mu_stag, Hx_stag, macro_cr, i, j, k, scomp);
+
+                Hx(i, j, k, PMLComp::xz) += 1._rt / mu_arrx * 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) -= 1._rt / mu_arrx * 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){
+
+                Real mu_arry = CoarsenIO::Interp( mu_arr, mu_stag, Hy_stag, macro_cr, i, j, k, scomp);
+
+                Hy(i, j, k, PMLComp::yx) += 1._rt / mu_arry * 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) -= 1._rt / mu_arry * 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){
+
+                Real mu_arrz = CoarsenIO::Interp( mu_arr, mu_stag, Hz_stag, macro_cr, i, j, k, scomp);
+
+                Hz(i, j, k, PMLComp::zy) += 1._rt / mu_arrz * 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) -= 1._rt / mu_arrz * 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) );
+            }
+
+        );
+
+    }
+
+}
+
+#endif // corresponds to ifndef WARPX_DIM_RZ
+
+#endif // #ifdef WARPX_MAG_LLG

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/Make.package

@@ -9,6 +9,7 @@ CEXE_sources += MacroscopicEvolveE.cpp
 CEXE_sources += MacroscopicEvolveHM.cpp
 CEXE_sources += MacroscopicEvolveHM_2nd.cpp
 CEXE_sources += EvolveHPML.cpp
+CEXE_sources += MacroscopicEvolveHPML.cpp
 #endif
 
 CEXE_sources += EvolveBPML.cpp

Source/FieldSolver/WarpXPushFieldsEM.cpp

@@ -414,11 +414,19 @@ WarpX::MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real a_dt) {
     // Evolve H field in PML cells
     if (do_pml && pml[lev]->ok()) {
         if (patch_type == PatchType::fine) {
-            m_fdtd_solver_fp[lev]->EvolveHPML(
-                pml[lev]->GetH_fp(), pml[lev]->GetE_fp(), a_dt );
+            m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML(
+                pml[lev]->GetH_fp(), 
+                pml[lev]->GetE_fp(), 
+                a_dt, 
+                m_macroscopic_properties, 
+                pml[lev]->Getmu_fp() );
         } else {
-            m_fdtd_solver_cp[lev]->EvolveHPML(
-                pml[lev]->GetH_cp(), pml[lev]->GetE_cp(), a_dt );
+            m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML(
+                pml[lev]->GetH_cp(), 
+                pml[lev]->GetE_cp(), 
+                a_dt, 
+                m_macroscopic_properties, 
+                pml[lev]->Getmu_cp() );
         }
     }
 }
@@ -457,11 +465,19 @@ WarpX::MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real a_dt)
     // Evolve H field in PML cells
     if (do_pml && pml[lev]->ok()) {
         if (patch_type == PatchType::fine) {
-            m_fdtd_solver_fp[lev]->EvolveHPML(
-                pml[lev]->GetH_fp(), pml[lev]->GetE_fp(), a_dt );
+            m_fdtd_solver_fp[lev]->MacroscopicEvolveHPML(
+                pml[lev]->GetH_fp(), 
+                pml[lev]->GetE_fp(), 
+                a_dt, 
+                m_macroscopic_properties, 
+                pml[lev]->Getmu_fp() );
         } else {
-            m_fdtd_solver_cp[lev]->EvolveHPML(
-                pml[lev]->GetH_cp(), pml[lev]->GetE_cp(), a_dt );
+            m_fdtd_solver_cp[lev]->MacroscopicEvolveHPML(
+                pml[lev]->GetH_cp(), 
+                pml[lev]->GetE_cp(), 
+                a_dt, 
+                m_macroscopic_properties, 
+                pml[lev]->Getmu_cp() );
         }
     }
 }