Adding support for pml in domain for multiple patches (ECP-WarpX#4632)

* this works only on single proc, but lifts the assumption of single patch for do pml in domain

* fix cbl

* building reduced ba and cba

* Removing print statements

* removing comment

* assert that box length must be greater than size of pml when its grown inwards

* same assertion for box length for coarse grid

* simplified_list and no need for assert in coarse, but if fine works, coarse should

* boxlist instead of simplified list

* comments

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* CI test for double patch with particles entering and leaving pml inside the patch

* checksum for new particles in pml 2d mr test with double patch

* renaming plt to diag1 again

* include level 0 maxfield to capture no j damping effect

* removing incorect comment about heavy paritcle in input

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py

@@ -0,0 +1,68 @@
+#!/usr/bin/env python3
+
+# Copyright 2019-2020 Luca Fedeli, Maxence Thevenet, Remi Lehe
+#
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+"""
+This script tests the absorption of particles in the PML.
+
+The input file inputs_2d/inputs is used: it features a positive and a
+negative particle, going in opposite direction and eventually
+leaving the box. This script tests that the field in the box
+is close to 0 once the particles have left. With regular
+PML, this test fails, since the particles leave a spurious
+charge, with associated fields, behind them.
+"""
+import os
+import sys
+
+import yt
+
+yt.funcs.mylog.setLevel(0)
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# Open plotfile specified in command line
+filename = sys.argv[1]
+ds = yt.load( filename )
+
+# Check that the field is low enough
+ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+Ex_array_lev0 = ad0[('mesh','Ex')].to_ndarray()
+Ey_array_lev0 = ad0[('mesh','Ey')].to_ndarray()
+Ez_array_lev0 = ad0[('mesh','Ez')].to_ndarray()
+max_Efield_lev0 = max(Ex_array_lev0.max(), Ey_array_lev0.max(), Ez_array_lev0.max())
+print( "max_Efield level0 = %s" %max_Efield_lev0 )
+
+ad1 = ds.covering_grid(level=1, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+Ex_array_lev1 = ad1[('mesh','Ex')].to_ndarray()
+Ey_array_lev1 = ad1[('mesh','Ey')].to_ndarray()
+Ez_array_lev1 = ad1[('mesh','Ez')].to_ndarray()
+max_Efield_lev1 = max(Ex_array_lev1.max(), Ey_array_lev1.max(), Ez_array_lev1.max())
+print( "max_Efield level1 = %s" %max_Efield_lev1 )
+
+# The field associated with the particle does not have
+# the same amplitude in 2d and 3d
+if ds.dimensionality == 2:
+    if ds.max_level == 0:
+        tolerance_abs = 0.0003
+    elif ds.max_level == 1:
+        tolerance_abs = 0.0006
+elif ds.dimensionality == 3:
+    if ds.max_level == 0:
+        tolerance_abs = 10
+    elif ds.max_level == 1:
+        tolerance_abs = 110
+else:
+    raise ValueError("Unknown dimensionality")
+
+print("tolerance_abs: " + str(tolerance_abs))
+assert max_Efield_lev0 < tolerance_abs
+assert max_Efield_lev1 < tolerance_abs
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename)

Examples/Tests/particles_in_pml/inputs_mr_2d

@@ -1,16 +1,15 @@
-max_step = 200
+max_step = 300
 amr.n_cell = 128 128
 
 amr.blocking_factor = 16
-amr.max_grid_size = 1024
+amr.max_grid_size = 64
 amr.max_level = 1
 
 # Geometry
 geometry.dims = 2
 geometry.prob_lo     = -32.e-6    -32.e-6      # physical domain
 geometry.prob_hi     =  32.e-6     32.e-6
-warpx.fine_tag_lo     = -8.e-6   -8.e-6      # physical domain
-warpx.fine_tag_hi     =  8.e-6    8.e-6
+warpx.ref_patch_function(x,y,z) = " ((x > -16.e-6)*(x<-10.e-6) + (x > 9.e-6)*(x<16.e-6))*(z>-8.e-6)*(z<8.e-6)"
 
 # Boundary condition
 boundary.field_lo = pml pml pml
@@ -38,14 +37,14 @@ particles.species_names = electron proton
 electron.charge = -q_e
 electron.mass = m_e
 electron.injection_style = "singleparticle"
-electron.single_particle_pos = 0. 0. 0.
+electron.single_particle_pos = -12.e-6 0. 0.
 electron.single_particle_u = 2. 0. 0.
 electron.single_particle_weight = 1.
 
 proton.charge = q_e
-proton.mass = m_p    # Very heavy ; should not move
+proton.mass = m_p
 proton.injection_style = "singleparticle"
-proton.single_particle_pos = 0. 0. 0.
+proton.single_particle_pos = -12.e-6 0. 0.
 proton.single_particle_u = -2. 0. 0.
 proton.single_particle_weight = 1.
 
@@ -54,6 +53,6 @@ algo.particle_shape = 1
 
 # Diagnostics
 diagnostics.diags_names = diag1
-diag1.intervals = 200
+diag1.intervals = 300
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz

Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json

@@ -1,22 +1,22 @@
 {
   "lev=0": {
     "Bx": 0.0,
-    "By": 3.578051588629885e-09,
+    "By": 1.959394057951299e-09,
     "Bz": 0.0,
-    "Ex": 1.9699822913484977,
+    "Ex": 1.2177330062543195,
     "Ey": 0.0,
-    "Ez": 0.5356004212513483,
+    "Ez": 0.28770171464461436,
     "jx": 0.0,
     "jy": 0.0,
     "jz": 0.0
   },
   "lev=1": {
     "Bx": 0.0,
-    "By": 2.7629151306453947e-09,
+    "By": 1.4999363537195907e-09,
     "Bz": 0.0,
-    "Ex": 2.4597363065789337,
+    "Ex": 1.5591272748392493,
     "Ey": 0.0,
-    "Ez": 0.45901250290130735,
+    "Ez": 0.29412053281248907,
     "jx": 0.0,
     "jy": 0.0,
     "jz": 0.0

Regression/WarpX-tests.ini

@@ -2525,7 +2525,7 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 0
-analysisRoutine = Examples/Tests/particles_in_pml/analysis_particles_in_pml.py
+analysisRoutine = Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py
 
 [particles_in_pml_3d_MR]
 buildDir = .

Source/BoundaryConditions/PML.cpp

@@ -560,40 +560,56 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
       m_geom(geom),
       m_cgeom(cgeom)
 {
-    // When `do_pml_in_domain` is true, the PML overlap with the last `ncell` of the physical domain
-    // (instead of extending `ncell` outside of the physical domain)
-    // In order to implement this, a reduced domain is created here (decreased by ncells in all direction)
-    // and passed to `MakeBoxArray`, which surrounds it by PML boxes
-    // (thus creating the PML boxes at the right position, where they overlap with the original domain)
-    // minimalBox provides the bounding box around grid_ba for level, lev.
-    // Note that this is okay to build pml inside domain for a single patch, or joint patches
-    // with same [min,max]. But it does not support multiple disjoint refinement patches.
-    Box domain0 = grid_ba.minimalBox();
+    // When `do_pml_in_domain` is true, the PML overlap with the last `ncell` of the physical domain or fine patch(es)
+    // (instead of extending `ncell` outside of the physical domain or fine patch(es))
+    // In order to implement this, we define a new reduced Box Array ensuring that it does not
+    // include ncells from the edges of the physical domain or fine patch.
+    // (thus creating the PML boxes at the right position, where they overlap with the original domain or fine patch(es))
+
+    BoxArray grid_ba_reduced = grid_ba;
     if (do_pml_in_domain) {
-        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-            if (do_pml_Lo[idim]){
-                domain0.growLo(idim, -ncell);
-            }
-            if (do_pml_Hi[idim]){
-                domain0.growHi(idim, -ncell);
+        BoxList bl = grid_ba.boxList();
+        // Here we loop over all the boxes in the original grid_ba BoxArray
+        // For each box, we find if its in the edge (or boundary), and the size of those boxes are decreased by ncell
+        for (auto& b : bl) {
+            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                if (do_pml_Lo[idim]) {
+                    // Get neighboring box on lower side in direction idim and check if it intersects with any of the boxes
+                    // in grid_ba. If no intersection, then the box, b, in the boxlist, is in the edge and we decrase
+                    // the size by ncells using growLo(idim,-ncell)
+                    Box const& bb = amrex::adjCellLo(b, idim);
+                    if ( ! grid_ba.intersects(bb) ) {
+                        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(b.length(idim) > ncell, " box length must be greater that pml size");
+                        b.growLo(idim, -ncell);
+                    }
+                }
+                if (do_pml_Hi[idim]) {
+                    // Get neighboring box on higher side in direction idim and check if it intersects with any of the boxes
+                    // in grid_ba. If no intersection, then the box, b, in the boxlist, is in the edge and we decrase
+                    // the size by ncells using growHi(idim,-ncell)
+                    Box const& bb = amrex::adjCellHi(b, idim);
+                    if ( ! grid_ba.intersects(bb) ) {
+                        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(b.length(idim) > ncell, " box length must be greater that pml size");
+                        b.growHi(idim, -ncell);
+                    }
+                }
             }
         }
+        grid_ba_reduced = BoxArray(std::move(bl));
     }
-    const BoxArray grid_ba_reduced = (do_pml_in_domain) ?
-        BoxArray(grid_ba.boxList().intersect(domain0)) : grid_ba;
-
+    Box const domain0 = grid_ba_reduced.minimalBox();
     const bool is_single_box_domain = domain0.numPts() == grid_ba_reduced.numPts();
     const BoxArray& ba = MakeBoxArray(is_single_box_domain, domain0, *geom, grid_ba_reduced,
                                       IntVect(ncell), do_pml_in_domain, do_pml_Lo, do_pml_Hi);
 
+
     if (ba.empty()) {
         m_ok = false;
         return;
     } else {
         m_ok = true;
     }
-
-    // Define the number of guard cells in each direction, for E, B, and F
+    // Define the number of guard cells in each di;rection, for E, B, and F
     auto nge = IntVect(AMREX_D_DECL(2, 2, 2));
     auto ngb = IntVect(AMREX_D_DECL(2, 2, 2));
     int ngf_int = 0;
@@ -760,24 +776,28 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         BoxArray grid_cba = grid_ba;
         grid_cba.coarsen(ref_ratio);
 
-        // assuming that the bounding box around grid_cba is a single patch, and not disjoint patches, similar to fine patch.
-        amrex::Box cdomain = grid_cba.minimalBox();
+        BoxArray grid_cba_reduced = grid_cba;
         if (do_pml_in_domain) {
-            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-                if (do_pml_Lo[idim]){
-                    // ncell is divided by refinement ratio to ensure that the
-                    // physical width of the PML region is equal in fine and coarse patch
-                    cdomain.growLo(idim, -ncell/ref_ratio[idim]);
-                }
-                if (do_pml_Hi[idim]){
-                    // ncell is divided by refinement ratio to ensure that the
-                    // physical width of the PML region is equal in fine and coarse patch
-                    cdomain.growHi(idim, -ncell/ref_ratio[idim]);
+            BoxList bl = grid_cba.boxList();
+            for (auto& b : bl) {
+                for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                    if (do_pml_Lo[idim]) {
+                        Box const& bb = amrex::adjCellLo(b, idim);
+                        if ( ! grid_cba.intersects(bb) ) {
+                            b.growLo(idim, -ncell/ref_ratio[idim]);
+                        }
+                    }
+                    if (do_pml_Hi[idim]) {
+                        Box const& bb = amrex::adjCellHi(b, idim);
+                        if ( ! grid_cba.intersects(bb) ) {
+                            b.growHi(idim, -ncell/ref_ratio[idim]);
+                        }
+                    }
                 }
             }
+            grid_cba_reduced = BoxArray(std::move(bl));
         }
-        const BoxArray grid_cba_reduced = (do_pml_in_domain) ?
-            BoxArray(grid_cba.boxList().intersect(cdomain)) : grid_cba;
+        Box const cdomain = grid_cba_reduced.minimalBox();
 
         const IntVect cncells = IntVect(ncell)/ref_ratio;
         const IntVect cdelta = IntVect(delta)/ref_ratio;
@@ -792,7 +812,6 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         } else {
             cdm.define(cba);
         }
-
         const amrex::BoxArray cba_Ex = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,0).ixType().toIntVect());
         const amrex::BoxArray cba_Ey = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,1).ixType().toIntVect());
         const amrex::BoxArray cba_Ez = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,2).ixType().toIntVect());