ERF_InitCustomPert_MultiSpeciesBubble.H

Go to the documentation of this file.

 
    // On a mesh other than ConstantDz, heights come from z_cc and the vertical levels,
    //    not from geomdata's prob_lo[2] and CellSize(2)
    const bool use_z_phys = use_physical_z();
 
    ParmParse pp_prob("prob");
    Real T_0 = amrex::Real(300.0); pp_prob.query("T_0", T_0);
 
    // center of thermal perturbation
    Real x_c = zero; pp_prob.query("x_c", x_c);
    Real y_c = zero; pp_prob.query("y_c", y_c);
    Real z_c = zero; pp_prob.query("z_c", z_c);
 
    // radial extent of thermal perturbation
    Real x_r = zero; pp_prob.query("x_r", x_r);
    Real y_r = zero; pp_prob.query("y_r", y_r);
    Real z_r = zero; pp_prob.query("z_r", z_r);
 
    // perturbation temperature
    Real T_pert            = -amrex::Real(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  = two; pp_prob.query("theta_pert", theta_pert);
    Real eq_pot_temp = amrex::Real(320.0); pp_prob.query("eq_pot_temp",eq_pot_temp);
    Real qt_init = amrex::Real(0.02); pp_prob.query("qt_init", qt_init);
    bool use_empirical = false; pp_prob.query("use_empirical_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 = RhoQ6_comp+1;  // Start after all 6 SuperDroplets components
 
    // 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,zero);
    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);
 
        if (use_z_phys) {
            HSEutils::init_isentropic_hse_no_terrain(h_t.data(), h_r.data(), h_p.data(),
                                                     h_q_v.data(), m_zlevels_stag[lev].data(), khi, qt_init, eq_pot_temp, use_empirical);
        } else {
            HSEutils::init_isentropic_hse_no_terrain(h_t.data(), h_r.data(), h_p.data(),
                                                     h_q_v.data(), dz, geomdata.ProbLo(2), 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 + myhalf) * dx[0];
            const auto y = prob_lo[1] + (j + myhalf) * dx[1];
            const auto z = use_z_phys ? z_cc(i,j,k) : prob_lo[2] + (k + myhalf) * 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 = zero;
            if (x_r > 0) rad += amrex::Math::powi<2>((x - x_c)/x_r);
            if (y_r > 0) rad += amrex::Math::powi<2>((y - y_c)/y_r);
            if (z_r > 0) rad += amrex::Math::powi<2>((z - z_c)/z_r);
            rad = std::sqrt(rad);
 
            if (rad <= one){
                delta_theta = theta_pert*amrex::Math::powi<2>(std::cos(PI*rad/two));
            } else {
                delta_theta = zero;
            }
 
            theta_total  = theta_back[k]*(delta_theta/amrex::Real(300.0) + one);
            Real T = getTgivenPandTh(p_back[k], theta_total, RdoCp);
            rho    = p_back[k]/(R_d*T*(one + RvoRd*q_v_back[k]));
            RH     = one; // 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], RdoCp);
            Real rho_back = p_back[k]/(R_d*T_back*(one + RvoRd*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]*(one + RvoRd*q_v_back[k]);
            state_pert(i, j, k, Rho_comp)      = rho - rho_back*(one + qt_init);
 
            // mean states
            state_pert(i, j, k, RhoQ1_comp) = rho*q_v_hot;
        });
 
        // The kernel above reads d_r, d_p, d_t and d_q_v; wait for it before they are freed
        Gpu::streamSynchronize();
    } 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 + myhalf) * dx[0];
            const Real y = prob_lo[1] + (j + myhalf) * dx[1];
            const Real z = use_z_phys ? z_cc(i,j,k) : prob_lo[2] + (k + myhalf) * 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 = zero;
            if (x_r > 0) L += amrex::Math::powi<2>((x - x_c)/x_r);
            if (y_r > 0) L += amrex::Math::powi<2>((y - y_c)/y_r);
            if (z_r > 0) L += amrex::Math::powi<2>((z - z_c)/z_r);
            L = std::sqrt(L);
            Real dT;
            if (L > one) {
                dT = zero;
            }
            else {
                dT = T_pert * amrex::Math::powi<2>(std::cos(PI*L/two));
            }
 
            // 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) = zero; // 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)      = zero; // 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 + myhalf) * dx[0];
        const Real y = prob_lo[1] + (j + myhalf) * dx[1];
        const Real z = use_z_phys ? z_cc(i,j,k) : prob_lo[2] + (k + myhalf) * dx[2];
 
        amrex::Real L = zero;
        if (x_r > 0) L += amrex::Math::powi<2>((x - x_c)/x_r);
        if (y_r > 0) L += amrex::Math::powi<2>((y - y_c)/y_r);
        if (z_r > 0) L += amrex::Math::powi<2>((z - z_c)/z_r);
        L = std::sqrt(L);
 
        if (L < one) {
            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) = zero;
            }
        }
    });
 
    // The kernel above reads qv_init_d; wait for it before it is freed
    Gpu::streamSynchronize();