diff --git a/Source/CAMR.cpp b/Source/CAMR.cpp index f359f98..7bb6799 100644 --- a/Source/CAMR.cpp +++ b/Source/CAMR.cpp @@ -1312,13 +1312,17 @@ CAMR::ZeroingOutForPlotting(amrex::MultiFab& S) amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { if (vfrac_arr(i,j,k) == 0.0) { - for (int n = 0; n < ncomp; ++n) { - Sarr(i, j, k, n) = 0.0; + for (int n = 0; n < ncomp; ++n) { + Sarr(i, j, k, n) = 0.0; } - } - { + } + else + { + // Clean up denormal-scale noise only. Note that we must test the + // magnitude here -- a signed test would zero every legitimately + // negative value, e.g. all negative momenta and velocities. for (int n = 0; n < ncomp; ++n) { - if(Sarr(i, j, k, n) < 1e-12){ + if (amrex::Math::abs(Sarr(i, j, k, n)) < 1.e-12) { Sarr(i, j, k, n) = 0.0; } } diff --git a/Source/Utils/CAMR_reset_internal_e.H b/Source/Utils/CAMR_reset_internal_e.H index 640e912..b8c9834 100644 --- a/Source/Utils/CAMR_reset_internal_e.H +++ b/Source/Utils/CAMR_reset_internal_e.H @@ -110,11 +110,16 @@ CAMR_reset_internal_e( S(i, j, k, UEDEN) = S(i,j,k,UEINT) + rho * ke; } else if (S(i,j,k,UEINT) <= 0.) { + // Note that we must take the delta against the (rho e) from + // before the floor is applied, so we hold on to it rather than + // relying on the order of the two writes below. + const amrex::Real rhoe_old = S(i, j, k, UEINT); + EOS::RTY2E(rho,small_T,massfrac,e_out); S(i, j, k, UEINT) = rho * e_out; if (l_dual_energy_update_E_from_e == 1) { - S(i, j, k, UEDEN) += (rho * e_out - S(i,j,k,UEINT)); + S(i, j, k, UEDEN) += (rho * e_out - rhoe_old); } } }