aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Docs/source/usage/parameters.rst4
-rwxr-xr-x[-rw-r--r--]Examples/Modules/embedded_boundary_cube/analysis_fields.py0
-rw-r--r--Examples/Modules/relativistic_space_charge_initialization/inputs_3d3
-rw-r--r--Examples/Modules/space_charge_initialization/inputs_3d3
-rwxr-xr-xExamples/Tests/ElectrostaticDirichletBC/analysis.py42
-rw-r--r--Examples/Tests/ElectrostaticDirichletBC/inputs_2d24
-rw-r--r--Python/pywarpx/Boundary.py9
-rw-r--r--Python/pywarpx/WarpX.py3
-rw-r--r--Python/pywarpx/__init__.py1
-rw-r--r--Python/pywarpx/picmi.py27
-rw-r--r--Regression/WarpX-tests.ini16
-rw-r--r--Source/FieldSolver/ElectrostaticSolver.cpp193
-rw-r--r--Source/WarpX.H8
13 files changed, 315 insertions, 18 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 8453d9f25..931fcfd54 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -80,7 +80,7 @@ Overall simulation parameters
The relative precision with which the electrostatic space-charge fields should
be calculated. More specifically, the space-charge fields are
computed with an iterative Multi-Level Multi-Grid (MLMG) solver.
- This solver can fail to reach the default precision within a reasonable
+ This solver can fail to reach the default precision within a reasonable time.
This only applies when warpx.do_electrostatic = labframe.
* ``warpx.self_fields_max_iters`` (`integer`, default: 200)
@@ -204,6 +204,7 @@ Domain Boundary Conditions
Options are:
* ``Periodic``: This option can be used to set periodic domain boundaries. Note that if the fields for lo in a certain dimension are set to periodic, then the corresponding upper boundary must also be set to periodic. If particle boundaries are not specified in the input file, then particles boundaries by default will be set to periodic. If particles boundaries are specified, then they must be set to periodic corresponding to the periodic field boundaries.
+ * ``pec``: This option can be used to add a Perfect Electric Conductor at the simulation boundary. If an electrostatic field solve is used the boundary potentials can also be set through ``boundary.potential_lo_x/y/z`` and ``boundary.potential_hi_x/y/z`` (default `0`).
* ``pml`` (default): This option can be used to add Perfectly Matched Layers (PML) around the simulation domain. It will override the user-defined value provided for ``warpx.do_pml``. See the :ref:`PML theory section <theory-bc>` for more details.
Additional pml algorithms can be explored using the parameters ``warpx.do_pml_in_domain``, ``warpx.do_particles_in_pml``, and ``warpx.do_pml_j_damping``.
@@ -1106,7 +1107,6 @@ following the algorithm given by `Perez et al. (Phys. Plasmas 19, 083104, 2012)
:math:`R\approx1.4A^{1/3}` is the effective Coulombic radius of the nucleus,
:math:`A` is the mass number.
If this is not provided, or if a non-positive value is provided,
- a Coulomb logarithm will be computed automatically according to the algorithm.
a Coulomb logarithm will be computed automatically according to the algorithm in
`Perez et al. (Phys. Plasmas 19, 083104, 2012) <https://doi.org/10.1063/1.4742167>`_.
diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py
index 58dad1bb7..58dad1bb7 100644..100755
--- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py
+++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py
diff --git a/Examples/Modules/relativistic_space_charge_initialization/inputs_3d b/Examples/Modules/relativistic_space_charge_initialization/inputs_3d
index 457fabcbd..d6758d160 100644
--- a/Examples/Modules/relativistic_space_charge_initialization/inputs_3d
+++ b/Examples/Modules/relativistic_space_charge_initialization/inputs_3d
@@ -4,7 +4,8 @@ amr.max_grid_size = 32
amr.max_level = 0
geometry.coord_sys = 0 # 0: Cartesian
-geometry.is_periodic = 0 0 0 # Is periodic?
+boundary.field_lo = pec pec pec
+boundary.field_hi = pec pec pec
geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6
geometry.prob_hi = 50.e-6 50.e-6 50.e-6
diff --git a/Examples/Modules/space_charge_initialization/inputs_3d b/Examples/Modules/space_charge_initialization/inputs_3d
index c6895630a..4fe2fbabb 100644
--- a/Examples/Modules/space_charge_initialization/inputs_3d
+++ b/Examples/Modules/space_charge_initialization/inputs_3d
@@ -4,7 +4,8 @@ amr.max_grid_size = 32
amr.max_level = 0
geometry.coord_sys = 0 # 0: Cartesian
-geometry.is_periodic = 0 0 0 # Is periodic?
+boundary.field_lo = pec pec pec
+boundary.field_hi = pec pec pec
geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6
geometry.prob_hi = 50.e-6 50.e-6 50.e-6
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)
diff --git a/Examples/Tests/ElectrostaticDirichletBC/inputs_2d b/Examples/Tests/ElectrostaticDirichletBC/inputs_2d
new file mode 100644
index 000000000..4b3327d9a
--- /dev/null
+++ b/Examples/Tests/ElectrostaticDirichletBC/inputs_2d
@@ -0,0 +1,24 @@
+max_step = 100
+warpx.verbose = 0
+warpx.const_dt = 7.5e-10
+warpx.do_electrostatic = labframe
+warpx.self_fields_required_precision = 1e-06
+amr.n_cell = 64 8
+amr.max_grid_size = 32
+amr.max_level = 0
+
+geometry.coord_sys = 0
+geometry.prob_lo = 0.0 0.0
+geometry.prob_hi = 0.032 0.004
+boundary.field_lo = pec periodic
+boundary.field_hi = pec periodic
+boundary.potential_lo_x = 150.0*sin(2*pi*6.78e+06*t)
+boundary.potential_hi_x = 450.0*sin(2*pi*13.56e+06*t)
+
+interpolation.nox = 1
+interpolation.noy = 1
+interpolation.noz = 1
+diagnostics.diags_names = diag1
+diag1.diag_type = Full
+diag1.intervals = ::4
+diag1.fields_to_plot = phi
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):
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index bbbaf5641..d0b79ca3f 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -2034,3 +2034,19 @@ numthreads = 1
compileTest = 0
doVis = 0
compareParticles = 0
+
+[dirichletbc]
+buildDir = .
+inputFile = Examples/Tests/ElectrostaticDirichletBC/inputs_2d
+runtime_params =
+dim = 2
+addToCompileString =
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 0
+numthreads = 0
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/ElectrostaticDirichletBC/analysis.py
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`.
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 8663fd9c4..7d66943d5 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -596,6 +596,14 @@ public:
amrex::Real const required_precision,
int const max_iters) const;
+ void setPhiBC (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi,
+ std::array<bool,AMREX_SPACEDIM> dirichlet_flag,
+ amrex::Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_lo,
+ amrex::Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_hi
+ ) const;
+ void getPhiBC (const int idim, amrex::Real &pot_lo,
+ amrex::Real &pot_hi) const;
+
void computeE (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E,
const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi,
std::array<amrex::Real, 3> const beta = {{0,0,0}} ) const;