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

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.

Most problems will see no difference

.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.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
 

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.