diff --git a/Source/ERF.H b/Source/ERF.H index e02d106d7..ea9c52dee 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -530,7 +530,7 @@ private: amrex::Vector rW_new; Microphysics micro; - amrex::Vector qmoist; // This has 6 components: qv, qc, qi, qr, qs, qg + amrex::Vector> qmoist; // This has up to 6 components: qv, qc, qi, qr, qs, qg #if defined(ERF_USE_RRTMGP) Radiation rad; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 1efc763a9..e2b51b376 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -548,7 +548,7 @@ ERF::InitData () if (solverChoice.moisture_type != MoistureType::None) { for (int lev = 0; lev <= finest_level; lev++) { - FillPatchMoistVars(lev, qmoist[lev]); + FillPatchMoistVars(lev, *(qmoist[lev])); } } } @@ -587,9 +587,9 @@ ERF::InitData () { // If not restarting we need to fill qmoist given qt and qp. if (restart_chkfile.empty()) { - micro.Init(vars_new[lev][Vars::cons], qmoist[lev], + micro.Init(vars_new[lev][Vars::cons], *(qmoist[lev]), grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose - micro.Update(vars_new[lev][Vars::cons], qmoist[lev]); + micro.Update(vars_new[lev][Vars::cons], *(qmoist[lev])); } } } @@ -856,6 +856,8 @@ ERF::init_only (int lev, Real time) lev_new[Vars::yvel].setVal(0.0); lev_old[Vars::yvel].setVal(0.0); lev_new[Vars::zvel].setVal(0.0); lev_old[Vars::zvel].setVal(0.0); + qmoist[lev]->setVal(0.); + // Initialize background flow (optional) if (init_type == "input_sounding") { // The base state is initialized by integrating vertically through the @@ -889,8 +891,6 @@ ERF::init_only (int lev, Real time) init_from_hse(lev); } - qmoist[lev].setVal(0.); - // Add problem-specific flow features // // Notes: @@ -1172,7 +1172,7 @@ ERF::MakeHorizontalAverages () if (use_moisture) { - MultiFab qv(qmoist[lev], make_alias, 0, 1); + MultiFab qv(*(qmoist[lev]), make_alias, 0, 1); for (MFIter mfi(mf); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index ff62dd201..b48b41b92 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -59,7 +59,8 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, //******************************************************************************************** // Microphysics // ******************************************************************************************* - qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg + int q_size = micro.Get_Qmoist_Size(); + qmoist[lev] = std::make_unique(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg // ******************************************************************************************** // Build the data structures for calculating diffusive/turbulent terms @@ -194,7 +195,8 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // Microphysics // ******************************************************************************************* int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; - qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg + int q_size = micro.Get_Qmoist_Size(); + qmoist[lev] = std::make_unique(ba, dm, q_size, ngrow_state); init_stuff(lev, ba, dm); @@ -296,7 +298,8 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp //******************************************************************************************** // Microphysics // ******************************************************************************************* - qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg + int q_size = micro.Get_Qmoist_Size(); + qmoist[lev]->define(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg init_stuff(lev,ba,dm); @@ -544,6 +547,9 @@ ERF::ClearLevel (int lev) rW_new[lev].clear(); rW_old[lev].clear(); + // Clears the qmoist memory + qmoist[lev].reset(); + // Clears the integrator memory mri_integrator_mem[lev].reset(); physbcs[lev].reset(); diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index 0fd1c7db5..effa60ba1 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -129,9 +129,9 @@ ERF::WriteCheckpointFile () const IntVect ng; if (solverChoice.moisture_type != MoistureType::None) { // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars - ng = qmoist[lev].nGrowVect(); - MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng); - MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),ng); + ng = qmoist[lev]->nGrowVect(); + MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng); + MultiFab::Copy(moist_vars,*(qmoist[lev]),0,0,qmoist[lev]->nComp(),ng); VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars")); } @@ -346,12 +346,12 @@ ERF::ReadCheckpointFile () IntVect ng; if (solverChoice.moisture_type == MoistureType::None) { - qmoist[lev].setVal(0.0); + qmoist[lev]->setVal(0.0); } else { - ng = qmoist[lev].nGrowVect(); - MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng); + ng = qmoist[lev]->nGrowVect(); + MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng); VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars")); - MultiFab::Copy(qmoist[lev],moist_vars,0,0,qmoist[lev].nComp(),ng); + MultiFab::Copy(*(qmoist[lev]),moist_vars,0,0,qmoist[lev]->nComp(),ng); } // Note that we read the ghost cells of the base state (unlike above) diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index cacc3d0b8..3b5c6ca24 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -246,7 +246,7 @@ ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVector plot_var_names) const Box& bx = mfi.tilebox(); const Array4& derdat = mf[lev].array(mfi); const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - const Array4& qv_arr = qmoist[0].const_array(mfi); + const Array4& qv_arr = qmoist[0]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -227,7 +227,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& derdat = mf[lev].array(mfi); const Array4& p0_arr = p_hse.const_array(mfi); const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - const Array4 & qv_arr = qmoist[0].const_array(mfi); + const Array4 & qv_arr = qmoist[0]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -593,12 +593,14 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } - MultiFab qv_mf(qmoist[lev], make_alias, 0, 1); - MultiFab qc_mf(qmoist[lev], make_alias, 1, 1); - MultiFab qi_mf(qmoist[lev], make_alias, 2, 1); - MultiFab qr_mf(qmoist[lev], make_alias, 3, 1); - MultiFab qs_mf(qmoist[lev], make_alias, 4, 1); - MultiFab qg_mf(qmoist[lev], make_alias, 5, 1); + // TODO: Protect against accessing non-existent data + int q_size = micro.Get_Qmoist_Size(); + MultiFab qv_mf(*(qmoist[lev]), make_alias, 0, 1); + MultiFab qc_mf(*(qmoist[lev]), make_alias, 1, 1); + MultiFab qi_mf(*(qmoist[lev]), make_alias, 2, 1); + MultiFab qr_mf(*(qmoist[lev]), make_alias, 3, 1); + MultiFab qs_mf(*(qmoist[lev]), make_alias, 4, 1); + MultiFab qg_mf(*(qmoist[lev]), make_alias, 5, 1); if (containerHasElement(plot_var_names, "qt")) { diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index fcc158133..2e2fb0586 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -25,7 +25,7 @@ void ERF::init_custom (int lev) { auto& lev_new = vars_new[lev]; - auto& qmoist_new = qmoist[lev]; + auto& qmoist_new = *(qmoist[lev]); MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component @@ -42,7 +42,7 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); - MultiFab qmoist_pert(qmoist[lev].boxArray(), qmoist[lev].DistributionMap(), 3, qmoist[lev].nGrow()); + MultiFab qmoist_pert(qmoist[lev]->boxArray(), qmoist[lev]->DistributionMap(), 3, qmoist[lev]->nGrow()); qmoist_pert.setVal(0.); MultiFab qv_pert(qmoist_pert, amrex::make_alias, 0, 1); diff --git a/Source/Initialization/InputSoundingData.H b/Source/Initialization/InputSoundingData.H index 0496fb2cb..1f18033da 100644 --- a/Source/Initialization/InputSoundingData.H +++ b/Source/Initialization/InputSoundingData.H @@ -79,7 +79,7 @@ public: iss_z >> z >> theta >> qv >> U >> V; if (z == 0) { AMREX_ALWAYS_ASSERT(theta == theta_inp_sound_tmp[0]); - AMREX_ALWAYS_ASSERT(qv == qv_inp_sound_tmp[0]); + AMREX_ALWAYS_ASSERT(qv*0.001 == qv_inp_sound_tmp[0]); U_inp_sound_tmp[0] = U; V_inp_sound_tmp[0] = V; } else { diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H index a60eb888d..e35354b62 100644 --- a/Source/Microphysics/FastEddy/FastEddy.H +++ b/Source/Microphysics/FastEddy/FastEddy.H @@ -34,109 +34,116 @@ namespace MicVar_FE { // class FastEddy : public NullMoist { - using FabPtr = std::shared_ptr; - - public: - // constructor - FastEddy () {} - - // destructor - virtual ~FastEddy () = default; - - // cloud physics - void AdvanceFE (); - - // diagnose - void Diagnose () override; - - // Set up for first time - void - define (SolverChoice& sc) override - { - docloud = sc.do_cloud; - doprecip = sc.do_precip; - m_fac_cond = lcond / sc.c_p; - m_fac_fus = lfus / sc.c_p; - m_fac_sub = lsub / sc.c_p; - m_gOcp = CONST_GRAV / sc.c_p; - m_axis = sc.ave_plane; - } - - // init - void - Init (const amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist, - const amrex::BoxArray& grids, - const amrex::Geometry& geom, - const amrex::Real& dt_advance) override; - - // update ERF variables - void - Update (amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist) override; - - // wrapper to do all the updating - void - Advance ( ) override - { - this->AdvanceFE(); - this->Diagnose(); - } + using FabPtr = std::shared_ptr; + +public: + // constructor + FastEddy () {} + + // destructor + virtual ~FastEddy () = default; + + // cloud physics + void AdvanceFE (); + + // diagnose + void Diagnose () override; + + // Set up for first time + void + define (SolverChoice& sc) override + { + docloud = sc.do_cloud; + doprecip = sc.do_precip; + m_fac_cond = lcond / sc.c_p; + m_fac_fus = lfus / sc.c_p; + m_fac_sub = lsub / sc.c_p; + m_gOcp = CONST_GRAV / sc.c_p; + m_axis = sc.ave_plane; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) override; + + // update ERF variables + void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist) override; + + // wrapper to do all the updating + void + Advance ( ) override + { + this->AdvanceFE(); + this->Diagnose(); + } + + int + Qmoist_Size ( ) override { return FastEddy::m_qmoist_size; } private: - // geometry - amrex::Geometry m_geom; - // valid boxes on which to evolve the solution - amrex::BoxArray m_gtoe; - - // timestep - amrex::Real dt; - - // number of vertical levels - int nlev, zlo, zhi; - - // plane average axis - int m_axis; - - // model options - bool docloud, doprecip; - - // constants - amrex::Real m_fac_cond; - amrex::Real m_fac_fus; - amrex::Real m_fac_sub; - amrex::Real m_gOcp; - - // microphysics parameters/coefficients - amrex::TableData accrrc; - amrex::TableData accrsi; - amrex::TableData accrsc; - amrex::TableData coefice; - amrex::TableData evaps1; - amrex::TableData evaps2; - amrex::TableData accrgi; - amrex::TableData accrgc; - amrex::TableData evapg1; - amrex::TableData evapg2; - amrex::TableData evapr1; - amrex::TableData evapr2; - - // vertical plane average data - amrex::TableData rho1d; - amrex::TableData pres1d; - amrex::TableData tabs1d; - amrex::TableData qt1d; - amrex::TableData qv1d; - amrex::TableData qn1d; - - // independent variables - amrex::Array mic_fab_vars; - - amrex::TableData gamaz; - amrex::TableData zmid; // mid value of vertical coordinate in physical domain - - // data (output) - amrex::TableData qifall; - amrex::TableData tlatqi; + // Number of qmoist variables + int m_qmoist_size = 2; + + // geometry + amrex::Geometry m_geom; + + // valid boxes on which to evolve the solution + amrex::BoxArray m_gtoe; + + // timestep + amrex::Real dt; + + // number of vertical levels + int nlev, zlo, zhi; + + // plane average axis + int m_axis; + + // model options + bool docloud, doprecip; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // microphysics parameters/coefficients + amrex::TableData accrrc; + amrex::TableData accrsi; + amrex::TableData accrsc; + amrex::TableData coefice; + amrex::TableData evaps1; + amrex::TableData evaps2; + amrex::TableData accrgi; + amrex::TableData accrgc; + amrex::TableData evapg1; + amrex::TableData evapg2; + amrex::TableData evapr1; + amrex::TableData evapr2; + + // vertical plane average data + amrex::TableData rho1d; + amrex::TableData pres1d; + amrex::TableData tabs1d; + amrex::TableData qt1d; + amrex::TableData qv1d; + amrex::TableData qn1d; + + // independent variables + amrex::Array mic_fab_vars; + + amrex::TableData gamaz; + amrex::TableData zmid; // mid value of vertical coordinate in physical domain + + // data (output) + amrex::TableData qifall; + amrex::TableData tlatqi; }; #endif diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp index 40f04132d..a07fdee14 100644 --- a/Source/Microphysics/FastEddy/FastEddy.cpp +++ b/Source/Microphysics/FastEddy/FastEddy.cpp @@ -11,61 +11,62 @@ using namespace amrex; /** * Compute Precipitation-related Microphysics quantities. */ -void FastEddy::AdvanceFE() { +void FastEddy::AdvanceFE () { - auto tabs = mic_fab_vars[MicVar_FE::tabs]; + auto tabs = mic_fab_vars[MicVar_FE::tabs]; - // get the temperature, dentisy, theta, qt and qc from input - for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); - auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + // get the temperature, dentisy, theta, qt and qc from input + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); - const auto& box3d = mfi.tilebox(); + const auto& box3d = mfi.tilebox(); - // Expose for GPU - Real d_fac_cond = m_fac_cond; + // Expose for GPU + Real d_fac_cond = m_fac_cond; - ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { - qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); + qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); - //------- Autoconversion/accretion - Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; + //------- Autoconversion/accretion + Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; - Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; - erf_qsatw(tabs_array(i,j,k), pressure, qsat); + Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; + erf_qsatw(tabs_array(i,j,k), pressure, qsat); - // If there is precipitating water (i.e. rain), and the cell is not saturated - // then the rain water can evaporate leading to extraction of latent heat, hence - // reducing temperature and creating negative buoyancy + // If there is precipitating water (i.e. rain), and the cell is not saturated + // then the rain water can evaporate leading to extraction of latent heat, hence + // reducing temperature and creating negative buoyancy - dq_vapor_to_clwater = 0.0; - dq_clwater_to_vapor = 0.0; + dq_vapor_to_clwater = 0.0; + dq_clwater_to_vapor = 0.0; - //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); - Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); + //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); + Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); - // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature - if(qv_array(i,j,k) > qsat){ - dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); - } - // If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and - // reducing temperature - if(qv_array(i,j,k) < qsat and qc_array(i,j,k) > 0.0){ - dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); - } + // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature + if(qv_array(i,j,k) > qsat){ + dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); + } + // If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and + // reducing temperature + if(qv_array(i,j,k) < qsat and qc_array(i,j,k) > 0.0){ + dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); + } - qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor; - qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor; + qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor; + qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor; - theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor); + theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor); - qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k)); - qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); + qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k)); + qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); - }); - } + }); + } } diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp index 7708a1243..164016fc9 100644 --- a/Source/Microphysics/FastEddy/Init_FE.cpp +++ b/Source/Microphysics/FastEddy/Init_FE.cpp @@ -23,54 +23,55 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist, const BoxArray& grids, const Geometry& geom, const Real& dt_advance) - { - m_geom = geom; - m_gtoe = grids; +{ + m_geom = geom; + m_gtoe = grids; - dt = dt_advance; + dt = dt_advance; - // initialize microphysics variables - for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) { - mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); - mic_fab_vars[ivar]->setVal(0.); - } + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) { + mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); + mic_fab_vars[ivar]->setVal(0.); + } - // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled - // The ghost cells of these arrays aren't filled in the boundary condition calls for the state + // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled + // The ghost cells of these arrays aren't filled in the boundary condition calls for the state - qmoist.setVal(0.); + qmoist.setVal(0.); - for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { + for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { - const auto& box3d = mfi.tilebox(); + const auto& box3d = mfi.tilebox(); - const auto& lo = amrex::lbound(box3d); - const auto& hi = amrex::ubound(box3d); + const auto& lo = amrex::lbound(box3d); + const auto& hi = amrex::ubound(box3d); - nlev = box3d.length(2); - zlo = lo.z; - zhi = hi.z; + nlev = box3d.length(2); + zlo = lo.z; + zhi = hi.z; } - // Get the temperature, density, theta, qt and qp from input - for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { - auto states_array = cons_in.array(mfi); + // Get the temperature, density, theta, qt and qp from input + for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { + auto states_array = cons_in.array(mfi); - auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); - const auto& box3d = mfi.tilebox(); + const auto& box3d = mfi.tilebox(); - // Get pressure, theta, temperature, density, and qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho_array(i,j,k) = states_array(i,j,k,Rho_comp); - theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); - temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); - }); - } + // Get pressure, theta, temperature, density, and qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + rho_array(i,j,k) = states_array(i,j,k,Rho_comp); + theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); + qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); + }); + } } diff --git a/Source/Microphysics/FastEddy/Update_FE.cpp b/Source/Microphysics/FastEddy/Update_FE.cpp index 07d9a86c3..eaedc3697 100644 --- a/Source/Microphysics/FastEddy/Update_FE.cpp +++ b/Source/Microphysics/FastEddy/Update_FE.cpp @@ -14,27 +14,28 @@ void FastEddy::Update (amrex::MultiFab& cons, amrex::MultiFab& qmoist) { - // Get the temperature, density, theta, qt and qp from input - for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto states_arr = cons.array(mfi); - - auto rho_arr = mic_fab_vars[MicVar_FE::rho]->array(mfi); - auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto qv_arr = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto qc_arr = mic_fab_vars[MicVar_FE::qc]->array(mfi); - const auto& box3d = mfi.tilebox(); - - // get potential total density, temperature, qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); - states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); - }); - } - - // Fill interior ghost cells and periodic boundaries - cons.FillBoundary(m_geom.periodicity()); - qmoist.FillBoundary(m_geom.periodicity()); + // Get the temperature, density, theta, qt and qp from input + for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto qv_arr = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto qc_arr = mic_fab_vars[MicVar_FE::qc]->array(mfi); + const auto& box3d = mfi.tilebox(); + + // get potential total density, temperature, qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); + qmoist.FillBoundary(m_geom.periodicity()); } diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp index 43210f501..c231053dd 100644 --- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -6,36 +6,37 @@ */ void Kessler::Diagnose () { - auto qt = mic_fab_vars[MicVar_Kess::qt]; - auto qp = mic_fab_vars[MicVar_Kess::qp]; - auto qv = mic_fab_vars[MicVar_Kess::qv]; - auto qn = mic_fab_vars[MicVar_Kess::qn]; - auto qcl = mic_fab_vars[MicVar_Kess::qcl]; - auto qci = mic_fab_vars[MicVar_Kess::qci]; - auto qpl = mic_fab_vars[MicVar_Kess::qpl]; - auto qpi = mic_fab_vars[MicVar_Kess::qpi]; - auto tabs = mic_fab_vars[MicVar_Kess::tabs]; + auto qt = mic_fab_vars[MicVar_Kess::qt]; + auto qp = mic_fab_vars[MicVar_Kess::qp]; + auto qv = mic_fab_vars[MicVar_Kess::qv]; + auto qn = mic_fab_vars[MicVar_Kess::qn]; + auto qcl = mic_fab_vars[MicVar_Kess::qcl]; + auto qci = mic_fab_vars[MicVar_Kess::qci]; + auto qpl = mic_fab_vars[MicVar_Kess::qpl]; + auto qpi = mic_fab_vars[MicVar_Kess::qpi]; + auto tabs = mic_fab_vars[MicVar_Kess::tabs]; - for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto qt_array = qt->array(mfi); - auto qp_array = qp->array(mfi); - auto qv_array = qv->array(mfi); - auto qn_array = qn->array(mfi); - auto qcl_array = qcl->array(mfi); - auto qci_array = qci->array(mfi); - auto qpl_array = qpl->array(mfi); - auto qpi_array = qpi->array(mfi); + for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qt_array = qt->array(mfi); + auto qp_array = qp->array(mfi); + auto qv_array = qv->array(mfi); + auto qn_array = qn->array(mfi); + auto qcl_array = qcl->array(mfi); + auto qci_array = qci->array(mfi); + auto qpl_array = qpl->array(mfi); + auto qpi_array = qpi->array(mfi); - const auto& box3d = mfi.tilebox(); + const auto& box3d = mfi.tilebox(); - amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - qv_array(i,j,k) = std::max(0.0,qt_array(i,j,k) - qn_array(i,j,k)); - amrex::Real omn = 1.0; - qcl_array(i,j,k) = qn_array(i,j,k)*omn; - qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn); - amrex::Real omp = 1.0;; - qpl_array(i,j,k) = qp_array(i,j,k)*omp; - qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp); - }); - } + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + qv_array(i,j,k) = std::max(0.0,qt_array(i,j,k) - qn_array(i,j,k)); + amrex::Real omn = 1.0; + qcl_array(i,j,k) = qn_array(i,j,k)*omn; + qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn); + amrex::Real omp = 1.0;; + qpl_array(i,j,k) = qp_array(i,j,k)*omp; + qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp); + }); + } } diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp index f8a3d835f..fcf469a2a 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -19,231 +19,236 @@ using namespace amrex; * @param[in] geom Geometry associated with these MultiFabs and grids * @param[in] dt_advance Timestep for the advance */ -void Kessler::Init (const MultiFab& cons_in, MultiFab& qmoist, - const BoxArray& grids, - const Geometry& geom, - const Real& dt_advance) +void Kessler::Init (const MultiFab& cons_in, + MultiFab& qmoist, + const BoxArray& grids, + const Geometry& geom, + const Real& dt_advance) { - m_geom = geom; - m_gtoe = grids; - - auto dz = m_geom.CellSize(2); - auto lowz = m_geom.ProbLo(2); - - dt = dt_advance; - - // initialize microphysics variables - for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) { - mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); - mic_fab_vars[ivar]->setVal(0.); - } - - // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled - // The ghost cells of these arrays aren't filled in the boundary condition calls for the state - - //qmoist.setVal(0.); - - for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { - - const auto& box3d = mfi.tilebox(); - - const auto& lo = amrex::lbound(box3d); - const auto& hi = amrex::ubound(box3d); - - nlev = box3d.length(2); - zlo = lo.z; - zhi = hi.z; - - // parameters - accrrc.resize({zlo}, {zhi}); - accrsi.resize({zlo}, {zhi}); - accrsc.resize({zlo}, {zhi}); - coefice.resize({zlo}, {zhi}); - evaps1.resize({zlo}, {zhi}); - evaps2.resize({zlo}, {zhi}); - accrgi.resize({zlo}, {zhi}); - accrgc.resize({zlo}, {zhi}); - evapg1.resize({zlo}, {zhi}); - evapg2.resize({zlo}, {zhi}); - evapr1.resize({zlo}, {zhi}); - evapr2.resize({zlo}, {zhi}); - - // data (input) - rho1d.resize({zlo}, {zhi}); - pres1d.resize({zlo}, {zhi}); - tabs1d.resize({zlo}, {zhi}); - gamaz.resize({zlo}, {zhi}); - zmid.resize({zlo}, {zhi}); - - // output - qifall.resize({zlo}, {zhi}); - tlatqi.resize({zlo}, {zhi}); - } - - auto accrrc_t = accrrc.table(); - auto accrsi_t = accrsi.table(); - auto accrsc_t = accrsc.table(); - auto coefice_t = coefice.table(); - auto evaps1_t = evaps1.table(); - auto evaps2_t = evaps2.table(); - auto accrgi_t = accrgi.table(); - auto accrgc_t = accrgc.table(); - auto evapg1_t = evapg1.table(); - auto evapg2_t = evapg2.table(); - auto evapr1_t = evapr1.table(); - auto evapr2_t = evapr2.table(); - - auto rho1d_t = rho1d.table(); - auto pres1d_t = pres1d.table(); - auto tabs1d_t = tabs1d.table(); - - auto gamaz_t = gamaz.table(); - auto zmid_t = zmid.table(); - - Real gam3 = erf_gammafff(3.0 ); - Real gamr1 = erf_gammafff(3.0+b_rain ); - Real gamr2 = erf_gammafff((5.0+b_rain)/2.0); - // Real gamr3 = erf_gammafff(4.0+b_rain ); - Real gams1 = erf_gammafff(3.0+b_snow ); - Real gams2 = erf_gammafff((5.0+b_snow)/2.0); - // Real gams3 = erf_gammafff(4.0+b_snow ); - Real gamg1 = erf_gammafff(3.0+b_grau ); - Real gamg2 = erf_gammafff((5.0+b_grau)/2.0); - // Real gamg3 = erf_gammafff(4.0+b_grau ); - - amrex::MultiFab qv(qmoist, amrex::make_alias, 0, 1); - amrex::MultiFab qc(qmoist, amrex::make_alias, 1, 1); - amrex::MultiFab qi(qmoist, amrex::make_alias, 2, 1); - - // Get the temperature, density, theta, qt and qp from input - for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { - auto states_array = cons_in.array(mfi); - auto qv_array_from_moist = qv.array(mfi); - auto qc_array = qc.array(mfi); - auto qi_array = qi.array(mfi); - - auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); - auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); - auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi); - - const auto& box3d = mfi.tilebox(); - - // Get pressure, theta, temperature, density, and qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho_array(i,j,k) = states_array(i,j,k,Rho_comp); - theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); - qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); - qv_array(i,j,k) = qv_array_from_moist(i,j,k); - temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); - pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; + m_geom = geom; + m_gtoe = grids; + + auto dz = m_geom.CellSize(2); + auto lowz = m_geom.ProbLo(2); + + dt = dt_advance; + + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) { + mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); + mic_fab_vars[ivar]->setVal(0.); + } + + // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled + // The ghost cells of these arrays aren't filled in the boundary condition calls for the state + + //qmoist.setVal(0.); + + for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { + + const auto& box3d = mfi.tilebox(); + + const auto& lo = amrex::lbound(box3d); + const auto& hi = amrex::ubound(box3d); + + nlev = box3d.length(2); + zlo = lo.z; + zhi = hi.z; + + // parameters + accrrc.resize({zlo}, {zhi}); + accrsi.resize({zlo}, {zhi}); + accrsc.resize({zlo}, {zhi}); + coefice.resize({zlo}, {zhi}); + evaps1.resize({zlo}, {zhi}); + evaps2.resize({zlo}, {zhi}); + accrgi.resize({zlo}, {zhi}); + accrgc.resize({zlo}, {zhi}); + evapg1.resize({zlo}, {zhi}); + evapg2.resize({zlo}, {zhi}); + evapr1.resize({zlo}, {zhi}); + evapr2.resize({zlo}, {zhi}); + + // data (input) + rho1d.resize({zlo}, {zhi}); + pres1d.resize({zlo}, {zhi}); + tabs1d.resize({zlo}, {zhi}); + gamaz.resize({zlo}, {zhi}); + zmid.resize({zlo}, {zhi}); + + // output + qifall.resize({zlo}, {zhi}); + tlatqi.resize({zlo}, {zhi}); + } + + auto accrrc_t = accrrc.table(); + auto accrsi_t = accrsi.table(); + auto accrsc_t = accrsc.table(); + auto coefice_t = coefice.table(); + auto evaps1_t = evaps1.table(); + auto evaps2_t = evaps2.table(); + auto accrgi_t = accrgi.table(); + auto accrgc_t = accrgc.table(); + auto evapg1_t = evapg1.table(); + auto evapg2_t = evapg2.table(); + auto evapr1_t = evapr1.table(); + auto evapr2_t = evapr2.table(); + + auto rho1d_t = rho1d.table(); + auto pres1d_t = pres1d.table(); + auto tabs1d_t = tabs1d.table(); + + auto gamaz_t = gamaz.table(); + auto zmid_t = zmid.table(); + + Real gam3 = erf_gammafff(3.0 ); + Real gamr1 = erf_gammafff(3.0+b_rain ); + Real gamr2 = erf_gammafff((5.0+b_rain)/2.0); + // Real gamr3 = erf_gammafff(4.0+b_rain ); + Real gams1 = erf_gammafff(3.0+b_snow ); + Real gams2 = erf_gammafff((5.0+b_snow)/2.0); + // Real gams3 = erf_gammafff(4.0+b_snow ); + Real gamg1 = erf_gammafff(3.0+b_grau ); + Real gamg2 = erf_gammafff((5.0+b_grau)/2.0); + // Real gamg3 = erf_gammafff(4.0+b_grau ); + + amrex::MultiFab qv(qmoist, amrex::make_alias, 0, 1); + amrex::MultiFab qc(qmoist, amrex::make_alias, 1, 1); + amrex::MultiFab qi(qmoist, amrex::make_alias, 2, 1); + + // Get the temperature, density, theta, qt and qp from input + for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { + auto states_array = cons_in.array(mfi); + auto qv_array_from_moist = qv.array(mfi); + auto qc_array = qc.array(mfi); + auto qi_array = qi.array(mfi); + + auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + auto temp_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + // Get pressure, theta, temperature, density, and qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + rho_array(i,j,k) = states_array(i,j,k,Rho_comp); + theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); + qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); + qv_array(i,j,k) = qv_array_from_moist(i,j,k); + temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; + }); + } + + // calculate the plane average variables + PlaneAverage cons_ave(&cons_in, m_geom, m_axis); + cons_ave.compute_averages(ZDir(), cons_ave.field()); + + // get host variable rho, and rhotheta + int ncell = cons_ave.ncell_line(); + + Gpu::HostVector rho_h(ncell), rhotheta_h(ncell); + cons_ave.line_average(Rho_comp, rho_h); + cons_ave.line_average(RhoTheta_comp, rhotheta_h); + + // copy data to device + Gpu::DeviceVector rho_d(ncell), rhotheta_d(ncell); + Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, rhotheta_h.begin(), rhotheta_h.end(), rhotheta_d.begin()); + Gpu::streamSynchronize(); + + Real* rho_dptr = rho_d.data(); + Real* rhotheta_dptr = rhotheta_d.data(); + + Real gOcp = m_gOcp; + + amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept + { + Real pressure = getPgivenRTh(rhotheta_dptr[k]); + rho1d_t(k) = rho_dptr[k]; + pres1d_t(k) = pressure/100.; + tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],rhotheta_dptr[k]); + zmid_t(k) = lowz + (k+0.5)*dz; + gamaz_t(k) = gOcp*zmid_t(k); }); - } - // calculate the plane average variables - PlaneAverage cons_ave(&cons_in, m_geom, m_axis); - cons_ave.compute_averages(ZDir(), cons_ave.field()); - - // get host variable rho, and rhotheta - int ncell = cons_ave.ncell_line(); - - Gpu::HostVector rho_h(ncell), rhotheta_h(ncell); - cons_ave.line_average(Rho_comp, rho_h); - cons_ave.line_average(RhoTheta_comp, rhotheta_h); - - // copy data to device - Gpu::DeviceVector rho_d(ncell), rhotheta_d(ncell); - Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, rhotheta_h.begin(), rhotheta_h.end(), rhotheta_d.begin()); - Gpu::streamSynchronize(); - - Real* rho_dptr = rho_d.data(); - Real* rhotheta_dptr = rhotheta_d.data(); - - Real gOcp = m_gOcp; - - amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { - Real pressure = getPgivenRTh(rhotheta_dptr[k]); - rho1d_t(k) = rho_dptr[k]; - pres1d_t(k) = pressure/100.; - tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],rhotheta_dptr[k]); - zmid_t(k) = lowz + (k+0.5)*dz; - gamaz_t(k) = gOcp*zmid_t(k); - }); - - // This fills qv - //Diagnose(); + // This fills qv + //Diagnose(); #if 0 - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) { - fluxbmk(l,j,i) = 0.0; - fluxtmk(l,j,i) = 0.0; - }); - - amrex::ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int k, int l, int i0) { - mkwle (l,k) = 0.0; - mkwsb (l,k) = 0.0; - mkadv (l,k) = 0.0; - mkdiff(l,k) = 0.0; - }); + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) + { + fluxbmk(l,j,i) = 0.0; + fluxtmk(l,j,i) = 0.0; + }); + + amrex::ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int k, int l, int i0) + { + mkwle (l,k) = 0.0; + mkwsb (l,k) = 0.0; + mkadv (l,k) = 0.0; + mkdiff(l,k) = 0.0; + }); #endif - if(round(gam3) != 2) { - std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; - std::exit(-1); - } - - amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { - - Real pratio = sqrt(1.29 / rho1d_t(k)); - Real rrr1=393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5); - Real rrr2=std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k)); - Real estw = 100.0*erf_esatw(tabs1d_t(k)); - Real esti = 100.0*erf_esati(tabs1d_t(k)); - - // accretion by snow: - Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0)); - Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); - accrsi_t(k) = coef1 * coef2 * esicoef; - accrsc_t(k) = coef1 * esccoef; - coefice_t(k) = coef2; - - // evaporation of snow: - coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); - evaps1_t(k) = 0.65*4.0*nzeros/sqrt(PI*rhos*nzeros)/(coef1+coef2)/sqrt(rho1d_t(k)); - evaps2_t(k) = 0.49*4.0*nzeros*gams2*sqrt(a_snow/(muelq*rrr1))/pow((PI*rhos*nzeros) , ((5.0+b_snow)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_snow)/8.0))*sqrt(pratio); - - // accretion by graupel: - coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0)); - coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); - accrgi_t(k) = coef1 * coef2 * egicoef; - accrgc_t(k) = coef1 * egccoef; - - // evaporation of graupel: - coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); - evapg1_t(k) = 0.65*4.0*nzerog/sqrt(PI*rhog*nzerog)/(coef1+coef2)/sqrt(rho1d_t(k)); - evapg2_t(k) = 0.49*4.0*nzerog*gamg2*sqrt(a_grau/(muelq*rrr1))/pow((PI * rhog * nzerog) , ((5.0+b_grau)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_grau)/8.0))*sqrt(pratio); - - // accretion by rain: - accrrc_t(k)= 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3+b_rain)/4.))* erccoef; - - // evaporation of rain: - coef1 =(lcond/(tabs1d_t(k)*rv)-1.0)*lcond/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq * rrr2 * estw); - evapr1_t(k) = 0.78 * 2.0 * PI * nzeror / - sqrt(PI * rhor * nzeror) / (coef1+coef2) / sqrt(rho1d_t(k)); - evapr2_t(k) = 0.31 * 2.0 * PI * nzeror * gamr2 * 0.89 * sqrt(a_rain/(muelq*rrr1))/ - pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); - }); + if(round(gam3) != 2) { + std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; + std::exit(-1); + } + + amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept + { + Real pratio = sqrt(1.29 / rho1d_t(k)); + Real rrr1=393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5); + Real rrr2=std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k)); + Real estw = 100.0*erf_esatw(tabs1d_t(k)); + Real esti = 100.0*erf_esati(tabs1d_t(k)); + + // accretion by snow: + Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0)); + Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); + accrsi_t(k) = coef1 * coef2 * esicoef; + accrsc_t(k) = coef1 * esccoef; + coefice_t(k) = coef2; + + // evaporation of snow: + coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); + evaps1_t(k) = 0.65*4.0*nzeros/sqrt(PI*rhos*nzeros)/(coef1+coef2)/sqrt(rho1d_t(k)); + evaps2_t(k) = 0.49*4.0*nzeros*gams2*sqrt(a_snow/(muelq*rrr1))/pow((PI*rhos*nzeros) , ((5.0+b_snow)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_snow)/8.0))*sqrt(pratio); + + // accretion by graupel: + coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0)); + coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); + accrgi_t(k) = coef1 * coef2 * egicoef; + accrgc_t(k) = coef1 * egccoef; + + // evaporation of graupel: + coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); + evapg1_t(k) = 0.65*4.0*nzerog/sqrt(PI*rhog*nzerog)/(coef1+coef2)/sqrt(rho1d_t(k)); + evapg2_t(k) = 0.49*4.0*nzerog*gamg2*sqrt(a_grau/(muelq*rrr1))/pow((PI * rhog * nzerog) , ((5.0+b_grau)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_grau)/8.0))*sqrt(pratio); + + // accretion by rain: + accrrc_t(k)= 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3+b_rain)/4.))* erccoef; + + // evaporation of rain: + coef1 =(lcond/(tabs1d_t(k)*rv)-1.0)*lcond/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq * rrr2 * estw); + evapr1_t(k) = 0.78 * 2.0 * PI * nzeror / + sqrt(PI * rhor * nzeror) / (coef1+coef2) / sqrt(rho1d_t(k)); + evapr2_t(k) = 0.31 * 2.0 * PI * nzeror * gamr2 * 0.89 * sqrt(a_rain/(muelq*rrr1))/ + pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); + }); } diff --git a/Source/Microphysics/Kessler/Kessler.H b/Source/Microphysics/Kessler/Kessler.H index 7ef58dd34..f771575e0 100644 --- a/Source/Microphysics/Kessler/Kessler.H +++ b/Source/Microphysics/Kessler/Kessler.H @@ -51,109 +51,115 @@ namespace MicVar_Kess { // class Kessler : public NullMoist { - using FabPtr = std::shared_ptr; - - public: - // constructor - Kessler () {} - - // destructor - virtual ~Kessler () = default; - - // cloud physics - void AdvanceKessler (); - - // diagnose - void Diagnose () override; - - // Set up for first time - void - define (SolverChoice& sc) override - { - docloud = sc.do_cloud; - doprecip = sc.do_precip; - m_fac_cond = lcond / sc.c_p; - m_fac_fus = lfus / sc.c_p; - m_fac_sub = lsub / sc.c_p; - m_gOcp = CONST_GRAV / sc.c_p; - m_axis = sc.ave_plane; - } - - // init - void - Init (const amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist, - const amrex::BoxArray& grids, - const amrex::Geometry& geom, - const amrex::Real& dt_advance) override; - - // update ERF variables - void - Update (amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist) override; - - // wrapper to do all the updating - void - Advance ( ) override - { - this->AdvanceKessler(); - this->Diagnose(); - } + using FabPtr = std::shared_ptr; + +public: + // constructor + Kessler () {} + + // destructor + virtual ~Kessler () = default; + + // cloud physics + void AdvanceKessler (); + + // diagnose + void Diagnose () override; + + // Set up for first time + void + define (SolverChoice& sc) override + { + docloud = sc.do_cloud; + doprecip = sc.do_precip; + m_fac_cond = lcond / sc.c_p; + m_fac_fus = lfus / sc.c_p; + m_fac_sub = lsub / sc.c_p; + m_gOcp = CONST_GRAV / sc.c_p; + m_axis = sc.ave_plane; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) override; + + // update ERF variables + void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist) override; + + // wrapper to do all the updating + void + Advance ( ) override + { + this->AdvanceKessler(); + this->Diagnose(); + } + + int + Qmoist_Size ( ) override { return Kessler::m_qmoist_size; } private: - // geometry - amrex::Geometry m_geom; - // valid boxes on which to evolve the solution - amrex::BoxArray m_gtoe; - - // timestep - amrex::Real dt; - - // number of vertical levels - int nlev, zlo, zhi; - - // plane average axis - int m_axis; - - // model options - bool docloud, doprecip; - - // constants - amrex::Real m_fac_cond; - amrex::Real m_fac_fus; - amrex::Real m_fac_sub; - amrex::Real m_gOcp; - - // microphysics parameters/coefficients - amrex::TableData accrrc; - amrex::TableData accrsi; - amrex::TableData accrsc; - amrex::TableData coefice; - amrex::TableData evaps1; - amrex::TableData evaps2; - amrex::TableData accrgi; - amrex::TableData accrgc; - amrex::TableData evapg1; - amrex::TableData evapg2; - amrex::TableData evapr1; - amrex::TableData evapr2; - - // vertical plane average data - amrex::TableData rho1d; - amrex::TableData pres1d; - amrex::TableData tabs1d; - amrex::TableData qt1d; - amrex::TableData qv1d; - amrex::TableData qn1d; - - // independent variables - amrex::Array mic_fab_vars; - - amrex::TableData gamaz; - amrex::TableData zmid; // mid value of vertical coordinate in physical domain - - // data (output) - amrex::TableData qifall; - amrex::TableData tlatqi; + // Number of qmoist variables + int m_qmoist_size = 6; + + // geometry + amrex::Geometry m_geom; + // valid boxes on which to evolve the solution + amrex::BoxArray m_gtoe; + + // timestep + amrex::Real dt; + + // number of vertical levels + int nlev, zlo, zhi; + + // plane average axis + int m_axis; + + // model options + bool docloud, doprecip; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // microphysics parameters/coefficients + amrex::TableData accrrc; + amrex::TableData accrsi; + amrex::TableData accrsc; + amrex::TableData coefice; + amrex::TableData evaps1; + amrex::TableData evaps2; + amrex::TableData accrgi; + amrex::TableData accrgc; + amrex::TableData evapg1; + amrex::TableData evapg2; + amrex::TableData evapr1; + amrex::TableData evapr2; + + // vertical plane average data + amrex::TableData rho1d; + amrex::TableData pres1d; + amrex::TableData tabs1d; + amrex::TableData qt1d; + amrex::TableData qv1d; + amrex::TableData qn1d; + + // independent variables + amrex::Array mic_fab_vars; + + amrex::TableData gamaz; + amrex::TableData zmid; // mid value of vertical coordinate in physical domain + + // data (output) + amrex::TableData qifall; + amrex::TableData tlatqi; }; #endif diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 64fe35639..580799d6c 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -11,36 +11,36 @@ using namespace amrex; /** * Compute Precipitation-related Microphysics quantities. */ -void Kessler::AdvanceKessler() { +void Kessler::AdvanceKessler () { - auto qt = mic_fab_vars[MicVar_Kess::qt]; - auto qp = mic_fab_vars[MicVar_Kess::qp]; - auto qn = mic_fab_vars[MicVar_Kess::qn]; - auto tabs = mic_fab_vars[MicVar_Kess::tabs]; + auto qt = mic_fab_vars[MicVar_Kess::qt]; + auto qp = mic_fab_vars[MicVar_Kess::qp]; + auto qn = mic_fab_vars[MicVar_Kess::qn]; + auto tabs = mic_fab_vars[MicVar_Kess::tabs]; - auto qcl = mic_fab_vars[MicVar_Kess::qcl]; - auto theta = mic_fab_vars[MicVar_Kess::theta]; - auto qv = mic_fab_vars[MicVar_Kess::qv]; - auto rho = mic_fab_vars[MicVar_Kess::rho]; + auto qcl = mic_fab_vars[MicVar_Kess::qcl]; + auto theta = mic_fab_vars[MicVar_Kess::theta]; + auto qv = mic_fab_vars[MicVar_Kess::qv]; + auto rho = mic_fab_vars[MicVar_Kess::rho]; - auto dz = m_geom.CellSize(2); - auto domain = m_geom.Domain(); - int k_lo = domain.smallEnd(2); - int k_hi = domain.bigEnd(2); + auto dz = m_geom.CellSize(2); + auto domain = m_geom.Domain(); + int k_lo = domain.smallEnd(2); + int k_hi = domain.bigEnd(2); - MultiFab fz; - auto ba = tabs->boxArray(); - auto dm = tabs->DistributionMap(); - fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells + MultiFab fz; + auto ba = tabs->boxArray(); + auto dm = tabs->DistributionMap(); + fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells - for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto fz_array = fz.array(mfi); - const Box& tbz = mfi.tilebox(); - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto fz_array = fz.array(mfi); + const Box& tbz = mfi.tilebox(); + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { Real rho_avg, qp_avg; if (k==k_lo) { @@ -61,12 +61,12 @@ void Kessler::AdvanceKessler() { fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; /*if(k==0){ - fz_array(i,j,k) = 0; - }*/ + fz_array(i,j,k) = 0; + }*/ }); } - Real dtn = dt; + Real dtn = dt; /*for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -107,56 +107,57 @@ void Kessler::AdvanceKessler() { - // get the temperature, dentisy, theta, qt and qp from input - for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); - auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); - auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + // get the temperature, dentisy, theta, qt and qp from input + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); + auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto theta_array = theta->array(mfi); - auto qv_array = qv->array(mfi); - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_array = theta->array(mfi); + auto qv_array = qv->array(mfi); + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - const auto& box3d = mfi.tilebox(); + const auto& box3d = mfi.tilebox(); - auto fz_array = fz.array(mfi); + auto fz_array = fz.array(mfi); - // Expose for GPU - Real d_fac_cond = m_fac_cond; + // Expose for GPU + Real d_fac_cond = m_fac_cond; - ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { - qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); - qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); - qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); + qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); + qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); - if(qt_array(i,j,k) == 0.0){ - qv_array(i,j,k) = 0.0; - qn_array(i,j,k) = 0.0; - } + if(qt_array(i,j,k) == 0.0){ + qv_array(i,j,k) = 0.0; + qn_array(i,j,k) = 0.0; + } - //------- Autoconversion/accretion - Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; + //------- Autoconversion/accretion + Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; - Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; - erf_qsatw(tabs_array(i,j,k), pressure, qsat); + Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; + erf_qsatw(tabs_array(i,j,k), pressure, qsat); - // If there is precipitating water (i.e. rain), and the cell is not saturated - // then the rain water can evaporate leading to extraction of latent heat, hence - // reducing temperature and creating negative buoyancy + // If there is precipitating water (i.e. rain), and the cell is not saturated + // then the rain water can evaporate leading to extraction of latent heat, hence + // reducing temperature and creating negative buoyancy - dq_clwater_to_rain = 0.0; - dq_rain_to_vapor = 0.0; - dq_vapor_to_clwater = 0.0; - dq_clwater_to_vapor = 0.0; + dq_clwater_to_rain = 0.0; + dq_rain_to_vapor = 0.0; + dq_vapor_to_clwater = 0.0; + dq_clwater_to_vapor = 0.0; - //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); - Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); + //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); + Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature - if(qv_array(i,j,k) > qsat){ + if(qv_array(i,j,k) > qsat){ dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); } @@ -167,52 +168,52 @@ void Kessler::AdvanceKessler() { dq_clwater_to_vapor = std::min(qn_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); } - if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) { + if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) { Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046); - dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/ - (5.4e5 + 2.55e6/(pressure*qsat))*dtn; - // The negative sign is to make this variable (vapor formed from evaporation) - // a poistive quantity (as qv/qs < 1) - dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor}); + dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/ + (5.4e5 + 2.55e6/(pressure*qsat))*dtn; + // The negative sign is to make this variable (vapor formed from evaporation) + // a poistive quantity (as qv/qs < 1) + dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor}); - // Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature - } + // Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature + } // If there is cloud water present then do accretion and autoconversion to rain - if (qn_array(i,j,k) > 0.0) { - qcc = qn_array(i,j,k); + if (qn_array(i,j,k) > 0.0) { + qcc = qn_array(i,j,k); - autor = 0.0; - if (qcc > qcw0) { - autor = 0.001; - } + autor = 0.0; + if (qcc > qcw0) { + autor = 0.001; + } - accrr = 0.0; - accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875); - dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001)); + accrr = 0.0; + accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875); + dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001)); - // If the amount of change is more than the amount of qc present, then dq = qc - dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k)); - } + // If the amount of change is more than the amount of qc present, then dq = qc + dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k)); + } - if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; - if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; - Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; - if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; - //dq_sed = 0.0; + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; + if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; + Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; + if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; + //dq_sed = 0.0; - qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; - qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; - qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; + qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; + qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; + qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; - theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor); + theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor); - qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); - qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); - qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); + qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); + qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); - }); - } + }); + } } diff --git a/Source/Microphysics/Kessler/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp index 304c812a3..470b839e3 100644 --- a/Source/Microphysics/Kessler/Update_Kessler.cpp +++ b/Source/Microphysics/Kessler/Update_Kessler.cpp @@ -11,48 +11,49 @@ * @param[out] qmoist: qv, qc, qi, qr, qs, qg */ void Kessler::Update (amrex::MultiFab& cons, - amrex::MultiFab& qmoist) + amrex::MultiFab& qmoist) { - // copy multifab data to qc, qv, and qi - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qv], 0, 0, 1, mic_fab_vars[MicVar_Kess::qv]->nGrowVect()); // vapor - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qcl], 0, 1, 1, mic_fab_vars[MicVar_Kess::qcl]->nGrowVect()); // cloud water - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qci], 0, 2, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // cloud ice - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpl], 0, 3, 1, mic_fab_vars[MicVar_Kess::qpl]->nGrowVect()); // rain - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 4, 1, mic_fab_vars[MicVar_Kess::qpi]->nGrowVect()); // snow - - // Don't need to copy this since it is filled below - // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 5, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // graupel - - amrex::MultiFab qgraup_mf(qmoist, amrex::make_alias, 5, 1); - - // Get the temperature, density, theta, qt and qp from input - for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto states_arr = cons.array(mfi); - - auto rho_arr = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi); - auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi); - auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - - auto qgraup_arr= qgraup_mf.array(mfi); - - const auto& box3d = mfi.tilebox(); - - // get potential total density, temperature, qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); - states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); - - // Graupel == precip total - rain - snow (but must be >= 0) - qgraup_arr(i,j,k) = 0.0;// - }); - } - - // Fill interior ghost cells and periodic boundaries - cons.FillBoundary(m_geom.periodicity()); - qmoist.FillBoundary(m_geom.periodicity()); + // copy multifab data to qc, qv, and qi + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qv], 0, 0, 1, mic_fab_vars[MicVar_Kess::qv]->nGrowVect()); // vapor + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qcl], 0, 1, 1, mic_fab_vars[MicVar_Kess::qcl]->nGrowVect()); // cloud water + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qci], 0, 2, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // cloud ice + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpl], 0, 3, 1, mic_fab_vars[MicVar_Kess::qpl]->nGrowVect()); // rain + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 4, 1, mic_fab_vars[MicVar_Kess::qpi]->nGrowVect()); // snow + + // Don't need to copy this since it is filled below + // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 5, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // graupel + + amrex::MultiFab qgraup_mf(qmoist, amrex::make_alias, 5, 1); + + // Get the temperature, density, theta, qt and qp from input + for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + + auto qgraup_arr= qgraup_mf.array(mfi); + + const auto& box3d = mfi.tilebox(); + + // get potential total density, temperature, qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + + // Graupel == precip total - rain - snow (but must be >= 0) + qgraup_arr(i,j,k) = 0.0;// + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); + qmoist.FillBoundary(m_geom.periodicity()); } diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index c384ad8a0..770b6cc80 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -56,6 +56,9 @@ public: m_moist_model->Update(cons_in, qmoist); } + int + Get_Qmoist_Size ( ) { return m_moist_model->Qmoist_Size(); } + private: std::unique_ptr m_moist_model; diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H index 1484058fa..150f62118 100644 --- a/Source/Microphysics/Null/NullMoist.H +++ b/Source/Microphysics/Null/NullMoist.H @@ -12,7 +12,8 @@ class NullMoist { virtual ~NullMoist () = default; - virtual void + virtual + void define (SolverChoice& /*sc*/) { } virtual @@ -22,17 +23,24 @@ class NullMoist { const amrex::Geometry& /*geom*/, const amrex::Real& /*dt_advance*/) { } - virtual void + virtual + void Advance ( ) { } - virtual void + virtual + void Diagnose ( ) { } - virtual void + virtual + void Update (amrex::MultiFab& /*cons_in*/, amrex::MultiFab& /*qmoist*/) { } - private: + virtual + int + Qmoist_Size ( ) { return NullMoist::m_qmoist_size; } + private: + int m_qmoist_size = 1; }; #endif diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index e9d263fe2..ad6434ecf 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -51,127 +51,133 @@ namespace MicVar { // class SAM : public NullMoist { - using FabPtr = std::shared_ptr; - - public: - // constructor - SAM () {} - - // destructor - virtual ~SAM () = default; - - // cloud physics - void Cloud (); - - // ice physics - void IceFall (); - - // precip - void Precip (); - - // precip fall - void PrecipFall (int hydro_type); - - // micro interface for precip fall - void MicroPrecipFall (); - - // process microphysics - void Proc (); - - // diagnose - void Diagnose () override; - - // Set up for first time - void - define (SolverChoice& sc) override - { - docloud = sc.do_cloud; - doprecip = sc.do_precip; - m_fac_cond = lcond / sc.c_p; - m_fac_fus = lfus / sc.c_p; - m_fac_sub = lsub / sc.c_p; - m_gOcp = CONST_GRAV / sc.c_p; - m_axis = sc.ave_plane; - } - - // init - void - Init (const amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist, - const amrex::BoxArray& grids, - const amrex::Geometry& geom, - const amrex::Real& dt_advance) override; - - // update ERF variables - void - Update (amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist) override; - - // wrapper to do all the updating - void - Advance ( ) override - { - this->Cloud(); - this->Diagnose(); - this->IceFall(); - this->Precip(); - this->MicroPrecipFall(); - } + using FabPtr = std::shared_ptr; + +public: + // constructor + SAM () {} + + // destructor + virtual ~SAM () = default; + + // cloud physics + void Cloud (); + + // ice physics + void IceFall (); + + // precip + void Precip (); + + // precip fall + void PrecipFall (int hydro_type); + + // micro interface for precip fall + void MicroPrecipFall (); + + // process microphysics + void Proc (); + + // diagnose + void Diagnose () override; + + // Set up for first time + void + define (SolverChoice& sc) override + { + docloud = sc.do_cloud; + doprecip = sc.do_precip; + m_fac_cond = lcond / sc.c_p; + m_fac_fus = lfus / sc.c_p; + m_fac_sub = lsub / sc.c_p; + m_gOcp = CONST_GRAV / sc.c_p; + m_axis = sc.ave_plane; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) override; + + // update ERF variables + void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist) override; + + // wrapper to do all the updating + void + Advance ( ) override + { + this->Cloud(); + this->Diagnose(); + this->IceFall(); + this->Precip(); + this->MicroPrecipFall(); + } + + int + Qmoist_Size ( ) override { return SAM::m_qmoist_size; } private: - // geometry - amrex::Geometry m_geom; - // valid boxes on which to evolve the solution - amrex::BoxArray m_gtoe; - - // timestep - amrex::Real dt; - - // number of vertical levels - int nlev, zlo, zhi; - - // plane average axis - int m_axis; - - // model options - bool docloud, doprecip; - - // constants - amrex::Real m_fac_cond; - amrex::Real m_fac_fus; - amrex::Real m_fac_sub; - amrex::Real m_gOcp; - - // microphysics parameters/coefficients - amrex::TableData accrrc; - amrex::TableData accrsi; - amrex::TableData accrsc; - amrex::TableData coefice; - amrex::TableData evaps1; - amrex::TableData evaps2; - amrex::TableData accrgi; - amrex::TableData accrgc; - amrex::TableData evapg1; - amrex::TableData evapg2; - amrex::TableData evapr1; - amrex::TableData evapr2; - - // vertical plane average data - amrex::TableData rho1d; - amrex::TableData pres1d; - amrex::TableData tabs1d; - amrex::TableData qt1d; - amrex::TableData qv1d; - amrex::TableData qn1d; - - // independent variables - amrex::Array mic_fab_vars; - - amrex::TableData gamaz; - amrex::TableData zmid; // mid value of vertical coordinate in physical domain - - // data (output) - amrex::TableData qifall; - amrex::TableData tlatqi; + // Number of qmoist variables + int m_qmoist_size = 6; + + // geometry + amrex::Geometry m_geom; + // valid boxes on which to evolve the solution + amrex::BoxArray m_gtoe; + + // timestep + amrex::Real dt; + + // number of vertical levels + int nlev, zlo, zhi; + + // plane average axis + int m_axis; + + // model options + bool docloud, doprecip; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // microphysics parameters/coefficients + amrex::TableData accrrc; + amrex::TableData accrsi; + amrex::TableData accrsc; + amrex::TableData coefice; + amrex::TableData evaps1; + amrex::TableData evaps2; + amrex::TableData accrgi; + amrex::TableData accrgc; + amrex::TableData evapg1; + amrex::TableData evapg2; + amrex::TableData evapr1; + amrex::TableData evapr2; + + // vertical plane average data + amrex::TableData rho1d; + amrex::TableData pres1d; + amrex::TableData tabs1d; + amrex::TableData qt1d; + amrex::TableData qv1d; + amrex::TableData qn1d; + + // independent variables + amrex::Array mic_fab_vars; + + amrex::TableData gamaz; + amrex::TableData zmid; // mid value of vertical coordinate in physical domain + + // data (output) + amrex::TableData qifall; + amrex::TableData tlatqi; }; #endif diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 7fec7dfc4..eba329f78 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -55,7 +55,7 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ FillPatch(lev, time, {&vars_old[lev][Vars::cons], &vars_old[lev][Vars::xvel], &vars_old[lev][Vars::yvel], &vars_old[lev][Vars::zvel]}); #if defined(ERF_USE_MOISTURE) - FillPatchMoistVars(lev, qmoist[lev]); + FillPatchMoistVars(lev, *(qmoist[lev])); #endif MultiFab* S_crse; diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index c85c1a5f0..0a1b81832 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -68,9 +68,11 @@ void ERF::advance_dycore(int level, MultiFab pi_hse(base_state[level], make_alias, 2, 1); // pi_0 is second component #if defined(ERF_USE_MOISTURE) - MultiFab qvapor (qmoist[level], make_alias, 0, 1); - MultiFab qcloud (qmoist[level], make_alias, 1, 1); - MultiFab qice (qmoist[level], make_alias, 2, 1); + // TODO: Protect from grabbing non-existing data + int q_size = micro.Get_Qmoist_Size(); + MultiFab qvapor (*(qmoist[level]), make_alias, 0, 1); + MultiFab qcloud (*(qmoist[level]), make_alias, 1, 1); + MultiFab qice (*(qmoist[level]), make_alias, 2, 1); #endif // These pointers are used in the MRI utility functions diff --git a/Source/TimeIntegration/ERF_advance_microphysics.cpp b/Source/TimeIntegration/ERF_advance_microphysics.cpp index 9c0e330c2..4ab7fde22 100644 --- a/Source/TimeIntegration/ERF_advance_microphysics.cpp +++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp @@ -7,11 +7,11 @@ void ERF::advance_microphysics (int lev, MultiFab& cons, const Real& dt_advance) { - micro.Init(cons, qmoist[lev], + micro.Init(cons, *(qmoist[lev]), grids[lev], Geom(lev), dt_advance); micro.Advance(); - micro.Update(cons, qmoist[lev]); + micro.Update(cons, *(qmoist[lev])); } #endif diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 3e28792fe..bbb08ea83 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -91,14 +91,14 @@ make_buoyancy(S_data, S_prim, buoyancy, #if defined(ERF_USE_MOISTURE) - qmoist[level], qv_d, qc_d, qi_d, + *(qmoist[level]), qv_d, qc_d, qi_d, #endif fine_geom, solverChoice, r0_new); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, #if defined(ERF_USE_MOISTURE) - qmoist[level], + *(qmoist[level]), #endif z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, @@ -188,14 +188,14 @@ // If not moving_terrain make_buoyancy(S_data, S_prim, buoyancy, #if defined(ERF_USE_MOISTURE) - qmoist[level], qv_d, qc_d, qi_d, + *(qmoist[level]), qv_d, qc_d, qi_d, #endif fine_geom, solverChoice, r0); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, #if defined(ERF_USE_MOISTURE) - qmoist[level], + *(qmoist[level]), #endif z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, @@ -362,7 +362,7 @@ // If not moving_terrain make_buoyancy(S_data, S_prim, buoyancy, #if defined(ERF_USE_MOISTURE) - qmoist[level], qv_d, qc_d, qi_d, + *(qmoist[level]), qv_d, qc_d, qi_d, #endif fine_geom, solverChoice, r0); @@ -370,7 +370,7 @@ S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, #if defined(ERF_USE_MOISTURE) - qmoist[level], + *(qmoist[level]), #endif z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs,