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>

4 files changed, 40 insertions, 0 deletions

diff --git a/Python/pywarpx/Boundary.py b/Python/pywarpx/Boundary.py
new file mode 100644
index 000000000..b7ebd17af
--- /dev/null
+++ b/Python/pywarpx/Boundary.py
@@ -0,0 +1,9 @@
+# Copyright 2021 Roelof Groenewald
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+from .Bucket import Bucket
+
+boundary = Bucket('boundary')
diff --git a/Python/pywarpx/WarpX.py b/Python/pywarpx/WarpX.py
index ae48bfcf6..dc09ba241 100644
--- a/Python/pywarpx/WarpX.py
+++ b/Python/pywarpx/WarpX.py
@@ -9,6 +9,7 @@ from .Bucket import Bucket
 from .Constants import my_constants
 from .Amr import amr
 from .Geometry import geometry
+from .Boundary import boundary
 from .Algo import algo
 from .Langmuirwave import langmuirwave
 from .Interpolation import interpolation
@@ -30,6 +31,7 @@ class WarpX(Bucket):
         argv += my_constants.attrlist()
         argv += amr.attrlist()
         argv += geometry.attrlist()
+        argv += boundary.attrlist()
         argv += algo.attrlist()
         argv += langmuirwave.attrlist()
         argv += interpolation.attrlist()
@@ -89,6 +91,7 @@ class WarpX(Bucket):
 
     def write_inputs(self, filename='inputs', **kw):
         argv = self.create_argv_list()
+
         with open(filename, 'w') as ff:
 
             for k, v in kw.items():
diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py
index ed2152c1d..918c300b2 100644
--- a/Python/pywarpx/__init__.py
+++ b/Python/pywarpx/__init__.py
@@ -7,6 +7,7 @@
 from .WarpX import warpx
 from .Constants import my_constants
 from .Amr import amr
+from .Boundary import boundary
 from .Geometry import geometry
 from .Algo import algo
 from .Langmuirwave import langmuirwave
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index d93396a82..6d1ace18e 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -355,6 +355,13 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid):
         self.max_grid_size = kw.pop('warpx_max_grid_size', 32)
         self.blocking_factor = kw.pop('warpx_blocking_factor', None)
 
+        self.potential_xmin = None
+        self.potential_xmax = None
+        self.potential_ymin = None
+        self.potential_ymax = None
+        self.potential_zmin = kw.pop('warpx_potential_lo_z', None)
+        self.potential_zmax = kw.pop('warpx_potential_hi_z', None)
+
     def initialize_inputs(self):
         pywarpx.amr.n_cell = self.number_of_cells
 
@@ -398,6 +405,13 @@ class Cartesian2DGrid(picmistandard.PICMI_Cartesian2DGrid):
         self.max_grid_size = kw.pop('warpx_max_grid_size', 32)
         self.blocking_factor = kw.pop('warpx_blocking_factor', None)
 
+        self.potential_xmin = kw.pop('warpx_potential_lo_x', None)
+        self.potential_xmax = kw.pop('warpx_potential_hi_x', None)
+        self.potential_ymin = None
+        self.potential_ymax = None
+        self.potential_zmin = kw.pop('warpx_potential_lo_z', None)
+        self.potential_zmax = kw.pop('warpx_potential_hi_z', None)
+
     def initialize_inputs(self):
         pywarpx.amr.n_cell = self.number_of_cells
 
@@ -437,6 +451,13 @@ class Cartesian3DGrid(picmistandard.PICMI_Cartesian3DGrid):
         self.max_grid_size = kw.pop('warpx_max_grid_size', 32)
         self.blocking_factor = kw.pop('warpx_blocking_factor', None)
 
+        self.potential_xmin = kw.pop('warpx_potential_lo_x', None)
+        self.potential_xmax = kw.pop('warpx_potential_hi_x', None)
+        self.potential_ymin = kw.pop('warpx_potential_lo_y', None)
+        self.potential_ymax = kw.pop('warpx_potential_hi_y', None)
+        self.potential_zmin = kw.pop('warpx_potential_lo_z', None)
+        self.potential_zmax = kw.pop('warpx_potential_hi_z', None)
+
     def initialize_inputs(self):
         pywarpx.amr.n_cell = self.number_of_cells
 
@@ -545,6 +566,12 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver):
             pywarpx.warpx.do_electrostatic = 'labframe'
             pywarpx.warpx.self_fields_required_precision = self.required_precision
             pywarpx.warpx.self_fields_max_iters = self.maximum_iterations
+            pywarpx.boundary.potential_lo_x = self.grid.potential_xmin
+            pywarpx.boundary.potential_lo_y = self.grid.potential_ymin
+            pywarpx.boundary.potential_lo_z = self.grid.potential_zmin
+            pywarpx.boundary.potential_hi_x = self.grid.potential_xmax
+            pywarpx.boundary.potential_hi_y = self.grid.potential_ymax
+            pywarpx.boundary.potential_hi_z = self.grid.potential_zmax
 
 
 class GaussianLaser(picmistandard.PICMI_GaussianLaser):