diff --git a/.github/workflows/wdmerger_collision-compare.yml b/.github/workflows/wdmerger_collision-compare.yml index 8f12f71c89..f8707b0ad6 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 plt00086 > 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 04727cea1d..4f6e0bf29a 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 = 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 diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 1b37fde218..3b19019039 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,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]; @@ -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.