ERF_InitCustomPert_MultiSpeciesBubble.H

Go to the documentation of this file.

 
    ParmParse pp_prob("prob");
    Real T_0 = 300.0; pp_prob.query("T_0", T_0);
 
    // center of thermal perturbation
    Real x_c = 0.0; pp_prob.query("x_c", x_c);
    Real y_c = 0.0; pp_prob.query("y_c", y_c);
    Real z_c = 0.0; pp_prob.query("z_c", z_c);
 
    // radial extent of thermal perturbation
    Real x_r = 0.0; pp_prob.query("x_r", x_r);
    Real y_r = 0.0; pp_prob.query("y_r", y_r);
    Real z_r = 0.0; pp_prob.query("z_r", z_r);
 
    // perturbation temperature
    Real T_pert            = -15.0; pp_prob.query("T_pert", T_pert);
    bool T_pert_is_airtemp = true ; pp_prob.query("T_pert_is_airtemp", T_pert_is_airtemp);
    bool perturb_rho = true       ; pp_prob.query("perturb_rho", perturb_rho);
 
    // Moist bubble params
    bool do_moist_bubble = false; pp_prob.query("do_moist_bubble", do_moist_bubble);
    Real theta_pert  = 2.0; pp_prob.query("theta_pert", theta_pert);
    Real eq_pot_temp = 320.0; pp_prob.query("eq_pot_temp",eq_pot_temp);
    Real qt_init = 0.02; pp_prob.query("qt_init", qt_init);
    bool use_empirical = false; pp_prob.query("use_empircal_psat",use_empirical);
 
    std::vector<std::string> m_species;
    int m_num_species = 0;
    const int m_ncomp_per_species = 2; // rho*qc, rho*qc
    const int m_nstart_sp = RhoQ3_comp+1;
 
    // check for SDM species
    m_species.clear();
    ParmParse pp_sdm("super_droplets_moisture");
    std::string species_input = "species";
    if (pp_sdm.contains(species_input.c_str())) {
        int num_species = pp_sdm.countval(species_input.c_str());
        std::string sp_name;
        for (int i = 0; i < num_species; i++) {
            pp_sdm.get(species_input.c_str(), sp_name, i);
            m_species.push_back(sp_name);
        }
    }
 
    m_num_species = m_species.size();
 
    std::vector<amrex::Real> m_qv_init_species = std::vector<amrex::Real>(m_num_species,0.0);
    for (int i = 0; i < m_num_species; i++) {
        std::string key_str = "qv_init_" + m_species[i];
        pp_prob.query(key_str.c_str(), m_qv_init_species[i]);
    }
 
    AMREX_ALWAYS_ASSERT(sc.moisture_type == MoistureType::SuperDroplets);
 
    const auto khi     = geomdata.Domain().bigEnd()[2];
    const Real dz      = geomdata.CellSize()[2];
 
#ifndef AMREX_USE_GPU
    if (T_0 <= 0)
    {
        amrex::Print() << "Ignoring T_0 = " << T_0
                       << ", background fields should have been initialized with erf.init_type"
                       << std::endl;
    }
