ERF_InitCustomPert_ParticleTests.H

Go to the documentation of this file.

  ParmParse pp("prob");
  Real x_c =    0.0; pp.query("x_c", x_c);
  Real z_c = 3000.0; pp.query("z_c", z_c);
  Real x_r = 4000.0; pp.query("x_r", x_r);
  Real z_r = 2000.0; pp.query("z_r", z_r);
  Real T_pert = -15.0; pp.query("T_pert", T_pert);
 
  // Overridden physical constants
  Real C_p = 1004.0;
 
  const int khi = geomdata.Domain().bigEnd()[2];
 
  AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
 
  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
  {
    // Geometry (note we must include these here to get the data on device)
    const auto prob_lo  = geomdata.ProbLo();
    const auto dx       = geomdata.CellSize();
    const Real x = prob_lo[0] + (i + 0.5) * dx[0];
    const Real z = z_cc(i,j,k);
 
    // Temperature that satisfies the EOS given the hydrostatically balanced (r,p)
    const Real Tbar_hse = p_hse(i,j,k) / (R_d * r_hse(i,j,k));
 
    Real L = std::sqrt(
        std::pow((x - x_c)/x_r, 2) +
        std::pow((z - z_c)/z_r, 2)
    );
    if (L <= 1.0) {
        Real dT = T_pert * (std::cos(PI*L) + 1.0)/2.0;
 
        // Note: dT is a temperature perturbation, theta_perturbed is base state + perturbation in theta
        Real theta_perturbed = (Tbar_hse+dT)*std::pow(p_0/p_hse(i,j,k), R_d/C_p);
 
        // This version perturbs rho but not p
        state_pert(i, j, k, Rho_comp) = getRhoThetagivenP(p_hse(i,j,k)) / theta_perturbed - r_hse(i,j,k);
    }
  });