From 52ce08841c9bb4727b90c7c330372e01c67a6a6d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 20 Apr 2024 08:13:10 -0400 Subject: [PATCH 1/7] fix the wdmerger case when the 2 stars are very close the initialization was overlapping and doing bad things now we simply look at which model provides the largest density at a location and use that. --- .../wdmerger/problem_initialize_state_data.H | 68 +++++++++++-------- 1 file changed, 40 insertions(+), 28 deletions(-) diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 1b37fde218..85b9975d54 100644 --- a/Exec/science/wdmerger/problem_initialize_state_data.H +++ b/Exec/science/wdmerger/problem_initialize_state_data.H @@ -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(); @@ -53,43 +56,52 @@ 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); - } -#ifdef AUX_THERMO - set_aux_comp_from_X(zone_state); -#endif + bool S_star_test = problem::mass_S > 0.0_rt && + (dist_S < problem::radius_S || + (problem::radius_S <= max_dx && dist_S < max_dx)); - eos(eos_input_rt, zone_state); + if (P_star_test || S_star_test) { - } - 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); + 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]}; + + auto rho_P = interpolate_3d(pos_P, dx, model::idens, problem::nsub, 0); + + 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]}; + + auto 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]; From 380a69d86c39fe84eb1eca6c540fae3434bd2670 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 20 Apr 2024 10:11:18 -0400 Subject: [PATCH 2/7] fix single star case --- .../science/wdmerger/problem_initialize_state_data.H | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 85b9975d54..747f8442c4 100644 --- a/Exec/science/wdmerger/problem_initialize_state_data.H +++ b/Exec/science/wdmerger/problem_initialize_state_data.H @@ -70,13 +70,21 @@ void problem_initialize_state_data (int i, int j, int k, loc[1] - problem::center_P_initial[1], loc[2] - problem::center_P_initial[2]}; - auto rho_P = interpolate_3d(pos_P, dx, model::idens, problem::nsub, 0); + double rho_P{0.0}; + + if (problem::mass_P > 0.0_rt) { + rho_P = interpolate_3d(pos_P, dx, model::idens, problem::nsub, 0); + } 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]}; - auto rho_S = interpolate_3d(pos_S, dx, model::idens, problem::nsub, 1); + double rho_S{0.0}; + + 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 From 8daa646f276c21f335a36b5b39979cdb8b25b413 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 20 Apr 2024 10:39:19 -0400 Subject: [PATCH 3/7] update benchmark --- .../workflows/wdmerger_collision-compare.yml | 2 +- .../ci-benchmarks/wdmerger_collision_2D.out | 48 +++++++++---------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/.github/workflows/wdmerger_collision-compare.yml b/.github/workflows/wdmerger_collision-compare.yml index 8f12f71c89..d760dc8fc0 100644 --- a/.github/workflows/wdmerger_collision-compare.yml +++ b/.github/workflows/wdmerger_collision-compare.yml @@ -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 plt00088 > fextrema.out diff fextrema.out ci-benchmarks/wdmerger_collision_2D.out diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 33103755cb..c23719c88d 100644 --- a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out +++ b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out @@ -1,29 +1,29 @@ - plotfile = plt00095 + plotfile = plt00088 time = 1.25 variables minimum value maximum value - density 8.7135211913e-05 13353524.435 - xmom -4.4634547423e+14 1.4982445308e+15 - ymom -1.8881098335e+15 1.9802028197e+15 + density 8.6729292612e-05 7845757.9641 + xmom -4.0633354972e+14 1.551400564e+15 + ymom -1.3597988369e+15 2.2475296277e+15 zmom 0 0 - rho_E 7.6866508899e+11 5.6480764484e+24 - rho_e 7.3273369729e+11 5.6257109543e+24 - Temp 100000 3970499251.5 - rho_He4 8.7135211913e-17 1472.557835 - rho_C12 3.4854084765e-05 5034202.6702 - rho_O16 5.2281127147e-05 7781357.4533 - rho_Ne20 8.7135211913e-17 643712.05721 - rho_Mg24 8.7135211913e-17 1137247.8977 - rho_Si28 8.7135211913e-17 4245283.1554 - rho_S32 8.7135211913e-17 2178470.1687 - rho_Ar36 8.7135211913e-17 497477.29422 - rho_Ca40 8.7135211913e-17 382010.10154 - rho_Ti44 8.7135211913e-17 1577.2314864 - rho_Cr48 8.7135211913e-17 1468.2939271 - rho_Fe52 8.7135211913e-17 14839.719434 - rho_Ni56 8.7135211913e-17 182888.36194 - phiGrav -4.614866944e+17 -2.2055818115e+16 - grav_x -461232534.93 -48603.543258 - grav_y -444759175.11 392332867.05 + rho_E 9.2398321093e+11 5.3915021975e+24 + rho_e 8.4936113732e+11 5.3523796964e+24 + Temp 100000 3949829619.9 + rho_He4 8.6729292612e-17 1288.8271801 + rho_C12 3.4691717045e-05 2795592.1766 + rho_O16 5.2037575567e-05 4714582.9885 + rho_Ne20 8.6729292612e-17 550048.75584 + rho_Mg24 8.6729292612e-17 996221.69275 + rho_Si28 8.6729292612e-17 4108826.6962 + rho_S32 8.6729292612e-17 2116229.6136 + rho_Ar36 8.6729292612e-17 509049.28375 + rho_Ca40 8.6729292612e-17 472504.19598 + rho_Ti44 8.6729292612e-17 1071.9217721 + rho_Cr48 8.6729292612e-17 5082.0773428 + rho_Fe52 8.6729292612e-17 70257.464505 + rho_Ni56 8.6729292612e-17 1284867.6299 + phiGrav -4.8152316901e+17 -2.3376725137e+16 + grav_x -526515781.34 -51483.706401 + grav_y -565949608.7 391647841.27 grav_z 0 0 - rho_enuc -7.502262752e+21 6.4140001008e+26 + rho_enuc -2.7110182763e+22 5.5221803639e+26 From 161093b22e7df20a02750eda25c08cdd9c18bc1d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 20 Apr 2024 20:14:10 -0400 Subject: [PATCH 4/7] also do velocities --- .../wdmerger/problem_initialize_state_data.H | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 747f8442c4..8afe3f185c 100644 --- a/Exec/science/wdmerger/problem_initialize_state_data.H +++ b/Exec/science/wdmerger/problem_initialize_state_data.H @@ -64,14 +64,15 @@ void problem_initialize_state_data (int i, int j, int k, (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]}; - double rho_P{0.0}; - if (problem::mass_P > 0.0_rt) { rho_P = interpolate_3d(pos_P, dx, model::idens, problem::nsub, 0); } @@ -80,8 +81,6 @@ void problem_initialize_state_data (int i, int j, int k, loc[1] - problem::center_S_initial[1], loc[2] - problem::center_S_initial[2]}; - double rho_S{0.0}; - if (problem::mass_S > 0.0_rt) { rho_S = interpolate_3d(pos_S, dx, model::idens, problem::nsub, 1); } @@ -151,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) { + 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 { + 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. From 20b8781d6f354e61853012d2357fcd035b32bca2 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 20 Apr 2024 20:31:30 -0400 Subject: [PATCH 5/7] fix benchmark --- .../ci-benchmarks/wdmerger_collision_2D.out | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index faca118c8f..4f6e0bf29a 100644 --- a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out +++ b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out @@ -1 +1,29 @@ - slicing along x-direction and output to plt00086.slice + plotfile = plt00086 + time = 1.25 + variables minimum value maximum value + 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 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 -8.1740856019e+12 7.6429034917e+23 + From 956f3630ea4720c7c5ec3a0d66a4b7a3de522649 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 21 Apr 2024 07:39:09 -0400 Subject: [PATCH 6/7] fix vel --- Exec/science/wdmerger/problem_initialize_state_data.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 8afe3f185c..6817c3d43c 100644 --- a/Exec/science/wdmerger/problem_initialize_state_data.H +++ b/Exec/science/wdmerger/problem_initialize_state_data.H @@ -151,11 +151,11 @@ void problem_initialize_state_data (int i, int j, int k, if (problem::problem != 1) { if (P_star_test || S_star_test) { - if (rho_P > rho_S) { + if (rho_P > rho_S && problem::mass_P > 0.0_rt) { 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 { + } else if (problem::mass_S > 0.0_rt) { 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); From 2d0d90f79e19a33b471359d85919c895bab1a52e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 21 Apr 2024 11:40:22 -0400 Subject: [PATCH 7/7] try again --- Exec/science/wdmerger/problem_initialize_state_data.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 6817c3d43c..3b19019039 100644 --- a/Exec/science/wdmerger/problem_initialize_state_data.H +++ b/Exec/science/wdmerger/problem_initialize_state_data.H @@ -151,11 +151,11 @@ void problem_initialize_state_data (int i, int j, int k, if (problem::problem != 1) { if (P_star_test || S_star_test) { - if (rho_P > rho_S && problem::mass_P > 0.0_rt) { + 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) { + } 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);