From 169109b023918f6197519ae3f6334eef2e4018bc Mon Sep 17 00:00:00 2001 From: David Grote <grote1@llnl.gov> Date: Fri, 13 Dec 2024 09:50:27 -0800 Subject: [PATCH] Fix the volume factor for overlapping boxes within the domain --- .../Diagnostics/ReducedDiags/FieldEnergy.cpp | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp index 30a0c7aa3c4..659b2ddb742 100644 --- a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp @@ -155,9 +155,6 @@ FieldEnergy::ComputeNorm2(amrex::MultiFab const& field, int lev) amrex::Geometry const & geom = warpx.Geom(lev); amrex::IntVect const is_nodal = field.ixType().toIntVect(); - amrex::Box const matched_domain_box = convert(geom.Domain(), is_nodal); - amrex::IntVect const domain_lo = matched_domain_box.smallEnd(); - amrex::IntVect const domain_hi = matched_domain_box.bigEnd(); amrex::ReduceOps<amrex::ReduceOpSum> reduce_ops; amrex::ReduceData<amrex::Real> reduce_data(reduce_ops); @@ -173,6 +170,8 @@ FieldEnergy::ComputeNorm2(amrex::MultiFab const& field, int lev) amrex::Box const tilebox = mfi.tilebox(); amrex::Box const tb = convert(tilebox, is_nodal); + amrex::IntVect const tb_lo = tb.smallEnd(); + amrex::IntVect const tb_hi = tb.bigEnd(); #if defined(WARPX_DIM_RZ) // Lower corner of tile box physical domain @@ -185,6 +184,12 @@ FieldEnergy::ComputeNorm2(amrex::MultiFab const& field, int lev) int const irmax = hi.x; #endif + // On the boundaries, if the grid is nodal, use half of the volume. + // This applies to all boundary conditions, and to the overlap of + // boxes within the domain. + // Previously, the code used MultiFab::norm2, but that does not do + // the half-volume scaling for the domain boundaries when not periodic. + auto volume_factor = [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { amrex::ignore_unused(i,j,k,n); #if defined WARPX_DIM_RZ @@ -195,21 +200,20 @@ FieldEnergy::ComputeNorm2(amrex::MultiFab const& field, int lev) } else if (rmin == 0._rt && i == irmax) { v_factor = r/2._rt - dr/8._rt; } - if (j == domain_lo[1] && is_nodal[1]) { v_factor *= 0.5_rt; } - if (j == domain_hi[1] && is_nodal[1]) { v_factor *= 0.5_rt; } + if (j == tb_lo[1] && is_nodal[1]) { v_factor *= 0.5_rt; } + if (j == tb_hi[1] && is_nodal[1]) { v_factor *= 0.5_rt; } amrex::Real const theta_integral = (n == 0 ? 2._rt : 1._rt); return MathConst::pi*v_factor*theta_integral; #else - // On the boundaries, if the grid is nodal, use half of the volume. amrex::Real v_factor = 1._rt; AMREX_D_TERM( - if (i == domain_lo[0] && is_nodal[0]) { v_factor *= 0.5_rt; }, - if (j == domain_lo[1] && is_nodal[1]) { v_factor *= 0.5_rt; }, - if (k == domain_lo[2] && is_nodal[2]) { v_factor *= 0.5_rt; }) + if (i == tb_lo[0] && is_nodal[0]) { v_factor *= 0.5_rt; }, + if (j == tb_lo[1] && is_nodal[1]) { v_factor *= 0.5_rt; }, + if (k == tb_lo[2] && is_nodal[2]) { v_factor *= 0.5_rt; }) AMREX_D_TERM( - if (i == domain_hi[0] && is_nodal[0]) { v_factor *= 0.5_rt; }, - if (j == domain_hi[1] && is_nodal[1]) { v_factor *= 0.5_rt; }, - if (k == domain_hi[2] && is_nodal[2]) { v_factor *= 0.5_rt; }) + if (i == tb_hi[0] && is_nodal[0]) { v_factor *= 0.5_rt; }, + if (j == tb_hi[1] && is_nodal[1]) { v_factor *= 0.5_rt; }, + if (k == tb_hi[2] && is_nodal[2]) { v_factor *= 0.5_rt; }) return v_factor; #endif };