aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xExamples/Tests/flux_injection/analysis_flux_injection_3d.py10
-rw-r--r--Examples/Tests/flux_injection/inputs_3d6
-rw-r--r--Regression/Checksum/benchmarks_json/FluxInjection3D.json8
-rw-r--r--Source/Initialization/InjectorPosition.H20
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp17
-rw-r--r--Source/Utils/ParticleUtils.H16
6 files changed, 59 insertions, 18 deletions
diff --git a/Examples/Tests/flux_injection/analysis_flux_injection_3d.py b/Examples/Tests/flux_injection/analysis_flux_injection_3d.py
index 804cf95eb..048ef70f9 100755
--- a/Examples/Tests/flux_injection/analysis_flux_injection_3d.py
+++ b/Examples/Tests/flux_injection/analysis_flux_injection_3d.py
@@ -57,8 +57,12 @@ def gaussian_dist(u, u_th):
return 1./((2*np.pi)**.5*u_th) * np.exp(-u**2/(2*u_th**2) )
def gaussian_flux_dist(u, u_th, u_m):
- normalization_factor = u_th**2 * np.exp(-u_m**2/(2*u_th**2)) + (np.pi/2)**.5*u_m*u_th * (1 + erf(u_m/(2**.5*u_th)))
- return 1./normalization_factor * np.where( u>0, u * np.exp(-(u-u_m)**2/(2*u_th**2)), 0 )
+ au_m = np.abs(u_m)
+ normalization_factor = u_th**2 * np.exp(-au_m**2/(2*u_th**2)) + (np.pi/2)**.5*au_m*u_th * (1 + erf(au_m/(2**.5*u_th)))
+ result = 1./normalization_factor * np.where( u>0, u * np.exp(-(u-au_m)**2/(2*u_th**2)), 0 )
+ if u_m < 0.:
+ result = result[::-1]
+ return result
def compare_gaussian(u, w, u_th, label=''):
du = (hist_range[1]-hist_range[0])/hist_bins
@@ -103,7 +107,7 @@ uy = ad['proton','particle_momentum_y'].to_ndarray()/(m_p*c)
uz = ad['proton','particle_momentum_z'].to_ndarray()/(m_p*c)
w = ad['proton', 'particle_weight'].to_ndarray()
-compare_gaussian_flux(ux, w, u_th=0.1, u_m=0.05, label='u_x')
+compare_gaussian_flux(ux, w, u_th=0.1, u_m=-0.05, label='u_x')
compare_gaussian(uy, w, u_th=0.1, label='u_y')
compare_gaussian(uz, w, u_th=0.1, label='u_z')
plt.legend(loc=0)
diff --git a/Examples/Tests/flux_injection/inputs_3d b/Examples/Tests/flux_injection/inputs_3d
index f7cb80716..e3de9dc36 100644
--- a/Examples/Tests/flux_injection/inputs_3d
+++ b/Examples/Tests/flux_injection/inputs_3d
@@ -35,7 +35,7 @@ electron.charge = -q_e
electron.mass = m_e
electron.injection_style = NFluxPerCell
electron.num_particles_per_cell = 100
-electron.surface_flux_pos = -1.
+electron.surface_flux_pos = -4.
electron.flux_normal_axis = y
electron.flux_direction = +1
electron.profile = constant
@@ -50,9 +50,9 @@ proton.charge = +q_e
proton.mass = m_p
proton.injection_style = NFluxPerCell
proton.num_particles_per_cell = 100
-proton.surface_flux_pos = 1.
+proton.surface_flux_pos = 4.
proton.flux_normal_axis = x
-proton.flux_direction = +1
+proton.flux_direction = -1
proton.profile = constant
proton.density = 1.
proton.momentum_distribution_type = gaussianflux
diff --git a/Regression/Checksum/benchmarks_json/FluxInjection3D.json b/Regression/Checksum/benchmarks_json/FluxInjection3D.json
index 953fdb644..b2b373373 100644
--- a/Regression/Checksum/benchmarks_json/FluxInjection3D.json
+++ b/Regression/Checksum/benchmarks_json/FluxInjection3D.json
@@ -2,9 +2,9 @@
"electron": {
"particle_momentum_x": 1.1192116199394354e-18,
"particle_momentum_y": 2.238114590066897e-18,
- "particle_momentum_z": 1.1156457989239728e-18,
+ "particle_momentum_z": 1.1156457989239732e-18,
"particle_position_x": 102495.14197173176,
- "particle_position_y": 34752.73800291744,
+ "particle_position_y": 188132.22608016344,
"particle_position_z": 102423.13701045913,
"particle_weight": 8.959999999999998e-07
},
@@ -13,9 +13,9 @@
"particle_momentum_x": 3.835423016604918e-15,
"particle_momentum_y": 2.0468371931479925e-15,
"particle_momentum_z": 2.055186547721331e-15,
- "particle_position_x": 66743.84539580689,
+ "particle_position_x": 189256.1546041931,
"particle_position_y": 102293.00576740496,
"particle_position_z": 102314.93877691089,
"particle_weight": 8.959999999999998e-07
}
-} \ No newline at end of file
+}
diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H
index 059590979..012ea0ffc 100644
--- a/Source/Initialization/InjectorPosition.H
+++ b/Source/Initialization/InjectorPosition.H
@@ -175,7 +175,11 @@ struct InjectorPosition
};
}
- // bool: whether position specified is within bounds.
+ /* \brief Flags whether the point (x, y, z) is inside the plasma region
+ * or on the lower boundary
+ * \param x, y, z the point to check
+ * \returns bool flag
+ */
AMREX_GPU_HOST_DEVICE
bool
insideBounds (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
@@ -185,6 +189,20 @@ struct InjectorPosition
z < zmax and z >= zmin);
}
+ /* \brief Flags whether the point (x, y, z) is inside the plasma region
+ * or on the lower or upper boundary
+ * \param x, y, z the point to check
+ * \returns bool flag
+ */
+ AMREX_GPU_HOST_DEVICE
+ bool
+ insideBoundsInclusive (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
+ {
+ return (x <= xmax and x >= xmin and
+ y <= ymax and y >= ymin and
+ z <= zmax and z >= zmin);
+ }
+
// bool: whether the region defined by lo and hi overaps with the plasma region
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index cbb0a6932..18cf17162 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -32,6 +32,7 @@
#include "Particles/SpeciesPhysicalProperties.H"
#include "Particles/WarpXParticleContainer.H"
#include "Utils/Parser/ParserUtils.H"
+#include "Utils/ParticleUtils.H"
#include "Utils/Physics/IonizationEnergiesTable.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpXAlgorithmSelection.H"
@@ -1682,29 +1683,31 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
pu.y *= PhysConst::c;
pu.z *= PhysConst::c;
+ // The containsInclusive is used to allow the case of the flux surface
+ // being on the boundary of the domain. After the UpdatePosition below,
+ // the particles will be within the domain.
#if defined(WARPX_DIM_3D)
- if (!tile_realbox.contains(XDim3{ppos.x,ppos.y,ppos.z})) {
+ if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.y,ppos.z})) {
p.id() = -1;
continue;
}
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
amrex::ignore_unused(k);
- if (!tile_realbox.contains(XDim3{ppos.x,ppos.z,0.0_prt})) {
+ if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.z,0.0_prt})) {
p.id() = -1;
continue;
}
#else
amrex::ignore_unused(j,k);
- if (!tile_realbox.contains(XDim3{ppos.z,0.0_prt,0.0_prt})) {
+ if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.z,0.0_prt,0.0_prt})) {
p.id() = -1;
continue;
}
#endif
// Lab-frame simulation
- // If the particle is not within the species's
- // xmin, xmax, ymin, ymax, zmin, zmax, go to
- // the next generated particle.
- if (!inj_pos->insideBounds(ppos.x, ppos.y, ppos.z)) {
+ // If the particle's initial position is not within or on the species's
+ // xmin, xmax, ymin, ymax, zmin, zmax, go to the next generated particle.
+ if (!inj_pos->insideBoundsInclusive(ppos.x, ppos.y, ppos.z)) {
p.id() = -1;
continue;
}
diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H
index b9d8aa0ec..a997412f6 100644
--- a/Source/Utils/ParticleUtils.H
+++ b/Source/Utils/ParticleUtils.H
@@ -169,6 +169,22 @@ namespace ParticleUtils {
uy = y * vp;
uz = z * vp;
}
+
+ /* \brief Determines whether the point is within the tilebox, inclusive of the boundaries.
+ * Note that this routine is needed since tilebox.contains excludes the boundaries.
+ * \param[in] tilebox The tilebox being checked
+ * \param[in] point The point being checked
+ * \result true if the point with within the boundary, otherwise false
+ */
+ AMREX_GPU_HOST_DEVICE AMREX_INLINE
+ bool containsInclusive (amrex::RealBox const& tilebox, amrex::XDim3 const point) {
+ const auto xlo = tilebox.lo();
+ const auto xhi = tilebox.hi();
+ return AMREX_D_TERM((xlo[0] <= point.x) && (point.x <= xhi[0]),
+ && (xlo[1] <= point.y) && (point.y <= xhi[1]),
+ && (xlo[2] <= point.z) && (point.z <= xhi[2]));
+ }
+
}
#endif // WARPX_PARTICLE_UTILS_H_
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-11 20:31:30 +0000