implement a second hotspot for subchandra (#2544)

this applies to 3d simulations

Exec/science/subchandra/_prob_params

@@ -18,3 +18,13 @@ mu_p real -5.0_rt y
 mu_n real -11.0_rt y
 
 ihe4 integer -1
+
+# second perturbation angle -- if it is negative, then it is disabled
+# this is only relevant in 3D
+second_pert_angle real -1 y
+
+second_R_pert real 1.e7_rt y
+
+second_pert_temp_factor real 10.0_rt y
+
+second_pert_rad_factor real 2.0_rt y

Exec/science/subchandra/problem_initialize_state_data.H

@@ -70,7 +70,7 @@ void problem_initialize_state_data (int i, int j, int k,
  state(i,j,k,UMY) = 0.0_rt;
  state(i,j,k,UMZ) = 0.0_rt;
 
- // add a perturbation
+ // add a perturbation at the north pole
 
  Real t0 = state(i,j,k,UTEMP);
 
@@ -96,6 +96,36 @@ void problem_initialize_state_data (int i, int j, int k,
  (0.150e0_rt * (1.0_rt + std::tanh(2.0_rt - r1))));
 
 
+#if AMREX_SPACEDIM == 3
+ // optional second perturbation
+
+ if (problem::second_pert_angle > 0) {
+
+ Real angle_rad = problem::second_pert_angle * M_PI / 180.0_rt;
+
+ Real second_pert_center = problem::second_R_pert + problem::R_He_base;
+
+ // these are measured with respect to the center
+
+ Real x_pert = second_pert_center * std::sin(angle_rad);
+ Real y_pert = 0.0_rt;
+ Real z_pert = second_pert_center * std::cos(angle_rad);
+
+ Real dx_pert = x - x_pert;
+ Real dy_pert = y - y_pert;
+ Real dz_pert = z - z_pert;
+
+ Real r2 = std::sqrt(dx_pert * dx_pert + dy_pert * dy_pert + dz_pert * dz_pert) /
+ (2.5e6_rt * problem::second_pert_rad_factor);
+
+ // we are assuming here that the 2 perturbations don't overlap
+
+ state(i,j,k,UTEMP) = t0 * (1.0_rt + X_he * problem::second_pert_temp_factor *
+ (0.150e0_rt * (1.0_rt + std::tanh(2.0_rt - r2))));
+
+ }
+#endif
+
  burn_state.rho = state(i,j,k,URHO);
  burn_state.T = state(i,j,k,UTEMP);
  // we don't need to refill xn, since it still holds unchanged from above