Add PECInsulator boundary condition (#4943)

This PR adds a mixed PEC and insulator boundary condition. This allows
an insulator to be placed on a portion of the boundary. The rest of that
boundary will be PEC. Within the insulator portion, the tangential
fields can be specified on the boundary (as functions of space and
time). The normal fields and fields not specified are extrapolated to
the guard cells from the valid cells. The fields are specified in pairs,
the two tangential electric fields, and the two tangential magnetic
fields. In each pair, if one is set, the other will be zeroed if not
set.

A use case is the simulation of a dynamic pinch, driven by an external
current, represented as a time dependent B field on the boundary.


[PECinsulatorBC_warpX_summaryOnly.pdf](https://github.com/user-attachments/files/17637695/PECinsulatorBC_warpX_summaryOnly.pdf)

---------

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

Docs/source/usage/parameters.rst

@@ -493,6 +493,37 @@ Domain Boundary Conditions
 
     * ``pec``: This option can be used to set a Perfect Electric Conductor at the simulation boundary. Please see the :ref:`PEC theory section <theory-bc-pec>` for more details. Note that PEC boundary is invalid at `r=0` for the RZ solver. Please use ``none`` option. This boundary condition does not work with the spectral solver.
 
+    * ``pec_insulator``: This option specifies a mixed perfect electric conductor and insulator boundary, where some part of the
+      boundary is PEC and some is insulator. In the insulator portion, the normal fields are extrapolated and the tangential fields
+      are either set to the specified value or extrapolated. The region that is insulator is specified using a spatially dependent expression with the insulator being in the area where the value of the expression is greater than zero.
+      The expressions are given for the low and high boundary on each axis, as listed below. The tangential fields are specified as
+      expressions that can depend on the location and time. The tangential fields are in two pairs, the electric fields and the
+      magnetic fields. In each pair, if one is specified, the other will be set to zero if not also specified.
+
+      * ``insulator.area_x_lo(y,z)``: For the lower x (or r) boundary, expression specifying the insulator location
+
+      * ``insulator.area_x_hi(y,z)``: For the upper x (or r) boundary, expression specifying the insulator location
+
+      * ``insulator.area_y_lo(x,z)``: For the lower y boundary, expression specifying the insulator location
+
+      * ``insulator.area_y_hi(x,z)``: For the upper y boundary, expression specifying the insulator location
+
+      * ``insulator.area_z_lo(x,y)``: For the lower z boundary, expression specifying the insulator location
+
+      * ``insulator.area_z_hi(x,y)``: For the upper z boundary, expression specifying the insulator location
+
+      * ``insulator.Ey_x_lo(y,z,t)``, ``insulator.Ez_x_lo(y,z,t)``, ``insulator.By_x_lo(y,z,t)``, ``insulator.Bz_x_lo(y,z,t)``: expressions of the tangential field values for the lower x (or r) boundary
+
+      * ``insulator.Ey_x_hi(y,z,t)``, ``insulator.Ez_x_hi(y,z,t)``, ``insulator.By_x_hi(y,z,t)``, ``insulator.Bz_x_hi(y,z,t)``: expressions of the tangential field values for the upper x (or r) boundary
+
+      * ``insulator.Ex_y_lo(x,z,t)``, ``insulator.Ez_y_lo(x,z,t)``, ``insulator.Bx_y_lo(x,z,t)``, ``insulator.Bz_y_lo(x,z,t)``: expressions of the tangential field values for the lower y boundary
+
+      * ``insulator.Ex_y_hi(x,z,t)``, ``insulator.Ez_y_hi(x,z,t)``, ``insulator.Bx_y_hi(x,z,t)``, ``insulator.Bz_y_hi(x,z,t)``: expressions of the tangential field values for the upper y boundary
+
+      * ``insulator.Ex_z_lo(x,y,t)``, ``insulator.Ey_z_lo(x,y,t)``, ``insulator.Bx_z_lo(x,y,t)``, ``insulator.By_z_lo(x,y,t)``: expressions of the tangential field values for the lower z boundary
+
+      * ``insulator.Ex_z_hi(x,y,t)``, ``insulator.Ey_z_hi(x,y,t)``, ``insulator.Bx_z_hi(x,y,t)``, ``insulator.By_z_hi(x,y,t)``: expressions of the tangential field values for the upper z boundary
+
     * ``none``: No boundary condition is applied to the fields with the electromagnetic solver. This option must be used for the RZ-solver at `r=0`.
 
     * ``neumann``: For the electrostatic multigrid solver, a Neumann boundary condition (with gradient of the potential equal to 0) will be applied on the specified boundary.

Examples/Tests/pec/CMakeLists.txt

@@ -30,3 +30,13 @@ add_warpx_test(
     diags/diag1000020  # output
     OFF  # dependency
 )
+
+add_warpx_test(
+    test_2d_pec_field_insulator  # name
+    2  # dims
+    2  # nprocs
+    inputs_test_2d_pec_field_insulator  # inputs
+    analysis_default_regression.py  # analysis
+    diags/diag1000010  # output
+    OFF  # dependency
+)

Examples/Tests/pec/inputs_test_2d_pec_field_insulator

@@ -0,0 +1,34 @@
+# Maximum number of time steps
+max_step = 10
+
+# number of grid points
+amr.n_cell = 32 32
+amr.blocking_factor = 16
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 2
+geometry.prob_lo = 0.     2.e-2 # physical domain
+geometry.prob_hi = 1.e-2  3.e-2
+
+# Boundary condition
+boundary.field_lo = neumann      periodic
+boundary.field_hi = PECInsulator periodic
+
+warpx.serialize_initial_conditions = 1
+
+# Verbosity
+warpx.verbose = 1
+
+# CFL
+warpx.cfl = 1.0
+
+insulator.area_x_hi(y,z) = (2.25e-2 <= z and z <= 2.75e-2)
+insulator.By_x_hi(y,z,t) = min(t/1.0e-12,1)*1.e1*3.3e-4
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 10
+diag1.diag_type = Full

Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator.json

@@ -0,0 +1,13 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.34938851065132936,
+    "Bz": 0.0,
+    "Ex": 31871402.236828588,
+    "Ey": 0.0,
+    "Ez": 104908439.18998256,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0
+  }
+}

