fix the value of g in the inputs to match the initial_models value (#…

…2883)

also report the HSE error after reading in the model

Exec/science/subch_planar/inputs_2d

@@ -36,7 +36,7 @@ castro.pslope_cutoff_density = 1.e4
 castro.react_rho_min = 100.0
 
 gravity.gravity_type = ConstantGrav
-gravity.const_grav   = -1.466e9
+gravity.const_grav   = -1.1742e9
 
 # TIME STEP CONTROL
 castro.cfl            = 0.7     # cfl number for hyperbolic system

Exec/science/subch_planar/problem_initialize.H

@@ -21,12 +21,31 @@ void problem_initialize ()
     // Read initial model
     read_model_file(problem::model_name);
 
-    if (ParallelDescriptor::IOProcessor()) {
-        for (int i = 0; i < model::npts; i++) {
-            std::cout << i << " " << model::profile(0).r(i) << " " << model::profile(0).state(i,model::idens) << std::endl;
+    Real max_hse_error = std::numeric_limits<Real>::lowest();
+
+    for (int i = 0; i < model::npts-2; ++i) {
+        Real dpdr = (model::profile(0).state(i+1, model::ipres) -
+                     model::profile(0).state(i, model::ipres)) /
+            (model::profile(0).r(i+1) - model::profile(0).r(i));
+
+        Real rhog = 0.5_rt * (model::profile(0).state(i+1, model::idens) +
+                              model::profile(0).state(i, model::idens)) * gravity::const_grav;
+
+        // we'll check that dpdr in the next zone up is also negative, to ensure we are
+        // not hitting the part where we stopped integrating HSE
+        Real p_check = model::profile(0).state(i+2, model::ipres) -
+                       model::profile(0).state(i+1, model::ipres);
+
+        if (dpdr < 0 && p_check < 0) {
+            max_hse_error = amrex::max(std::abs(dpdr - rhog) / std::abs(dpdr),
+                                       max_hse_error);
         }
     }
 
+    amrex::Print() << Font::Bold << FGColor::Yellow
+                   << "HSE error in initial model = " << max_hse_error
+                   << ResetDisplay << std::endl;
+
     // set the ambient state for the upper boundary condition
 
     ambient::ambient_state[URHO] = model::profile(0).state(model::npts-1, model::idens);