Merge branch 'development' into new_shock_flag_call

.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

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.7135211913e-05 13353524.435
- xmom -4.4634547423e+14 1.4982445308e+15
- ymom -1.8881098335e+15 1.9802028197e+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 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 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  -7.502262752e+21 6.4140001008e+26
+ rho_enuc -8.1740856019e+12 7.6429034917e+23
 

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.