aboutsummaryrefslogtreecommitdiff
path: root/Source/FieldSolver
diff options
context:
space:
mode:
Diffstat (limited to 'Source/FieldSolver')
-rw-r--r--Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp80
-rw-r--r--Source/FieldSolver/WarpXPushFieldsEM.cpp71
2 files changed, 113 insertions, 38 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp
index f62eaaa16..a87bfdb54 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp
@@ -545,7 +545,7 @@ SpectralFieldDataRZ::BackwardTransform (const int lev,
}
amrex::Real wt = amrex::second();
- amrex::Box const& realspace_bx = tempHTransformed[mfi].box();
+ amrex::Box realspace_bx = tempHTransformed[mfi].box();
FABZBackwardTransform(mfi, realspace_bx, field_index, tempHTransformedSplit, is_nodal_z);
@@ -553,7 +553,43 @@ SpectralFieldDataRZ::BackwardTransform (const int lev,
// tempHTransformedSplit includes the imaginary component of mode 0.
// field_mf does not.
multi_spectral_hankel_transformer[mfi].SpectralToPhysical_Scalar(tempHTransformedSplit[mfi], field_mf_copy[mfi]);
- field_mf[mfi].copy<amrex::RunOn::Device>(field_mf_copy[mfi], 0, i_comp*ncomp, ncomp);
+
+ amrex::Array4<amrex::Real> const & field_mf_array = field_mf[mfi].array();
+ amrex::Array4<amrex::Real> const & field_mf_copy_array = field_mf_copy[mfi].array();
+
+ // The box will be extended to include the guards cells below the axis
+ // so that they can be filled in. This will not be a simple copy of the
+ // fields since the signs will change when there is anti-symmetry.
+ amrex::Box const& realspace_bx_with_guards = field_mf[mfi].box();
+ const int* lo_with_guards = realspace_bx_with_guards.loVect();
+
+ // Grow the lower side of realspace_bx by the number of guard cells.
+ // This assumes that the box extends over the full extent radially, so
+ // lo_with_guards[0] will be equal to minus the number of guard cells radially.
+ const int nguard_r = -lo_with_guards[0];
+ realspace_bx.growLo(0, nguard_r);
+
+ // Get the intersection of the two boxes in case the field_mf has fewer z-guard cells
+ realspace_bx &= realspace_bx_with_guards;
+
+ ParallelFor(realspace_bx, ncomp,
+ [=] AMREX_GPU_DEVICE(int i, int j, int k, int icomp) noexcept {
+ int ii = i;
+ amrex::Real sign = +1._rt;
+ if (i < 0) {
+ ii = -i - 1;
+ if (icomp == 0) {
+ // Mode zero is symmetric
+ sign = +1._rt;
+ } else {
+ // Odd modes are anti-symmetric
+ int imode = (icomp + 1)/2;
+ sign = std::pow(-1._rt, imode);
+ }
+ }
+ int ic = icomp + i_comp;
+ field_mf_array(i,j,k,ic) = sign*field_mf_copy_array(ii,j,k,icomp);
+ });
if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
{
@@ -593,7 +629,7 @@ SpectralFieldDataRZ::BackwardTransform (const int lev,
}
amrex::Real wt = amrex::second();
- amrex::Box const& realspace_bx = tempHTransformed[mfi].box();
+ amrex::Box realspace_bx = tempHTransformed[mfi].box();
FABZBackwardTransform(mfi, realspace_bx, field_index_r, tempHTransformedSplit_p, is_nodal_z);
FABZBackwardTransform(mfi, realspace_bx, field_index_t, tempHTransformedSplit_m, is_nodal_z);
@@ -610,15 +646,43 @@ SpectralFieldDataRZ::BackwardTransform (const int lev,
amrex::Array4<amrex::Real> const & field_mf_r_copy_array = field_mf_r_copy[mfi].array();
amrex::Array4<amrex::Real> const & field_mf_t_copy_array = field_mf_t_copy[mfi].array();
+ // The box will be extended to include the guards cells below the axis
+ // so that they can be filled in. This will not be a simple copy of the
+ // fields since the signs will change when there is anti-symmetry.
+ amrex::Box const& realspace_bx_with_guards = field_mf_r[mfi].box();
+ const int* lo_with_guards = realspace_bx_with_guards.loVect();
+
+ // Grow the lower side of realspace_bx by the number of guard cells.
+ // This assumes that the box extends over the full extent radially, so
+ // lo_with_guards[0] will be equal to minus the number of guard cells radially.
+ const int nguard_r = -lo_with_guards[0];
+ realspace_bx.growLo(0, nguard_r);
+
+ // Get the intersection of the two boxes in case the field_mf has fewer z-guard cells
+ realspace_bx &= realspace_bx_with_guards;
+
ParallelFor(realspace_bx, 2*n_rz_azimuthal_modes-1,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int icomp) noexcept {
+ int ii = i;
+ amrex::Real sign = +1._rt;
+ if (i < 0) {
+ ii = -i - 1;
+ if (icomp == 0) {
+ // Mode zero is anti-symmetric
+ sign = -1._rt;
+ } else {
+ // Even modes are anti-symmetric
+ int imode = (icomp + 1)/2;
+ sign = std::pow(-1._rt, imode+1);
+ }
+ }
if (icomp == 0) {
- // mode 0
- field_mf_r_array(i,j,k,icomp) = field_mf_r_copy_array(i,j,k,icomp);
- field_mf_t_array(i,j,k,icomp) = field_mf_t_copy_array(i,j,k,icomp);
+ // mode zero
+ field_mf_r_array(i,j,k,icomp) = sign*field_mf_r_copy_array(ii,j,k,icomp);
+ field_mf_t_array(i,j,k,icomp) = sign*field_mf_t_copy_array(ii,j,k,icomp);
} else {
- field_mf_r_array(i,j,k,icomp) = field_mf_r_copy_array(i,j,k,icomp+1);
- field_mf_t_array(i,j,k,icomp) = field_mf_t_copy_array(i,j,k,icomp+1);
+ field_mf_r_array(i,j,k,icomp) = sign*field_mf_r_copy_array(ii,j,k,icomp+1);
+ field_mf_t_array(i,j,k,icomp) = sign*field_mf_t_copy_array(ii,j,k,icomp+1);
}
});
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 1e375f0eb..15ce940d6 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -751,7 +751,7 @@ WarpX::DampFieldsInGuards(std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield
void
WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, MultiFab* Jz, int lev)
{
- const long ngJ = Jx->nGrow();
+ const amrex::IntVect ngJ = Jx->nGrowVect();
const std::array<Real,3>& dx = WarpX::CellSize(lev);
const Real dr = dx[0];
@@ -793,16 +793,16 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
// Grow the tileboxes to include the guard cells, except for the
// guard cells at negative radius.
if (rmin > 0.) {
- tbr.growLo(0, ngJ);
- tbt.growLo(0, ngJ);
- tbz.growLo(0, ngJ);
+ tbr.growLo(0, ngJ[0]);
+ tbt.growLo(0, ngJ[0]);
+ tbz.growLo(0, ngJ[0]);
}
- tbr.growHi(0, ngJ);
- tbt.growHi(0, ngJ);
- tbz.growHi(0, ngJ);
- tbr.grow(1, ngJ);
- tbt.grow(1, ngJ);
- tbz.grow(1, ngJ);
+ tbr.growHi(0, ngJ[0]);
+ tbt.growHi(0, ngJ[0]);
+ tbz.growHi(0, ngJ[0]);
+ tbr.grow(1, ngJ[1]);
+ tbt.grow(1, ngJ[1]);
+ tbz.grow(1, ngJ[1]);
// Rescale current in r-z mode since the inverse volume factor was not
// included in the current deposition.
@@ -812,7 +812,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
// Wrap the current density deposited in the guard cells around
// to the cells above the axis.
// Note that Jr(i==0) is at 1/2 dr.
- if (rmin == 0. && 0 <= i && i < ngJ) {
+ if (rmin == 0. && 0 <= i && i < ngJ[0]) {
Jr_arr(i,j,0,0) -= Jr_arr(-1-i,j,0,0);
}
// Apply the inverse volume scaling
@@ -825,9 +825,9 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
// Wrap the current density deposited in the guard cells around
// to the cells above the axis.
// Note that Jr(i==0) is at 1/2 dr.
- if (rmin == 0. && 0 <= i && i < ngJ) {
- Jr_arr(i,j,0,2*imode-1) -= Jr_arr(-1-i,j,0,2*imode-1);
- Jr_arr(i,j,0,2*imode) -= Jr_arr(-1-i,j,0,2*imode);
+ if (rmin == 0. && 0 <= i && i < ngJ[0]) {
+ Jr_arr(i,j,0,2*imode-1) += std::pow(-1, imode+1)*Jr_arr(-1-i,j,0,2*imode-1);
+ Jr_arr(i,j,0,2*imode) += std::pow(-1, imode+1)*Jr_arr(-1-i,j,0,2*imode);
}
// Apply the inverse volume scaling
// Since Jr is never node centered in r, no need for distinction
@@ -842,7 +842,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
// to the cells above the axis.
// If Jt is node centered, Jt[0] is located on the boundary.
// If Jt is cell centered, Jt[0] is at 1/2 dr.
- if (rmin == 0. && 1-ishift_t <= i && i <= ngJ-ishift_t) {
+ if (rmin == 0. && 1-ishift_t <= i && i <= ngJ[0]-ishift_t) {
Jt_arr(i,j,0,0) -= Jt_arr(-ishift_t-i,j,0,0);
}
@@ -858,9 +858,9 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
for (int imode=1 ; imode < nmodes ; imode++) {
// Wrap the current density deposited in the guard cells around
// to the cells above the axis.
- if (rmin == 0. && 1-ishift_t <= i && i <= ngJ-ishift_t) {
- Jt_arr(i,j,0,2*imode-1) -= Jt_arr(-ishift_t-i,j,0,2*imode-1);
- Jt_arr(i,j,0,2*imode) -= Jt_arr(-ishift_t-i,j,0,2*imode);
+ if (rmin == 0. && 1-ishift_t <= i && i <= ngJ[0]-ishift_t) {
+ Jt_arr(i,j,0,2*imode-1) += std::pow(-1, imode+1)*Jt_arr(-ishift_t-i,j,0,2*imode-1);
+ Jt_arr(i,j,0,2*imode) += std::pow(-1, imode+1)*Jt_arr(-ishift_t-i,j,0,2*imode);
}
// Apply the inverse volume scaling
@@ -880,8 +880,8 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
// to the cells above the axis.
// If Jz is node centered, Jt[0] is located on the boundary.
// If Jz is cell centered, Jt[0] is at 1/2 dr.
- if (rmin == 0. && 1-ishift_z <= i && i <= ngJ-ishift_z) {
- Jz_arr(i,j,0,0) -= Jz_arr(-ishift_z-i,j,0,0);
+ if (rmin == 0. && 1-ishift_z <= i && i <= ngJ[0]-ishift_z) {
+ Jz_arr(i,j,0,0) += Jz_arr(-ishift_z-i,j,0,0);
}
// Apply the inverse volume scaling
@@ -896,9 +896,9 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
for (int imode=1 ; imode < nmodes ; imode++) {
// Wrap the current density deposited in the guard cells around
// to the cells above the axis.
- if (rmin == 0. && 1-ishift_z <= i && i <= ngJ-ishift_z) {
- Jz_arr(i,j,0,2*imode-1) -= Jz_arr(-ishift_z-i,j,0,2*imode-1);
- Jz_arr(i,j,0,2*imode) -= Jz_arr(-ishift_z-i,j,0,2*imode);
+ if (rmin == 0. && 1-ishift_z <= i && i <= ngJ[0]-ishift_z) {
+ Jz_arr(i,j,0,2*imode-1) -= std::pow(-1, imode+1)*Jz_arr(-ishift_z-i,j,0,2*imode-1);
+ Jz_arr(i,j,0,2*imode) -= std::pow(-1, imode+1)*Jz_arr(-ishift_z-i,j,0,2*imode);
}
// Apply the inverse volume scaling
@@ -919,7 +919,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
void
WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev)
{
- const long ngRho = Rho->nGrow();
+ const amrex::IntVect ngRho = Rho->nGrowVect();
const std::array<Real,3>& dx = WarpX::CellSize(lev);
const Real dr = dx[0];
@@ -949,23 +949,34 @@ WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev)
// Grow the tilebox to include the guard cells, except for the
// guard cells at negative radius.
if (rmin > 0.) {
- tb.growLo(0, ngRho);
+ tb.growLo(0, ngRho[0]);
}
- tb.growHi(0, ngRho);
- tb.grow(1, ngRho);
+ tb.growHi(0, ngRho[0]);
+ tb.grow(1, ngRho[1]);
// Rescale charge in r-z mode since the inverse volume factor was not
// included in the charge deposition.
// Note that the loop is also over ncomps, which takes care of the RZ modes,
// as well as the old and new rho.
- amrex::ParallelFor(tb, Rho->nComp(),
+ int const ncomp = Rho->nComp();
+ amrex::ParallelFor(tb, ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int /*k*/, int icomp)
{
// Wrap the charge density deposited in the guard cells around
// to the cells above the axis.
// Rho is located on the boundary
- if (rmin == 0. && 1-ishift <= i && i <= ngRho-ishift) {
- Rho_arr(i,j,0,icomp) -= Rho_arr(-ishift-i,j,0,icomp);
+ if (rmin == 0. && 1-ishift <= i && i <= ngRho[0]-ishift) {
+ int imode;
+ if (icomp == 0 || icomp == ncomp/2) {
+ imode = 0;
+ }
+ else if (icomp < ncomp/2) {
+ imode = (icomp+1)/2;
+ }
+ else {
+ imode = (icomp - ncomp/2 + 1)/2;
+ }
+ Rho_arr(i,j,0,icomp) -= std::pow(-1, imode+1)*Rho_arr(-ishift-i,j,0,icomp);
}
// Apply the inverse volume scaling
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-12 05:16:56 +0000