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, 42 insertions, 0 deletions

diff --git a/Examples/Tests/ElectrostaticDirichletBC/analysis.py b/Examples/Tests/ElectrostaticDirichletBC/analysis.py
new file mode 100755
index 000000000..b08c993c9
--- /dev/null
+++ b/Examples/Tests/ElectrostaticDirichletBC/analysis.py
@@ -0,0 +1,42 @@
+#! /usr/bin/env python
+
+# Copyright 2021 Roelof Groenewald
+
+# This script tests the time dependent Dirichlet boundary
+# conditions in a 2D electrostatic simulation. An empty
+# domain of 64 x 8 cells is simulated with periodic boundary
+# conditions in the x directions and Dirichlet boundary
+# conditions in the y direction with specified potentials
+# of sine waves with different periods on the lo and hi side.
+# One period of the hi side sine wave is simulated and the
+# potentials at the boundaries compared to expectation.
+
+# Possible running time: ~ 19 s
+
+import numpy as np
+
+import yt
+import glob
+
+files = sorted(glob.glob('dirichletbc_plt*'))[1:]
+
+times = np.ones(len(files))
+potentials_lo = np.zeros(len(files))
+potentials_hi = np.zeros(len(files))
+
+for ii, file in enumerate(files):
+    ds = yt.load( file )
+    times[ii] = (
+        ds.current_time.item() - float(ds.parameters.get('warpx.const_dt'))
+    )
+    data = ds.covering_grid(
+        level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions
+    )
+    potentials_lo[ii] = np.mean(data['phi'].to_ndarray()[0])
+    potentials_hi[ii] = np.mean(data['phi'].to_ndarray()[-1])
+
+expected_potentials_lo = 150.0 * np.sin(2.0 * np.pi * 6.78e6 * times)
+expected_potentials_hi = 450.0 * np.sin(2.0 * np.pi * 13.56e6 * times)
+
+assert np.allclose(potentials_lo, expected_potentials_lo, rtol=0.1)
+assert np.allclose(potentials_hi, expected_potentials_hi, rtol=0.1)