From 6be87d6fa38290fb44756b9679846db38b615385 Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com> Date: Wed, 28 Feb 2024 15:46:44 -0800 Subject: [PATCH] Adding support for pml in domain for multiple patches (#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> --- .../analysis_particles_in_pml_2dmr.py | 68 ++++++++++++++ Examples/Tests/particles_in_pml/inputs_mr_2d | 15 ++-- .../particles_in_pml_2d_MR.json | 12 +-- Regression/WarpX-tests.ini | 2 +- Source/BoundaryConditions/PML.cpp | 89 +++++++++++-------- 5 files changed, 136 insertions(+), 50 deletions(-) create mode 100755 Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py diff --git a/Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py b/Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py new file mode 100755 index 00000000000..81a981e70a6 --- /dev/null +++ b/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) diff --git a/Examples/Tests/particles_in_pml/inputs_mr_2d b/Examples/Tests/particles_in_pml/inputs_mr_2d index a6db42a3897..40e869a5e16 100644 --- a/Examples/Tests/particles_in_pml/inputs_mr_2d +++ b/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 diff --git a/Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json b/Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json index 11c0a2a6108..dfe63ddc1c2 100644 --- a/Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json +++ b/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 diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 128495195da..b25d41db5d2 100644 --- a/Regression/WarpX-tests.ini +++ b/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 = . diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 1f0b6f09a33..805b4fec181 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/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());