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

fix the wdmerger case when the 2 stars are very close #2829

Merged
merged 8 commits into from
Apr 22, 2024
Merged

fix the wdmerger case when the 2 stars are very close #2829

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
52ce088
fix the wdmerger case when the 2 stars are very close
zingale Apr 20, 2024
380a69d
fix single star case
zingale Apr 20, 2024
8daa646
update benchmark
zingale Apr 20, 2024
161093b
also do velocities
zingale Apr 21, 2024
31c6057
Merge branch 'development' into fix_close_stars
zingale Apr 21, 2024
20b8781
fix benchmark
zingale Apr 21, 2024
956f363
fix vel
zingale Apr 21, 2024
2d0d90f
try again
zingale Apr 21, 2024
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
2 changes: 1 addition & 1 deletion .github/workflows/wdmerger_collision-compare.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,5 +48,5 @@ jobs:
- name: Check the extrema
run: |
cd Exec/science/wdmerger
../../../external/amrex/Tools/Plotfile/fextrema.gnu.ex plt00095 > fextrema.out
../../../external/amrex/Tools/Plotfile/fextrema.gnu.ex plt00086 > fextrema.out
diff fextrema.out ci-benchmarks/wdmerger_collision_2D.out
48 changes: 24 additions & 24 deletions Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,29 +1,29 @@
plotfile = plt00095
plotfile = plt00086
time = 1.25
variables minimum value maximum value
density 8.6596894178e-05 14342742.997
xmom -5.3795431399e+14 1.5534999809e+15
ymom -1.8900763787e+15 1.9125216773e+15
density 8.693611703e-05 19565534.299
xmom -5.4964100651e+14 1.3559128302e+14
ymom -2.5530096328e+15 2.5530122744e+15
zmom 0 0
rho_E 1.0478130927e+12 1.1816191282e+25
rho_e 9.9978448556e+11 1.1805863756e+25
Temp 100000 4617810497.4
rho_He4 8.6596894178e-17 36171.490906
rho_C12 3.4638757671e-05 5329342.135
rho_O16 5.1958136506e-05 8013661.8633
rho_Ne20 8.6596894178e-17 641339.21986
rho_Mg24 8.6596894178e-17 885671.32868
rho_Si28 8.6596894178e-17 7306845.0691
rho_S32 8.6596894178e-17 4315180.4037
rho_Ar36 8.6596894178e-17 1186734.5033
rho_Ca40 8.6596894178e-17 860488.91828
rho_Ti44 8.6596894178e-17 4978.4502978
rho_Cr48 8.6596894178e-17 17110.643453
rho_Fe52 8.6596894178e-17 110923.97673
rho_Ni56 8.6596894178e-17 827475.12913
phiGrav -4.6592040725e+17 -2.205539637e+16
grav_x -466332248.65 -48558.246229
grav_y -470363988.33 399464334.41
rho_E 7.4982062146e+11 5.0669247218e+24
rho_e 7.1077581849e+11 5.0640768325e+24
Temp 242288.68588 1409652233.6
rho_He4 8.693611703e-17 3.5999033031
rho_C12 3.4774446812e-05 7825956.6934
rho_O16 5.2161670217e-05 11739149.75
rho_Ne20 8.693611703e-17 181951.05725
rho_Mg24 8.693611703e-17 1192.7969771
rho_Si28 8.693611703e-17 6.6913703161
rho_S32 8.693611703e-17 0.00019493291757
rho_Ar36 8.693611703e-17 1.9565534609e-05
rho_Ca40 8.693611703e-17 1.9565534331e-05
rho_Ti44 8.693611703e-17 1.9565534308e-05
rho_Cr48 8.693611703e-17 1.9565534308e-05
rho_Fe52 8.693611703e-17 1.9565534308e-05
rho_Ni56 8.693611703e-17 1.9565534308e-05
phiGrav -5.8708033902e+17 -2.3375498341e+16
grav_x -684991644 -51428.243166
grav_y -739606241.84 739606820.44
grav_z 0 0
rho_enuc -4.0195280523e+22 2.0312948109e+26
rho_enuc -8.1740856019e+12 7.6429034917e+23

