Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add the ability to have C/O in the He envelope #57

Merged
merged 3 commits into from
Sep 15, 2023
Merged

add the ability to have C/O in the He envelope #57

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
355f3d4
allow for C/O in the He envelope
zingale Sep 15, 2023
8276839
this is an initial model that matches Rivas et al. mostly
zingale Sep 15, 2023
4d7770f
this seems to work
zingale Sep 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions sub_chandra/GNUmakefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,7 @@ EBASE := initialmodel
# EOS and network
EOS_DIR := helmholtz

NETWORK_DIR := subch
#NETWORK_INPUTS := triple_alpha_plus_o.net
NETWORK_DIR := subch_simple

Bpack := ./Make.package
Blocs := . ..
Expand Down
6 changes: 6 additions & 0 deletions sub_chandra/_parameters
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,12 @@ M_He real 0.2
# mass fraction of N14 to put in the accreted envelope
X_N14 real 0.00

# mass fraction of C12 to put in the accreted envelope
X_C12 real 0.00

# mass fraction of O16 to put in the accreted envelope
X_O16 real 0.00

delta real 1.e-6

xmin real 0.0
Expand Down
81 changes: 55 additions & 26 deletions sub_chandra/init_1d.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ AMREX_INLINE void init_1d() {

const Real TOL_HE_MASS = 1.e-3_rt;

const int MAX_ITER = 250;
const int MAX_ITER = 500;


// convert the envelope and WD mass into CGS
Expand Down Expand Up @@ -78,11 +78,18 @@ AMREX_INLINE void init_1d() {

if (problem_rp::X_N14 > 0.0_rt) {
xn_he[in14] = problem_rp::X_N14;
xn_he[ihe4] = 1.0_rt - (NumSpec - 2) * network_rp::small_x - problem_rp::X_N14;
} else {
xn_he[ihe4] = 1.0_rt - (NumSpec - 1) * network_rp::small_x;
}

if (problem_rp::X_C12 > 0.0_rt) {
xn_he[ic12] = problem_rp::X_C12;
}

if (problem_rp::X_O16 > 0.0_rt) {
xn_he[io16] = problem_rp::X_O16;
}

xn_he[ihe4] = 1.0_rt - xn_he[ic12] - xn_he[in14] - xn_he[io16] - (NumSpec-4) * network_rp::small_x;

// Create a 1-d uniform grid that is identical to the mesh that we
// are mapping onto, and then we want to force it into HSE on that
// mesh.
Expand Down Expand Up @@ -457,36 +464,44 @@ AMREX_INLINE void init_1d() {
} // end loop over zones


// compute the total mass of the He layer and C/O WD
mass_he = 0.0;
mass_wd = 0.0;

Real zone_X = model_hse(0, model::ispec+ihe4);
if (in14 > 0) {
zone_X += model_hse(0, model::ispec+in14);
}
for (int i = 0; i < icutoff; ++i) {

mass_he = (4.0_rt / 3.0_rt) * M_PI *
(std::pow(xznr(0), 3) - std::pow(xznl(0), 3)) *
model_hse(0, model::idens) * zone_X;
Real vol{0.0};
if (i == 0) {
vol = (4.0_rt / 3.0_rt) * M_PI * (std::pow(xznr(0), 3) - std::pow(xznl(0), 3));
} else {
vol = (4.0_rt / 3.0_rt) * M_PI *
(xznr(i) - xznl(i)) * (std::pow(xznr(i), 2) + xznl(i) * xznr(i) + std::pow(xznl(i), 2));
}

mass_wd = (4.0_rt / 3.0_rt) * M_PI *
(std::pow(xznr(0), 3) - std::pow(xznl(0), 3)) *
model_hse(0, model::idens) * (model_hse(0, model::ispec+ic12) + model_hse(0, model::ispec+io16));
// compute the total mass of the He layer and C/O WD
// note: only count the C and O if we are in the He layer

for (int i = 1; i < icutoff; ++i) {
zone_X = model_hse(i, model::ispec+ihe4);
if (in14 > 0) {
Real zone_X = model_hse(i, model::ispec+ihe4);
if (problem_rp::X_N14 > 0.0) {
zone_X += model_hse(i, model::ispec+in14);
}

mass_he += (4.0_rt / 3.0_rt) * M_PI *
(xznr(i) - xznl(i)) *
(std::pow(xznr(i), 2) + xznl(i) * xznr(i) + std::pow(xznl(i), 2)) *
model_hse(i,model::idens) * zone_X;
if (problem_rp::X_C12 > 0.0 && zone_X > 10.0 * network_rp::small_x) {
zone_X += model_hse(i, model::ispec+ic12);
}

if (problem_rp::X_O16 > 0.0 && zone_X > 10.0 * network_rp::small_x) {
zone_X += model_hse(i, model::ispec+io16);
}

mass_wd += (4.0_rt / 3.0_rt) * M_PI *
(xznr(i) - xznl(i)) *
(std::pow(xznr(i), 2) + xznl(i) * xznr(i) + std::pow(xznl(i), 2)) *
model_hse(i, model::idens) * (model_hse(i, model::ispec+ic12) + model_hse(i, model::ispec+io16));
Real core_X{0.0};
if (problem_rp::X_C12 == 0.0 && problem_rp::X_O16 == 0.0) {
core_X += model_hse(i, model::ispec+ic12) + model_hse(i, model::ispec+io16);
} else if (model_hse(i, model::ispec+ihe4) <= network_rp::small_x) {
core_X += model_hse(i, model::ispec+ic12) + model_hse(i, model::ispec+io16);
}

mass_he += vol * model_hse(i, model::idens) * zone_X;
mass_wd += vol * model_hse(i, model::idens) * core_X;
}

if (rho_c_old < 0.0_rt) {
Expand Down Expand Up @@ -571,10 +586,18 @@ AMREX_INLINE void init_1d() {
outfile += ".C";
}

if (problem_rp::X_C12 > 0.0_rt) {
outfile += ".C12";
}

if (problem_rp::X_N14 > 0.0_rt) {
outfile += ".N14";
}

if (problem_rp::X_O16 > 0.0_rt) {
outfile += ".O16";
}

write_model(outfile, xzn_hse, model_hse);

// extra info
Expand All @@ -587,9 +610,15 @@ AMREX_INLINE void init_1d() {
outfile += ".C";
}

if (problem_rp::X_C12 > 0.0_rt) {
outfile += ".C12";
}
if (problem_rp::X_N14 > 0.0_rt) {
outfile += ".N14";
}
if (problem_rp::X_O16 > 0.0_rt) {
outfile += ".O16";
}

outfile += "." + dxstr;

Expand Down
24 changes: 24 additions & 0 deletions sub_chandra/inputs.M_WD-0.981.M_He-0.05.CO.CO
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
problem.nx = 4096

problem.xmin = 0.0
problem.xmax = 4.096e9

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.low_density_cutoff = 1.e-4
problem.temp_fluff = 7.5e7

problem.small_temp = 1.e6