From edbc077759997fce819e72ec91329a41ea97a60b Mon Sep 17 00:00:00 2001 From: cgilet Date: Tue, 28 Nov 2023 14:54:26 -0500 Subject: [PATCH] Better conservation check --- src/incflo_advance.cpp | 33 ++++++++++++++----- .../hydro_state_redistribute.cpp | 4 ++- .../test_conservation/prob_init_fluid.cpp | 8 +++-- 3 files changed, 33 insertions(+), 12 deletions(-) diff --git a/src/incflo_advance.cpp b/src/incflo_advance.cpp index f760ab2c..431ba8ca 100644 --- a/src/incflo_advance.cpp +++ b/src/incflo_advance.cpp @@ -93,6 +93,11 @@ void incflo::Advance(Real orig_mass, Real& prev_mass, Real& prev_vol) #endif auto const dx = geom[my_lev].CellSize(); + Box const& domain = geom[my_lev].Domain(); + const auto dlo = amrex::lbound(domain); + const auto dhi = amrex::ubound(domain); + const int ncelly = dhi.y-dlo.y+1; + #if (AMREX_SPACEDIM == 2) sum *= dx[0] * dx[1]; vol *= dx[0] * dx[1]; @@ -103,15 +108,15 @@ void incflo::Advance(Real orig_mass, Real& prev_mass, Real& prev_vol) amrex::Print() << "Pred:Change in domain volume " << (vol - prev_vol) << std::endl; // For plane geoms: - amrex::Print() << "Pred:Change in domain mass " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl; - amrex::Print() << "Pred:Change speed " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl; - amrex::Print() << "Pred:Change error " << (sum - prev_mass - vol + prev_vol)/ (vol-prev_vol) << std::endl; + // amrex::Print() << "Pred:Change in domain mass " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl; + // amrex::Print() << "Pred:Change speed " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl; + // amrex::Print() << "Pred:Change error " << (sum - prev_mass - vol + prev_vol)/ (vol-prev_vol) << std::endl; // amrex::Print() << "Pred:Change over time in last time step divided by time and area " << (sum - prev_mass) / m_dt / 1.6 << // " using " << sum << " " << prev_mass << " " << m_dt << std::endl; // For cyl : - // amrex::Print() << "Pred:Change in domain mass 2D " << (sum - (prev_mass + (2.45-1.)*dx[1]*64*m_dt*1.)) << std::endl; + amrex::Print() << "Pred:Change in domain mass 2D " << (sum - (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.)) << std::endl; // amrex::Print() << "Pred:Change speed " << (sum - prev_mass - (vol-prev_vol)*(2.45+1.625)/2.) / ((2.45-1.)* dx[1]*64 *m_dt) << std::endl; // For sphere : @@ -150,12 +155,22 @@ void incflo::Advance(Real orig_mass, Real& prev_mass, Real& prev_vol) #endif amrex::Print() << "Corr:Change in domain volume " << (vol - prev_vol) < 1e-14 ) Abort("Vol diff and dens diff are off"); // amrex::Print() << "Corr:Sum of mass at time = " << m_cur_time+m_dt << " " << sum << " " << std::endl; // amrex::Print() << "Change over time divided by time and area " << (sum - orig_mass) / (m_cur_time+m_dt) / 1.6 << diff --git a/src/redistribution/hydro_state_redistribute.cpp b/src/redistribution/hydro_state_redistribute.cpp index 7c044f55..81415d71 100644 --- a/src/redistribution/hydro_state_redistribute.cpp +++ b/src/redistribution/hydro_state_redistribute.cpp @@ -224,7 +224,9 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp, { // Add NU advective term to Qhat of it's merging partner // FIXME - this won't work for when SRD slopes are used; can't - // garantee the MP slope_soln_hat will be correct... + // garantee the MP slope_soln_hat will be correct because maybe this comes + // first and then the loop hits MP and slope_soln_hat gets reset, or maybe + // the order is opposite... soln_hat(r,s,t,n) += alpha(r,s,t,1)*U_in(i,j,k,n)/nrs(i,j,k)/nbhd_vol(r,s,t); slope_soln_hat(i,j,k,n) = soln_hat(r,s,t,n); diff --git a/test_2d_moving/test_conservation/prob_init_fluid.cpp b/test_2d_moving/test_conservation/prob_init_fluid.cpp index 4e5d70a3..66c69ee5 100644 --- a/test_2d_moving/test_conservation/prob_init_fluid.cpp +++ b/test_2d_moving/test_conservation/prob_init_fluid.cpp @@ -27,6 +27,10 @@ void incflo::prob_init_fluid (int lev) amrex::Print() <<" TYPE " << m_probtype << std::endl; + const auto dlo = amrex::lbound(domain); + const auto dhi = amrex::ubound(domain); + const int ncell = dhi.y-dlo.y+1; + for (MFIter mfi(ld.density); mfi.isValid(); ++mfi) { const Box& vbx = mfi.validbox(); @@ -35,11 +39,11 @@ void incflo::prob_init_fluid (int lev) Array4 den_arr = ld.density.array(mfi); amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - den_arr(i,j,k) = 1.0 + 0.1 * static_cast(23-(i/Real(4.))); + den_arr(i,j,k) = 1.0 + 0.1 * static_cast(23-(i/(ncell/16.))); if ( den_arr(i,j,k) < 1.0 ){ den_arr(i,j,k) = 1.0; } - if ( i < 9*4 ) { + if ( i < 9*(ncell/16) ) { den_arr(i,j,k) = 2.45; //1.0 + 0.1 * static_cast(23-9); } if (j == 3) amrex::Print() << "DEN " << i << " " << den_arr(i,j,k) << std::endl;