#endif
 
    const Real rdOcp   = sc.rdOcp;
 
    if (do_moist_bubble) {
        Vector<Real> h_r(khi+2);
        Vector<Real> h_p(khi+2);
        Vector<Real> h_t(khi+2);
        Vector<Real> h_q_v(khi+2);
 
        Gpu::DeviceVector<Real> d_r(khi+2);
        Gpu::DeviceVector<Real> d_p(khi+2);
        Gpu::DeviceVector<Real> d_t(khi+2);
        Gpu::DeviceVector<Real> d_q_v(khi+2);
 
        HSEutils::init_isentropic_hse_no_terrain(h_t.data(), h_r.data(), h_p.data(),
                                                 h_q_v.data(), dz, khi, qt_init, eq_pot_temp, use_empirical);
 
        Gpu::copy(Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin());
        Gpu::copy(Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin());
        Gpu::copy(Gpu::hostToDevice, h_t.begin(), h_t.end(), d_t.begin());
        Gpu::copy(Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin());
 
        Real* theta_back   = d_t.data();
        Real* p_back   = d_p.data();
        Real* q_v_back = d_q_v.data();
 
        int which_zone = -1;
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
        {
            // 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 auto x = prob_lo[0] + (i + 0.5) * dx[0];
            const auto y = prob_lo[1] + (j + 0.5) * dx[1];
            const auto z = prob_lo[2] + (k + 0.5) * dx[2];
 
            Real rad, delta_theta, theta_total, rho, RH;
 
            // Introduce the warm bubble.
            // Assume that the bubble is pressure matched with the background
            rad = 0.0;
            if (x_r > 0) rad += std::pow((x - x_c)/x_r, 2);
            if (y_r > 0) rad += std::pow((y - y_c)/y_r, 2);
            if (z_r > 0) rad += std::pow((z - z_c)/z_r, 2);
            rad = std::sqrt(rad);
 
            if (rad <= 1.0){
                delta_theta = theta_pert*std::pow(cos(PI*rad/2.0),2);
            } else {
                delta_theta = 0.0;
            }
 
            theta_total  = theta_back[k]*(delta_theta/300.0 + 1);
            Real T = getTgivenPandTh(p_back[k], theta_total, (R_d/Cp_d));
            rho    = p_back[k]/(R_d*T*(1.0 + (R_v/R_d)*q_v_back[k]));
            RH     = 1.0; // compute_relative_humidity();
            Real q_v_hot = HSEutils::vapor_mixing_ratio(p_back[k], T, RH, use_empirical, which_zone);
 
            // Compute background quantities
            Real T_back   = getTgivenPandTh(p_back[k], theta_back[k], (R_d/Cp_d));
            Real rho_back = p_back[k]/(R_d*T_back*(1.0 + (R_v/R_d)*q_v_back[k]));
 
            // This version perturbs rho but not p
            state_pert(i, j, k, RhoTheta_comp) = rho*theta_total - rho_back*theta_back[k]*(1.0 + (R_v/R_d)*q_v_back[k]);
            state_pert(i, j, k, Rho_comp)      = rho - rho_back*(1.0 + qt_init);
 
            // mean states
            state_pert(i, j, k, RhoQ1_comp) = rho*q_v_hot;
        });
    } else {
        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 y = prob_lo[1] + (j + 0.5) * dx[1];
            const Real z = prob_lo[2] + (k + 0.5) * dx[2];
 
            // perturb_rho_theta(x, y, z, p_hse(i,j,k), r_hse(i,j,k), rdOcp,
            //                   state_pert(i, j, k, Rho_comp),
            //                   state_pert(i, j, k, RhoTheta_comp));
 
            // Perturbation air temperature
            // - The bubble is either cylindrical (for 2-D problems, if two
            //   radial extents are specified) or an ellipsoid (if all three
            //   radial extents are specified).
            Real L = 0.0;
            if (x_r > 0) L += std::pow((x - x_c)/x_r, 2);
            if (y_r > 0) L += std::pow((y - y_c)/y_r, 2);
            if (z_r > 0) L += std::pow((z - z_c)/z_r, 2);
            L = std::sqrt(L);
            Real dT;
            if (L > 1.0) {
                dT = 0.0;
            }
            else {
                dT = T_pert * std::pow(cos(PI*L/2.0),2);
            }
 
            // 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));
 
            // Note: theta_perturbed is theta PLUS perturbation in theta
            Real theta_perturbed;
            if (T_pert_is_airtemp) {
                // dT is air temperature
                theta_perturbed = (Tbar_hse + dT)*std::pow(p_0/p_hse(i,j,k), rdOcp);
            } else {
                // dT is potential temperature
                theta_perturbed = Tbar_hse*std::pow(p_0/p_hse(i,j,k), rdOcp) + dT;
            }
 
            if (perturb_rho)
            {
                // this version perturbs rho but not p (i.e., rho*theta)
                // - hydrostatic rebalance is needed (TODO: is this true?)
                // - this is the approach taken in the density current problem
                state_pert(i,j,k,RhoTheta_comp) = 0.0; // i.e., hydrostatically balanced pressure stays const
                state_pert(i,j,k,Rho_comp     ) = getRhoThetagivenP(p_hse(i,j,k)) / theta_perturbed - r_hse(i,j,k);
            }
            else
            {
                // this version perturbs rho*theta (i.e., p) but not rho
                state_pert(i,j,k,Rho_comp)      = 0.0; // i.e., hydrostatically balanced density stays const
                state_pert(i,j,k,RhoTheta_comp) = r_hse(i,j,k) * theta_perturbed - getRhoThetagivenP(p_hse(i,j,k));
            }
 
        });
    } // do_moist_bubble
 
    // if other species are present
    auto n_sp = m_num_species;
    auto nstart = m_nstart_sp;
    int ncomp_species = state_pert.nComp() - nstart;
    AMREX_ALWAYS_ASSERT(n_sp*m_ncomp_per_species == ncomp_species);
    Gpu::DeviceVector<Real> qv_init_d(n_sp);
    Gpu::copy( Gpu::hostToDevice,
               m_qv_init_species.begin(),
               m_qv_init_species.end(),
               qv_init_d.begin() );
    auto qv_arr = qv_init_d.data();
 
    ParallelFor(bx, [=]
                AMREX_GPU_DEVICE(int i, int j, int k) noexcept
    {
        const auto prob_lo         = geomdata.ProbLo();
        const auto dx              = geomdata.CellSize();
 
        const Real x = prob_lo[0] + (i + 0.5) * dx[0];
        const Real y = prob_lo[1] + (j + 0.5) * dx[1];
        const Real z = prob_lo[2] + (k + 0.5) * dx[2];
 
        amrex::Real L = 0.0;
        if (x_r > 0) L += std::pow((x - x_c)/x_r, 2);
        if (y_r > 0) L += std::pow((y - y_c)/y_r, 2);
        if (z_r > 0) L += std::pow((z - z_c)/z_r, 2);
        L = std::sqrt(L);
 
        if (L < 1.0) {
            auto rho = state(i,j,k,Rho_comp) + state_pert(i,j,k,Rho_comp);
            for (int ns = 0; ns < n_sp; ns++) {
                state_pert(i, j, k, nstart+2*ns)   = rho*qv_arr[ns];
                state_pert(i, j, k, nstart+2*ns+1) = 0.0;
            }
        }
    });