add an isothermal model

sub_chandra/init_1d.H

@@ -466,7 +466,13 @@ AMREX_INLINE void init_1d() {
         mass_he = 0.0;
         mass_wd = 0.0;
 
-        for (int i = 0; i < icutoff; ++i) {
+        // it might be that we never reach the cutoff density in our
+        // domain.  This is especially the case if we do an isothermal
+        // model.  Make sure we integrate over everything in that
+        // case.
+        int max_index = icutoff == -1 ? problem_rp::nx : icutoff;
+
+        for (int i = 0; i < max_index; ++i) {
 
             Real vol{0.0};
             if (i == 0) {
@@ -503,6 +509,12 @@ AMREX_INLINE void init_1d() {
             mass_wd += vol * model_hse(i, model::idens) * core_X;
         }
 
+        if (mass_wd == 0.0) {
+            std::string err_file = "zero_mass";
+            write_model(err_file, xzn_hse, model_hse);
+            amrex::Error("zero mass");
+        }
+
         if (rho_c_old < 0.0_rt) {
             // not enough iterations yet -- store the old central
             // density and mass and pick a new value
@@ -548,7 +560,6 @@ AMREX_INLINE void init_1d() {
             rho_he = amrex::min(1.5_rt * rho_he_old,
                                 amrex::max((rho_he + drho_he), 0.75_rt * rho_he_old));
 
-            std::cout << "current mass = " << mass_wd / C::M_solar << " " << mass_he / C::M_solar << std::endl;
         }
 
     } // end mass constraint loop

sub_chandra/inputs.M_WD-0.981.M_He-0.05.CO.CO.isothermal

@@ -0,0 +1,28 @@
+problem.nx = 6144
+
+problem.xmin = 0.0
+problem.xmax = 6.144e9
+
+problem.M_tot = 0.981
+problem.M_He  = 0.05
+
+problem.delta = 5.e6
+
+problem.temp_core = 5.e7
+problem.temp_base = 2.e8
+
+problem.mixed_co_wd = 1
+
+problem.X_C12 = 0.05
+problem.X_O16 = 0.05
+
+problem.isothermal_layer = 1
+
+problem.tol_hse = 1.e-9
+
+problem.low_density_cutoff = 1.e-4
+problem.temp_fluff = 7.5e7
+
+problem.small_temp = 1.e6
+
+