diff --git a/Examples/Tests/PML/analysis_pml_psatd_rz.py b/Examples/Tests/PML/analysis_pml_psatd_rz.py new file mode 100755 index 00000000000..b905d938bb8 --- /dev/null +++ b/Examples/Tests/PML/analysis_pml_psatd_rz.py @@ -0,0 +1,55 @@ +#! /usr/bin/env python + +# Copyright 2021 David Grote +# +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +""" +This script tests the absorption of fields in the PML in RZ geometry. + +The input file inputs_particle_rz is used: it features an electron +moving radially that launches a pulse. This scripts runs until +most of the pulse escapes the radial boundary. If the PML fails, +the pulse will remain with in the domain. +""" +import os +import sys + +import numpy as np +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 ) + +# yt 4.0+ has rounding issues with our domain data: +# RuntimeError: yt attempted to read outside the boundaries +# of a non-periodic domain along dimension 0. +if 'force_periodicity' in dir(ds): ds.force_periodicity() + +# 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 = ad0['boxlib', 'Ex'].to_ndarray() +Ez_array = ad0['boxlib', 'Ez'].to_ndarray() +max_Ex = np.abs(Ex_array).max() +max_Ez = np.abs(Ez_array).max() +print( f'max Ex = {max_Ex}' ) +print( f'max Ez = {max_Ez}' ) +max_Efield = max(max_Ex, max_Ez) + +# This tolerance was obtained empirically. As the simulation progresses, the field energy is leaking +# out through PML so that the max field diminishes with time. When the PML is working properly, +# the field level falls below 2 at the end of the simulation. +tolerance_abs = 2. +print('tolerance_abs: ' + str(tolerance_abs)) +assert max_Efield < tolerance_abs + +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, filename) diff --git a/Examples/Tests/PML/inputs_rz b/Examples/Tests/PML/inputs_rz new file mode 100644 index 00000000000..bd798c23ba1 --- /dev/null +++ b/Examples/Tests/PML/inputs_rz @@ -0,0 +1,49 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 500 +amr.n_cell = 32 128 +amr.max_level = 0 + +geometry.dims = RZ +geometry.prob_lo = 0. -100.e-6 +geometry.prob_hi = 50.e-6 +100.e-6 + +warpx.n_rz_azimuthal_modes = 2 + +################################# +######## Boundary condition ##### +################################# +boundary.field_lo = none periodic +boundary.field_hi = pml periodic + +# PML +warpx.pml_ncell = 10 +warpx.do_pml_in_domain = 0 + +################################# +############ NUMERICS ########### +################################# +algo.maxwell_solver = psatd +warpx.use_filter = 0 +algo.particle_shape = 1 + +################################# +############ PARTICLE ########### +################################# +particles.species_names = electron + +electron.charge = -q_e +electron.mass = m_e +electron.injection_style = "singleparticle" +electron.single_particle_pos = 0. 0. 0. +electron.single_particle_vel = 10. 0. 0. +electron.single_particle_weight = 1. + +################################# +########## DIAGNOSTICS ########## +################################# +diagnostics.diags_names = diag1 +diag1.intervals = 500 +diag1.fields_to_plot = Bx By Bz Ex Ey Ez jx jy jz rho +diag1.diag_type = Full diff --git a/Regression/Checksum/benchmarks_json/pml_psatd_rz.json b/Regression/Checksum/benchmarks_json/pml_psatd_rz.json new file mode 100644 index 00000000000..b228b11f522 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/pml_psatd_rz.json @@ -0,0 +1,14 @@ +{ + "lev=0": { + "Ex": 1317.2337418101674, + "Ey": 7.056477765087081e-12, + "Ez": 1156.1377910176022, + "Bx": 2.594232124519114e-20, + "By": 3.2688462949057505e-06, + "Bz": 4.292128402662449e-21, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0, + "rho": 0.0 + } +} diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 16ce07bbce5..b987c662c4a 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -130,6 +130,23 @@ compileTest = 0 doVis = 0 analysisRoutine = Examples/analysis_default_regression.py +[pml_psatd_rz] +buildDir = . +inputFile = Examples/Tests/PML/inputs_rz +runtime_params = +dim = 2 +addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE +cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_PSATD=ON +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +analysisRoutine = Examples/Tests/PML/analysis_pml_psatd_rz.py +tolerance = 1.e-14 + [silver_mueller_2d_x] buildDir = . inputFile = Examples/Tests/SilverMueller/inputs_2d_x diff --git a/Source/BoundaryConditions/CMakeLists.txt b/Source/BoundaryConditions/CMakeLists.txt index 3fa90d1f180..a8a24ca6dda 100644 --- a/Source/BoundaryConditions/CMakeLists.txt +++ b/Source/BoundaryConditions/CMakeLists.txt @@ -5,3 +5,10 @@ target_sources(WarpX WarpXFieldBoundaries.cpp WarpX_PEC.cpp ) + +if(WarpX_DIMS STREQUAL RZ) + target_sources(WarpX + PRIVATE + PML_RZ.cpp + ) +endif() diff --git a/Source/BoundaryConditions/Make.package b/Source/BoundaryConditions/Make.package index 10fbdf08051..43d18425ffc 100644 --- a/Source/BoundaryConditions/Make.package +++ b/Source/BoundaryConditions/Make.package @@ -1,5 +1,8 @@ CEXE_sources += PML.cpp WarpXEvolvePML.cpp CEXE_sources += WarpXFieldBoundaries.cpp WarpX_PEC.cpp +ifeq ($(USE_RZ),TRUE) + CEXE_sources += PML_RZ.cpp +endif VPATH_LOCATIONS += $(WARPX_HOME)/Source/BoundaryConditions diff --git a/Source/BoundaryConditions/PML_RZ.H b/Source/BoundaryConditions/PML_RZ.H new file mode 100644 index 00000000000..d6e6c991ea4 --- /dev/null +++ b/Source/BoundaryConditions/PML_RZ.H @@ -0,0 +1,75 @@ +/* Copyright 2021 David Grote + * + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PML_RZ_H_ +#define WARPX_PML_RZ_H_ + +#include "PML_RZ_fwd.H" + +#ifdef WARPX_USE_PSATD +# include "FieldSolver/SpectralSolver/SpectralSolverRZ.H" +#endif + +#include +#include +#include +#include + +#include + +#include +#include + +enum struct PatchType : int; + +class PML_RZ +{ +public: + PML_RZ (const int lev, const amrex::BoxArray& grid_ba, const amrex::DistributionMapping& grid_dm, + const amrex::Geometry* geom, const int ncell, const int do_pml_in_domain); + + void ApplyDamping(amrex::MultiFab* Et_fp, amrex::MultiFab* Ez_fp, + amrex::MultiFab* Bt_fp, amrex::MultiFab* Bz_fp, + amrex::Real dt); + + std::array GetE_fp (); + std::array GetB_fp (); + +#ifdef WARPX_USE_PSATD + void PushPSATD (const int lev); +#endif + + void FillBoundaryE (); + void FillBoundaryB (); + void FillBoundaryE (PatchType patch_type); + void FillBoundaryB (PatchType patch_type); + + void CheckPoint (const std::string& dir) const; + void Restart (const std::string& dir); + + ~PML_RZ () = default; + +private: + + const int m_ncell; + const int m_do_pml_in_domain; + const amrex::Geometry* m_geom; + + // Only contains Er and Et, and Br and Bt + std::array,2> pml_E_fp; + std::array,2> pml_B_fp; + +#ifdef WARPX_USE_PSATD + void PushPMLPSATDSinglePatchRZ ( const int lev, + SpectralSolverRZ& solver, + std::array,2>& pml_E, + std::array,2>& pml_B); +#endif + +}; + +#endif diff --git a/Source/BoundaryConditions/PML_RZ.cpp b/Source/BoundaryConditions/PML_RZ.cpp new file mode 100644 index 00000000000..422747d07ff --- /dev/null +++ b/Source/BoundaryConditions/PML_RZ.cpp @@ -0,0 +1,227 @@ +/* Copyright 2021 David Grote + * + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PML_RZ.H" + +#include "BoundaryConditions/PML_RZ.H" +#ifdef WARPX_USE_PSATD +# include "FieldSolver/SpectralSolver/SpectralFieldDataRZ.H" +#endif +#include "Utils/WarpXConst.H" +#include "WarpX.H" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace amrex; + +PML_RZ::PML_RZ (const int lev, const amrex::BoxArray& grid_ba, const amrex::DistributionMapping& grid_dm, + const amrex::Geometry* geom, const int ncell, const int do_pml_in_domain) + : m_ncell(ncell), + m_do_pml_in_domain(do_pml_in_domain), + m_geom(geom) +{ + + const amrex::MultiFab & Er_fp = WarpX::GetInstance().getEfield_fp(lev,0); + const amrex::MultiFab & Et_fp = WarpX::GetInstance().getEfield_fp(lev,1); + pml_E_fp[0] = std::make_unique( amrex::convert( grid_ba, Er_fp.ixType().toIntVect() ), grid_dm, + Er_fp.nComp(), Er_fp.nGrow() ); + pml_E_fp[1] = std::make_unique( amrex::convert( grid_ba, Et_fp.ixType().toIntVect() ), grid_dm, + Et_fp.nComp(), Et_fp.nGrow() ); + + const amrex::MultiFab & Br_fp = WarpX::GetInstance().getBfield_fp(lev,0); + const amrex::MultiFab & Bt_fp = WarpX::GetInstance().getBfield_fp(lev,1); + pml_B_fp[0] = std::make_unique( amrex::convert( grid_ba, Br_fp.ixType().toIntVect() ), grid_dm, + Br_fp.nComp(), Br_fp.nGrow() ); + pml_B_fp[1] = std::make_unique( amrex::convert( grid_ba, Bt_fp.ixType().toIntVect() ), grid_dm, + Bt_fp.nComp(), Bt_fp.nGrow() ); + + pml_E_fp[0]->setVal(0.0); + pml_E_fp[1]->setVal(0.0); + pml_B_fp[0]->setVal(0.0); + pml_B_fp[1]->setVal(0.0); + +} + +void +PML_RZ::ApplyDamping (amrex::MultiFab* Et_fp, amrex::MultiFab* Ez_fp, + amrex::MultiFab* Bt_fp, amrex::MultiFab* Bz_fp, + amrex::Real dt) +{ + + const amrex::Real dr = m_geom->CellSize(0); + const amrex::Real cdt_over_dr = PhysConst::c*dt/dr; + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( amrex::MFIter mfi(*Et_fp, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + amrex::Array4 const& Et_arr = Et_fp->array(mfi); + amrex::Array4 const& Ez_arr = Ez_fp->array(mfi); + amrex::Array4 const& Bt_arr = Bt_fp->array(mfi); + amrex::Array4 const& Bz_arr = Bz_fp->array(mfi); + + amrex::Array4 const& pml_Et_arr = pml_E_fp[1]->array(mfi); + amrex::Array4 const& pml_Bt_arr = pml_B_fp[1]->array(mfi); + + // Get the tileboxes from Efield and Bfield so that they include the guard cells + // They are all the same, cell centered + amrex::Box tilebox = amrex::convert((*Et_fp)[mfi].box(), Et_fp->ixType().toIntVect()); + + // Box for the whole simulation domain + amrex::Box const& domain = m_geom->Domain(); + const int nr_domain = domain.bigEnd(0); + + // Set tilebox to only include the upper radial cells + const int nr_damp = m_ncell; + int nr_damp_min; + if (m_do_pml_in_domain) { + nr_damp_min = nr_domain - nr_damp; + } else { + nr_damp_min = nr_domain; + } + tilebox.setSmall(0, nr_damp_min + 1); + + amrex::ParallelFor( tilebox, Et_fp->nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + { + const amrex::Real rr = static_cast(i - nr_damp_min); + const amrex::Real wr = rr/nr_damp; + const amrex::Real damp_factor = std::exp( -4._rt * cdt_over_dr * wr*wr ); + + // Substract the theta PML fields from the regular theta fields + Et_arr(i,j,k,icomp) -= pml_Et_arr(i,j,k,icomp); + Bt_arr(i,j,k,icomp) -= pml_Bt_arr(i,j,k,icomp); + // Damp the theta PML fields + pml_Et_arr(i,j,k,icomp) *= damp_factor; + pml_Bt_arr(i,j,k,icomp) *= damp_factor; + // Add the theta PML fields back to the regular theta fields + Et_arr(i,j,k,icomp) += pml_Et_arr(i,j,k,icomp); + Bt_arr(i,j,k,icomp) += pml_Bt_arr(i,j,k,icomp); + // Damp the z fields + Bz_arr(i,j,k,icomp) *= damp_factor; + Ez_arr(i,j,k,icomp) *= damp_factor; + }); + } +} + +std::array +PML_RZ::GetE_fp () +{ + return {pml_E_fp[0].get(), pml_E_fp[1].get()}; +} + +std::array +PML_RZ::GetB_fp () +{ + return {pml_B_fp[0].get(), pml_B_fp[1].get()}; +} + +void +PML_RZ::FillBoundaryE () +{ + FillBoundaryE(PatchType::fine); +} + +void +PML_RZ::FillBoundaryE (PatchType patch_type) +{ + if (patch_type == PatchType::fine && pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0) + { + const amrex::Periodicity& period = m_geom->periodicity(); + Vector mf{pml_E_fp[0].get(),pml_E_fp[1].get()}; + amrex::FillBoundary(mf, period); + } +} + +void +PML_RZ::FillBoundaryB () +{ + FillBoundaryB(PatchType::fine); +} + +void +PML_RZ::FillBoundaryB (PatchType patch_type) +{ + if (patch_type == PatchType::fine && pml_B_fp[0]) + { + const amrex::Periodicity& period = m_geom->periodicity(); + Vector mf{pml_B_fp[0].get(),pml_B_fp[1].get()}; + amrex::FillBoundary(mf, period); + } +} + +void +PML_RZ::CheckPoint (const std::string& dir) const +{ + if (pml_E_fp[0]) + { + VisMF::AsyncWrite(*pml_E_fp[0], dir+"_Er_fp"); + VisMF::AsyncWrite(*pml_E_fp[1], dir+"_Et_fp"); + VisMF::AsyncWrite(*pml_B_fp[0], dir+"_Br_fp"); + VisMF::AsyncWrite(*pml_B_fp[1], dir+"_Bt_fp"); + } +} + +void +PML_RZ::Restart (const std::string& dir) +{ + if (pml_E_fp[0]) + { + VisMF::Read(*pml_E_fp[0], dir+"_Er_fp"); + VisMF::Read(*pml_E_fp[1], dir+"_Et_fp"); + VisMF::Read(*pml_B_fp[0], dir+"_Br_fp"); + VisMF::Read(*pml_B_fp[1], dir+"_Bt_fp"); + } +} + +#ifdef WARPX_USE_PSATD +void +PML_RZ::PushPSATD (const int lev) +{ + // Update the fields on the fine and coarse patch + WarpX& warpx = WarpX::GetInstance(); + SpectralSolverRZ& solver = warpx.get_spectral_solver_fp(lev); + PushPMLPSATDSinglePatchRZ(lev, solver, pml_E_fp, pml_B_fp); +} + +void +PML_RZ::PushPMLPSATDSinglePatchRZ ( + const int lev, + SpectralSolverRZ& solver, + std::array,2>& pml_E, + std::array,2>& pml_B) +{ + const SpectralFieldIndex& Idx = solver.m_spectral_index; + + // Perform forward Fourier transforms + solver.ForwardTransform(lev, *pml_E[0], Idx.Er_pml, *pml_E[1], Idx.Et_pml); + solver.ForwardTransform(lev, *pml_B[0], Idx.Br_pml, *pml_B[1], Idx.Bt_pml); + + // Advance fields in spectral space + const bool doing_pml = true; + solver.pushSpectralFields(doing_pml); + + // Perform backward Fourier transforms + solver.BackwardTransform(lev, *pml_E[0], Idx.Er_pml, *pml_E[1], Idx.Et_pml); + solver.BackwardTransform(lev, *pml_B[0], Idx.Br_pml, *pml_B[1], Idx.Bt_pml); +} +#endif diff --git a/Source/BoundaryConditions/PML_RZ_fwd.H b/Source/BoundaryConditions/PML_RZ_fwd.H new file mode 100644 index 00000000000..36111a22480 --- /dev/null +++ b/Source/BoundaryConditions/PML_RZ_fwd.H @@ -0,0 +1,8 @@ +/* Copyright 2021 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +class PML_RZ; diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index fce75d21e7e..852e83e1ce0 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -8,6 +8,9 @@ #include "WarpX.H" #include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif #include "PML_current.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX_PML_kernels.H" @@ -45,19 +48,33 @@ WarpX::DampPML () } void -WarpX::DampPML (int lev) +WarpX::DampPML (const int lev) { DampPML(lev, PatchType::fine); if (lev > 0) DampPML(lev, PatchType::coarse); } void -WarpX::DampPML (int lev, PatchType patch_type) +WarpX::DampPML (const int lev, PatchType patch_type) { if (!do_pml) return; WARPX_PROFILE("WarpX::DampPML()"); +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (pml_rz[lev]) { + pml_rz[lev]->ApplyDamping(Efield_fp[lev][1].get(), Efield_fp[lev][2].get(), + Bfield_fp[lev][1].get(), Bfield_fp[lev][2].get(), + dt[lev]); + } +#endif + if (pml[lev]) { + DampPML_Cartesian (lev, patch_type); + } +} +void +WarpX::DampPML_Cartesian (const int lev, PatchType patch_type) +{ const bool dive_cleaning = WarpX::do_pml_dive_cleaning; const bool divb_cleaning = WarpX::do_pml_divb_cleaning; @@ -245,6 +262,7 @@ WarpX::DampJPML (int lev, PatchType patch_type) { if (!do_pml) return; if (!do_pml_j_damping) return; + if (!pml[lev]) return; WARPX_PROFILE("WarpX::DampJPML()"); @@ -349,7 +367,7 @@ WarpX::CopyJPML () { for (int lev = 0; lev <= finest_level; ++lev) { - if (pml[lev]->ok()){ + if (pml[lev] && pml[lev]->ok()){ pml[lev]->CopyJtoPMLs({ current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get() }, diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 890d4ceeceb..8e512221963 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp @@ -1,6 +1,9 @@ #include "FlushFormatCheckpoint.H" #include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif #include "Diagnostics/ParticleDiag/ParticleDiag.H" #include "Particles/WarpXParticleContainer.H" #include "Utils/WarpXProfilerWrapper.H" @@ -138,9 +141,17 @@ FlushFormatCheckpoint::WriteToFile ( } } - if (warpx.DoPML() && warpx.GetPML(lev)) { - warpx.GetPML(lev)->CheckPoint( - amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pml")); + if (warpx.DoPML()) { + if (warpx.GetPML(lev)) { + warpx.GetPML(lev)->CheckPoint( + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pml")); + } +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (warpx.GetPML_RZ(lev)) { + warpx.GetPML_RZ(lev)->CheckPoint( + amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pml_rz")); + } +#endif } } diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 2eb392f489c..1347c660855 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -8,6 +8,9 @@ * License: BSD-3-Clause-LBNL */ #include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif #include "FieldIO.H" #include "Particles/MultiParticleContainer.H" #include "Utils/CoarsenIO.H" @@ -330,7 +333,12 @@ WarpX::InitFromCheckpoint () if (do_pml) { for (int lev = 0; lev < nlevs; ++lev) { - pml[lev]->Restart(amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "pml")); + if (pml[lev]) + pml[lev]->Restart(amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "pml")); +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (pml_rz[lev]) + pml_rz[lev]->Restart(amrex::MultiFabFileFullPrefix(lev, restart_chkfile, level_prefix, "pml_rz")); +#endif } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index f46fd9eb99f..1ce24961d3a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -12,5 +12,6 @@ if(WarpX_DIMS STREQUAL RZ) SpectralBaseAlgorithmRZ.cpp PsatdAlgorithmRZ.cpp GalileanPsatdAlgorithmRZ.cpp + PMLPsatdAlgorithmRZ.cpp ) endif() diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index 11a10bb91e2..ab32564521e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -7,6 +7,7 @@ ifeq ($(USE_RZ),TRUE) CEXE_sources += SpectralBaseAlgorithmRZ.cpp CEXE_sources += PsatdAlgorithmRZ.cpp CEXE_sources += GalileanPsatdAlgorithmRZ.cpp + CEXE_sources += PMLPsatdAlgorithmRZ.cpp endif VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H new file mode 100644 index 00000000000..e0956a514d2 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.H @@ -0,0 +1,72 @@ +/* Copyright 2021 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PMLPSATD_ALGORITHM_RZ_H_ +#define WARPX_PMLPSATD_ALGORITHM_RZ_H_ + +#include "SpectralBaseAlgorithmRZ.H" + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PMLPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ +{ + + public: + PMLPsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, + amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, + int const n_rz_azimuthal_modes, int const norder_z, + bool const nodal, amrex::Real const dt_step); + + // Redefine functions from base class + virtual void pushSpectralFields (SpectralFieldDataRZ & f) override final; + + void InitializeSpectralCoefficients (SpectralFieldDataRZ const & f); + + /** + * \brief Virtual function for current correction in Fourier space + * ( Vay et al, 2013). + * This function overrides the virtual function \c CurrentCorrection in the + * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density + */ + virtual void CurrentCorrection (const int lev, + SpectralFieldDataRZ& field_data, + std::array,3>& current, + const std::unique_ptr& rho) override final; + + /** + * \brief Virtual function for Vay current deposition in Fourier space + * ( Vay et al, 2013). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (const int lev, + SpectralFieldDataRZ& field_data, + std::array,3>& current) override final; + + private: + + SpectralFieldIndex m_spectral_index; + + bool coefficients_initialized; + // Note that dt is saved to use in InitializeSpectralCoefficients + amrex::Real m_dt; + SpectralRealCoefficients C_coef, S_ck_coef; +}; + +#endif // WARPX_PMLPSATD_ALGORITHM_RZ_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp new file mode 100644 index 00000000000..892c542c745 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithmRZ.cpp @@ -0,0 +1,177 @@ +/* Copyright 2021 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PMLPsatdAlgorithmRZ.H" +#include "FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include + +using amrex::operator""_rt; + + +/* \brief Initialize coefficients for the update equation */ +PMLPsatdAlgorithmRZ::PMLPsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, + amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, + int const n_rz_azimuthal_modes, int const norder_z, + bool const nodal, amrex::Real const dt) + // Initialize members of base class + : SpectralBaseAlgorithmRZ(spectral_kspace, dm, spectral_index, norder_z, nodal), + m_spectral_index(spectral_index), + m_dt(dt) +{ + // Allocate the arrays of coefficients + amrex::BoxArray const & ba = spectral_kspace.spectralspace_ba; + C_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + S_ck_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); + + coefficients_initialized = false; +} + +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ +void +PMLPsatdAlgorithmRZ::pushSpectralFields (SpectralFieldDataRZ & f) +{ + + if (not coefficients_initialized) { + // This is called from here since it needs the kr values + // which can be obtained from the SpectralFieldDataRZ + InitializeSpectralCoefficients(f); + coefficients_initialized = true; + } + + const SpectralFieldIndex& Idx = m_spectral_index; + + // Loop over boxes + for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){ + + amrex::Box const & bx = f.fields[mfi].box(); + + // Extract arrays for the fields to be updated + amrex::Array4 const& fields = f.fields[mfi].array(); + // Extract arrays for the coefficients + amrex::Array4 const& C_arr = C_coef[mfi].array(); + amrex::Array4 const& S_ck_arr = S_ck_coef[mfi].array(); + + // Extract pointers for the k vectors + HankelTransform::RealVector const & kr_modes = f.getKrArray(mfi); + amrex::Real const* kr_arr = kr_modes.dataPtr(); + int const nr = bx.length(0); + + // Loop over indices within one box + // Note that k = 0 + int const modes = f.n_rz_azimuthal_modes; + amrex::ParallelFor(bx, modes, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept + { + + // All of the fields of each mode are grouped together + int const Ep_m_pml = Idx.Er_pml + Idx.n_fields*mode; + int const Em_m_pml = Idx.Et_pml + Idx.n_fields*mode; + int const Bp_m_pml = Idx.Br_pml + Idx.n_fields*mode; + int const Bm_m_pml = Idx.Bt_pml + Idx.n_fields*mode; + int const Ez_m = Idx.Ez + Idx.n_fields*mode; + int const Bz_m = Idx.Bz + Idx.n_fields*mode; + + // Record old values of the fields to be updated + Complex const Ep_old_pml = fields(i,j,k,Ep_m_pml); + Complex const Em_old_pml = fields(i,j,k,Em_m_pml); + Complex const Bp_old_pml = fields(i,j,k,Bp_m_pml); + Complex const Bm_old_pml = fields(i,j,k,Bm_m_pml); + Complex const Ez_old = fields(i,j,k,Ez_m); + Complex const Bz_old = fields(i,j,k,Bz_m); + + // k vector values, and coefficients + // The k values for each mode are grouped together + int const ir = i + nr*mode; + amrex::Real const kr = kr_arr[ir]; + + constexpr amrex::Real c2 = PhysConst::c*PhysConst::c; + Complex const I = Complex{0._rt,1._rt}; + amrex::Real const C = C_arr(i,j,k,mode); + amrex::Real const S_ck = S_ck_arr(i,j,k,mode); + + // Update E (see WarpX online documentation: theory section) + fields(i,j,k,Ep_m_pml) = C*Ep_old_pml + + S_ck*(-c2*I*kr/2._rt*Bz_old); + fields(i,j,k,Em_m_pml) = C*Em_old_pml + + S_ck*(-c2*I*kr/2._rt*Bz_old); + // Update B (see WarpX online documentation: theory section) + fields(i,j,k,Bp_m_pml) = C*Bp_old_pml + - S_ck*(-I*kr/2._rt*Ez_old); + fields(i,j,k,Bm_m_pml) = C*Bm_old_pml + - S_ck*(-I*kr/2._rt*Ez_old); + + }); + } +} + +void PMLPsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f) +{ + + // Fill them with the right values: + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc + for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){ + + amrex::Box const & bx = f.fields[mfi].box(); + + // Extract pointers for the k vectors + amrex::Real const* const modified_kz = modified_kz_vec[mfi].dataPtr(); + + // Extract arrays for the coefficients + amrex::Array4 const& C = C_coef[mfi].array(); + amrex::Array4 const& S_ck = S_ck_coef[mfi].array(); + + HankelTransform::RealVector const & kr_modes = f.getKrArray(mfi); + amrex::Real const* kr_arr = kr_modes.dataPtr(); + int const nr = bx.length(0); + amrex::Real const dt = m_dt; + + // Loop over indices within one box + int const modes = f.n_rz_azimuthal_modes; + amrex::ParallelFor(bx, modes, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept + { + // Calculate norm of vector + int const ir = i + nr*mode; + amrex::Real const kr = kr_arr[ir]; + amrex::Real const kz = modified_kz[j]; + amrex::Real const k_norm = std::sqrt(kr*kr + kz*kz); + + // Calculate coefficients + constexpr amrex::Real c = PhysConst::c; + if (k_norm != 0._rt){ + C(i,j,k,mode) = std::cos(c*k_norm*dt); + S_ck(i,j,k,mode) = std::sin(c*k_norm*dt)/(c*k_norm); + } else { // Handle k_norm = 0, by using the analytical limit + C(i,j,k,mode) = 1._rt; + S_ck(i,j,k,mode) = dt; + } + }); + } +} + +void +PMLPsatdAlgorithmRZ::CurrentCorrection (const int /* lev */, + SpectralFieldDataRZ& /* field_data */, + std::array,3>& /* current */, + const std::unique_ptr& /* rho */) +{ + amrex::Abort("Current correction not implemented in RZ geometry PML"); +} + +void +PMLPsatdAlgorithmRZ::VayDeposition (const int /* lev */, + SpectralFieldDataRZ& /*field_data*/, + std::array,3>& /*current*/) +{ + amrex::Abort("Vay deposition not implemented in RZ geometry PML"); +} diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 9b748a048fc..4f0ad1c2054 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -50,13 +50,16 @@ class SpectralFieldIndex * div(B) = 0 law (new field G in the update equations) * \param[in] pml whether the indices are used to access spectral data * for the PML spectral solver + * \param[in] pml_rz whether the indices are used to access spectral data + * for the RZ PML spectral solver */ SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning, - const bool pml); + const bool pml, + const bool pml_rz = false); /** * \brief Default constructor @@ -97,6 +100,9 @@ class SpectralFieldIndex // PML with div(E) and/or div(B) cleaning int Exx = -1, Eyy = -1, Ezz = -1, Bxx = -1, Byy = -1, Bzz = -1; int Fx = -1, Fy = -1, Fz = -1, Gx = -1, Gy = -1, Gz = -1; + + // PML RZ + int Er_pml = -1, Et_pml = -1, Br_pml = -1, Bt_pml = -1; }; /** \brief Class that stores the fields in spectral space, and performs the diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 86fe8f7bdd8..23498214a41 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -37,7 +37,8 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning, - const bool pml) + const bool pml, + const bool pml_rz) { // TODO Use these to allocate rho_old, rho_new, F, and G only when needed amrex::ignore_unused(update_with_rho); @@ -72,6 +73,13 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, Jy_new = c++; Jz_new = c++; } + + if (pml_rz) + { + Er_pml = c++; Et_pml = c++; + Br_pml = c++; Bt_pml = c++; + } + } else // PML { diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 827d226c2c2..9f9b8dbefa1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -34,6 +34,7 @@ class SpectralSolverRZ int const norder_z, bool const nodal, const amrex::Vector& v_galilean, amrex::RealVect const dx, amrex::Real const dt, + bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, const bool do_multi_J, @@ -63,7 +64,7 @@ class SpectralSolverRZ amrex::MultiFab& field_mf2, int const field_index2); /* \brief Update the fields in spectral space, over one timestep */ - void pushSpectralFields (); + void pushSpectralFields (const bool doing_pml=false); /* \brief Initialize K space filtering arrays */ void InitFilter (amrex::IntVect const & filter_npass_each_dir, @@ -160,6 +161,7 @@ class SpectralSolverRZ SpectralFieldDataRZ field_data; // Store field in spectral space // and perform the Fourier transforms std::unique_ptr algorithm; + std::unique_ptr PML_algorithm; // Defines field update equation in spectral space, // and the associated coefficients. // SpectralBaseAlgorithmRZ is a base class ; this pointer is meant diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index e187b9d27be..558d4a78a39 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -6,6 +6,7 @@ */ #include "SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H" #include "SpectralAlgorithms/PsatdAlgorithmRZ.H" +#include "SpectralAlgorithms/PMLPsatdAlgorithmRZ.H" #include "SpectralKSpaceRZ.H" #include "SpectralSolverRZ.H" #include "Utils/WarpXProfilerWrapper.H" @@ -22,8 +23,7 @@ * \param nodal Whether the solver is applied to a nodal or staggered grid * \param dx Cell size along each dimension * \param dt Time step - * \param pml Whether the boxes in which the solver is applied are PML boxes - * PML is not supported. + * \param with_pml Whether PML boundary will be used */ SpectralSolverRZ::SpectralSolverRZ (const int lev, amrex::BoxArray const & realspace_ba, @@ -32,6 +32,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, int const norder_z, bool const nodal, const amrex::Vector& v_galilean, amrex::RealVect const dx, amrex::Real const dt, + bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, const bool do_multi_J, @@ -45,13 +46,16 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, // the spectral space corresponding to each box in `realspace_ba`, // as well as the value of the corresponding k coordinates. - const bool pml = false; + const bool is_pml = false; m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - do_multi_J, dive_cleaning, divb_cleaning, pml); + do_multi_J, dive_cleaning, divb_cleaning, is_pml, with_pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space - // PML is not supported. + if (with_pml) { + PML_algorithm = std::make_unique( + k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, dt); + } if (v_galilean[2] == 0) { // v_galilean is 0: use standard PSATD algorithm algorithm = std::make_unique( @@ -117,12 +121,16 @@ SpectralSolverRZ::BackwardTransform (const int lev, /* \brief Update the fields in spectral space, over one timestep */ void -SpectralSolverRZ::pushSpectralFields () { +SpectralSolverRZ::pushSpectralFields (const bool doing_pml) { WARPX_PROFILE("SpectralSolverRZ::pushSpectralFields"); // Virtual function: the actual function used here depends // on the sub-class of `SpectralBaseAlgorithm` that was // initialized in the constructor of `SpectralSolverRZ` - algorithm->pushSpectralFields(field_data); + if (doing_pml) { + PML_algorithm->pushSpectralFields(field_data); + } else { + algorithm->pushSpectralFields(field_data); + } } /** diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index a096401bc9d..87bad1b85f4 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -15,6 +15,7 @@ # include "FieldSolver/SpectralSolver/SpectralFieldData.H" # ifdef WARPX_DIM_RZ # include "FieldSolver/SpectralSolver/SpectralSolverRZ.H" +# include "BoundaryConditions/PML_RZ.H" # else # include "FieldSolver/SpectralSolver/SpectralSolver.H" # endif @@ -414,12 +415,18 @@ WarpX::PushPSATD () PSATDForwardTransformEB(); PSATDForwardTransformJ(); + // Do rho FFTs only if needed if (WarpX::update_with_rho || WarpX::current_correction || WarpX::do_dive_cleaning) { PSATDForwardTransformRho(0,0); // rho old PSATDForwardTransformRho(1,1); // rho new } + +#ifdef WARPX_DIM_RZ + if (pml_rz[0]) pml_rz[0]->PushPSATD(0); +#endif + if (WarpX::do_dive_cleaning) PSATDForwardTransformF(); if (WarpX::do_divb_cleaning) PSATDForwardTransformG(); PSATDPushSpectralFields(); @@ -431,7 +438,7 @@ WarpX::PushPSATD () // Evolve the fields in the PML boxes for (int lev = 0; lev <= finest_level; ++lev) { - if (do_pml && pml[lev]->ok()) + if (pml[lev] && pml[lev]->ok()) { pml[lev]->PushPSATD(lev); } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index be810d21e2f..f0ffbd04b39 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -11,6 +11,9 @@ #include "WarpX.H" #include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif #include "Diagnostics/BackTransformedDiagnostic.H" #include "Diagnostics/MultiDiagnostics.H" #include "Diagnostics/ReducedDiags/MultiReducedDiags.H" @@ -236,9 +239,10 @@ WarpX::InitPML () { amrex::IntVect do_pml_Lo_corrected = do_pml_Lo; -#ifdef WARPX_DIM_RZ +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) do_pml_Lo_corrected[0] = 0; // no PML at r=0, in cylindrical geometry -#endif + pml_rz[0] = std::make_unique(0, boxArray(0), DistributionMap(0), &Geom(0), pml_ncell, do_pml_in_domain); +#else pml[0] = std::make_unique(0, boxArray(0), DistributionMap(0), &Geom(0), nullptr, pml_ncell, pml_delta, amrex::IntVect::TheZeroVector(), dt[0], nox_fft, noy_fft, noz_fft, do_nodal, @@ -247,6 +251,7 @@ WarpX::InitPML () do_pml_dive_cleaning, do_pml_divb_cleaning, guard_cells.ng_FieldSolver.max(), do_pml_Lo_corrected, do_pml_Hi); +#endif for (int lev = 1; lev <= finest_level; ++lev) { @@ -288,7 +293,8 @@ WarpX::ComputePMLFactors () { for (int lev = 0; lev <= finest_level; ++lev) { - pml[lev]->ComputePMLFactors(dt[lev]); + if (pml[lev]) + pml[lev]->ComputePMLFactors(dt[lev]); } } } diff --git a/Source/Parallelization/GuardCellManager.H b/Source/Parallelization/GuardCellManager.H index 182b4aa80be..05d99de4df6 100644 --- a/Source/Parallelization/GuardCellManager.H +++ b/Source/Parallelization/GuardCellManager.H @@ -44,6 +44,9 @@ public: * \param do_electrostatic Whether to run in electrostatic mode i.e. solving the Poisson equation instead of the Maxwell equations. * \param do_multi_J Whether to use the multi-J PSATD scheme * \param fft_do_time_averaging Whether to average the E and B field in time (with PSATD) before interpolating them onto the macro-particles + * \param do_pml whether pml is turned on (only used by RZ PSATD) + * \param do_pml_in_domain whether pml is done in the domain (only used by RZ PSATD) + * \param pml_ncell number of cells on the pml layer (only used by RZ PSATD) * \param ref_ratios mesh refinement ratios between mesh-refinement levels */ void Init( @@ -65,6 +68,9 @@ public: const int do_electrostatic, const int do_multi_J, const bool fft_do_time_averaging, + const bool do_pml, + const int do_pml_in_domain, + const int pml_ncell, const amrex::Vector& ref_ratios); // Guard cells allocated for MultiFabs E and B diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp index 2fa0c4c4a71..d09e469519b 100644 --- a/Source/Parallelization/GuardCellManager.cpp +++ b/Source/Parallelization/GuardCellManager.cpp @@ -49,6 +49,9 @@ guardCellManager::Init ( const int do_electrostatic, const int do_multi_J, const bool fft_do_time_averaging, + const bool do_pml, + const int do_pml_in_domain, + const int pml_ncell, const amrex::Vector& ref_ratios) { // When using subcycling, the particles on the fine level perform two pushes @@ -202,6 +205,16 @@ guardCellManager::Init ( IntVect ngFFT = IntVect(ngFFt_z); #endif +#ifdef WARPX_DIM_RZ + if (do_pml) { + if (!do_pml_in_domain) { + ngFFT[0] = std::max(ngFFT[0], pml_ncell); + } + } +#else + amrex::ignore_unused(do_pml, do_pml_in_domain, pml_ncell); +#endif + // All boxes should have the same number of guard cells, to avoid temporary parallel copies: // thus we take the maximum of the required number of guard cells over all available fields. for (int i_dim = 0; i_dim < AMREX_SPACEDIM; i_dim++) { diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index d0f7e352a75..9e4a7de884a 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -9,6 +9,9 @@ #include "WarpX.H" #include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif #include "Filter/BilinearFilter.H" #include "Utils/CoarsenMR.H" #include "Utils/IntervalsParser.H" @@ -548,14 +551,22 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type, IntVect ng) { if (patch_type == PatchType::fine) { - if (do_pml && pml[lev]->ok()) - { - pml[lev]->ExchangeE(patch_type, - { Efield_fp[lev][0].get(), - Efield_fp[lev][1].get(), - Efield_fp[lev][2].get() }, - do_pml_in_domain); - pml[lev]->FillBoundaryE(patch_type); + if (do_pml) { + if (pml[lev] && pml[lev]->ok()) + { + pml[lev]->ExchangeE(patch_type, + { Efield_fp[lev][0].get(), + Efield_fp[lev][1].get(), + Efield_fp[lev][2].get() }, + do_pml_in_domain); + pml[lev]->FillBoundaryE(patch_type); + } +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (pml_rz[lev]) + { + pml_rz[lev]->FillBoundaryE(patch_type); + } +#endif } const amrex::Periodicity& period = Geom(lev).periodicity(); @@ -610,15 +621,25 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type, IntVect ng) { if (patch_type == PatchType::fine) { - if (do_pml && pml[lev]->ok()) + if (do_pml) { - pml[lev]->ExchangeB(patch_type, - { Bfield_fp[lev][0].get(), - Bfield_fp[lev][1].get(), - Bfield_fp[lev][2].get() }, - do_pml_in_domain); - pml[lev]->FillBoundaryB(patch_type); + if (pml[lev] && pml[lev]->ok()) + { + pml[lev]->ExchangeB(patch_type, + { Bfield_fp[lev][0].get(), + Bfield_fp[lev][1].get(), + Bfield_fp[lev][2].get() }, + do_pml_in_domain); + pml[lev]->FillBoundaryB(patch_type); + } +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (pml_rz[lev]) + { + pml_rz[lev]->FillBoundaryB(patch_type); + } +#endif } + const amrex::Periodicity& period = Geom(lev).periodicity(); if ( safe_guard_cells ) { Vector mf{Bfield_fp[lev][0].get(),Bfield_fp[lev][1].get(),Bfield_fp[lev][2].get()}; @@ -777,7 +798,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type, IntVect ng) { if (patch_type == PatchType::fine) { - if (do_pml && pml[lev]->ok()) + if (do_pml && pml[lev] && pml[lev]->ok()) { if (F_fp[lev]) pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); @@ -792,7 +813,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type, IntVect ng) } else if (patch_type == PatchType::coarse) { - if (do_pml && pml[lev]->ok()) + if (do_pml && pml[lev] && pml[lev]->ok()) { if (F_cp[lev]) pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); @@ -821,7 +842,7 @@ void WarpX::FillBoundaryG (int lev, PatchType patch_type, IntVect ng) { if (patch_type == PatchType::fine) { - if (do_pml && pml[lev]->ok()) + if (do_pml && pml[lev] && pml[lev]->ok()) { if (G_fp[lev]) pml[lev]->ExchangeG(patch_type, G_fp[lev].get(), do_pml_in_domain); pml[lev]->FillBoundaryG(patch_type); @@ -836,7 +857,7 @@ void WarpX::FillBoundaryG (int lev, PatchType patch_type, IntVect ng) } else if (patch_type == PatchType::coarse) { - if (do_pml && pml[lev]->ok()) + if (do_pml && pml[lev] && pml[lev]->ok()) { if (G_cp[lev]) pml[lev]->ExchangeG(patch_type, G_cp[lev].get(), do_pml_in_domain); pml[lev]->FillBoundaryG(patch_type); @@ -1284,7 +1305,7 @@ void WarpX::NodalSyncPML (int lev) void WarpX::NodalSyncPML (int lev, PatchType patch_type) { - if (pml[lev]->ok()) + if (pml[lev] && pml[lev]->ok()) { const std::array& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); @@ -1308,6 +1329,24 @@ void WarpX::NodalSyncPML (int lev, PatchType patch_type) WarpXCommUtil::OverrideSync(*pml_G, period); } } + +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (pml_rz[lev]) + { + // This is not actually needed with RZ PSATD since the + // arrays are always cell centered. Keep for now since + // it may be useful if the PML is used with FDTD + const std::array pml_rz_E = pml_rz[lev]->GetE_fp(); + const std::array pml_rz_B = pml_rz[lev]->GetB_fp(); + + // Always synchronize nodal points + const amrex::Periodicity& period = Geom(lev).periodicity(); + pml_rz_E[0]->OverrideSync(period); + pml_rz_E[1]->OverrideSync(period); + pml_rz_B[0]->OverrideSync(period); + pml_rz_B[1]->OverrideSync(period); + } +#endif } void WarpX::NodalSync (amrex::Vector,3>>& mf_fp, diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index d844c23170f..34087ab2706 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -9,6 +9,9 @@ #include "WarpX.H" #include "BoundaryConditions/PML.H" +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) +# include "BoundaryConditions/PML_RZ.H" +#endif #include "Particles/MultiParticleContainer.H" #include "Parallelization/WarpXCommUtil.H" #include "Utils/WarpXConst.H" @@ -179,13 +182,20 @@ WarpX::MoveWindow (const int step, bool move_j) if (move_j) { shiftMF(*current_fp[lev][dim], geom[lev], num_shift, dir); } - if (do_pml && pml[lev]->ok()) { + if (pml[lev] && pml[lev]->ok()) { const std::array& pml_B = pml[lev]->GetB_fp(); const std::array& pml_E = pml[lev]->GetE_fp(); shiftMF(*pml_B[dim], geom[lev], num_shift, dir); shiftMF(*pml_E[dim], geom[lev], num_shift, dir); } - +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (pml_rz[lev] && dim < 2) { + const std::array& pml_rz_B = pml_rz[lev]->GetB_fp(); + const std::array& pml_rz_E = pml_rz[lev]->GetE_fp(); + shiftMF(*pml_rz_B[dim], geom[lev], num_shift, dir); + shiftMF(*pml_rz_E[dim], geom[lev], num_shift, dir); + } +#endif if (lev > 0) { // coarse grid shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, B_external_grid[dim], use_Bparser, Bfield_parser); @@ -407,6 +417,33 @@ WarpX::shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, dstfab(i,j,k,n) = srcfab(i+shift.x,j+shift.y,k+shift.z,n); }) } + +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + if (WarpX::GetInstance().getPMLRZ()) { + // This does the exchange of data in the corner guard cells, the cells that are in the + // guard region both radially and longitudinally. These are the PML cells in the overlapping + // longitudinal region. FillBoundary normally does not update these cells. + // This update is needed so that the cells at the end of the FABs are updated appropriately + // with the data shifted from the nieghboring FAB. Without this update, the RZ PML becomes + // unstable with the moving grid. + // This code creates a temporary MultiFab using a BoxList where the radial size of all of + // its boxes is increased so that the radial guard cells are included in the boxes valid domain. + // The temporary MultiFab is setup to refer to the data of the original Multifab (this can + // be done since the shape of the data is all the same, just the indexing is different). + amrex::BoxList bl; + for (int i = 0, N=ba.size(); i < N; ++i) { + bl.push_back(amrex::grow(ba[i], 0, mf.nGrowVect()[0])); + } + amrex::BoxArray rba(std::move(bl)); + amrex::MultiFab rmf(rba, dm, mf.nComp(), IntVect(0,mf.nGrowVect()[1]), MFInfo().SetAlloc(false)); + + for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) { + rmf.setFab(mfi, FArrayBox(mf[mfi], amrex::make_alias, 0, mf.nComp())); + } + rmf.FillBoundary(false); + } +#endif + } void diff --git a/Source/WarpX.H b/Source/WarpX.H index 8160a428b1e..55ce039aa97 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -25,6 +25,7 @@ #ifdef WARPX_USE_PSATD # ifdef WARPX_DIM_RZ # include "FieldSolver/SpectralSolver/SpectralSolverRZ_fwd.H" +# include "BoundaryConditions/PML_RZ_fwd.H" # else # include "FieldSolver/SpectralSolver/SpectralSolver_fwd.H" # endif @@ -344,6 +345,9 @@ public: const amrex::MultiFab& getBfield_avg_cp (int lev, int direction) {return *Bfield_avg_cp[lev][direction];} bool DoPML () const {return do_pml;} +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + const PML_RZ* getPMLRZ() {return pml_rz[0].get();} +#endif /** get low-high-low-high-... vector for each direction indicating if mother grid PMLs are enabled */ std::vector getPMLdirections() const; @@ -506,8 +510,9 @@ public: void ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType dt_type); void DampPML (); - void DampPML (int lev); - void DampPML (int lev, PatchType patch_type); + void DampPML (const int lev); + void DampPML (const int lev, PatchType patch_type); + void DampPML_Cartesian (const int lev, PatchType patch_type); void DampJPML (); void DampJPML (int lev); @@ -532,6 +537,9 @@ public: void NodalSyncPML (int lev, PatchType patch_type); PML* GetPML (int lev); +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + PML_RZ* GetPML_RZ (int lev); +#endif /** Run the ionization module on all species */ void doFieldIonization (); @@ -1123,6 +1131,9 @@ private: amrex::IntVect do_pml_Lo = amrex::IntVect::TheZeroVector(); amrex::IntVect do_pml_Hi = amrex::IntVect::TheZeroVector(); amrex::Vector > pml; +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + amrex::Vector > pml_rz; +#endif amrex::Real moving_window_x = std::numeric_limits::max(); amrex::Real current_injection_position = 0; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index b2f6480a715..65657341858 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -21,6 +21,7 @@ # include "FieldSolver/SpectralSolver/SpectralKSpace.H" # ifdef WARPX_DIM_RZ # include "FieldSolver/SpectralSolver/SpectralSolverRZ.H" +# include "BoundaryConditions/PML_RZ.H" # else # include "FieldSolver/SpectralSolver/SpectralSolver.H" # endif // RZ ifdef @@ -320,6 +321,9 @@ WarpX::WarpX () charge_buf.resize(nlevs_max); pml.resize(nlevs_max); +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) + pml_rz.resize(nlevs_max); +#endif costs.resize(nlevs_max); load_balance_efficiency.resize(nlevs_max); @@ -807,8 +811,10 @@ WarpX::ReadParameters () } #ifdef WARPX_DIM_RZ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( isAnyBoundaryPML() == false, - "PML are not implemented in RZ geometry; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( isAnyBoundaryPML() == false || maxwell_solver_id == MaxwellSolverAlgo::PSATD, + "PML are not implemented in RZ geometry with FDTD; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( field_boundary_lo[1] != FieldBoundaryType::PML && field_boundary_hi[1] != FieldBoundaryType::PML, + "PML are not implemented in RZ geometry along z; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); #endif if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){ @@ -1481,6 +1487,9 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d WarpX::do_electrostatic, WarpX::do_multi_J, WarpX::fft_do_time_averaging, + WarpX::isAnyBoundaryPML(), + WarpX::do_pml_in_domain, + WarpX::pml_ncell, this->refRatio()); @@ -1759,6 +1768,11 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if ( fft_periodic_single_box == false ) { realspace_ba.grow(1, ngE[1]); // add guard cells only in z } + if (field_boundary_hi[0] == FieldBoundaryType::PML && !do_pml_in_domain) { + // Extend region that is solved for to include the guard cells + // which is where the PML boundary is applied. + realspace_ba.growHi(0, pml_ncell); + } AllocLevelSpectralSolverRZ(spectral_solver_fp, lev, realspace_ba, @@ -1899,6 +1913,11 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // Define spectral solver #ifdef WARPX_DIM_RZ c_realspace_ba.grow(1, ngE[1]); // add guard cells only in z + if (field_boundary_hi[0] == FieldBoundaryType::PML && !do_pml_in_domain) { + // Extend region that is solved for to include the guard cells + // which is where the PML boundary is applied. + c_realspace_ba.growHi(0, pml_ncell); + } AllocLevelSpectralSolverRZ(spectral_solver_cp, lev, c_realspace_ba, @@ -2009,6 +2028,7 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector