diff options
33 files changed, 1136 insertions, 162 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 50e4d5984..10cd401d2 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -205,13 +205,16 @@ Domain Boundary Conditions * ``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``. * ``damped``: This is the recommended option in the moving direction when using the spectral solver with moving window (currently only supported along z). This boundary condition applies a damping factor to the electric and magnetic fields in the outer half of the guard cells, using a sine squared profile. As the spectral solver is by nature periodic, the damping prevents fields from wrapping around to the other end of the domain when the periodicity is not desired. This boundary condition is only valid when using the spectral solver. + * ``pec``: This option can be used to set a Perfect Electric Conductor at the simulation boundary. For the electromagnetic solve, at PEC, the tangential electric field and the normal magnetic field are set to 0. This boundary can be used to model a dielectric or metallic surface. In the guard-cell region, the tangential electric field is set equal and opposite to the respective field component in the mirror location across the PEC boundary, and the normal electric field is set equal to the field component in the mirror location in the domain across the PEC boundary. Similarly, the tangential (and normal) magnetic field components are set equal (and opposite) to the respective magnetic field components in the mirror locations across the PEC boundary. 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. +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`). + + * ``none``: No boundary condition is applied to the fields. This option must be used for the RZ-solver at `r=0`. + * ``boundary.particle_lo`` and ``boundary.particle_hi`` (`2 strings` for 2D, `3 strings` for 3D) Options are: * ``Periodic``: Particles leaving the boundary will re-enter from the opposite boundary. The field boundary condition must be consistenly set to periodic and both lower and upper boundaries must be periodic. diff --git a/Examples/Tests/Langmuir/inputs_2d_multi_rz_rt b/Examples/Tests/Langmuir/inputs_2d_multi_rz_rt index e86e7ef9d..78e0a3694 100644 --- a/Examples/Tests/Langmuir/inputs_2d_multi_rz_rt +++ b/Examples/Tests/Langmuir/inputs_2d_multi_rz_rt @@ -13,9 +13,10 @@ amr.max_level = 0 # Geometry geometry.coord_sys = 1 # 0: Cartesian -geometry.is_periodic = 0 1 # Is periodic? geometry.prob_lo = 0.e-6 -20.e-6 # physical domain geometry.prob_hi = 20.e-6 20.e-6 +boundary.field_lo = none periodic +boundary.field_hi = none periodic warpx.serialize_ics = 1 diff --git a/Examples/Tests/PEC/analysis_pec.py b/Examples/Tests/PEC/analysis_pec.py new file mode 100755 index 000000000..3257d8e4e --- /dev/null +++ b/Examples/Tests/PEC/analysis_pec.py @@ -0,0 +1,84 @@ +#! /usr/bin/env python + +# +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + + +# This is a script that analyses the simulation results from +# the script `inputs_field_PEC_3d`. This simulates a sinusoindal wave. +# The electric field (Ey) is a standing wave due to the PEC boundary condition, +# and as a result, the minimum and maximum value after reflection would be two times the value at initialization due to constructive interference. +# Additionally, the value of Ey at the boundary must be equal to zero. +import sys +import re +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import yt +yt.funcs.mylog.setLevel(50) +import numpy as np +from scipy.constants import e, m_e, epsilon_0, c +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# this will be the name of the plot file +fn = sys.argv[1] + +# Parameters (these parameters must match the parameters in `inputs.multi.rt`) +Ey_in = 1.e5 +E_th = 2.0*Ey_in + +# Read the file +ds = yt.load(fn) + +t0 = ds.current_time.to_value() +data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, + dims=ds.domain_dimensions) +Ey_array = data['Ey'].to_ndarray() +max_Ey_sim = Ey_array.max() +min_Ey_sim = Ey_array.min() +max_Ey_error = abs(max_Ey_sim-E_th)/abs(E_th) +min_Ey_error = abs(min_Ey_sim - (-E_th))/abs(E_th) +print('Ey_max is %s: Max Eth is %s : Max error for Ey_max: %.2e' %(max_Ey_sim,E_th,max_Ey_error)) +print('Ey_min is %s: Min Eth is %s : Max error for Ey_min: %.2e' %(min_Ey_sim,(-E_th), min_Ey_error)) +error_rel = 0. +max_Ey_error_rel = max( error_rel, max_Ey_error ) +min_Ey_error_rel = max( error_rel, min_Ey_error ) + +# Plot the last field from the loop (Ez at iteration 40) +field = 'Ey' +xCell = ds.domain_dimensions[0] +yCell = ds.domain_dimensions[1] +zCell = ds.domain_dimensions[2] +z_array = data['z'].to_ndarray() + +plt.figure(figsize=(8,4)) +plt.plot(z_array[int(xCell/2),int(yCell/2),:]-z_array[int(xCell/2),int(yCell/2),int(zCell/2)],Ey_array[int(xCell/2),int(yCell/2),:]) +plt.ylim(-1.2e-3, 1.2e-3) +plt.ylim(-2.2e5, 2.2e5) +plt.xlim(-4.0e-6, 4.000001e-6) +plt.xticks(np.arange(-4.e-6,4.000001e-6, step=1e-6)) +plt.xlabel('z (m)') +plt.ylabel(field +' (V/m)') +plt.tight_layout() +plt.savefig('Ey_pec_analysis.png') + + +tolerance_rel = 0.01 + +print("min_Ey_error_rel : " + str(min_Ey_error_rel)) +print("max_Ey_error_rel : " + str(max_Ey_error_rel)) +print("tolerance_rel: " + str(tolerance_rel)) + +assert( max_Ey_error_rel < tolerance_rel ) +assert( min_Ey_error_rel < tolerance_rel ) + +test_name = fn[:-9] # Could also be os.path.split(os.getcwd())[1] + +if re.search( 'single_precision', fn ): + checksumAPI.evaluate_checksum(test_name, fn, rtol=1.e-3) +else: + checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Tests/PEC/analysis_pec_mr.py b/Examples/Tests/PEC/analysis_pec_mr.py new file mode 100755 index 000000000..067f2e658 --- /dev/null +++ b/Examples/Tests/PEC/analysis_pec_mr.py @@ -0,0 +1,84 @@ +#! /usr/bin/env python + +# +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + + +# This is a script that analyses the simulation results from +# the script `inputs_field_PEC_mr_3d`. This simulates a sinusoindal wave. +# The electric field (Ey) is a standing wave due to the PEC boundary condition, +# and as a result, the minimum and maximum value after reflection would be two times the value at initialization due to constructive interference. +# Additionally, the value of Ey at the boundary must be equal to zero. +import sys +import re +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import yt +yt.funcs.mylog.setLevel(50) +import numpy as np +from scipy.constants import e, m_e, epsilon_0, c +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# this will be the name of the plot file +fn = sys.argv[1] + +# Parameters (these parameters must match the parameters in `inputs.multi.rt`) +Ey_in = 1.e5 +E_th = 2.0*Ey_in + +# Read the file +ds = yt.load(fn) + +t0 = ds.current_time.to_value() +data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, + dims=ds.domain_dimensions) +Ey_array = data['Ey'].to_ndarray() +max_Ey_sim = Ey_array.max() +min_Ey_sim = Ey_array.min() +max_Ey_error = abs(max_Ey_sim-E_th)/abs(E_th) +min_Ey_error = abs(min_Ey_sim - (-E_th))/abs(E_th) +print('Ey_max is %s: Max Eth is %s : Max error for Ey_max: %.2e' %(max_Ey_sim,E_th,max_Ey_error)) +print('Ey_min is %s: Min Eth is %s : Max error for Ey_min: %.2e' %(min_Ey_sim,(-E_th), min_Ey_error)) +error_rel = 0. +max_Ey_error_rel = max( error_rel, max_Ey_error ) +min_Ey_error_rel = max( error_rel, min_Ey_error ) + +# Plot the last field from the loop (Ez at iteration 40) +# Plot the last field from the loop (Ez at iteration 40) +field = 'Ey' +xCell = ds.domain_dimensions[0] +yCell = ds.domain_dimensions[1] +zCell = ds.domain_dimensions[2] +z_array = data['z'].to_ndarray() + +plt.figure(figsize=(8,4)) +plt.plot(z_array[int(xCell/2),int(yCell/2),:]-z_array[int(xCell/2),int(yCell/2),int(zCell/2)],Ey_array[int(xCell/2),int(yCell/2),:]) +plt.ylim(-1.2e-3, 1.2e-3) +plt.ylim(-2.2e5, 2.2e5) +plt.xlim(-4.0e-6, 4.000001e-6) +plt.xticks(np.arange(-4.e-6,4.000001e-6, step=1e-6)) +plt.xlabel('z (m)') +plt.ylabel(field +' (V/m)') +plt.tight_layout() +plt.savefig('Ey_pec_analysis_mr.png') + +tolerance_rel = 0.05 + +print("min_Ey_error_rel : " + str(min_Ey_error_rel)) +print("max_Ey_error_rel : " + str(max_Ey_error_rel)) +print("tolerance_rel: " + str(tolerance_rel)) + +assert( max_Ey_error_rel < tolerance_rel ) +assert( min_Ey_error_rel < tolerance_rel ) + +test_name = fn[:-9] # Could also be os.path.split(os.getcwd())[1] + +if re.search( 'single_precision', fn ): + checksumAPI.evaluate_checksum(test_name, fn, rtol=1.e-3) +else: + checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Tests/PEC/inputs_field_PEC_3d b/Examples/Tests/PEC/inputs_field_PEC_3d new file mode 100644 index 000000000..ebe5dc423 --- /dev/null +++ b/Examples/Tests/PEC/inputs_field_PEC_3d @@ -0,0 +1,54 @@ +# Set-up to test the PEC Boundary condition for the fields +# Constructive interference between the incident and reflected wave result in a +# standing wave. + +# max step +max_step = 125 + +# number of grid points +amr.n_cell = 32 32 256 + +# Maximum allowable size of each subdomain +amr.max_grid_size = 1024 +amr.blocking_factor = 32 + +amr.max_level = 0 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.prob_lo = -8.e-6 -8.e-6 -4.e-6 +geometry.prob_hi = 8.e-6 8.e-6 4.e-6 + +# Boundary condition +boundary.field_lo = periodic periodic pec +boundary.field_hi = periodic periodic pec + +warpx.serialize_ics = 1 + +# Verbosity +warpx.verbose = 1 + +# Algorithms +algo.current_deposition = esirkepov +# CFL +warpx.cfl = 0.9 + + +my_constants.z1 = -2.e-6 +my_constants.z2 = 2.e-6 +my_constants.wavelength = 1.e-6 +warpx.E_ext_grid_init_style = parse_E_ext_grid_function +warpx.Ez_external_grid_function(x,y,z) = "0." +warpx.Ex_external_grid_function(x,y,z) = "0." +warpx.Ey_external_grid_function(x,y,z) = "((1.e5*sin(2*pi*(z)/wavelength)) * (z<z2) * (z>z1))" + +warpx.B_ext_grid_init_style = parse_B_ext_grid_function +warpx.Bx_external_grid_function(x,y,z)= "(((-1.e5*sin(2*pi*(z)/wavelength))/clight))*(z<z2) * (z>z1) " +warpx.By_external_grid_function(x,y,z)= "0." +warpx.Bz_external_grid_function(x,y,z) = "0." + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 125 +diag1.diag_type = Full +diag1.fields_to_plot = Ey Bx diff --git a/Examples/Tests/PEC/inputs_field_PEC_mr_3d b/Examples/Tests/PEC/inputs_field_PEC_mr_3d new file mode 100644 index 000000000..5d39f11b3 --- /dev/null +++ b/Examples/Tests/PEC/inputs_field_PEC_mr_3d @@ -0,0 +1,54 @@ +# Set-up to test the PEC Boundary condition for the fields + +# max step +max_step = 125 + +# number of grid points +amr.n_cell = 16 16 128 + +# Maximum allowable size of each subdomain +amr.max_grid_size = 1024 +amr.blocking_factor = 16 + +amr.max_level = 1 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.prob_lo = -8.e-6 -8.e-6 -4.e-6 +geometry.prob_hi = 8.e-6 8.e-6 4.e-6 +warpx.fine_tag_lo = -4.e-6 -4.e-6 -2.e-6 # physical domain +warpx.fine_tag_hi = 4.e-6 4.e-6 2.e-6 # physical domain + +# Boundary condition +boundary.field_lo = periodic periodic pec +boundary.field_hi = periodic periodic pec + +warpx.serialize_ics = 1 + +# Verbosity +warpx.verbose = 1 + +# Algorithms +algo.current_deposition = esirkepov +# CFL +warpx.cfl = 0.9 + + +my_constants.z1 = -2.e-6 +my_constants.z2 = 2.e-6 +my_constants.wavelength = 1.e-6 +warpx.E_ext_grid_init_style = parse_E_ext_grid_function +warpx.Ez_external_grid_function(x,y,z) = "0." +warpx.Ex_external_grid_function(x,y,z) = "0." +warpx.Ey_external_grid_function(x,y,z) = "((1.e5*sin(2*pi*(z)/wavelength)) * (z<z2) * (z>z1))" + +warpx.B_ext_grid_init_style = parse_B_ext_grid_function +warpx.Bx_external_grid_function(x,y,z)= "(((-1.e5*sin(2*pi*(z)/wavelength))/clight))*(z<z2) * (z>z1) " +warpx.By_external_grid_function(x,y,z)= "0." +warpx.Bz_external_grid_function(x,y,z) = "0." + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 125 +diag1.diag_type = Full +diag1.fields_to_plot = Ey Bx diff --git a/Examples/Tests/PEC/inputs_particle_PEC_3d b/Examples/Tests/PEC/inputs_particle_PEC_3d new file mode 100644 index 000000000..128c4abf5 --- /dev/null +++ b/Examples/Tests/PEC/inputs_particle_PEC_3d @@ -0,0 +1,57 @@ +# Set-up to test PEC boundary condition with particles. +# A proton and negatively charged particle with equal mass as a proton are initialized +# near the boundary. The proton is given a finite bulk velocity in y and the negatively +# charged particle is stationary, initially. Due to the PEC boundary setting zero Ey +# tangential to the PEC surface at x_min and x_max, the Ey field gathered by the particles +# is two orders of magnitude smaller than the periodic boundary. + +max_step = 20 +amr.n_cell = 128 64 64 + +amr.blocking_factor = 16 +amr.max_grid_size = 1024 +amr.max_level = 0 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.prob_lo = -32.e-6 -32.e-6 -32.e-6 # physical domain +geometry.prob_hi = 32.e-6 32.e-6 32.e-6 + +# boundary condition +boundary.field_lo = pec periodic periodic +boundary.field_hi = pec periodic periodic + + +# Algorithms +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.particle_pusher = vay +algo.maxwell_solver = yee +warpx.cfl = 0.9 +warpx.use_filter = 1 + +# Particle species +particles.species_names = electron proton + +electron.charge = -q_e +electron.mass = m_p +electron.injection_style = "singleparticle" +electron.single_particle_pos = 31.998e-6 0. 0. +electron.single_particle_vel = 0. 0. 0. +electron.single_particle_weight = 1.0 + +proton.charge = q_e +proton.mass = m_p # Very heavy ; should not move +proton.injection_style = "singleparticle" +proton.single_particle_pos = 31.998e-6 0.e-6 0. +proton.single_particle_vel = 0. -2. 0. +proton.single_particle_weight = 1.0 + +# Particle shape factor in each direction +algo.particle_shape = 3 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 20 +diag1.diag_type = Full +diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz diff --git a/Examples/Tests/galilean/inputs_rz b/Examples/Tests/galilean/inputs_rz index e20c37f49..2704e00fe 100644 --- a/Examples/Tests/galilean/inputs_rz +++ b/Examples/Tests/galilean/inputs_rz @@ -10,9 +10,10 @@ amr.max_level = 0 psatd.v_galilean = 0. 0. -0.99498743710662 geometry.coord_sys = 1 -geometry.is_periodic = 0 1 geometry.prob_lo = 0. -38.68 geometry.prob_hi = 38.68 38.68 +boundary.field_lo = none periodic +boundary.field_hi = none periodic ################################# ############ NUMERICS ########### diff --git a/Examples/Tests/spectral_staggered/inputs_hybrid_2d b/Examples/Tests/spectral_staggered/inputs_hybrid_2d index 08c92efd4..284a16521 100644 --- a/Examples/Tests/spectral_staggered/inputs_hybrid_2d +++ b/Examples/Tests/spectral_staggered/inputs_hybrid_2d @@ -5,7 +5,7 @@ amr.n_cell = 96 704 warpx.numprocs = 1 2 geometry.coord_sys = 0 -geometry.is_periodic = 0 0 +geometry.is_periodic = 1 0 geometry.prob_lo = -90.e-6 -70.e-6 geometry.prob_hi = 90.e-6 0.e-6 diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json index b6477c2a0..6ef636f75 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json @@ -2,30 +2,30 @@ "electrons": { "particle_cpu": 29440.0, "particle_id": 913434880.0, - "particle_momentum_x": 3.3220167752726847e-21, - "particle_momentum_y": 3.3049663196342287e-21, - "particle_momentum_z": 7.085667539965317e-21, - "particle_position_x": 0.5290008445390879, - "particle_position_y": 0.5888, + "particle_momentum_x": 3.3220167756236166e-21, + "particle_momentum_y": 3.30496631998347e-21, + "particle_momentum_z": 7.085667540713589e-21, + "particle_position_x": 0.5290008445390881, + "particle_position_y": 0.5888000000000001, "particle_theta": 92506.29869360026, "particle_weight": 81147583679.15044 }, "ions": { "particle_cpu": 29440.0, "particle_id": 2743896320.0, - "particle_momentum_x": 1.6915725245148433e-20, - "particle_momentum_y": 1.6994902822356307e-20, - "particle_momentum_z": 3.6847066227487176e-20, + "particle_momentum_x": 1.691572524693474e-20, + "particle_momentum_y": 1.6994902824150732e-20, + "particle_momentum_z": 3.684706623137782e-20, "particle_position_x": 0.5290000037284651, "particle_position_y": 0.5888, "particle_theta": 92089.14050622113, "particle_weight": 81147583679.15044 }, "lev=0": { - "By": 11.642563221840676, - "Ex": 1290320281471.9312, - "Ez": 1798383207735.9521, - "jx": 252864052657309.56, - "jz": 347280901333537.25 + "By": 11.64256322418941, + "Ex": 1290320281608.1353, + "Ez": 1798383207925.8325, + "jx": 252864052684005.03, + "jz": 347280901370212.9 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd.json index 20c859fc5..b2d9348d0 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd.json @@ -2,10 +2,10 @@ "electrons": { "particle_cpu": 29440.0, "particle_id": 913434880.0, - "particle_momentum_x": 1.2708268121853588e-20, - "particle_momentum_y": 1.2649567906287055e-20, - "particle_momentum_z": 2.7840518005688094e-20, - "particle_position_x": 0.5289998224216842, + "particle_momentum_x": 1.2708268123195493e-20, + "particle_momentum_y": 1.2649567907622783e-20, + "particle_momentum_z": 2.784051800862764e-20, + "particle_position_x": 0.5289998224216841, "particle_position_y": 0.5887999999999999, "particle_theta": 92506.29869360024, "particle_weight": 81147583679.15044 @@ -13,19 +13,19 @@ "ions": { "particle_cpu": 29440.0, "particle_id": 2743896320.0, - "particle_momentum_x": 1.0076728520983276e-21, - "particle_momentum_y": 1.0107440645323043e-21, - "particle_momentum_z": 1.9638565442671807e-21, - "particle_position_x": 0.5290000097838686, - "particle_position_y": 0.5888000000000001, + "particle_momentum_x": 1.0076728522150053e-21, + "particle_momentum_y": 1.0107440646493333e-21, + "particle_momentum_z": 1.9638565444745213e-21, + "particle_position_x": 0.5290000097838687, + "particle_position_y": 0.5888, "particle_theta": 92089.14050622113, "particle_weight": 81147583679.15044 }, "lev=0": { - "By": 3.408082904349099, - "Ex": 470238522779.2839, - "Ez": 660998748155.4016, - "jx": 892696391771030.2, - "jz": 1240343452920434.2 + "By": 3.4080829049007804, + "Ex": 470238522828.9618, + "Ez": 660998748225.2086, + "jx": 892696391865291.1, + "jz": 1240343453051394.5 } -} +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd_current_correction.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd_current_correction.json index 755ccdb6d..54b32c9a6 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd_current_correction.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz_psatd_current_correction.json @@ -2,32 +2,32 @@ "electrons": { "particle_cpu": 0.0, "particle_id": 1815005440.0, - "particle_momentum_x": 1.2741450985804339e-20, - "particle_momentum_y": 1.2747912182450875e-20, - "particle_momentum_z": 2.788652064780249e-20, - "particle_position_x": 0.5290000916698099, - "particle_position_y": 0.5888, - "particle_theta": 92606.42261144849, + "particle_momentum_x": 1.2741450987149708e-20, + "particle_momentum_y": 1.274791218379693e-20, + "particle_momentum_z": 2.788652065074695e-20, + "particle_position_x": 0.5290000916698098, + "particle_position_y": 0.5888000000000001, + "particle_theta": 92606.42261144848, "particle_weight": 81147583679.15044 }, "ions": { "particle_cpu": 0.0, "particle_id": 5475928320.0, - "particle_momentum_x": 9.468760984312313e-22, - "particle_momentum_y": 9.505742271609918e-22, - "particle_momentum_z": 1.9177141146735318e-21, + "particle_momentum_x": 9.468760985422032e-22, + "particle_momentum_y": 9.505742272721782e-22, + "particle_momentum_z": 1.9177141148761584e-21, "particle_position_x": 0.5290000096367489, "particle_position_y": 0.5888, "particle_theta": 92621.92887164818, "particle_weight": 81147583679.15044 }, "lev=0": { - "By": 3.3876510906847885, - "Ex": 462288551669.3938, - "Ez": 654508949022.0504, - "divE": 4.0082894251404314e+17, - "jx": 893199738249705.1, - "jz": 1241183345696335.2, - "rho": 3549014.737825297 + "By": 3.3876510912177427, + "Ex": 462288551718.21625, + "Ez": 654508949091.1531, + "divE": 4.008289425563525e+17, + "jx": 893199738344011.5, + "jz": 1241183345827358.0, + "rho": 3549014.7381998966 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json b/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json index 158a77310..218df2c8a 100644 --- a/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json +++ b/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json @@ -13,24 +13,24 @@ "electrons": { "particle_cpu": 4128.0, "particle_id": 10445216.0, - "particle_momentum_x": 3.221375738272786e-24, - "particle_momentum_y": 3.906701694384273e-22, - "particle_momentum_z": 1.371753970917068e-23, + "particle_momentum_x": 3.221375736217804e-24, + "particle_momentum_y": 3.9067016943869266e-22, + "particle_momentum_z": 1.371753970900497e-23, "particle_position_x": 0.041602500045359654, "particle_position_y": 0.04789125041898898, "particle_theta": 6754.424205218055, "particle_weight": 813672305.532158 }, "lev=0": { - "Bx": 101568.56451607398, + "Bx": 101568.56443034408, "By": 1224.3046117091226, - "Bz": 4142.575901037649, - "Ex": 204754846120.36102, - "Ey": 37715481947555.95, - "Ez": 548241772266.403, - "jx": 913428495166.8202, - "jy": 2.1552131562722282e+17, + "Bz": 4142.575962368743, + "Ex": 204754846120.36096, + "Ey": 37715481920582.85, + "Ez": 548241772266.4032, + "jx": 913428495169.505, + "jy": 2.1552131562722285e+17, "jz": 1772965958612206.2, - "rho": 38578930.47726793 + "rho": 38578930.47726787 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/LaserInjection_2d.json b/Regression/Checksum/benchmarks_json/LaserInjection_2d.json index 4f0b5adb4..7048a5a9b 100644 --- a/Regression/Checksum/benchmarks_json/LaserInjection_2d.json +++ b/Regression/Checksum/benchmarks_json/LaserInjection_2d.json @@ -1,13 +1,13 @@ { "lev=0": { - "Bx": 19661649.2843472324311733245849609375000000000000, - "By": 101030587.3703132718801498413085937500000000000000, - "Bz": 39719434.3724947720766067504882812500000000000000, - "Ex": 13838472033041280.0000000000000000000000000000000000000000, - "Ey": 13195219839199646.0000000000000000000000000000000000000000, - "Ez": 26772619800801276.0000000000000000000000000000000000000000, - "jx": 40346610039036864.0000000000000000000000000000000000000000, - "jy": 40346609403930968.0000000000000000000000000000000000000000, - "jz": 80693222329514336.0000000000000000000000000000000000000000 + "Bx": 19661649.262577523, + "By": 101030587.90280046, + "Bz": 39719434.779972926, + "Ex": 1.3838472028306398e+16, + "Ey": 1.319521974773849e+16, + "Ez": 2.6772619620846296e+16, + "jx": 4.034661003911076e+16, + "jy": 4.034660940393097e+16, + "jz": 8.06932223295597e+16 } -} +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/PEC_field.json b/Regression/Checksum/benchmarks_json/PEC_field.json new file mode 100644 index 000000000..c85a451c2 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/PEC_field.json @@ -0,0 +1,6 @@ +{ + "lev=0": { + "Bx": 7.325639113482142, + "Ey": 8374867179.431679 + } +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/PEC_field_mr.json b/Regression/Checksum/benchmarks_json/PEC_field_mr.json new file mode 100644 index 000000000..c6d84b011 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/PEC_field_mr.json @@ -0,0 +1,10 @@ +{ + "lev=0": { + "Bx": 1.3137691514773655, + "Ey": 1052208601.7204809 + }, + "lev=1": { + "Bx": 0.36174107941906025, + "Ey": 94466060.63059175 + } +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/PEC_particle.json b/Regression/Checksum/benchmarks_json/PEC_particle.json new file mode 100644 index 000000000..f6c1bf3a4 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/PEC_particle.json @@ -0,0 +1,35 @@ +{ + "electron": { + "particle_cpu": 0.0, + "particle_id": 1.0, + "particle_momentum_x": 5.957203757427428e-31, + "particle_momentum_y": 1.9967937404466168e-32, + "particle_momentum_z": 1.98439211397183e-47, + "particle_position_x": 3.199800000000332e-05, + "particle_position_y": 2.779673004428267e-19, + "particle_position_z": 1.0949416650810793e-34, + "particle_weight": 1.0 + }, + "lev=0": { + "Bx": 5.66137047937551e-05, + "By": 1.8816537559788902e-16, + "Bz": 0.00011100031731970804, + "Ex": 26731.847337613115, + "Ey": 29057.33399043876, + "Ez": 16060.200852115468, + "jx": 5.7823931096313506e-05, + "jy": 43090052.266495675, + "jz": 0.0 + }, + "proton": { + "particle_cpu": 0.0, + "particle_id": 2.0, + "particle_momentum_x": 5.652705594934373e-32, + "particle_momentum_y": 1.0028788756153835e-18, + "particle_momentum_z": 1.9377019463348945e-48, + "particle_position_x": 3.199799999999942e-05, + "particle_position_y": 6.572670690061942e-06, + "particle_position_z": 5.250319610855047e-36, + "particle_weight": 1.0 + } +} diff --git a/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json b/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json index 06c4b8367..fa936c3c7 100644 --- a/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json +++ b/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json @@ -1,44 +1,44 @@ { "beam": { "particle_cpu": 0.0, - "particle_id": 466643.0, - "particle_momentum_x": 5.935248102846229e-19, - "particle_momentum_y": 4.2079455139390086e-19, - "particle_momentum_z": 9.069157096209629e-18, - "particle_position_x": 0.0013709782662344476, - "particle_position_y": 0.14863704108231895, - "particle_weight": 2911663983235.947 + "particle_id": 500500.0, + "particle_momentum_x": 7.587450161877821e-19, + "particle_momentum_y": 4.515786832153623e-19, + "particle_momentum_z": 9.228375405679396e-18, + "particle_position_x": 0.005388318043043357, + "particle_position_y": 0.16843460833541193, + "particle_weight": 3120754537230.3823 }, "electrons": { "particle_cpu": 155246.0, "particle_id": 26701298988.0, - "particle_momentum_x": 9.169991045389105e-19, - "particle_momentum_y": 8.134032712971527e-19, - "particle_momentum_z": 5.716870930487574e-16, - "particle_position_x": 6.634385105820577, - "particle_position_y": 58.08979694901561, - "particle_weight": 3.212795210423339e+18 + "particle_momentum_x": 9.19759655761835e-19, + "particle_momentum_y": 8.111395741629531e-19, + "particle_momentum_z": 5.717012342804088e-16, + "particle_position_x": 6.634386766502452, + "particle_position_y": 58.08979821925373, + "particle_weight": 3.212795210423338e+18 }, "ions": { "particle_cpu": 155246.0, "particle_id": 26830093836.0, - "particle_momentum_x": 2.0028127909408492e-18, - "particle_momentum_y": 8.233828374969269e-19, - "particle_momentum_z": 1.0089821936116652e-12, - "particle_position_x": 6.6222971327574776, - "particle_position_y": 58.088776020644396, + "particle_momentum_x": 2.0052369737454595e-18, + "particle_momentum_y": 8.210971641417303e-19, + "particle_momentum_z": 1.0089821790887407e-12, + "particle_position_x": 6.622297130412207, + "particle_position_y": 58.08877602017003, "particle_weight": 3.212794492131793e+18 }, "lev=0": { - "Bx": 1305305.431740205, - "By": 4004955.3295348613, - "Bz": 203986.21157566475, - "Ex": 1202527548976529.5, - "Ey": 408729463466964.56, - "Ez": 217606643896681.66, - "jx": 3.5156420344831892e+16, - "jy": 2.330006936826398e+16, - "jz": 4.1305946927278195e+17, - "rho": 1350581221.909034 + "Bx": 1305297.4828241316, + "By": 4009300.510439389, + "Bz": 203715.33013685, + "Ex": 1203889010242205.5, + "Ey": 408608837447717.06, + "Ez": 217844734969419.84, + "jx": 3.5197389872605124e+16, + "jy": 2.314621849516091e+16, + "jz": 4.137802990284829e+17, + "rho": 1352926282.578562 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json b/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json index 4dc35e93f..36159b1cb 100644 --- a/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json +++ b/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json @@ -1,44 +1,44 @@ { "beam": { "particle_cpu": 0.0, - "particle_id": 447573.0, - "particle_momentum_x": 5.384525633353412e-19, - "particle_momentum_y": 4.011357939321231e-19, - "particle_momentum_z": 8.08848605215067e-18, - "particle_position_x": 0.0014094341950741979, - "particle_position_y": 0.1412193837490303, - "particle_weight": 2780592292672.2705 + "particle_id": 500500.0, + "particle_momentum_x": 7.781818973078852e-19, + "particle_momentum_y": 4.515677165595182e-19, + "particle_momentum_z": 8.375527032802002e-18, + "particle_position_x": 0.006743036325434187, + "particle_position_y": 0.1716091883553772, + "particle_weight": 3120754537230.3823 }, "electrons": { "particle_cpu": 155383.0, "particle_id": 26721507984.0, - "particle_momentum_x": 9.809438398217426e-19, - "particle_momentum_y": 5.687153095622055e-19, - "particle_momentum_z": 5.697506964559476e-16, - "particle_position_x": 6.639908741842812, - "particle_position_y": 58.210506112500646, + "particle_momentum_x": 7.343166622433748e-19, + "particle_momentum_y": 5.669478712783679e-19, + "particle_momentum_z": 5.697759351248577e-16, + "particle_position_x": 6.639863939125528, + "particle_position_y": 58.2105156240904, "particle_weight": 3.215903744897337e+18 }, "ions": { "particle_cpu": 155428.0, "particle_id": 26814790002.0, - "particle_momentum_x": 1.906616117322184e-18, - "particle_momentum_y": 5.75419216528982e-19, - "particle_momentum_z": 1.0101677927464587e-12, - "particle_position_x": 6.630049949959092, - "particle_position_y": 58.24914762101138, + "particle_momentum_x": 1.6709338576817887e-18, + "particle_momentum_y": 5.73611349305925e-19, + "particle_momentum_z": 1.010167756546284e-12, + "particle_position_x": 6.63004997278591, + "particle_position_y": 58.249147619440286, "particle_weight": 3.21662600232034e+18 }, "lev=0": { - "Bx": 1225109.2628475595, - "By": 3748536.8460838096, - "Bz": 195216.2190878017, - "Ex": 1092766647326914.6, - "Ey": 384575347595082.3, - "Ez": 343130879909148.75, - "jx": 3.51564554699967e+16, - "jy": 2.2405130712647624e+16, - "jz": 3.0976785868961254e+17, - "rho": 1054402097.6298534 + "Bx": 1225047.3704144938, + "By": 3461610.5122063416, + "Bz": 194622.52078763963, + "Ex": 1039200462527705.4, + "Ey": 384410222541170.8, + "Ez": 200388928565719.44, + "jx": 3.100345182047955e+16, + "jy": 2.2227570738047532e+16, + "jz": 2.9254388838934234e+17, + "rho": 995460563.3392936 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json b/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json index e40fa7112..190f062f8 100644 --- a/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json +++ b/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json @@ -2,10 +2,10 @@ "electrons": { "particle_cpu": 0.0, "particle_id": 540033024.0, - "particle_momentum_x": 6.999660495409049e-22, - "particle_momentum_y": 2.7286257042710144e-22, - "particle_momentum_z": 8.903628651318955e-17, - "particle_position_x": 633733.7284462163, + "particle_momentum_x": 6.999660495445657e-22, + "particle_momentum_y": 2.7286257042962753e-22, + "particle_momentum_z": 8.903628651318945e-17, + "particle_position_x": 633733.7284462162, "particle_position_y": 7862151.373823494, "particle_theta": 51362.08760999037, "particle_weight": 1.0261080645329302e+20 @@ -13,8 +13,8 @@ "ions": { "particle_cpu": 0.0, "particle_id": 1630552064.0, - "particle_momentum_x": 1.3125868827082553e-18, - "particle_momentum_y": 2.7234814546935865e-22, + "particle_momentum_x": 1.312586882708253e-18, + "particle_momentum_y": 2.72348145473723e-22, "particle_momentum_z": 1.634880473424221e-13, "particle_position_x": 633733.7424979067, "particle_position_y": 7862151.99748233, @@ -22,13 +22,13 @@ "particle_weight": 1.0261080645329302e+20 }, "lev=0": { - "By": 0.0024307114063147504, - "Ex": 732693.5768885817, - "Ey": 130003.11944446711, - "Ez": 36772.52480157818, - "divE": 1348744.2569353955, - "jx": 163.0127498331044, - "jz": 50246.47997764914, - "rho": 1.1999059634758576e-05 + "By": 0.0024307114063306778, + "Ex": 732693.5769027277, + "Ey": 130003.11944446625, + "Ez": 36772.524796743244, + "divE": 1348744.2569426426, + "jx": 163.01274983651132, + "jz": 50246.47997764428, + "rho": 1.1999059634818442e-05 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json b/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json index 580e7f19c..ab57d3f93 100644 --- a/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json +++ b/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json @@ -2,9 +2,9 @@ "electrons": { "particle_cpu": 0.0, "particle_id": 540033024.0, - "particle_momentum_x": 7.043225585491615e-22, - "particle_momentum_y": 2.8095488450457388e-22, - "particle_momentum_z": 8.903838826207317e-17, + "particle_momentum_x": 7.043225585505337e-22, + "particle_momentum_y": 2.809548845143752e-22, + "particle_momentum_z": 8.903838826207275e-17, "particle_position_x": 633734.0562112605, "particle_position_y": 7862151.998116621, "particle_theta": 51362.087830016804, @@ -13,22 +13,22 @@ "ions": { "particle_cpu": 0.0, "particle_id": 1630552064.0, - "particle_momentum_x": 1.3125871812988232e-18, - "particle_momentum_y": 2.8053964233138836e-22, - "particle_momentum_z": 1.634880452407025e-13, - "particle_position_x": 633733.7423183098, + "particle_momentum_x": 1.3125871812988242e-18, + "particle_momentum_y": 2.8053964233878157e-22, + "particle_momentum_z": 1.6348804524070253e-13, + "particle_position_x": 633733.7423183097, "particle_position_y": 7862151.997142417, "particle_theta": 51470.93289806633, "particle_weight": 1.0261080645329302e+20 }, "lev=0": { - "By": 0.0018845897904908619, - "Ex": 568097.4835035336, - "Ey": 130003.2620179485, - "Ez": 40677.99022297442, - "divE": 1279972.8318518812, - "jx": 166.90827028769894, - "jz": 50430.778236016544, - "rho": 1.1333119848498335e-05 + "By": 0.0018845897905120875, + "Ex": 568097.483503956, + "Ey": 130003.26201794606, + "Ez": 40677.99022514904, + "divE": 1279972.8318505273, + "jx": 166.90827030470052, + "jz": 50430.77823601886, + "rho": 1.1333119848485656e-05 } }
\ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 22e012993..7be3e2210 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2050,3 +2050,54 @@ compileTest = 0 doVis = 0 compareParticles = 0 analysisRoutine = Examples/Tests/ElectrostaticDirichletBC/analysis.py + +[PEC_field] +buildDir = . +inputFile = Examples/Tests/PEC/inputs_field_PEC_3d +runtime_params = +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Tests/PEC/analysis_pec.py + +[PEC_field_mr] +buildDir = . +inputFile = Examples/Tests/PEC/inputs_field_PEC_mr_3d +runtime_params = +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Tests/PEC/analysis_pec_mr.py + +[PEC_particle] +buildDir = . +inputFile = Examples/Tests/PEC/inputs_particle_PEC_3d +runtime_params = +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electron proton +analysisRoutine = Examples/analysis_default_regression.py +tolerance = 1.e-14 + diff --git a/Source/BoundaryConditions/CMakeLists.txt b/Source/BoundaryConditions/CMakeLists.txt index 265f5c210..3fa90d1f1 100644 --- a/Source/BoundaryConditions/CMakeLists.txt +++ b/Source/BoundaryConditions/CMakeLists.txt @@ -2,4 +2,6 @@ target_sources(WarpX PRIVATE PML.cpp WarpXEvolvePML.cpp + WarpXFieldBoundaries.cpp + WarpX_PEC.cpp ) diff --git a/Source/BoundaryConditions/Make.package b/Source/BoundaryConditions/Make.package index aab9158a4..10fbdf080 100644 --- a/Source/BoundaryConditions/Make.package +++ b/Source/BoundaryConditions/Make.package @@ -1,3 +1,5 @@ CEXE_sources += PML.cpp WarpXEvolvePML.cpp +CEXE_sources += WarpXFieldBoundaries.cpp WarpX_PEC.cpp + VPATH_LOCATIONS += $(WARPX_HOME)/Source/BoundaryConditions diff --git a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp new file mode 100644 index 000000000..7e99dbc8e --- /dev/null +++ b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp @@ -0,0 +1,29 @@ +#include "WarpX.H" +#include "WarpX_PEC.H" +#include <AMReX.H> +#include <AMReX_Vector.H> +#include <AMReX_MultiFab.H> +using namespace amrex::literals; + +void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type) +{ + if (PEC::isAnyBoundaryPEC()) { + if (patch_type == PatchType::fine) { + PEC::ApplyPECtoEfield( Efield_fp[lev], lev, patch_type); + } else { + PEC::ApplyPECtoEfield( Efield_cp[lev], lev, patch_type); + } + } +} + +void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type) +{ + if (PEC::isAnyBoundaryPEC()) { + if (patch_type == PatchType::fine) { + PEC::ApplyPECtoBfield( Bfield_fp[lev], lev, patch_type); + } else { + PEC::ApplyPECtoBfield( Bfield_cp[lev], lev, patch_type); + } + } +} + diff --git a/Source/BoundaryConditions/WarpX_PEC.H b/Source/BoundaryConditions/WarpX_PEC.H new file mode 100644 index 000000000..f196d8bf8 --- /dev/null +++ b/Source/BoundaryConditions/WarpX_PEC.H @@ -0,0 +1,286 @@ +#ifndef WARPX_PEC_KERNELS_H_ +#define WARPX_PEC_KERNELS_H_ + +#include "WarpX.H" +#include "Utils/WarpXAlgorithmSelection.H" + + +namespace PEC { +using namespace amrex; + /** + * \brief Determines if the field boundary condition stored in fboundary is PEC + * in direction, dir, is PEC. + * + * \param[in] fboundary Value containing boundary type + * \param[in] dir direction + * + * \param[out] 1 if the boundary type is PEC else 0 + */ + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + bool is_boundary_PEC (amrex::GpuArray<int, 3> const& fboundary, int dir) { + return ( fboundary[dir] == FieldBoundaryType::PEC); + } + + /** + * \brief Sets the electric field value tangential to the PEC boundary to zero. The + * tangential Efield components in the guard cells outside the + * domain boundary are set equal and opposite to the field in the valid cells + * at their mirrored locations. The normal Efield components in the guard cells + * are set equal to the field in the valid cells at their mirrored locations. + * The number or depth of guard cells updated is equal to the shape factor of + * particles in each dimension. + * For corner cells with mixed boundaries, the mirror location could be outside + * valid region, while still ensuring PEC condition is maintained across the + * PEC boundary, and the necessary sign change is accounted for depending on + * if the component, icomp, is tangential or normal to the PEC boundary. + * + * For 3D : + * x component is tangential to the y-boundary and z-boundary + * y component is tangential to the x-boundary and z-boundary + * z component is tangential to the x-boundary and y-boundary + * x component is normal to the x-boundary + * y component is normal to the y-boundary + * z component is normal to the z-boundary + * where, x-boundary is the yz-plane at x=xmin and x=xmax + * y-boundary is the xz-plane at y=ymin and y=ymax + * z-boundary is the xy-plane at z=zmin and z=zmax + * + * For 2D : WarpX uses X-Z as the two dimensions + * x component is tangential to the z-boundary + * y component is tangential to the x-boundary and z-boundary + * z component is tangential to the x-boundary + * x component is normal to the x-boundary + * y component is not normal to any boundary (Only xz dimensions in 2D) + * z component is normal to the z-boundary + * where, x-boundary is along the line z at x=xmin and x=xmax + * z-boundary is along the line x at z=zmin and z=zmax + * + * For RZ : WarpX uses R-Z as the two dimensions + * r component is tangential to the z-boundary + * theta_component is tangential to the r-boundary and z-boundary + * z component is tangential to the r-boundary + * r component is normal to the r-boundary + * theta_component is not normal to any boundary (on RZ dimensions are modeled) + * z component is normal to the z-boundary + * where, r-boundary is along the line z at r=rmin and r=rmax + * z-boundary is along the line r at z=zmin and z=zmax + * + * + * \param[in] icomp component of the Efield being updated + * (0=x, 1=y, 2=z in Cartesian) + * (0=r, 1=theta, 2=z in RZ) + * \param[in] dom_lo index value of the lower domain boundary (cell-centered) + * \param[in] dom_hi index value of the higher domain boundary (cell-centered) + * \param[in] ijk_vec indices along the x(i), y(j), z(k) of Efield Array4 + * \param[in] n index of the MultiFab component being updated + * \param[in] Efield field data to be updated if (ijk) is at the boundary or a guard cell + * \param[in] is_nodal staggering of the field data being updated. + * \param[in] fbndry_lo Field boundary type at the lower boundaries + * \param[in] fbndry_hi Field boundary type at the upper boundaries + */ + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void SetEfieldOnPEC (const int icomp, const amrex::IntVect & dom_lo, + const amrex::IntVect &dom_hi, + const amrex::IntVect &ijk_vec, const int n, + amrex::Array4<amrex::Real> const& Efield, + const amrex::IntVect& is_nodal, + amrex::GpuArray<int, 3> const& fbndry_lo, + amrex::GpuArray<int, 3> const& fbndry_hi ) + { + // Tangential Efield componentes in guard cells set equal and opposite to cells + // in the mirror locations across the PEC boundary, whereas normal E-field + // components are set equal to values in the mirror locations across the PEC + // boundary. Here we just initialize it. + amrex::IntVect ijk_mirror = ijk_vec; + bool OnPECBoundary = false; + bool GuardCell = false; + amrex::Real sign = 1._rt; + // Loop over all the dimensions + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + // Loop over sides, iside = 0 (lo), iside = 1 (hi) + // Loop over sides, iside = 0 (lo), iside = 1 (hi) + for (int iside = 0; iside < 2; ++iside) { + const bool isPECBoundary = ( (iside == 0) + ? is_boundary_PEC(fbndry_lo, idim) + : is_boundary_PEC(fbndry_hi, idim) ); +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + // For 2D : for icomp==1, (Ey in XZ, Etheta in RZ), + // icomp=1 is tangential to both x and z boundaries + // The logic below ensures that the flags are set right for 2D + const bool is_tangent_to_PEC = ( (icomp == AMREX_SPACEDIM*idim) + ? false : true ); +#else + const bool is_tangent_to_PEC = ( ( icomp == idim) ? false : true ); +#endif + if (isPECBoundary == true) { + // guard cell outside the domain by "ig" number cells, in direction, idim + const int ig = ( (iside == 0) + ? (dom_lo[idim] - ijk_vec[idim]) + : (ijk_vec[idim] - (dom_hi[idim] + is_nodal[idim]) ) ); + if (ig == 0) { + if (is_tangent_to_PEC == true and is_nodal[idim] == 1) { + OnPECBoundary = true; + } + } else if (ig > 0) { + // Find mirror location across PEC boundary + ijk_mirror[idim] = ( ( iside == 0) + ? (dom_lo[idim] + ig) + : (dom_hi[idim] + is_nodal[idim] - ig)); + GuardCell = true; + // tangential components are inverted across PEC boundary + if (is_tangent_to_PEC == true) sign *= -1._rt; + } + } // is PEC boundary + } // loop over iside + } // loop over dimensions + if (OnPECBoundary == true) { + // if ijk_vec is on a PEC boundary in any direction, set Etangential to 0. + Efield(ijk_vec,n) = 0._rt; + } else if (GuardCell == true) { + Efield(ijk_vec,n) = sign * Efield(ijk_mirror,n); + } + } + + /** + * \brief Sets the magnetic field value normal to the PEC boundary to zero. The + * tangential (and normal) field value of the guard cells outside the + * domain boundary are set equal (and opposite) to the respective field components + * in the valid cells at their mirrored locations. + * The number or depth of guard cells updated is equal to the shape factor of + * particles in each dimension. + * + * For 3D : + * x component is tangential to the y-boundary and z-boundary + * y component is tangential to the x-boundary and z-boundary + * z component is tangential to the x-boundary and y-boundary + * x component is normal to the x-boundary + * y component is normal to the y-boundary + * z component is normal to the z-boundary + * where, x-boundary is the yz-plane at x=xmin and x=xmax + * y-boundary is the xz-plane at y=ymin and y=ymax + * z-boundary is the xy-plane at z=zmin and z=zmax + * + * For 2D : WarpX uses X-Z as the two dimensions + * x component is tangential to the z-boundary + * y component is tangential to the x-boundary and z-boundary + * z component is tangential to the x-boundary + * x component is normal to the x-boundary + * y component is not normal to any boundary (Only xz dimensions in 2D) + * z component is normal to the z-boundary + * where, x-boundary is along the line z at x=xmin and x=xmax + * z-boundary is along the line x at z=zmin and z=zmax + * + * For RZ : WarpX uses R-Z as the two dimensions + * r component is tangential to the z-boundary + * theta_component is tangential to the r-boundary and z-boundary + * z component is tangential to the r-boundary + * r component is normal to the r-boundary + * theta_component is not normal to any boundary (on RZ dimensions are modeled) + * z component is normal to the z-boundary + * where, r-boundary is along the line z at r=rmin and r=rmax + * z-boundary is along the line r at z=zmin and z=zmax + * + * + * \param[in] icomp component of the Bfield being updated + * (0=x, 1=y, 2=z in Cartesian) + * (0=r, 1=theta, 2=z in RZ) + * \param[in] dom_lo index value of the lower domain boundary (cell-centered) + * \param[in] dom_hi index value of the higher domain boundary (cell-centered) + * \param[in] ijk_vec indices along the x(i), y(j), z(k) of Efield Array4 + * \param[in] n index of the MultiFab component being updated + * \param[in] Bfield field data to be updated if (ijk) is at the boundary + or a guard cell + * \param[in] is_nodal staggering of the field data being updated. + * \param[in] fbndry_lo Field boundary type at the lower boundaries + * \param[in] fbndry_hi Field boundary type at the upper boundaries + */ + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void SetBfieldOnPEC (const int icomp, const amrex::IntVect & dom_lo, + const amrex::IntVect & dom_hi, + const amrex::IntVect & ijk_vec, const int n, + amrex::Array4<amrex::Real> const& Bfield, + const amrex::IntVect & is_nodal, + amrex::GpuArray<int, 3> const& fbndry_lo, + amrex::GpuArray<int, 3> const& fbndry_hi ) + { + amrex::IntVect ijk_mirror = ijk_vec; + bool OnPECBoundary = false; + bool GuardCell = false; + amrex::Real sign = 1._rt; + // Loop over all dimensions + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + // Loop over sides, iside = 0 (lo), iside = 1 (hi) + for (int iside = 0; iside < 2; ++iside) { + const int isPECBoundary = ( (iside == 0 ) + ? is_boundary_PEC(fbndry_lo, idim) + : is_boundary_PEC(fbndry_hi, idim) ); + if (isPECBoundary == true) { +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + // For 2D : for icomp==1, (By in XZ, Btheta in RZ), + // icomp=1 is not normal to x or z boundary + // The logic below ensures that the flags are set right for 2D + const bool is_normal_to_PEC = ( (icomp == AMREX_SPACEDIM*idim) + ? true : false ); +#else + const bool is_normal_to_PEC = ( ( icomp == idim) ? true : false ); +#endif + // guard cell outside the domain by "ig" number cells, in direction, idim + const int ig = ( (iside == 0) + ? (dom_lo[idim] - ijk_vec[idim]) + : (ijk_vec[idim] - (dom_hi[idim] + is_nodal[idim]) ) ); + if (ig == 0) { + // Only normal component is set to 0 + if (is_normal_to_PEC == true and is_nodal[idim]==1) { + OnPECBoundary = true; + } + } else if ( ig > 0) { + // Mirror location inside the domain by "ig" number of cells + // across PEC boundary in direction, idim, and side, iside + ijk_mirror[idim] = ( (iside == 0) + ? (dom_lo[idim] + ig) + : (dom_hi[idim] + is_nodal[idim] - ig)); + GuardCell = true; + // Sign of the normal component in guard cell is inverted + if (is_normal_to_PEC == true) sign *= -1._rt; + } + } // if PEC Boundary + } // loop over sides + } // loop of dimensions + + if (OnPECBoundary == true) { + // if ijk_vec is on a PEC boundary in any direction, set Bnormal to 0. + Bfield(ijk_vec,n) = 0._rt; + } else if (GuardCell == true) { + // Bnormal and Btangential is set opposite and equal to the value + // in the mirror location, respectively. + Bfield(ijk_vec,n) = sign * Bfield(ijk_mirror,n); + } + } + + /** Returns 1 if any domain boundary is set to PEC, else returns 0.*/ + bool isAnyBoundaryPEC(); + /** + * \brief Sets the tangential electric field at the PEC boundary to zero. + * The guard cell values are set equal and opposite to the valid cell + * field value at the respective mirror locations. + * + * \param[in,out] Efield Boundary values of tangential Efield are set to zero. + * \param[in] lev level of the Multifab + * \param[in] patch_type coarse or fine + */ + void ApplyPECtoEfield ( std::array<std::unique_ptr<amrex::MultiFab>, 3>& Efield, + const int lev, PatchType patch_type); + /** + * \brief Sets the normal component of the magnetic field at the PEC boundary to zero. + * The guard cell values are set equal and opposite to the valid cell + * field value at the respective mirror locations. + * + * \param[in,out] Bfield Boundary values of normal Bfield are set to zero. + * \param[in] lev level of the Multifab + * \param[in] patch_type coarse or fine + */ + void ApplyPECtoBfield ( std::array<std::unique_ptr<amrex::MultiFab>, 3>& Bfield, + const int lev, PatchType patch_type); +} + +#endif // WarpX_PEC_KERNELS_H_ diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp new file mode 100644 index 000000000..29eb8d60a --- /dev/null +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -0,0 +1,179 @@ +#include "BoundaryConditions/WarpX_PEC.H" +#include "WarpX.H" +#include <AMReX.H> +#include <AMReX_Vector.H> +#include <AMReX_MultiFab.H> +using namespace amrex::literals; + + +bool +PEC::isAnyBoundaryPEC() { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC) return true; + if ( WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC) return true; + } + return false; +} + +void +PEC::ApplyPECtoEfield (std::array<std::unique_ptr<amrex::MultiFab>, 3>& Efield, const int lev, + PatchType patch_type) +{ + auto& warpx = WarpX::GetInstance(); + amrex::Box domain_box = warpx.Geom(lev).Domain(); + if (patch_type == PatchType::coarse) { + amrex::IntVect ref_ratio = ( (lev > 0) ? WarpX::RefRatio(lev-1) : amrex::IntVect(1) ); + domain_box.coarsen(ref_ratio); + } + amrex::IntVect domain_lo = domain_box.smallEnd(); + amrex::IntVect domain_hi = domain_box.bigEnd(); + amrex::IntVect shape_factor(AMREX_D_DECL(WarpX::nox, WarpX::noy, WarpX::noz)); +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + shape_factor[1] = WarpX::noz; +#endif + shape_factor.max( amrex::IntVect::TheUnitVector() ); + amrex::GpuArray<int, 3> fbndry_lo; + amrex::GpuArray<int, 3> fbndry_hi; + for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { + fbndry_lo[idim] = WarpX::field_boundary_lo[idim]; + fbndry_hi[idim] = WarpX::field_boundary_hi[idim]; + } + amrex::IntVect Ex_nodal = Efield[0]->ixType().toIntVect(); + amrex::IntVect Ey_nodal = Efield[1]->ixType().toIntVect(); + amrex::IntVect Ez_nodal = Efield[2]->ixType().toIntVect(); + int nComp_x = Efield[0]->nComp(); + int nComp_y = Efield[1]->nComp(); + int nComp_z = Efield[2]->nComp(); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(*Efield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + // Extract field data + amrex::Array4<amrex::Real> const& Ex = Efield[0]->array(mfi); + amrex::Array4<amrex::Real> const& Ey = Efield[1]->array(mfi); + amrex::Array4<amrex::Real> const& Ez = Efield[2]->array(mfi); + + // Extract tileboxes for which to loop + amrex::Box const& tex = mfi.tilebox(Efield[0]->ixType().toIntVect(), shape_factor); + amrex::Box const& tey = mfi.tilebox(Efield[1]->ixType().toIntVect(), shape_factor); + amrex::Box const& tez = mfi.tilebox(Efield[2]->ixType().toIntVect(), shape_factor); + + // loop over cells and update fields + amrex::ParallelFor( + tex, nComp_x, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + amrex::ignore_unused(k); +#endif + amrex::IntVect iv(AMREX_D_DECL(i,j,k)); + const int icomp = 0; + PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, + Ex, Ex_nodal, fbndry_lo, fbndry_hi); + }, + tey, nComp_y, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + amrex::ignore_unused(k); +#endif + amrex::IntVect iv(AMREX_D_DECL(i,j,k)); + const int icomp = 1; + PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, + Ey, Ey_nodal, fbndry_lo, fbndry_hi); + }, + tez, nComp_z, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + amrex::ignore_unused(k); +#endif + amrex::IntVect iv(AMREX_D_DECL(i,j,k)); + const int icomp = 2; + PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, + Ez, Ez_nodal, fbndry_lo, fbndry_hi); + } + ); + + } +} + + +void +PEC::ApplyPECtoBfield (std::array<std::unique_ptr<amrex::MultiFab>, 3>& Bfield, const int lev, + PatchType patch_type) +{ + auto& warpx = WarpX::GetInstance(); + amrex::Box domain_box = warpx.Geom(lev).Domain(); + if (patch_type == PatchType::coarse) { + amrex::IntVect ref_ratio = ( (lev > 0) ? WarpX::RefRatio(lev-1) : amrex::IntVect(1) ); + domain_box.coarsen(ref_ratio); + } + amrex::IntVect domain_lo = domain_box.smallEnd(); + amrex::IntVect domain_hi = domain_box.bigEnd(); + amrex::IntVect shape_factor(AMREX_D_DECL(WarpX::nox, WarpX::noy, WarpX::noz)); +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + shape_factor[1] = WarpX::noz; +#endif + amrex::GpuArray<int, 3> fbndry_lo; + amrex::GpuArray<int, 3> fbndry_hi; + for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { + fbndry_lo[idim] = WarpX::field_boundary_lo[idim]; + fbndry_hi[idim] = WarpX::field_boundary_hi[idim]; + } + amrex::IntVect Bx_nodal = Bfield[0]->ixType().toIntVect(); + amrex::IntVect By_nodal = Bfield[1]->ixType().toIntVect(); + amrex::IntVect Bz_nodal = Bfield[2]->ixType().toIntVect(); + const int nComp_x = Bfield[0]->nComp(); + const int nComp_y = Bfield[1]->nComp(); + const int nComp_z = Bfield[2]->nComp(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + // Extract field data + amrex::Array4<amrex::Real> const& Bx = Bfield[0]->array(mfi); + amrex::Array4<amrex::Real> const& By = Bfield[1]->array(mfi); + amrex::Array4<amrex::Real> const& Bz = Bfield[2]->array(mfi); + + // Extract tileboxes for which to loop + amrex::Box const& tbx = mfi.tilebox(Bfield[0]->ixType().toIntVect(), shape_factor); + amrex::Box const& tby = mfi.tilebox(Bfield[1]->ixType().toIntVect(), shape_factor); + amrex::Box const& tbz = mfi.tilebox(Bfield[2]->ixType().toIntVect(), shape_factor); + + // loop over cells and update fields + amrex::ParallelFor( + tbx, nComp_x, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + amrex::ignore_unused(k); +#endif + amrex::IntVect iv(AMREX_D_DECL(i,j,k)); + const int icomp = 0; + PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, + Bx, Bx_nodal, fbndry_lo, fbndry_hi); + }, + tby, nComp_y, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + amrex::ignore_unused(k); +#endif + amrex::IntVect iv(AMREX_D_DECL(i,j,k)); + const int icomp = 1; + PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, + By, By_nodal, fbndry_lo, fbndry_hi); + }, + tbz, nComp_z, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + amrex::ignore_unused(k); +#endif + amrex::IntVect iv(AMREX_D_DECL(i,j,k)); + const int icomp = 2; + PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, + Bz, Bz_nodal, fbndry_lo, fbndry_hi); + } + ); + + } + +} diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 011c88b4f..2329854f2 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -140,6 +140,10 @@ WarpX::PushPSATD (amrex::Real a_dt) if (do_pml && pml[lev]->ok()) { pml[lev]->PushPSATD(lev); } + ApplyEfieldBoundary(lev,PatchType::fine); + if (lev > 0) ApplyEfieldBoundary(lev,PatchType::coarse); + ApplyBfieldBoundary(lev,PatchType::fine); + if (lev > 0) ApplyBfieldBoundary(lev,PatchType::coarse); } #endif } @@ -221,6 +225,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) } } + ApplyBfieldBoundary(lev, patch_type); } void @@ -280,6 +285,9 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) a_dt, pml_has_particles ); } } + + ApplyEfieldBoundary(lev, patch_type); + } @@ -420,6 +428,8 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { a_dt, pml_has_particles ); } } + + ApplyEfieldBoundary(lev, patch_type); } void diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 107b0f224..07d7c89a2 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -100,8 +100,10 @@ struct FieldBoundaryType { Periodic = 1, PEC = 2, //!< perfect electric conductor (PEC) with E_tangential=0 PMC = 3, //!< perfect magnetic conductor (PMC) with B_tangential=0 - Damped = 4 // Fields in the guard cells are damped for PSATD + Damped = 4, // Fields in the guard cells are damped for PSATD //in the moving window direction + None = 5 // The fields values at the boundary are not updated. This is + // useful for RZ simulations, at r=0. }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index a0b90efec..cdf896eb2 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -81,6 +81,7 @@ const std::map<std::string, int> FieldBCType_algo_to_int = { {"pec", FieldBoundaryType::PEC}, {"pmc", FieldBoundaryType::PMC}, {"damped", FieldBoundaryType::Damped}, + {"none", FieldBoundaryType::None}, {"default", FieldBoundaryType::PML} }; diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index d85cd17a4..e9cb6e34f 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -446,6 +446,7 @@ void ReadBCParams () ParmParse pp_geometry("geometry"); ParmParse pp_warpx("warpx"); ParmParse pp_algo("algo"); + int maxwell_solver_id = GetAlgorithmInteger(pp_algo, "maxwell_solver"); if (pp_geometry.queryarr("is_periodic", geom_periodicity)) { // set default field and particle boundary appropriately for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { @@ -458,15 +459,24 @@ void ReadBCParams () } else { // if non-periodic and do_pml=0, then set default boundary to PEC int pml_input = 1; + int silverMueller_input = 0; pp_warpx.query("do_pml", pml_input); - if (pml_input == 0) { - WarpX::field_boundary_lo[idim] = FieldBoundaryType::PEC; - WarpX::field_boundary_hi[idim] = FieldBoundaryType::PEC; + pp_warpx.query("do_silver_mueller", silverMueller_input); + if (pml_input == 0 and silverMueller_input == 0) { + if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) { + WarpX::field_boundary_lo[idim] = FieldBoundaryType::None; + WarpX::field_boundary_hi[idim] = FieldBoundaryType::None; + } else { + WarpX::field_boundary_lo[idim] = FieldBoundaryType::PEC; + WarpX::field_boundary_hi[idim] = FieldBoundaryType::PEC; + } +#ifdef WARPX_DIM_RZ + if (idim == 0) WarpX::field_boundary_lo[idim] = FieldBoundaryType::None; +#endif } } } - // Temporarily setting default boundary to Damped until PEC Boundary Type is enabled - int maxwell_solver_id = GetAlgorithmInteger(pp_algo, "maxwell_solver"); + // Temporarily setting default boundary to Damped until new boundary interface is introduced if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) { ParmParse pp_psatd("psatd"); int do_moving_window = 0; @@ -528,9 +538,19 @@ void ReadBCParams () WarpX::particle_boundary_lo[idim] = ParticleBoundaryType::Periodic; WarpX::particle_boundary_hi[idim] = ParticleBoundaryType::Periodic; } - + } + if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) { + if (WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC || + WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC) { + amrex::Abort(" PEC boundary not implemented for PSATD, yet!"); + } } } +#ifdef WARPX_DIM_RZ + // Ensure code aborts if PEC is specified at r=0 for RZ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::field_boundary_lo[0] == FieldBoundaryType::None, + "Error : Field boundary at r=0 must be ``none``. \n"); +#endif pp_geometry.addarr("is_periodic", geom_periodicity); } diff --git a/Source/WarpX.H b/Source/WarpX.H index 7b86251bd..26f6045d0 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -406,6 +406,9 @@ public: int lev); #endif + void ApplyEfieldBoundary (const int lev, PatchType patch_type); + void ApplyBfieldBoundary (const int lev, PatchType patch_type); + void DampPML (); void DampPML (int lev); void DampPML (int lev, PatchType patch_type); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 9c9091cb7..6339a5697 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -992,7 +992,7 @@ WarpX::ReadParameters () WarpX::field_boundary_hi[zdir] == FieldBoundaryType::Damped ) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::field_boundary_lo[zdir] == WarpX::field_boundary_hi[zdir], - "field boundary in both lo and high must be set to Damped for PSATD" + "field boundary in both lo and hi must be set to Damped for PSATD" ); } } |