93 changes: 56 additions & 37 deletions Exec/science/wdmerger/problem_initialize_state_data.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,10 @@ void problem_initialize_state_data (int i, int j, int k,
// things right on the coarse levels. So we can still use the
// interpolation scheme, because it handles this special case
// for us by simply using the central zone of the model; we
// just need to make sure we catch it.
// just need to make sure we catch it. Finally, sometimes
// the stars are so close that the interpolation will overlap.
// in this case, look at which model has the highest density
// at that location and use that model.

const Real* dx = geomdata.CellSize();

Expand All @@ -53,43 +56,59 @@ void problem_initialize_state_data (int i, int j, int k,

eos_t zone_state;

if (problem::mass_P > 0.0_rt && (dist_P < problem::radius_P ||
(problem::radius_P <= max_dx && dist_P < max_dx))) {
Real pos[3] = {loc[0] - problem::center_P_initial[0],
loc[1] - problem::center_P_initial[1],
loc[2] - problem::center_P_initial[2]};
bool P_star_test = problem::mass_P > 0.0_rt &&
(dist_P < problem::radius_P ||
(problem::radius_P <= max_dx && dist_P < max_dx));

zone_state.rho = interpolate_3d(pos, dx, model::idens, problem::nsub, 0);
zone_state.T = interpolate_3d(pos, dx, model::itemp, problem::nsub, 0);
for (int n = 0; n < NumSpec; ++n) {
zone_state.xn[n] = interpolate_3d(pos, dx, model::ispec + n, problem::nsub, 0);
bool S_star_test = problem::mass_S > 0.0_rt &&
(dist_S < problem::radius_S ||
(problem::radius_S <= max_dx && dist_S < max_dx));

double rho_P{0.0};
double rho_S{0.0};

if (P_star_test || S_star_test) {

Real pos_P[3] = {loc[0] - problem::center_P_initial[0],
loc[1] - problem::center_P_initial[1],
loc[2] - problem::center_P_initial[2]};

if (problem::mass_P > 0.0_rt) {
rho_P = interpolate_3d(pos_P, dx, model::idens, problem::nsub, 0);
}
#ifdef AUX_THERMO
set_aux_comp_from_X(zone_state);
#endif

eos(eos_input_rt, zone_state);
Real pos_S[3] = {loc[0] - problem::center_S_initial[0],
loc[1] - problem::center_S_initial[1],
loc[2] - problem::center_S_initial[2]};

}
else if (problem::mass_S > 0.0_rt && (dist_S < problem::radius_S ||
(problem::radius_S <= max_dx && dist_S < max_dx))) {
Real pos[3] = {loc[0] - problem::center_S_initial[0],
loc[1] - problem::center_S_initial[1],
loc[2] - problem::center_S_initial[2]};

zone_state.rho = interpolate_3d(pos, dx, model::idens, problem::nsub, 1);
zone_state.T = interpolate_3d(pos, dx, model::itemp, problem::nsub, 1);
for (int n = 0; n < NumSpec; ++n) {
zone_state.xn[n] = interpolate_3d(pos, dx, model::ispec + n, problem::nsub, 1);
if (problem::mass_S > 0.0_rt) {
rho_S = interpolate_3d(pos_S, dx, model::idens, problem::nsub, 1);
}

if (rho_P > rho_S) {
// use the primary star initialization
zone_state.rho = rho_P;
zone_state.T = interpolate_3d(pos_P, dx, model::itemp, problem::nsub, 0);
for (int n = 0; n < NumSpec; ++n) {
zone_state.xn[n] = interpolate_3d(pos_P, dx, model::ispec + n, problem::nsub, 0);
}

} else {
// use the secondary star initialization
zone_state.rho = rho_S;
zone_state.T = interpolate_3d(pos_S, dx, model::itemp, problem::nsub, 1);
for (int n = 0; n < NumSpec; ++n) {
zone_state.xn[n] = interpolate_3d(pos_S, dx, model::ispec + n, problem::nsub, 1);
}
}

#ifdef AUX_THERMO
set_aux_comp_from_X(zone_state);
#endif

eos(eos_input_rt, zone_state);

}
else {
} else {

zone_state.rho = ambient::ambient_state[URHO];
zone_state.T = ambient::ambient_state[UTEMP];
Expand Down Expand Up @@ -131,17 +150,17 @@ void problem_initialize_state_data (int i, int j, int k,

if (problem::problem != 1) {

if (dist_P < problem::radius_P) {
state(i,j,k,UMX) += problem::vel_P[0] * state(i,j,k,URHO);
state(i,j,k,UMY) += problem::vel_P[1] * state(i,j,k,URHO);
state(i,j,k,UMZ) += problem::vel_P[2] * state(i,j,k,URHO);
}
else if (dist_S < problem::radius_S) {
state(i,j,k,UMX) += problem::vel_S[0] * state(i,j,k,URHO);
state(i,j,k,UMY) += problem::vel_S[1] * state(i,j,k,URHO);
state(i,j,k,UMZ) += problem::vel_S[2] * state(i,j,k,URHO);
if (P_star_test || S_star_test) {
if (rho_P > rho_S && problem::mass_P > 0.0_rt && dist_P < problem::radius_P) {
state(i,j,k,UMX) += problem::vel_P[0] * state(i,j,k,URHO);
state(i,j,k,UMY) += problem::vel_P[1] * state(i,j,k,URHO);
state(i,j,k,UMZ) += problem::vel_P[2] * state(i,j,k,URHO);
} else if (problem::mass_S > 0.0_rt && dist_S < problem::radius_S) {
state(i,j,k,UMX) += problem::vel_S[0] * state(i,j,k,URHO);
state(i,j,k,UMY) += problem::vel_S[1] * state(i,j,k,URHO);
state(i,j,k,UMZ) += problem::vel_S[2] * state(i,j,k,URHO);
}
}

}

// If we're in the inertial reference frame, use rigid body rotation with velocity omega x r.
Expand Down