sub_chandra: allow for an isothermal envelope

.github/workflows/sub_chandra.yml

@@ -50,7 +50,7 @@ jobs:
       - name: Compare to stored output
         run: |
           cd sub_chandra
-          diff sub_chandra.M_WD-1.10.M_He-0.050.hse.CO.N14.10.00km ci-benchmarks/sub_chandra.M_WD-1.10.M_He-0.050.hse.CO.N14.10.00km
+          diff sub_chandra.M_WD-1.10.M_He-0.050.delta50.00km.hse.CO.N14.10.00km ci-benchmarks/sub_chandra.M_WD-1.10.M_He-0.050.delta50.00km.hse.CO.N14.10.00km
 
 
 

sub_chandra/_parameters

@@ -31,3 +31,11 @@ temp_fluff           real         1.e5
 
 # not used, but needed for the reader
 model_file           character    ""
+
+isothermal_layer     int          0
+
+# tol_hse is the tolerance used when iterating over a zone to
+# force it into HSE by adjusting the current density (and
+# possibly temperature).  TOL_HSE should be very small (~
+# 1.e-10).
+tol_hse              real         1.e-10

sub_chandra/init_1d.H

@@ -21,13 +21,6 @@ using namespace amrex;
 
 AMREX_INLINE void init_1d() {
 
-    // TOL_HSE is the tolerance used when iterating over a zone to
-    // force it into HSE by adjusting the current density (and
-    // possibly temperature).  TOL_HSE should be very small (~
-    // 1.e-10).
-
-    const Real TOL_HSE = 1.e-10_rt;
-
     // TOL_WD_MASS is tolerance used for getting the total WD mass
     // equal to M_tot (defined below).  It can be reasonably small,
     // since there will always be a central density value that can
@@ -42,7 +35,7 @@ AMREX_INLINE void init_1d() {
 
     const Real TOL_HE_MASS = 1.e-3_rt;
 
-    const int MAX_ITER = 500;
+    const int MAX_ITER = 750;
 
 
     // convert the envelope and WD mass into CGS
@@ -245,15 +238,21 @@ AMREX_INLINE void init_1d() {
 
                 } else {
 
-                    // fully isentropic
-
-                    if (ihe_entropy == -1) {
-                        ihe_entropy = i;
+                    if (problem_rp::isothermal_layer) {
+                        // isothermal He layer no matter what
                         temp_zone = problem_rp::temp_base;
                         isentropic = false;
                     } else {
-                        temp_zone = model_hse(i-1, model::itemp);
-                        isentropic = true;
+                        // fully isentropic
+
+                        if (ihe_entropy == -1) {
+                            ihe_entropy = i;
+                            temp_zone = problem_rp::temp_base;
+                            isentropic = false;
+                        } else {
+                            temp_zone = model_hse(i-1, model::itemp);
+                            isentropic = true;
+                        }
                     }
 
                     for (int n = 0; n < NumSpec; ++n) {
@@ -343,8 +342,8 @@ AMREX_INLINE void init_1d() {
 
                         }
 
-                        if (std::abs(drho) < TOL_HSE * dens_zone &&
-                            std::abs(dtemp) < TOL_HSE * temp_zone) {
+                        if (std::abs(drho) < problem_rp::tol_hse * dens_zone &&
+                            std::abs(dtemp) < problem_rp::tol_hse * temp_zone) {
                             converged_hse = true;
                             break;
                         }
@@ -383,7 +382,7 @@ AMREX_INLINE void init_1d() {
                             amrex::max(0.9_rt * dens_zone,
                                        amrex::min(dens_zone + drho, 1.1_rt * dens_zone));
 
-                        if (std::abs(drho) < TOL_HSE * dens_zone) {
+                        if (std::abs(drho) < problem_rp::tol_hse * dens_zone) {
                             converged_hse = true;
                             break;
                         }
@@ -411,8 +410,8 @@ AMREX_INLINE void init_1d() {
 
                     std::cout << "Error zone " << i <<  " did not converge in init_1d" << std::endl;
                     std::cout << dens_zone << " " << temp_zone << std::endl;
-                    std::cout << p_want;
-                    std::cout << drho;
+                    std::cout << p_want << std::endl;
+                    std::cout << drho << std::endl;
                     amrex::Error("Error: HSE non-convergence");
                 }
 
@@ -467,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) {
@@ -504,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
@@ -549,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
+
+