Create alpha and beta that depend on the state. This gets the 3D

sphere const vel, rho case working.

For the nbhd in question, the MEB state redist equations create a
quadratic equation for beta. The code checks for this type of nbhd
and then solved the quadratic.

It does not treat another possible challenging nbhd that could
exist but doesn't seem to arise in our 3D sphere set up (at least
for a while). This more complex nbhd would require solving a cubic
equation.

Code needs a lot of clean up.

src/incflo_advance.cpp

@@ -108,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]*ncelly*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 :
@@ -157,30 +157,31 @@ void incflo::Advance(Real orig_mass, Real& prev_mass, Real& prev_vol)
  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: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 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 <<
-// " using " << sum << " " << orig_mass << " " << m_cur_time+m_dt << std::endl;
- // amrex::Print() << "Corr: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;
+ 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 <<
+// // " using " << sum << " " << orig_mass << " " << m_cur_time+m_dt << std::endl;
+// // amrex::Print() << "Corr: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;
+
+#endif
 
  prev_mass = sum;
  prev_vol = vol;
-#endif
 
  // Stop timing current time step
  Real end_step = static_cast<Real>(ParallelDescriptor::second()) - strt_step;