diff --git a/.github/workflows/sub_chandra.yml b/.github/workflows/sub_chandra.yml index 7d7bb8a..5dfbb44 100644 --- a/.github/workflows/sub_chandra.yml +++ b/.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 diff --git a/sub_chandra/_parameters b/sub_chandra/_parameters index afeaa56..94bc348 100644 --- a/sub_chandra/_parameters +++ b/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 diff --git a/sub_chandra/ci-benchmarks/sub_chandra.M_WD-1.10.M_He-0.050.hse.CO.N14.10.00km b/sub_chandra/ci-benchmarks/sub_chandra.M_WD-1.10.M_He-0.050.delta50.00km.hse.CO.N14.10.00km similarity index 100% rename from sub_chandra/ci-benchmarks/sub_chandra.M_WD-1.10.M_He-0.050.hse.CO.N14.10.00km rename to sub_chandra/ci-benchmarks/sub_chandra.M_WD-1.10.M_He-0.050.delta50.00km.hse.CO.N14.10.00km diff --git a/sub_chandra/init_1d.H b/sub_chandra/init_1d.H index fd0e0ba..30ffc6e 100644 --- a/sub_chandra/init_1d.H +++ b/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 diff --git a/sub_chandra/inputs.M_WD-0.981.M_He-0.05.CO.CO.isothermal b/sub_chandra/inputs.M_WD-0.981.M_He-0.05.CO.CO.isothermal new file mode 100644 index 0000000..a1940aa --- /dev/null +++ b/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 + +