Feature - Time dependent Dirichlet boundary conditions for electrostatic simulations (#1761)

* Update copyright notices

* allow specification of boundary potentials at runtime when using Dirichlet boundary conditions in the electrostatic solver (labframe)

* added parsing to boundary potentials specified at runtime to allow time dependence through a mathematical expression with t (time)

* updated to picmistandard 0.0.14 in order to set the electrostatic solver convergence threshold

* update docs

* various changes requested during PR review

* fixed issue causing old tests to break and added a new test for time varying boundary potentials

* possibly a fix for the failed time varying boundary condition test

* changed permission on the analysis file for the time varying BCs test

* switched to using yt for data analysis since h5py is not available

* made changes compatible with PR#1730; changed potential boundary setting routine to use the ParallelFor construct and set all boundaries in a single call

* fixed typo in computePhiRZ

* updated docs and fixed other minor typos

* fixed bug in returning from loop over dimensions when setting boundary potentials rather than continuing

* changed to setting potentials on domain boundaries rather than tilebox boundaries and changed picmi.py to accept boundary potentials

* now using domain.surroundingNodes() to get the proper boundary cells for the physical domain

* fixed typo in variable name specifying z-boundary potential

* changed boundary value parameter for Dirichlet BC to boundary.field_lo/hi and changed setPhiBC() to only loop over the grid points when a boundary value has changed

* switched specifying potential boundary values though individual inputs of the form boundary.potential_lo/hi_x/y/z and incorporated the new BC formalism through FieldBoundaryType::Periodic and FieldBoundaryType::PEC rather than Geom(0).isPeriodic(idim)

* removed incorrect check of whether the potential boundary values are already correct, also had to change the input to test space_charge_initialization_2d to comply with the new boundary condition input parameters and finally changed permissions to analysis_fields.py file for the embedded boundary test since it was failing

* remove line from WarpX-tests.ini that was incorrectly added during upstream merge

* changed input file for relativistic space charge initialization to new boundary condition specification

* fixed outdated comment and updated documentation to reflect that the Dirichlet BCs can also be specified when using the relativistic electrostatic field solver

* moved call to get domain boundaries inside the loop over levels

* cleaned up the code some by using domain.smallEnd and domain.bigEnd rather than lbound and ubound

* added check that a box contains boundary cells before launching a loop over that box cells to set the boundary conditions

Co-authored-by: Peter Scherpelz <peter.scherpelz@modernelectron.com>

1 files changed, 179 insertions, 14 deletions

diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp
index ab2fc2b25..7aea1d4e4 100644
--- a/Source/FieldSolver/ElectrostaticSolver.cpp
+++ b/Source/FieldSolver/ElectrostaticSolver.cpp
@@ -266,15 +266,35 @@ WarpX::computePhiRZ (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho
     Array<LinOpBCType,AMREX_SPACEDIM> lobc, hibc;
     lobc[0] = LinOpBCType::Neumann;
     hibc[0] = LinOpBCType::Dirichlet;
-    if ( Geom(0).isPeriodic(1) ) {
+    std::array<bool,AMREX_SPACEDIM> dirichlet_flag;
+    dirichlet_flag[0] = false;
+    Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_lo, phi_bc_values_hi;
+    if ( WarpX::field_boundary_lo[1] == FieldBoundaryType::Periodic
+         && WarpX::field_boundary_hi[1] == FieldBoundaryType::Periodic ) {
         lobc[1] = LinOpBCType::Periodic;
         hibc[1] = LinOpBCType::Periodic;
-    } else {
+        dirichlet_flag[1] = false;
+    } else if ( WarpX::field_boundary_lo[1] == FieldBoundaryType::PEC
+         && WarpX::field_boundary_hi[1] == FieldBoundaryType::PEC ) {
         // Use Dirichlet boundary condition by default.
         // Ideally, we would often want open boundary conditions here.
         lobc[1] = LinOpBCType::Dirichlet;
         hibc[1] = LinOpBCType::Dirichlet;
+
+        // set flag so we know which dimensions to fix the potential for
+        dirichlet_flag[1] = true;
+        // parse the input file for the potential at the current time
+        getPhiBC(1, phi_bc_values_lo[1], phi_bc_values_hi[1]);
     }
+    else {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+            "Field boundary conditions have to be either periodic or PEC "
+            "when using the electrostatic solver"
+        );
+    }
+
+    // set the boundary potential values if needed
+    setPhiBC(phi, dirichlet_flag, phi_bc_values_lo, phi_bc_values_hi);
 
     // Define the linear operator (Poisson operator)
     MLNodeLaplacian linop( geom_scaled, boxArray(), dmap );
@@ -282,23 +302,22 @@ WarpX::computePhiRZ (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho
         linop.setSigma( lev, *sigma[lev] );
     }
 
+    for (int lev=0; lev < rho.size(); lev++){
+        rho[lev]->mult(-1._rt/PhysConst::ep0);
+    }
+
     // Solve the Poisson equation
     linop.setDomainBC( lobc, hibc );
     MLMG mlmg(linop);
     mlmg.setVerbose(2);
     mlmg.setMaxIter(max_iters);
     mlmg.solve( GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), required_precision, 0.0);
-
-    // Normalize by the correct physical constant
-    for (int lev=0; lev < rho.size(); lev++){
-        phi[lev]->mult(-1._rt/PhysConst::ep0);
-    }
 }
 
 #else
 
 /* Compute the potential `phi` in Cartesian geometry by solving the Poisson equation
-   hwith `rho` as a source, assuming that the source moves at a constant
+   with `rho` as a source, assuming that the source moves at a constant
    speed \f$\vec{\beta}\f$.
    This uses the amrex solver.
 
@@ -321,20 +340,39 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab>
 
     // Define the boundary conditions
     Array<LinOpBCType,AMREX_SPACEDIM> lobc, hibc;
+    std::array<bool,AMREX_SPACEDIM> dirichlet_flag;
+    Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_lo, phi_bc_values_hi;
     for (int idim=0; idim<AMREX_SPACEDIM; idim++){
-        if ( Geom(0).isPeriodic(idim) ) {
+        if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::Periodic
+             && WarpX::field_boundary_hi[idim] == FieldBoundaryType::Periodic ) {
             lobc[idim] = LinOpBCType::Periodic;
             hibc[idim] = LinOpBCType::Periodic;
-        } else {
-            // Use Dirichlet boundary condition by default.
+            dirichlet_flag[idim] = false;
+        } else if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC
+             && WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC ) {
             // Ideally, we would often want open boundary conditions here.
             lobc[idim] = LinOpBCType::Dirichlet;
             hibc[idim] = LinOpBCType::Dirichlet;
+
+            // set flag so we know which dimensions to fix the potential for
+            dirichlet_flag[idim] = true;
+            // parse the input file for the potential at the current time
+            getPhiBC(idim, phi_bc_values_lo[idim], phi_bc_values_hi[idim]);
+        }
+        else {
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                "Field boundary conditions have to be either periodic or PEC "
+                "when using the electrostatic solver"
+            );
         }
     }
 
+    // set the boundary potential values if needed
+    setPhiBC(phi, dirichlet_flag, phi_bc_values_lo, phi_bc_values_hi);
+
     // Define the linear operator (Poisson operator)
     MLNodeTensorLaplacian linop( Geom(), boxArray(), DistributionMap() );
+
     // Set the value of beta
     amrex::Array<amrex::Real,AMREX_SPACEDIM> beta_solver =
 #if (AMREX_SPACEDIM==2)
@@ -346,18 +384,145 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab>
 
     // Solve the Poisson equation
     linop.setDomainBC( lobc, hibc );
+
+    for (int lev=0; lev < rho.size(); lev++){
+        rho[lev]->mult(-1._rt/PhysConst::ep0);
+    }
+
     MLMG mlmg(linop);
     mlmg.setVerbose(2);
     mlmg.setMaxIter(max_iters);
     mlmg.solve( GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), required_precision, 0.0);
+}
+#endif
 
-    // Normalize by the correct physical constant
-    for (int lev=0; lev < rho.size(); lev++){
-        phi[lev]->mult(-1._rt/PhysConst::ep0);
+/* \bried Set Dirichlet boundary conditions for the electrostatic solver.
+
+    The given potential's values are fixed on the boundaries of the given
+    dimension according to the desired values from the simulation input file,
+    boundary.potential_lo and boundary.potential_hi.
+
+   \param[inout] phi The electrostatic potential
+   \param[in] idim The dimension for which the Dirichlet boundary condition is set
+*/
+void
+WarpX::setPhiBC( amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi,
+                 std::array<bool,AMREX_SPACEDIM> dirichlet_flag,
+                 Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_lo,
+                 Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_hi ) const
+{
+    // check if any dimension has Dirichlet boundary conditions
+    bool has_Dirichlet = false;
+    for (int idim=0; idim<AMREX_SPACEDIM; idim++){
+        if (dirichlet_flag[idim]) {
+            has_Dirichlet = true;
+        }
     }
+    if (!has_Dirichlet) return;
+
+    // loop over all mesh refinement levels and set the boundary values
+    for (int lev=0; lev <= max_level; lev++) {
+
+        amrex::Box domain = Geom(lev).Domain();
+        domain.surroundingNodes();
+
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for ( MFIter mfi(*phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+            // Extract the potential
+            auto phi_arr = phi[lev]->array(mfi);
+            // Extract tileboxes for which to loop
+            const Box& tb  = mfi.tilebox( phi[lev]->ixType().toIntVect() );
+
+            // loop over dimensions
+            for (int idim=0; idim<AMREX_SPACEDIM; idim++){
+                // check if the boundary in this dimension should be set
+                if (!dirichlet_flag[idim]) continue;
+
+                // a check can be added below to test if the boundary values
+                // are already correct, in which case the ParallelFor over the
+                // cells can be skipped
+
+                if (!domain.strictly_contains(tb)) {
+                    amrex::ParallelFor( tb,
+                        [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+
+                            IntVect iv(AMREX_D_DECL(i,j,k));
+
+                            if (iv[idim] == domain.smallEnd(idim)){
+                                phi_arr(i,j,k) = phi_bc_values_lo[idim];
+                            }
+                            if (iv[idim] == domain.bigEnd(idim)) {
+                                phi_arr(i,j,k) = phi_bc_values_hi[idim];
+                            }
+
+                        } // loop ijk
+                    );
+                }
+            } // idim
+    }} // lev & MFIter
 }
+
+/* \bried Utility function to parse input file for boundary potentials.
+
+    The input values are parsed to allow math expressions for the potentials
+    that specify time dependence.
+
+   \param[in] idim The dimension for which the potential is queried
+   \param[inout] pot_lo The specified value of `phi` on the lower boundary.
+   \param[inout] pot_hi The specified value of `phi` on the upper boundary.
+*/
+void
+WarpX::getPhiBC( const int idim, amrex::Real &pot_lo, amrex::Real &pot_hi ) const
+{
+    // set default potentials to zero in order for current tests to pass
+    // but forcing the user to specify a potential might be better
+    std::string potential_lo_str = "0";
+    std::string potential_hi_str = "0";
+
+    // Get the boundary potentials specified in the simulation input file
+    // first as strings and then parse those for possible math expressions
+    ParmParse pp_boundary("boundary");
+
+#ifdef WARPX_DIM_RZ
+    if (idim == 1) {
+        pp_boundary.query("potential_lo_z", potential_lo_str);
+        pp_boundary.query("potential_hi_z", potential_hi_str);
+    }
+    else {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+            "Field boundary condition values can currently only be specified "
+            "for z when using RZ geometry."
+        );
+    }
+#else
+    if (idim == 0) {
+        pp_boundary.query("potential_lo_x", potential_lo_str);
+        pp_boundary.query("potential_hi_x", potential_hi_str);
+    }
+    else if (idim == 1){
+        if (AMREX_SPACEDIM == 2){
+            pp_boundary.query("potential_lo_z", potential_lo_str);
+            pp_boundary.query("potential_hi_z", potential_hi_str);
+        }
+        else {
+            pp_boundary.query("potential_lo_y", potential_lo_str);
+            pp_boundary.query("potential_hi_y", potential_hi_str);
+        }
+    }
+    else {
+        pp_boundary.query("potential_lo_z", potential_lo_str);
+        pp_boundary.query("potential_hi_z", potential_hi_str);
+    }
 #endif
 
+    auto parser_lo = makeParser(potential_lo_str, {"t"});
+    pot_lo = parser_lo.eval(gett_new(0));
+    auto parser_hi = makeParser(potential_hi_str, {"t"});
+    pot_hi = parser_hi.eval(gett_new(0));
+}
+
 /* \bried Compute the electric field that corresponds to `phi`, and
           add it to the set of MultiFab `E`.