two layer implementation from GDS

Examples/Tests/circuits/Candice/inputs_candice20

@@ -0,0 +1,152 @@
+# See parameters in https://github.com/ECP-WarpX/artemis/pull/43
+
+################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+
+amr.n_cell = n_cellx n_celly n_cellz
+amr.max_grid_size = max_grid_sizex max_grid_sizey max_grid_sizez
+amr.blocking_factor = blocking_factor
+amr.refine_grid_layout = 1  # if n_MPI > n_grids, the grids will be successively divided in half until n_MPI <= n_grids
+
+geometry.dims = 3
+geometry.prob_lo = prob_lox prob_loy prob_loz
+geometry.prob_hi = prob_hix prob_hiy prob_hiz
+
+amr.max_level = 0
+
+# use pec instead of pml overlaying current source so you don't get a reflection
+boundary.field_lo = pml pml pml
+boundary.field_hi = pml pml pml
+
+my_constants.offset_y = -1100.e-6
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.verbose = 1
+
+warpx.cfl = 0.8
+
+# vacuum or macroscopic
+algo.em_solver_medium = macroscopic
+
+# laxwendroff or backwardeuler
+algo.macroscopic_sigma_method = laxwendroff
+
+macroscopic.sigma_function(x,y,z) = "sigma_0 + (sigma_si - sigma_0) * (z <= h_si)"
+macroscopic.epsilon_function(x,y,z) = "eps_0 + eps_0 * (eps_r_si - 1.) * (z <= h_si)"
+macroscopic.mu_function(x,y,z) = "mu_0 + mu_0 * (mu_r_si - 1.) * (z <= h_si)"
+
+macroscopic.npy_k_index = 0
+macroscopic.npy_k_index2 = 3
+
+# comment for PEC, uncomment these for GDS-PEC
+macroscopic.sigma_npy_file = "voyageurs_layer1.npy"
+macroscopic.sigma_npy_file2 = "voyageurs_layer2.npy"
+macroscopic.sigma_npy_value = 1.e11
+algo.use_PEC_mask = 1
+
+macroscopic.mu_npy_file = "voyageurs_layer1.npy"
+macroscopic.mu_npy_file2 = "voyageurs_layer2.npy"
+macroscopic.mu_npy_value = 1.25663706e-06
+
+#################################
+############ FIELDS #############
+#################################
+
+###############
+# domain size
+# n_cellx/y/z and Lx/y/z are needed to calculate dx/dy/dz
+###############
+my_constants.n_cellx = 7544
+my_constants.n_celly = 7544
+my_constants.n_cellz = 4
+
+my_constants.max_grid_sizex = 92
+my_constants.max_grid_sizey = 92
+my_constants.max_grid_sizez = 4
+my_constants.blocking_factor = 2
+
+my_constants.prob_lox = 0.
+my_constants.prob_loy = 0.
+my_constants.prob_loz = 0.
+
+my_constants.prob_hix = 9430.e-6
+my_constants.prob_hiy = 9430.e-6
+my_constants.prob_hiz = 300.e-6
+
+my_constants.Lx = prob_hix - prob_lox
+my_constants.Ly = prob_hiy - prob_loy
+my_constants.Lz = prob_hiz - prob_loz
+
+###############
+# material properties
+###############
+my_constants.sigma_0 = 0.0
+my_constants.sigma_si = 0.
+
+my_constants.eps_0 = 8.8541878128e-12
+my_constants.eps_r_si = 11.7
+
+my_constants.mu_0 = 1.25663706212e-06
+my_constants.mu_r_si = 1.0
+
+###############
+# silicon and palladium cross section
+###############
+my_constants.h_si = 150.e-6
+
+###############
+# derived quantities and fundamental constants - don't touch these
+###############
+
+my_constants.pi = 3.14159265358979
+
+my_constants.freq = 8.6e9
+my_constants.TP = 1./freq
+
+# grid spacing
+my_constants.dx = Lx / n_cellx
+my_constants.dy = Ly / n_celly
+my_constants.dz = Lz / n_cellz
+
+my_constants.ddx = dx/1.e6
+my_constants.ddy = dy/1.e6
+my_constants.ddz = dz/1.e6
+
+my_constants.flag_none = 0
+my_constants.flag_hs = 1
+my_constants.flag_ss = 2
+
+###############
+# excitation
+###############
+
+warpx.E_excitation_on_grid_style = parse_E_excitation_grid_function
+
+warpx.Ex_excitation_flag_function(x,y,z) = "flag_none"
+warpx.Ey_excitation_flag_function(x,y,z) = "flag_ss * 0"
+warpx.Ez_excitation_flag_function(x,y,z) = "flag_none"
+
+
+# This is a source on a qubit control line
+warpx.Ex_excitation_grid_function(x,y,z,t) = "0."
+warpx.Ey_excitation_grid_function(x,y,z,t) = "0."
+warpx.Ez_excitation_grid_function(x,y,z,t) = "0."
+
+###############
+# diagnostics
+###############
+
+warpx.field_io_nfiles=64
+warpx.particle_io_nfiles=64
+
+diagnostics.diags_names = plt
+###############
+# full plotfiles
+plt.intervals = 1
+plt.fields_to_plot = Ex Ey Ez Bx By Bz mu
+plt.diag_type = Full
+plt.file_min_digits = 7

Source/BoundaryConditions/PML.cpp

@@ -723,6 +723,13 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
                                                                      lev,
                                                                      macroscopic_properties->m_npy_k_index,
                                                                      macroscopic_properties->m_sigma_npy_value);
+            if (!macroscopic_properties->m_sigma_npy_filename2.empty()) {
+                macroscopic_properties->InitializeMacroMultiFabFromNumpy(pml_sigma_fp.get(),
+                                                                         macroscopic_properties->m_sigma_npy_filename2,
+                                                                         lev,
+                                                                         macroscopic_properties->m_npy_k_index2,
+                                                                         macroscopic_properties->m_sigma_npy_value);
+            }
         }
 
         if (warpx.use_PEC_mask) {
@@ -755,6 +762,13 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
                                                                      lev,
                                                                      macroscopic_properties->m_npy_k_index,
                                                                      macroscopic_properties->m_epsilon_npy_value);
+            if (!macroscopic_properties->m_epsilon_npy_filename2.empty()) {
+                macroscopic_properties->InitializeMacroMultiFabFromNumpy(pml_eps_fp.get(),
+                                                                         macroscopic_properties->m_epsilon_npy_filename2,
+                                                                         lev,
+                                                                         macroscopic_properties->m_npy_k_index2,
+                                                                         macroscopic_properties->m_epsilon_npy_value);
+            }
         }
 
         // Initialize mu, permeability
@@ -775,6 +789,13 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
                                                                      lev,
                                                                      macroscopic_properties->m_npy_k_index,
                                                                      macroscopic_properties->m_mu_npy_value);
+            if (!macroscopic_properties->m_mu_npy_filename2.empty()) {
+                macroscopic_properties->InitializeMacroMultiFabFromNumpy(pml_mu_fp.get(),
+                                                                         macroscopic_properties->m_mu_npy_filename2,
+                                                                         lev,
+                                                                         macroscopic_properties->m_npy_k_index2,
+                                                                         macroscopic_properties->m_mu_npy_value);
+            }
         }
 
     }

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H

@@ -100,10 +100,14 @@ public:
      amrex::Real m_mu = PhysConst::mu0;
      /** index to insert the mask */
      int m_npy_k_index = 0;
+     int m_npy_k_index2 = 0;
      /** numpy arrays for layer */
      std::string m_sigma_npy_filename;
+     std::string m_sigma_npy_filename2;
      std::string m_epsilon_npy_filename;
+     std::string m_epsilon_npy_filename2;
      std::string m_mu_npy_filename;