Source/BoundaryConditions/CMakeLists.txt

@@ -2,6 +2,7 @@ foreach(D IN LISTS WarpX_DIMS)
     warpx_set_suffix_dims(SD ${D})
     target_sources(lib_${SD}
       PRIVATE
+        PEC_Insulator.cpp
         PML.cpp
         WarpXEvolvePML.cpp
         WarpXFieldBoundaries.cpp

Source/BoundaryConditions/Make.package

@@ -1,3 +1,4 @@
+CEXE_sources += PEC_Insulator.cpp
 CEXE_sources += PML.cpp WarpXEvolvePML.cpp
 CEXE_sources += WarpXFieldBoundaries.cpp WarpX_PEC.cpp
 

Source/BoundaryConditions/PEC_Insulator.H

@@ -0,0 +1,180 @@
+#ifndef PEC_INSULATOR_H_
+#define PEC_INSULATOR_H_
+
+#include "Utils/WarpXAlgorithmSelection.H"
+
+#include <AMReX_Array.H>
+#include <AMReX_Geometry.H>
+#include <AMReX_Vector.H>
+
+#include <AMReX_BaseFwd.H>
+
+#include <array>
+#include <memory>
+
+class PEC_Insulator
+{
+public:
+
+    PEC_Insulator();
+
+    /**
+     * \brief Apply either the PEC or insulator boundary condition on the boundary and in the
+     *        guard cells.
+     *        In the PEC, the nodal fields (in a Yee mesh) are made even relative to the boundary,
+     *        the non-nodal fields are made odd.
+     *        In the insulator, the tangential fields are set to the value if specified, otherwise unchanged,
+     *        and the normal fields extrapolated from the valid cells.
+     *
+     * \param[in,out] Efield
+     * \param[in]     field_boundary_lo   lower field boundary conditions
+     * \param[in]     field_boundary_hi   upper field boundary conditions
+     * \param[in]     ng_fieldgather      number of guard cells used by field gather
+     * \param[in]     geom                geometry object of level "lev"
+     * \param[in]     lev                 level of the Multifab
+     * \param[in]     patch_type          coarse or fine
+     * \param[in]     ref_ratios          vector containing the refinement ratios of the refinement levels
+     * \param[in]     time                current time of the simulation
+     * \param[in]     split_pml_field     whether pml the multifab is the regular Efield or
+     *                                    split pml field
+     */
+    void ApplyPEC_InsulatortoEfield (std::array<amrex::MultiFab*, 3> Efield,
+                                     amrex::Array<FieldBoundaryType,AMREX_SPACEDIM> const & field_boundary_lo,
+                                     amrex::Array<FieldBoundaryType,AMREX_SPACEDIM> const & field_boundary_hi,
+                                     amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom,
+                                     int lev, PatchType patch_type, amrex::Vector<amrex::IntVect> const & ref_ratios,
+                                     amrex::Real time,
+                                     bool split_pml_field = false);
+    /**
+     * \brief Apply either the PEC or insulator boundary condition on the boundary and in the
+     *        guard cells.
+     *        In the PEC, the nodal fields (in a Yee mesh) are made even relative to the boundary,
+     *        the non-nodal fields are made odd.
+     *        In the insulator, the tangential fields are set to the value if specified, otherwise unchanged,
+     *        and the normal fields extrapolated from the valid cells.
+     *
+     * \param[in,out] Bfield
+     * \param[in]     field_boundary_lo   lower field boundary conditions
+     * \param[in]     field_boundary_hi   upper field boundary conditions
+     * \param[in]     ng_fieldgather      number of guard cells used by field gather
+     * \param[in]     geom                geometry object of level "lev"
+     * \param[in]     lev                 level of the Multifab
+     * \param[in]     patch_type          coarse or fine
+     * \param[in]     ref_ratios          vector containing the refinement ratios of the refinement levels
+     * \param[in]     time                current time of the simulation
+     */
+    void ApplyPEC_InsulatortoBfield (std::array<amrex::MultiFab*, 3> Bfield,
+                                     amrex::Array<FieldBoundaryType,AMREX_SPACEDIM> const & field_boundary_lo,
+                                     amrex::Array<FieldBoundaryType,AMREX_SPACEDIM> const & field_boundary_hi,
+                                     amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom,
+                                     int lev, PatchType patch_type, amrex::Vector<amrex::IntVect> const & ref_ratios,
+                                     amrex::Real time);
+
+    /**
+     * \brief The work routine applying the boundary condition
+     *
+     * \param[in,out] field
+     * \param[in]     field_boundary_lo   lower field boundary conditions
+     * \param[in]     field_boundary_hi   upper field boundary conditions
+     * \param[in]     ng_fieldgather      number of guard cells used by field gather
+     * \param[in]     geom                geometry object of level "lev"
+     * \param[in]     lev                 level of the Multifab
+     * \param[in]     patch_type          coarse or fine
+     * \param[in]     ref_ratios          vector containing the refinement ratios of the refinement levels
+     * \param[in]     time                current time of the simulation
+     * \param[in]     split_pml_field     whether pml the multifab is the regular Efield or
+     *                                    split pml field
+     * \param[in]     E_like whether the field is E like or B like
+     * \param[in]     set_F_x_lo whether the tangential field at the boundary was specified
+     * \param[in]     set_F_x_hi whether the tangential field at the boundary was specified
+     * \param[in]     a_Fy_x_lo the parser for the tangential field at the boundary
+     * \param[in]     a_Fz_x_lo the parser for the tangential field at the boundary
+     * \param[in]     a_Fy_x_hi the parser for the tangential field at the boundary
+     * \param[in]     a_Fz_x_hi the parser for the tangential field at the boundary
+     * \param[in]     set_F_y_lo whether the tangential field at the boundary was specified
+     * \param[in]     set_F_y_hi whether the tangential field at the boundary was specified
+     * \param[in]     a_Fx_y_lo the parser for the tangential field at the boundary
+     * \param[in]     a_Fz_y_lo the parser for the tangential field at the boundary
+     * \param[in]     a_Fx_y_hi the parser for the tangential field at the boundary
+     * \param[in]     a_Fz_y_hi the parser for the tangential field at the boundary
+     * \param[in]     set_F_z_lo whether the tangential field at the boundary was specified
+     * \param[in]     set_F_z_hi whether the tangential field at the boundary was specified
+     * \param[in]     a_Fx_z_lo the parser for the tangential field at the boundary
+     * \param[in]     a_Fy_z_lo the parser for the tangential field at the boundary
+     * \param[in]     a_Fx_z_hi the parser for the tangential field at the boundary
+     * \param[in]     a_Fy_z_hi the parser for the tangential field at the boundary
+     */
+    void
+    ApplyPEC_InsulatortoField (std::array<amrex::MultiFab*, 3> field,
+                               amrex::Array<FieldBoundaryType,AMREX_SPACEDIM> const & field_boundary_lo,
+                               amrex::Array<FieldBoundaryType,AMREX_SPACEDIM> const & field_boundary_hi,
+                               amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom,
+                               int lev, PatchType patch_type, amrex::Vector<amrex::IntVect> const & ref_ratios,
+                               amrex::Real time,
+                               bool split_pml_field,
+                               bool E_like,
+#if (AMREX_SPACEDIM > 1)
+                               bool set_F_x_lo, bool set_F_x_hi,
+                               std::unique_ptr<amrex::Parser> const & a_Fy_x_lo, std::unique_ptr<amrex::Parser> const & a_Fz_x_lo,
+                               std::unique_ptr<amrex::Parser> const & a_Fy_x_hi, std::unique_ptr<amrex::Parser> const & a_Fz_x_hi,
+#endif
+#if defined(WARPX_DIM_3D)
+                               bool set_F_y_lo, bool set_F_y_hi,
+                               std::unique_ptr<amrex::Parser> const & a_Fx_y_lo, std::unique_ptr<amrex::Parser> const & a_Fz_y_lo,
+                               std::unique_ptr<amrex::Parser> const & a_Fx_y_hi, std::unique_ptr<amrex::Parser> const & a_Fz_y_hi,
+#endif
+                               bool set_F_z_lo, bool set_F_z_hi,
+                               std::unique_ptr<amrex::Parser> const & a_Fx_z_lo, std::unique_ptr<amrex::Parser> const & a_Fy_z_lo,
+                               std::unique_ptr<amrex::Parser> const & a_Fx_z_hi, std::unique_ptr<amrex::Parser> const & a_Fy_z_hi);
+
+private:
+
+    /* \brief Reads in the parsers for the tangential fields, returning whether
+     *        the input parameter was specified.
+     * \param[in]  pp_insulator ParmParse instance
+     * \param[out] parser the parser generated from the input
+     * \param[in]  input_name the name of the input parameter
+     * \param[in]  coord1 the first coordinate in the plane
+     * \param[in]  coord2 the second coordinate in the plane
+     */
+    bool ReadTangentialFieldParser (amrex::ParmParse const & pp_insulator,
+                                    std::unique_ptr<amrex::Parser> & parser,
+                                    std::string const & input_name,
+                                    std::string const & coord1,
+                                    std::string const & coord2);
+
+    std::vector<std::unique_ptr<amrex::Parser>> m_insulator_area_lo;
+    std::vector<std::unique_ptr<amrex::Parser>> m_insulator_area_hi;
+
+#if (AMREX_SPACEDIM > 1)
+    bool m_set_B_x_lo = false, m_set_B_x_hi = false;
+    std::unique_ptr<amrex::Parser> m_By_x_lo, m_Bz_x_lo;
+    std::unique_ptr<amrex::Parser> m_By_x_hi, m_Bz_x_hi;
+#endif
+#if defined(WARPX_DIM_3D)
+    bool m_set_B_y_lo = false, m_set_B_y_hi = false;
+    std::unique_ptr<amrex::Parser> m_Bx_y_lo, m_Bz_y_lo;
+    std::unique_ptr<amrex::Parser> m_Bx_y_hi, m_Bz_y_hi;
+#endif
+    bool m_set_B_z_lo = false, m_set_B_z_hi = false;
+    std::unique_ptr<amrex::Parser> m_Bx_z_lo, m_By_z_lo;
+    std::unique_ptr<amrex::Parser> m_Bx_z_hi, m_By_z_hi;
+
+
+#if (AMREX_SPACEDIM > 1)
+    bool m_set_E_x_lo = false, m_set_E_x_hi = false;
+    std::unique_ptr<amrex::Parser> m_Ey_x_lo, m_Ez_x_lo;
+    std::unique_ptr<amrex::Parser> m_Ey_x_hi, m_Ez_x_hi;
+#endif
+#if defined(WARPX_DIM_3D)
+    bool m_set_E_y_lo = false, m_set_E_y_hi = false;
+    std::unique_ptr<amrex::Parser> m_Ex_y_lo, m_Ez_y_lo;
+    std::unique_ptr<amrex::Parser> m_Ex_y_hi, m_Ez_y_hi;
+#endif
+    bool m_set_E_z_lo = false, m_set_E_z_hi = false;
+    std::unique_ptr<amrex::Parser> m_Ex_z_lo, m_Ey_z_lo;
+    std::unique_ptr<amrex::Parser> m_Ex_z_hi, m_Ey_z_hi;
+
+
+};
+#endif // PEC_INSULATOR_H_