Better conservation check

AMReX-Fluids · Nov 28, 2023 · edbc077 · edbc077
1 parent e8b5f24
commit edbc077
Show file tree

Hide file tree

Showing 3 changed files with 33 additions and 12 deletions.
diff --git 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) <<std::endl;
+ static Real verr = 0.;
+ verr += (vol-prev_vol);
+ amrex::Print() << "Corr:Sum volume error " << verr << std::endl;
  //amrex::Print() << "Corr:Change in domain mass " << (sum - prev_mass) << std::endl;
- amrex::Print() << "Corr:Change in domain mass " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl;
- amrex::Print() << "Corr:Change speed " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl;
-
-
- //amrex::Print() << "Corr:Change in domain mass " << (sum - prev_mass - (2.45*.1*.001*16 - .1*.001*16) )/3. << std::endl;
+ // amrex::Print() << "Corr:Change in domain mass " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl;
+ // amrex::Print() << "Corr:Change speed " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl;
+
+
+ amrex::Print() << "Corr:Change in domain mass 2D " << (sum - (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.)) << std::endl;
+ // This only works if we start at time 0
+ // Real m_calc = (orig_mass + (2.45-1.)*dx[1]*ncelly*m_t_new[0]*1.);
+ // amrex::Print() << "Corr:Mass conservation error 2D " << (sum - m_calc) << std::endl;
+ Real m_calc = (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.);
+ static Real err = 0.;
+ err += (sum - m_calc);
+ amrex::Print() << "Corr:Sum conservation error " << err << std::endl;
  //if ( std::abs( (vol - prev_vol) - (sum - prev_mass) ) > 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
@@ -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
@@ -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<Real> 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<Real>(23-(i/Real(4.)));
+ den_arr(i,j,k) = 1.0 + 0.1 * static_cast<Real>(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<Real>(23-9);
  }
  if (j == 3) amrex::Print() << "DEN " << i << " " << den_arr(i,j,k) << std::endl;