+     std::string m_mu_npy_filename2;
      /** numerical value of sigma, epsilon, and mu for gds numpy array layer */
      amrex::Real m_sigma_npy_value = 0.;
      amrex::Real m_epsilon_npy_value = 0.;

Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp

@@ -44,6 +44,7 @@ MacroscopicProperties::ReadParameters ()
 
     // Query mask index
     pp_macroscopic.query("npy_k_index", m_npy_k_index);
+    pp_macroscopic.query("npy_k_index2", m_npy_k_index2);
 
     // Query input for material conductivity, sigma.
     pp_macroscopic.query("sigma_npy_value", m_sigma_npy_value);
@@ -63,6 +64,8 @@ MacroscopicProperties::ReadParameters ()
         m_sigma_s = "parse_sigma_npy_file";
         sigma_specified = true;
         sigma_npy_specified = true;
+        // parse second layer
+        pp_macroscopic.query("sigma_npy_file2", m_sigma_npy_filename2);
     }
     if (sigma_func_specified && sigma_npy_specified) {
          // initialize both later in InitData
@@ -101,6 +104,7 @@ MacroscopicProperties::ReadParameters ()
         m_epsilon_s = "parse_epsilon_npy_file";
         epsilon_specified = true;
         epsilon_npy_specified = true;
+        pp_macroscopic.query("epsilon_npy_file2", m_epsilon_npy_filename2);
     }
     if (epsilon_func_specified && epsilon_npy_specified) {
          // initialize both later in InitData
@@ -139,6 +143,7 @@ MacroscopicProperties::ReadParameters ()
         m_mu_s = "parse_mu_npy_file";
         mu_specified = true;
         mu_npy_specified = true;
+        pp_macroscopic.query("mu_npy_file2", m_mu_npy_filename2);
     }
     if (mu_func_specified && mu_npy_specified) {
          // initialize both later in InitData
@@ -269,6 +274,9 @@ MacroscopicProperties::InitData ()
     // Step 2: Overwrite with numpy mask in valid region if provided
     if (m_sigma_s == "parse_sigma_npy_file" || m_sigma_s == "parse_sigma_both") {
         InitializeMacroMultiFabFromNumpy(m_sigma_mf.get(), m_sigma_npy_filename, lev, m_npy_k_index, m_sigma_npy_value);
+        if (!m_sigma_npy_filename2.empty()) {
+            InitializeMacroMultiFabFromNumpy(m_sigma_mf.get(), m_sigma_npy_filename2, lev, m_npy_k_index2, m_sigma_npy_value);
+        }
     }
 
     if (warpx.use_PEC_mask) {
@@ -301,6 +309,9 @@ MacroscopicProperties::InitData ()
     // Step 2: Overwrite with numpy mask in valid region if provided
     if (m_epsilon_s == "parse_epsilon_npy_file" || m_epsilon_s == "parse_epsilon_both") {
         InitializeMacroMultiFabFromNumpy(m_eps_mf.get(), m_epsilon_npy_filename, lev, m_npy_k_index, m_epsilon_npy_value);
+        if (!m_epsilon_npy_filename2.empty()) {
+            InitializeMacroMultiFabFromNumpy(m_eps_mf.get(), m_epsilon_npy_filename2, lev, m_npy_k_index2, m_epsilon_npy_value);
+        }
     }
 
     ////////////////////////
@@ -318,6 +329,9 @@ MacroscopicProperties::InitData ()
     // Step 2: Overwrite with numpy mask in valid region if provided
     if (m_mu_s == "parse_mu_npy_file" || m_mu_s == "parse_mu_both") {
         InitializeMacroMultiFabFromNumpy(m_mu_mf.get(), m_mu_npy_filename, lev, m_npy_k_index, m_mu_npy_value);
+        if (!m_mu_npy_filename2.empty()) {
+            InitializeMacroMultiFabFromNumpy(m_mu_mf.get(), m_mu_npy_filename2, lev, m_npy_k_index2, m_mu_npy_value);
+        }
     }
 
 #ifdef WARPX_MAG_LLG