From c78b9ab3a6cbf390e0ca32e80e00ff1d589cc0a4 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 12 Jul 2021 12:17:11 -0700 Subject: Fixes to RZ PSATD (#1945) * For RZ, changed the sign of the density corrections near the axis * Further fixes for deposition correction near axis * Yet one more sign fix for charge density * For RZ spectral solver, filled in the guard cells below the radial axis * Fix white space at end of line * In RZ spectral backtransform, ensure box is valid * For RZ inverse volume scaling, fixed use of nGrow to use nGrowVect * Temporary fix adding damped cells in the domain interior * Bug fix for RZ PSATD scalar backward transform * Fixes for damping of the fields in the z-guards * Bug fix in DampFieldsInGuards * Bug fix in DampFieldsInGuards (for tiling) * Added warpx_amr_check_input input parameter * Removed unneeded damp and zero_in_domain input * Removed damping related code from picmi * Improved some comments in code copying field to the radial guard cells * Update Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp Simplify the expression for the sign Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Updated benchmarks * Updated tolerance for Langmuir analysis script * Updated CI test galilean_rz_psatd_current_correction Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> --- .../SpectralSolver/SpectralFieldDataRZ.cpp | 80 +++++++++++++++++++--- 1 file changed, 72 insertions(+), 8 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp') 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(field_mf_copy[mfi], 0, i_comp*ncomp, ncomp); + + amrex::Array4 const & field_mf_array = field_mf[mfi].array(); + amrex::Array4 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 const & field_mf_r_copy_array = field_mf_r_copy[mfi].array(); amrex::Array4 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); } }); -- cgit v1.2.3