diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index dc32173d9..68270c55a 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -15,7 +15,6 @@ PhysBCFunctNoOp null_bc; * @param[in] time time at which the data should be filled * @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data */ - void ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) { @@ -100,6 +99,10 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) // We call this even if init_type == real because this routine will fill the vertical bcs (*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time); + + // Update vars in the micro model + if (solverChoice.moisture_type != MoistureType::None) + micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]); } /* @@ -109,7 +112,6 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) * @param[in] time time at which the data should be filled * @param[out] mf MultiFab to be filled (qmoist[lev]) */ - void ERF::FillPatchMoistVars (int lev, MultiFab& mf) { @@ -151,7 +153,6 @@ ERF::FillPatchMoistVars (int lev, MultiFab& mf) * @param[in] eddyDiffs diffusion coefficients for LES turbulence models * @param[in] allow_most_bcs if true then use MOST bcs at the low boundary */ - void ERF::FillIntermediatePatch (int lev, Real time, const Vector& mfs, @@ -260,9 +261,11 @@ ERF::FillIntermediatePatch (int lev, Real time, // MOST boundary conditions if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs) m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get()); -} -// + // Update vars in the micro model + if (solverChoice.moisture_type != MoistureType::None) + micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]); +} /* * Fill valid and ghost data. @@ -274,7 +277,6 @@ ERF::FillIntermediatePatch (int lev, Real time, * @param[in] time time at which the data should be filled * @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data */ - void ERF::FillCoarsePatch (int lev, Real time, const Vector& mfs) { diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 01276a7b8..a3b31199a 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -87,7 +87,11 @@ struct SolverChoice { // Set a different default for moist vs dry if (moisture_type != MoistureType::None) { - buoyancy_type = 2; // uses Tprime + if (moisture_type == MoistureType::FastEddy) { + buoyancy_type = 1; // asserted in make buoyancy + } else { + buoyancy_type = 2; // uses Tprime + } } // Is the terrain static or moving? diff --git a/Source/ERF.H b/Source/ERF.H index 122cb0116..c6ff0ebe7 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -525,7 +525,7 @@ private: amrex::Vector rW_new; Microphysics micro; - amrex::Vector> qmoist; // This has up to 6 components: qv, qc, qi, qr, qs, qg + amrex::Vector> qmoist; // (lev,ncomp) 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 8df292a95..53a1f65c1 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -109,6 +109,12 @@ ERF::ERF () amrex::Print() << "\n"; } + int nlevs_max = max_level + 1; + + // NOTE: size micro before readparams (chooses the model at all levels) + micro.ReSize(nlevs_max); + qmoist.resize(nlevs_max); + ReadParameters(); const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1); const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2); @@ -128,8 +134,6 @@ ERF::ERF () // No valid BoxArray and DistributionMapping have been defined. // But the arrays for them have been resized. - int nlevs_max = max_level + 1; - istep.resize(nlevs_max, 0); nsubsteps.resize(nlevs_max, 1); for (int lev = 1; lev <= max_level; ++lev) { @@ -157,8 +161,6 @@ ERF::ERF () vars_old[lev].resize(Vars::NumTypes); } - qmoist.resize(nlevs_max); - mri_integrator_mem.resize(nlevs_max); physbcs.resize(nlevs_max); @@ -217,8 +219,7 @@ ERF::ERF () } } -ERF::~ERF () -= default; +ERF::~ERF () = default; // advance solution to final time void @@ -529,7 +530,6 @@ ERF::InitData () } } - if (solverChoice.coupling_type == CouplingType::TwoWay) { int src_comp_reflux = 0; int num_comp_reflux = vars_new[0][Vars::cons].nComp(); @@ -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][0])); // qv component } } } @@ -578,23 +578,6 @@ ERF::InitData () } } - // Initialize microphysics here - micro.define(solverChoice); - - // Call Init which will call Diagnose to fill qmoist - if (solverChoice.moisture_type != MoistureType::None) - { - for (int lev = 0; lev <= finest_level; ++lev) - { - // 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]), - grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose - micro.Update(vars_new[lev][Vars::cons], *(qmoist[lev])); - } - } - } - // Configure ABLMost params if used MostWall boundary condition // NOTE: we must set up the MOST routine before calling WritePlotFile because // WritePlotFile calls FillPatch in order to compute gradients @@ -817,6 +800,7 @@ ERF::InitData () void ERF::restart () { + // TODO: This could be deleted since ba/dm are not created yet? for (int lev = 0; lev <= finest_level; ++lev) { auto& lev_new = vars_new[lev]; @@ -859,8 +843,6 @@ 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 @@ -1177,13 +1159,11 @@ ERF::MakeHorizontalAverages () if (use_moisture) { - MultiFab qv(*(qmoist[lev]), make_alias, 0, 1); - for (MFIter mfi(mf); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); auto fab_arr = mf.array(mfi); auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi); - auto const qv_arr = qv.const_array(mfi); + auto const qv_arr = qmoist[lev][0]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real dens = cons_arr(i, j, k, Rho_comp); @@ -1464,6 +1444,12 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, amrex::Print() << "\n"; } + int nlevs_max = max_level + 1; + + // NOTE: size micro before readparams (chooses the model at all levels) + micro.ReSize(nlevs_max); + qmoist.resize(nlevs_max); + ReadParameters(); const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1); const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2); @@ -1475,8 +1461,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, // No valid BoxArray and DistributionMapping have been defined. // But the arrays for them have been resized. - int nlevs_max = max_level + 1; - istep.resize(nlevs_max, 0); nsubsteps.resize(nlevs_max, 1); for (int lev = 1; lev <= max_level; ++lev) { @@ -1504,8 +1488,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, rV_old.resize(nlevs_max); rW_old.resize(nlevs_max); - qmoist.resize(nlevs_max); - mri_integrator_mem.resize(nlevs_max); physbcs.resize(nlevs_max); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index cbd27688e..7d9c77e7a 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -62,7 +62,15 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, // Microphysics // ******************************************************************************************* int q_size = micro.Get_Qmoist_Size(); - qmoist[lev] = std::make_unique(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg + qmoist[lev].resize(q_size); + micro.Define(lev, solverChoice); + if (solverChoice.moisture_type != MoistureType::None) + { + micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value + } + for (int mvar(0); mvar(ba, dm, q_size, ngrow_state); + qmoist[lev].resize(q_size); + micro.Define(lev, solverChoice); + if (solverChoice.moisture_type != MoistureType::None) + { + micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value + } + for (int mvar(0); mvardefine(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg + qmoist[lev].resize(q_size); + micro.Define(lev, solverChoice); + if (solverChoice.moisture_type != MoistureType::None) + { + micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value + } + for (int mvar(0); mvarnGrowVect(); - 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")); + for (int mvar(0); mvarnGrowVect(); + int nvar = qmoist[lev][mvar]->nComp(); + MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng); + MultiFab::Copy(moist_vars,*(qmoist[lev][mvar]),0,0,nvar,ng); + VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars")); + } } // Note that we write the ghost cells of the base state (unlike above) @@ -300,7 +303,6 @@ ERF::ReadCheckpointFile () } for (int lev = 0; lev <= finest_level; ++lev) { - // read in level 'lev' BoxArray from Header BoxArray ba; ba.readFrom(is); @@ -338,13 +340,14 @@ ERF::ReadCheckpointFile () MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0); IntVect ng; - if (solverChoice.moisture_type == MoistureType::None) { - qmoist[lev]->setVal(0.0); - } else { - 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); + if (solverChoice.moisture_type != MoistureType::None) { + for (int mvar(0); mvarnGrowVect(); + int nvar = qmoist[lev][mvar]->nComp(); + MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng); + VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars")); + MultiFab::Copy(*(qmoist[lev][mvar]),moist_vars,0,0,nvar,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 3b5c6ca24..bd06f78f2 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -246,18 +246,17 @@ ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVector& fab_arr = mf_out.array(mfi); + const Array4& fab_arr = mf_out.array(mfi); const Array4& cons_arr = mf_cons.array(mfi); const Array4& u_cc_arr = u_cc.array(mfi); const Array4& v_cc_arr = v_cc.array(mfi); const Array4& w_cc_arr = w_cc.array(mfi); const Array4& p0_arr = p_hse.array(mfi); - const Array4& qv_arr = qv.array(mfi); + const Array4& qv_arr = qmoist[0][0]->array(mfi); // TODO: Is this written only on lev 0? ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 3b1347b8d..b237edd0e 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -119,6 +119,14 @@ ERF::WritePlotFile (int which, Vector plot_var_names) fillset); } + // Get qmoist pointers if using moisture + bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None); + for (int lev = 0; lev <= finest_level; ++lev) { + for (int mvar(0); mvar mf(finest_level+1); for (int lev = 0; lev <= finest_level; ++lev) { mf[lev].define(grids[lev], dmap[lev], ncomp_mf, 0); @@ -211,11 +219,12 @@ ERF::WritePlotFile (int which, Vector 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[lev][0]) ? qmoist[lev][0]->const_array(mfi) : + Array4 {}; ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real qv_for_p = qv_arr(i,j,k); + Real qv_for_p = (qv_arr) ? qv_arr(i,j,k) : 0; const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p); }); @@ -233,11 +242,12 @@ 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[lev][0]) ? qmoist[lev][0]->const_array(mfi) : + Array4 {}; ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real qv_for_p = qv_arr(i,j,k); + Real qv_for_p = (qv_arr) ? qv_arr(i,j,k) : 0; const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k); }); @@ -600,69 +610,71 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } // NOTE: Protect against accessing non-existent data - int q_size = qmoist[lev]->nComp(); - MultiFab qv_mf(*(qmoist[lev]), make_alias, 0, 1); + if (l_use_moisture) { + int q_size = qmoist[lev].size(); + MultiFab qv_mf(*(qmoist[lev][0]), make_alias, 0, 1); - if (containerHasElement(plot_var_names, "qt") && (q_size >= 3)) - { - MultiFab qc_mf(*(qmoist[lev]), make_alias, 1, 1); - MultiFab qi_mf(*(qmoist[lev]), make_alias, 2, 1); - MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qc_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qi_mf,0,mf_comp,1,0); - mf_comp += 1; - } + if (containerHasElement(plot_var_names, "qt") && (q_size >= 3)) + { + MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1); + MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qc_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qi_mf,0,mf_comp,1,0); + mf_comp += 1; + } - if (containerHasElement(plot_var_names, "qp") && (q_size >= 6)) - { - 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); - MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qs_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qg_mf,0,mf_comp,1,0); - mf_comp += 1; - } + if (containerHasElement(plot_var_names, "qp") && (q_size >= 6)) + { + MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1); + MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1); + MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qs_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qg_mf,0,mf_comp,1,0); + mf_comp += 1; + } - if (containerHasElement(plot_var_names, "qv")) - { - MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); - mf_comp += 1; - } + if (containerHasElement(plot_var_names, "qv")) + { + MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); + mf_comp += 1; + } - if (containerHasElement(plot_var_names, "qc") && (q_size >= 2)) - { - MultiFab qc_mf(*(qmoist[lev]), make_alias, 1, 1); - MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); - mf_comp += 1; - } + if (containerHasElement(plot_var_names, "qc") && (q_size >= 2)) + { + MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); + mf_comp += 1; + } - if (containerHasElement(plot_var_names, "qi") && (q_size >= 3)) - { - MultiFab qi_mf(*(qmoist[lev]), make_alias, 2, 1); - MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); - mf_comp += 1; - } + if (containerHasElement(plot_var_names, "qi") && (q_size >= 3)) + { + MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); + mf_comp += 1; + } - if (containerHasElement(plot_var_names, "qrain") && (q_size >= 4)) - { - MultiFab qr_mf(*(qmoist[lev]), make_alias, 3, 1); - MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); - mf_comp += 1; - } + if (containerHasElement(plot_var_names, "qrain") && (q_size >= 4)) + { + MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); + mf_comp += 1; + } - if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 5)) - { - MultiFab qs_mf(*(qmoist[lev]), make_alias, 4, 1); - MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0); - mf_comp += 1; - } + if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 5)) + { + MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0); + mf_comp += 1; + } - if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 6)) - { - MultiFab qg_mf(*(qmoist[lev]), make_alias, 5, 1); - MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); - mf_comp += 1; + if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 6)) + { + MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); + mf_comp += 1; + } } #ifdef ERF_USE_PARTICLES diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 2e2fb0586..94a832c7f 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -24,8 +24,7 @@ using namespace amrex; void ERF::init_custom (int lev) { - auto& lev_new = vars_new[lev]; - auto& qmoist_new = *(qmoist[lev]); + auto& lev_new = vars_new[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,12 +41,12 @@ 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()); - qmoist_pert.setVal(0.); - - MultiFab qv_pert(qmoist_pert, amrex::make_alias, 0, 1); - MultiFab qc_pert(qmoist_pert, amrex::make_alias, 1, 1); - MultiFab qi_pert(qmoist_pert, amrex::make_alias, 2, 1); + MultiFab qv_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow()); + MultiFab qc_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow()); + MultiFab qi_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow()); + qv_pert.setVal(0.); + qc_pert.setVal(0.); + qi_pert.setVal(0.); int fix_random_seed = 0; ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); @@ -108,11 +107,11 @@ ERF::init_custom (int lev) if (solverChoice.moisture_type != MoistureType::None) { MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ1_comp, RhoQ1_comp, 1, cons_pert.nGrow()); MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ2_comp, RhoQ2_comp, 1, cons_pert.nGrow()); - MultiFab::Add( qmoist_new, qmoist_pert, 0, 0, 3, qmoist_pert.nGrow()); // qv, qc, qi -#if defined(ERF_USE_WARM_NO_PRECIP) - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQv_comp, RhoQv_comp, 1, cons_pert.nGrow()); - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQc_comp, RhoQc_comp, 1, cons_pert.nGrow()); -#endif + + MultiFab::Add(*(qmoist[lev][0]) , qv_pert , 0, 0, 1, qv_pert.nGrow()); + MultiFab::Add(*(qmoist[lev][1]) , qc_pert , 0, 0, 1, qc_pert.nGrow()); + if (qmoist[lev].size() > 2) + MultiFab::Add(*(qmoist[lev][2]) , qi_pert , 0, 0, 1, qi_pert.nGrow()); } MultiFab::Add(lev_new[Vars::xvel], xvel_pert, 0, 0, 1, xvel_pert.nGrowVect()); diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H index e35354b62..11a59a1d7 100644 --- a/Source/Microphysics/FastEddy/FastEddy.H +++ b/Source/Microphysics/FastEddy/FastEddy.H @@ -51,7 +51,7 @@ public: // Set up for first time void - define (SolverChoice& sc) override + Define (SolverChoice& sc) override { docloud = sc.do_cloud; doprecip = sc.do_precip; @@ -65,31 +65,63 @@ public: // 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; + // Copy state into micro vars + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + + // Copy micro into state vars + void + Copy_Micro_to_State (amrex::MultiFab& cons_in) override; + // update ERF variables void Update (amrex::MultiFab& cons_in, amrex::MultiFab& qmoist) override; + void + Update_Micro_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_State_to_Micro(cons_in); + this->Diagnose(); + } + + void + Update_State_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_Micro_to_State(cons_in); + } + // wrapper to do all the updating void - Advance ( ) override + Advance (const amrex::Real& dt_advance) override { + dt = dt_advance; + this->AdvanceFE(); this->Diagnose(); } + amrex::MultiFab* + Qmoist_Ptr (const int& varIdx) override + { + AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size); + return mic_fab_vars[MicVarMap[varIdx]].get(); + } + int - Qmoist_Size ( ) override { return FastEddy::m_qmoist_size; } + Qmoist_Size () override { return FastEddy::m_qmoist_size; } private: // Number of qmoist variables int m_qmoist_size = 2; + // MicVar map (Qmoist indices -> MicVar enum) + amrex::Vector MicVarMap; + // geometry amrex::Geometry m_geom; diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp index a07fdee14..0eb69f9b1 100644 --- a/Source/Microphysics/FastEddy/FastEddy.cpp +++ b/Source/Microphysics/FastEddy/FastEddy.cpp @@ -46,8 +46,8 @@ void FastEddy::AdvanceFE () { 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){ diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp index 164016fc9..1e27d4e49 100644 --- a/Source/Microphysics/FastEddy/Init_FE.cpp +++ b/Source/Microphysics/FastEddy/Init_FE.cpp @@ -8,6 +8,93 @@ using namespace amrex; +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @param[in] geom Geometry associated with these MultiFabs and grids + * @param[in] dt_advance Timestep for the advance + */ +void FastEddy::Init (const MultiFab& cons_in, + const BoxArray& grids, + const Geometry& geom, + const Real& dt_advance) +{ + dt = dt_advance; + m_geom = geom; + m_gtoe = grids; + + MicVarMap.resize(m_qmoist_size); + MicVarMap = {MicVar_FE::qv, MicVar_FE::qc}; + + // 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.); + } + + // Set class data members + 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; + } +} + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + */ +void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in) +{ + // Get the temperature, density, theta, qt and qp from input + for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); + + 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); + + // 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)); + }); + } +} + + + +#if 0 +#include +#include "Microphysics.H" +#include "IndexDefines.H" +#include "PlaneAverage.H" +#include "EOS.H" +#include "TileNoZ.H" + +using namespace amrex; + /** * Initializes the Microphysics module. * @@ -75,3 +162,4 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist, }); } } +#endif diff --git a/Source/Microphysics/FastEddy/Update_FE.cpp b/Source/Microphysics/FastEddy/Update_FE.cpp index eaedc3697..d6155ece7 100644 --- a/Source/Microphysics/FastEddy/Update_FE.cpp +++ b/Source/Microphysics/FastEddy/Update_FE.cpp @@ -38,4 +38,38 @@ void FastEddy::Update (amrex::MultiFab& cons, qmoist.FillBoundary(m_geom.periodicity()); } +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc, qi, qr, qs, qg + */ +void FastEddy::Copy_Micro_to_State (amrex::MultiFab& cons) +{ + + // Get the temperature, density, theta, qt and qp from input + for (amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); + + 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); + + // 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()); +} + diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp index c231053dd..55b613ca3 100644 --- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -1,11 +1,12 @@ #include "Microphysics.H" +#include "EOS.H" /** * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid * from the existing Microphysics variables. */ -void Kessler::Diagnose () { - +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]; @@ -40,3 +41,34 @@ void Kessler::Diagnose () { }); } } + + +/** + * Computes diagnostic quantities like temperature and pressure + * from the existing Microphysics variables. + */ +void Kessler::Compute_Tabs_and_Pres () +{ + auto rho = mic_fab_vars[MicVar_Kess::rho]; + auto theta = mic_fab_vars[MicVar_Kess::theta]; + auto qv = mic_fab_vars[MicVar_Kess::qv]; + auto tabs = mic_fab_vars[MicVar_Kess::tabs]; + auto pres = mic_fab_vars[MicVar_Kess::pres]; + + for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); + + auto rho_array = rho->array(mfi); + auto theta_array = theta->array(mfi); + auto qv_array = qv->array(mfi); + auto tabs_array = tabs->array(mfi); + auto pres_array = pres->array(mfi); + + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + amrex::Real RhoTheta = rho_array(i,j,k)*theta_array(i,j,k); + tabs_array(i,j,k) = getTgivenRandRTh(rho_array(i,j,k), RhoTheta, qv_array(i,j,k)); + pres_array(i,j,k) = getPgivenRTh(RhoTheta, qv_array(i,j,k))/100.; + }); + } +} diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp index 3c96dac8b..56f8c6be4 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -8,6 +8,88 @@ using namespace amrex; + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @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, + const BoxArray& grids, + const Geometry& geom, + const Real& dt_advance) +{ + dt = dt_advance; + m_geom = geom; + m_gtoe = grids; + + MicVarMap.resize(m_qmoist_size); + MicVarMap = {MicVar_Kess::qv, MicVar_Kess::qcl, MicVar_Kess::qci, + MicVar_Kess::qr, MicVar_Kess::qpi, MicVar_Kess::qg}; + + // 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.); + } + + // Set class data members + 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; + } +} + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + */ +void Kessler::Copy_State_to_Micro (const MultiFab& cons_in) +{ + // Get the temperature, density, theta, qt and qp from input + for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); + + auto states_array = cons_in.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 qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); + auto qi_array = mic_fab_vars[MicVar_Kess::qci]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + + // 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.; + }); + } +} + +#if 0 /** * Initializes the Microphysics module. * @@ -55,6 +137,7 @@ void Kessler::Init (const MultiFab& cons_in, zlo = lo.z; zhi = hi.z; + // parameters accrrc.resize({zlo}, {zhi}); accrsi.resize({zlo}, {zhi}); @@ -79,8 +162,10 @@ void Kessler::Init (const MultiFab& cons_in, // output qifall.resize({zlo}, {zhi}); tlatqi.resize({zlo}, {zhi}); + } + auto accrrc_t = accrrc.table(); auto accrsi_t = accrsi.table(); auto accrsc_t = accrsc.table(); @@ -112,6 +197,7 @@ void Kessler::Init (const MultiFab& cons_in, 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); @@ -148,6 +234,7 @@ void Kessler::Init (const MultiFab& cons_in, }); } + // calculate the plane average variables PlaneAverage cons_ave(&cons_in, m_geom, m_axis); cons_ave.compute_averages(ZDir(), cons_ave.field()); @@ -251,4 +338,6 @@ void Kessler::Init (const MultiFab& cons_in, pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) / (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); }); + } +#endif diff --git a/Source/Microphysics/Kessler/Kessler.H b/Source/Microphysics/Kessler/Kessler.H index f771575e0..00ade80ef 100644 --- a/Source/Microphysics/Kessler/Kessler.H +++ b/Source/Microphysics/Kessler/Kessler.H @@ -40,6 +40,7 @@ namespace MicVar_Kess { qcl, // cloud water qpl, // precip rain qpi, // precip ice + qg, // graupel // temporary variable omega, NumVars @@ -68,7 +69,7 @@ public: // Set up for first time void - define (SolverChoice& sc) override + Define (SolverChoice& sc) override { docloud = sc.do_cloud; doprecip = sc.do_precip; @@ -82,33 +83,73 @@ public: // 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; + // Copy state into micro vars + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + + // Copy state into micro vars + void + Copy_Micro_to_State (amrex::MultiFab& cons_in) override; + // update ERF variables void Update (amrex::MultiFab& cons_in, amrex::MultiFab& qmoist) override; - // wrapper to do all the updating + // update micro vars + void + Update_Micro_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_State_to_Micro(cons_in); + this->Diagnose(); + this->Compute_Tabs_and_Pres(); + } + + // update state vars + void + Update_State_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_Micro_to_State(cons_in); + } + + // wrapper to advance micro vars void - Advance ( ) override + Advance (const amrex::Real& dt_advance) override { + dt = dt_advance; + this->AdvanceKessler(); this->Diagnose(); + this->Compute_Tabs_and_Pres(); + } + + amrex::MultiFab* + Qmoist_Ptr (const int& varIdx) override + { + AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size); + return mic_fab_vars[MicVarMap[varIdx]].get(); } + void + Compute_Tabs_and_Pres (); + int - Qmoist_Size ( ) override { return Kessler::m_qmoist_size; } + Qmoist_Size () override { return Kessler::m_qmoist_size; } - private: +private: // Number of qmoist variables int m_qmoist_size = 6; + // MicVar map (Qmoist indices -> MicVar enum) + amrex::Vector MicVarMap; + // geometry amrex::Geometry m_geom; + // valid boxes on which to evolve the solution amrex::BoxArray m_gtoe; @@ -130,6 +171,10 @@ public: amrex::Real m_fac_sub; amrex::Real m_gOcp; + // independent variables + amrex::Array mic_fab_vars; + +#if 0 // microphysics parameters/coefficients amrex::TableData accrrc; amrex::TableData accrsi; @@ -161,5 +206,7 @@ public: // data (output) amrex::TableData qifall; amrex::TableData tlatqi; +#endif + }; #endif diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 580799d6c..cb5ca7dfe 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -11,15 +11,16 @@ 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 pres = mic_fab_vars[MicVar_Kess::pres]; auto qcl = mic_fab_vars[MicVar_Kess::qcl]; - auto theta = mic_fab_vars[MicVar_Kess::theta]; + auto theta = mic_fab_vars[MicVar_Kess::theta]; auto qv = mic_fab_vars[MicVar_Kess::qv]; auto rho = mic_fab_vars[MicVar_Kess::rho]; @@ -34,89 +35,51 @@ void Kessler::AdvanceKessler () { 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 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) { - rho_avg = rho_array(i,j,k); - qp_avg = qp_array(i,j,k); - } else if (k==k_hi+1) { - rho_avg = rho_array(i,j,k-1); - qp_avg = qp_array(i,j,k-1); - } else { - rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3 - qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k)); - } - - qp_avg = std::max(0.0, qp_avg); - - Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s - - fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; - - /*if(k==0){ - fz_array(i,j,k) = 0; - }*/ - }); - } - - Real dtn = dt; - - /*for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++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 qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); - - const auto& box3d = mfi.tilebox(); - - 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)); - - if(qt_array(i,j,k) == 0.0){ - qv_array(i,j,k) = 0.0; - qn_array(i,j,k) = 0.0; - } - }); - } + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real rho_avg, qp_avg; + + if (k==k_lo) { + rho_avg = rho_array(i,j,k); + qp_avg = qp_array(i,j,k); + } else if (k==k_hi+1) { + rho_avg = rho_array(i,j,k-1); + qp_avg = qp_array(i,j,k-1); + } else { + rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3 + qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k)); + } + qp_avg = std::max(0.0, qp_avg); - for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s - const auto& box3d = mfi.tilebox(); - auto fz_array = fz.array(mfi); + fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; - ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; - qp_array(i,j,k) = qp_array(i,j,k) + dq_sed; + /*if(k==0){ + fz_array(i,j,k) = 0; + }*/ }); - }*/ - - + } + Real dtn = dt; // 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 pres_array = mic_fab_vars[MicVar_Kess::pres]->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 rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); const auto& box3d = mfi.tilebox(); @@ -140,7 +103,7 @@ void Kessler::AdvanceKessler () { //------- 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; + Real pressure = pres_array(i,j,k); //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 diff --git a/Source/Microphysics/Kessler/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp index 470b839e3..cc085f9e0 100644 --- a/Source/Microphysics/Kessler/Update_Kessler.cpp +++ b/Source/Microphysics/Kessler/Update_Kessler.cpp @@ -57,3 +57,37 @@ void Kessler::Update (amrex::MultiFab& cons, } +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc, qi, qr, qs, qg + */ +void Kessler::Copy_Micro_to_State (amrex::MultiFab& cons) +{ + // Get the temperature, density, theta, qt and qp from input + for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); + + 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); + + // 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); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); +} + diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index 770b6cc80..4113fe8db 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -12,55 +12,76 @@ public: Microphysics () { } - ~Microphysics () = default; + ~Microphysics () = default; + + void + ReSize (const int& nlev) { m_moist_model.resize(nlev); } template void SetModel () { - m_moist_model = std::make_unique(); + for (int lev(0); lev(); + } } void - define (SolverChoice& sc) + Define (const int& lev, + SolverChoice& sc) { - m_moist_model->define(sc); + m_moist_model[lev]->Define(sc); } void - Init (const amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist, + Init (const int& lev, + const amrex::MultiFab& cons_in, const amrex::BoxArray& grids, const amrex::Geometry& geom, const amrex::Real& dt_advance) { - m_moist_model->Init(cons_in, qmoist, grids, geom, dt_advance); + m_moist_model[lev]->Init(cons_in, grids, geom, dt_advance); } void - Advance ( ) + Advance (const int& lev, const amrex::Real& dt_advance) { - m_moist_model->Advance(); + m_moist_model[lev]->Advance(dt_advance); } void - Diagnose ( ) + Diagnose (const int& lev) { - m_moist_model->Diagnose(); + m_moist_model[lev]->Diagnose(); } void - Update (amrex::MultiFab& cons_in, + Update (const int& lev, + amrex::MultiFab& cons_in, amrex::MultiFab& qmoist) { - m_moist_model->Update(cons_in, qmoist); + m_moist_model[lev]->Update(cons_in, qmoist); } - int - Get_Qmoist_Size ( ) { return m_moist_model->Qmoist_Size(); } + void + Update_Micro_Vars_Lev (const int& lev, amrex::MultiFab& cons_in) + { + m_moist_model[lev]->Update_Micro_Vars(cons_in); + } + void + Update_State_Vars_Lev (const int& lev, amrex::MultiFab& cons_in) + { + m_moist_model[lev]->Update_State_Vars(cons_in); + } + + amrex::MultiFab* + Get_Qmoist_Ptr (const int& lev, const int& varIdx) { return m_moist_model[lev]->Qmoist_Ptr(varIdx); } + + int + Get_Qmoist_Size () { return m_moist_model[0]->Qmoist_Size(); } private: - std::unique_ptr m_moist_model; + amrex::Vector> m_moist_model; }; #endif diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H index 150f62118..5dfb0836e 100644 --- a/Source/Microphysics/Null/NullMoist.H +++ b/Source/Microphysics/Null/NullMoist.H @@ -14,31 +14,50 @@ class NullMoist { virtual void - define (SolverChoice& /*sc*/) { } + Define (SolverChoice& /*sc*/) { } virtual void Init (const amrex::MultiFab& /*cons_in*/, - amrex::MultiFab& /*qmoist*/, const amrex::BoxArray& /*grids*/, const amrex::Geometry& /*geom*/, const amrex::Real& /*dt_advance*/) { } virtual void - Advance ( ) { } + Advance (const amrex::Real& /*dt_advance*/) { } virtual void - Diagnose ( ) { } + Update (amrex::MultiFab& /*cons_in*/, + amrex::MultiFab& /*qmoist*/) { } virtual void - Update (amrex::MultiFab& /*cons_in*/, - amrex::MultiFab& /*qmoist*/) { } + Update_Micro_Vars (amrex::MultiFab& /*cons_in*/) { } + + virtual + void + Update_State_Vars (amrex::MultiFab& /*cons_in*/) { } + + virtual + void + Diagnose () { } + + virtual + void + Copy_State_to_Micro (const amrex::MultiFab& /*cons_in*/) { } + + virtual + void + Copy_Micro_to_State (amrex::MultiFab& /*cons_in*/) { } + + virtual + amrex::MultiFab* + Qmoist_Ptr (const int& /*varIdx*/ ) { return nullptr; } virtual int - Qmoist_Size ( ) { return NullMoist::m_qmoist_size; } + Qmoist_Size () { return NullMoist::m_qmoist_size; } private: int m_qmoist_size = 1; diff --git a/Source/Microphysics/SAM/Diagnose_SAM.cpp b/Source/Microphysics/SAM/Diagnose_SAM.cpp index 784f21a8e..544a202f6 100644 --- a/Source/Microphysics/SAM/Diagnose_SAM.cpp +++ b/Source/Microphysics/SAM/Diagnose_SAM.cpp @@ -1,4 +1,5 @@ #include "Microphysics.H" +#include "EOS.H" /** * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid @@ -14,6 +15,7 @@ void SAM::Diagnose () { auto qci = mic_fab_vars[MicVar::qci]; auto qpl = mic_fab_vars[MicVar::qpl]; auto qpi = mic_fab_vars[MicVar::qpi]; + auto qg = mic_fab_vars[MicVar::qg]; auto tabs = mic_fab_vars[MicVar::tabs]; for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -25,35 +27,55 @@ void SAM::Diagnose () { auto qci_array = qci->array(mfi); auto qpl_array = qpl->array(mfi); auto qpi_array = qpi->array(mfi); + auto qg_array = qg->array(mfi); auto tabs_array = tabs->array(mfi); const auto& box3d = mfi.tilebox(); - amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - qv_array(i,j,k) = qt_array(i,j,k) - qn_array(i,j,k); - amrex::Real omn = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); - 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 = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); - 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) = qt_array(i,j,k) - qn_array(i,j,k); + amrex::Real omn = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); + 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 = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + qpl_array(i,j,k) = qp_array(i,j,k)*omp; + qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp); + + // TODO: If omp above is bound by [0, 1], shouldn't this always be 0? + // Graupel == precip total - rain - snow (but must be >= 0) + qg_array(i,j,k) = std::max(0.0, qp_array(i,j,k)-qpl_array(i,j,k)-qpi_array(i,j,k)); }); } } + /** - * Wrapper for the PrecipFall, Cloud, Precipitation, and Diagnostics routines. + * Computes diagnostic quantities like temperature and pressure + * from the existing Microphysics variables. */ -void SAM::Proc () { +void SAM::Compute_Tabs_and_Pres () +{ + auto rho = mic_fab_vars[MicVar::rho]; + auto theta = mic_fab_vars[MicVar::theta]; + auto qv = mic_fab_vars[MicVar::qv]; + auto tabs = mic_fab_vars[MicVar::tabs]; + auto pres = mic_fab_vars[MicVar::pres]; - MicroPrecipFall(); + for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); - if (docloud) { - Cloud(); - if (doprecip) { - Precip(); - } - } + auto rho_array = rho->array(mfi); + auto theta_array = theta->array(mfi); + auto qv_array = qv->array(mfi); + auto tabs_array = tabs->array(mfi); + auto pres_array = pres->array(mfi); - Diagnose(); + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + amrex::Real RhoTheta = rho_array(i,j,k)*theta_array(i,j,k); + tabs_array(i,j,k) = getTgivenRandRTh(rho_array(i,j,k), RhoTheta, qv_array(i,j,k)); + pres_array(i,j,k) = getPgivenRTh(RhoTheta, qv_array(i,j,k))/100.; + }); + } } diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp index ad5f0c484..19970f4f3 100644 --- a/Source/Microphysics/SAM/Init_SAM.cpp +++ b/Source/Microphysics/SAM/Init_SAM.cpp @@ -8,6 +8,240 @@ using namespace amrex; +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @param[in] geom Geometry associated with these MultiFabs and grids + * @param[in] dt_advance Timestep for the advance + */ +void SAM::Init (const MultiFab& cons_in, + const BoxArray& grids, + const Geometry& geom, + const Real& dt_advance) +{ + dt = dt_advance; + m_geom = geom; + m_gtoe = grids; + + MicVarMap.resize(m_qmoist_size); + MicVarMap = {MicVar::qv, MicVar::qcl, MicVar::qci, + MicVar::qr, MicVar::qpi, MicVar::qg}; + + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar::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.); + } + + // Set class data members + 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}); + } +} + + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + */ +void SAM::Copy_State_to_Micro (const MultiFab& cons_in) +{ + // Get the temperature, density, theta, qt and qp from input + for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); + + auto states_array = cons_in.array(mfi); + + auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi); + auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi); + auto qc_array = mic_fab_vars[MicVar::qcl]->array(mfi); + auto qi_array = mic_fab_vars[MicVar::qci]->array(mfi); + auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi); + + // 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.; + }); + } +} + +void SAM::Compute_Coefficients () +{ + auto dz = m_geom.CellSize(2); + auto lowz = m_geom.ProbLo(2); + + 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 ); + + // calculate the plane average variables + PlaneAverage rho_ave(mic_fab_vars[MicVar::rho].get(), m_geom, m_axis); + PlaneAverage theta_ave(mic_fab_vars[MicVar::theta].get(), m_geom, m_axis); + rho_ave.compute_averages(ZDir(), rho_ave.field()); + theta_ave.compute_averages(ZDir(), theta_ave.field()); + + // get host variable rho, and rhotheta + int ncell = rho_ave.ncell_line(); + + Gpu::HostVector rho_h(ncell), theta_h(ncell); + rho_ave.line_average(0, rho_h); + theta_ave.line_average(0, theta_h); + + // copy data to device + Gpu::DeviceVector rho_d(ncell), theta_d(ncell); + Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); + Gpu::streamSynchronize(); + + Real* rho_dptr = rho_d.data(); + Real* theta_dptr = theta_d.data(); + + Real gOcp = m_gOcp; + + // TODO: None of this is qv aware...? + amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept + { + Real RhoTheta = rho_dptr[k]*theta_dptr[k]; + Real pressure = getPgivenRTh(RhoTheta); + rho1d_t(k) = rho_dptr[k]; + pres1d_t(k) = pressure/100.; + tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],RhoTheta); + zmid_t(k) = lowz + (k+0.5)*dz; + gamaz_t(k) = gOcp*zmid_t(k); + }); + + if(round(gam3) != 2) { + std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; + std::exit(-1); + } + + // Populate all the coefficients + 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 0 /** * Initializes the Microphysics module. * @@ -243,3 +477,4 @@ void SAM::Init (const MultiFab& cons_in, MultiFab& qmoist, (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); }); } +#endif diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index ad6434ecf..4cccf0995 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -40,6 +40,7 @@ namespace MicVar { qcl, // cloud water qpl, // precip rain qpi, // precip ice + qg, // graupel // temporary variable omega, NumVars @@ -75,15 +76,12 @@ public: // 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 + Define (SolverChoice& sc) override { docloud = sc.do_cloud; doprecip = sc.do_precip; @@ -97,20 +95,44 @@ public: // 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; + // Copy state into micro vars + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + + // Copy state into micro vars + void + Copy_Micro_to_State (amrex::MultiFab& cons_in) override; + // update ERF variables void Update (amrex::MultiFab& cons_in, amrex::MultiFab& qmoist) override; + void + Update_Micro_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_State_to_Micro(cons_in); + this->Diagnose(); + this->Compute_Tabs_and_Pres(); + this->Compute_Coefficients(); + } + + void + Update_State_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_Micro_to_State(cons_in); + } + // wrapper to do all the updating void - Advance ( ) override + Advance (const amrex::Real& dt_advance) override { + dt = dt_advance; + this->Cloud(); this->Diagnose(); this->IceFall(); @@ -118,15 +140,35 @@ public: this->MicroPrecipFall(); } + amrex::MultiFab* + Qmoist_Ptr (const int& varIdx) override + { + AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size); + return mic_fab_vars[MicVarMap[varIdx]].get(); + } + + void + Compute_Tabs_and_Pres ( ); + + void + Compute_Coefficients ( ); + int Qmoist_Size ( ) override { return SAM::m_qmoist_size; } - private: + // independent variables + amrex::Array mic_fab_vars; + +private: // Number of qmoist variables int m_qmoist_size = 6; + // MicVar map (Qmoist indices -> MicVar enum) + amrex::Vector MicVarMap; + // geometry amrex::Geometry m_geom; + // valid boxes on which to evolve the solution amrex::BoxArray m_gtoe; @@ -170,9 +212,6 @@ public: 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 diff --git a/Source/Microphysics/SAM/Update_SAM.cpp b/Source/Microphysics/SAM/Update_SAM.cpp index 7e391d1dd..bd7531d7f 100644 --- a/Source/Microphysics/SAM/Update_SAM.cpp +++ b/Source/Microphysics/SAM/Update_SAM.cpp @@ -57,4 +57,38 @@ void SAM::Update (amrex::MultiFab& cons, qmoist.FillBoundary(m_geom.periodicity()); } +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc, qi, qr, qs, qg + */ +void SAM::Copy_Micro_to_State (amrex::MultiFab& cons) +{ + // Get the temperature, density, theta, qt and qp from input + for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); + + 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); + + // 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); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); +} + diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H index 0eb4e1beb..0c072ff5f 100644 --- a/Source/Radiation/Radiation.H +++ b/Source/Radiation/Radiation.H @@ -34,56 +34,57 @@ // Radiation code interface class class Radiation { public: - Radiation() { + Radiation () { // First, make sure yakl has been initialized if (!yakl::isInitialized()) yakl::init(); } - ~Radiation() = default; + ~Radiation () = default; // init - void initialize(const amrex::MultiFab& cons_in, amrex::MultiFab& qmoist, - const amrex::BoxArray& grids, - const amrex::Geometry& geom, - const amrex::Real& dt_advance, - const bool& do_sw_rad, - const bool& do_lw_rad, - const bool& do_aero_rad, - const bool& do_snow_opt, - const bool& is_cmip6_volcano); + void initialize (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance, + const bool& do_sw_rad, + const bool& do_lw_rad, + const bool& do_aero_rad, + const bool& do_snow_opt, + const bool& is_cmip6_volcano); // run radiation model - void run(); + void run (); // call back - void on_complete(); - - void radiation_driver_lw(int ncol, int nlev, - const real3d& gas_vmr, - const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint, - const real3d& cld_tau_gpt, const real3d& aer_tau_bnd, - FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, - const real2d& qrl, const real2d& qrlc); - - void radiation_driver_sw(int ncol, - const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid, - const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs, - const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt, - const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd, - FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, + void on_complete (); + + void radiation_driver_lw (int ncol, int nlev, + const real3d& gas_vmr, + const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint, + const real3d& cld_tau_gpt, const real3d& aer_tau_bnd, + FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, + const real2d& qrl, const real2d& qrlc); + + void radiation_driver_sw (int ncol, + const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid, + const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs, + const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt, + const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd, + FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, const real2d& qrs, const real2d& qrsc); - void set_daynight_indices(const real1d& coszrs, - const int1d& day_indices, - const int1d& night_indices); + void set_daynight_indices (const real1d& coszrs, + const int1d& day_indices, + const int1d& night_indices); - void get_gas_vmr(const std::vector& gas_names, - const real3d& gas_vmr); + void get_gas_vmr (const std::vector& gas_names, + const real3d& gas_vmr); - void calculate_heating_rate(const real2d& flux_up, - const real2d& flux_dn, - const real2d& pint, - const real2d& heating_rate); + void calculate_heating_rate (const real2d& flux_up, + const real2d& flux_dn, + const real2d& pint, + const real2d& heating_rate); private: // geometry amrex::Geometry m_geom; diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index 8767ea9a2..56bee93c3 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -74,7 +74,8 @@ namespace internal { } // init -void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, +void Radiation::initialize(const MultiFab& cons_in, + MultiFab& qmoist, const BoxArray& grids, const Geometry& geom, const Real& dt_advance, diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 362464d16..b7f987cbb 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -56,7 +56,8 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ &vars_old[lev][Vars::yvel], &vars_old[lev][Vars::zvel]}); if (solverChoice.moisture_type != MoistureType::None) { - FillPatchMoistVars(lev, *qmoist[lev]); + // TODO: This is only qv + FillPatchMoistVars(lev, *(qmoist[lev][0])); } MultiFab* S_crse; diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 68327b12d..f6525cc5c 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -256,38 +256,40 @@ void ERF::advance_dycore(int level, Gpu::DeviceVector qv_d(ncell), qc_d(ncell), qi_d(ncell); - if (q_size >= 1) { - MultiFab qvapor (*(qmoist[level]), make_alias, 0, 1); - PlaneAverage qv_ave(&qvapor, geom[level], solverChoice.ave_plane); - qv_ave.compute_averages(ZDir(), qv_ave.field()); + if (solverChoice.moisture_type != MoistureType::None) { + if (q_size >= 1) { + MultiFab qvapor (*(qmoist[level][0]), make_alias, 0, 1); + PlaneAverage qv_ave(&qvapor, geom[level], solverChoice.ave_plane); + qv_ave.compute_averages(ZDir(), qv_ave.field()); - Gpu::HostVector qv_h(ncell); + Gpu::HostVector qv_h(ncell); - qv_ave.line_average(0, qv_h); - Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); - } + qv_ave.line_average(0, qv_h); + Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); + } - if (q_size >= 2) { - MultiFab qcloud (*(qmoist[level]), make_alias, 1, 1); - PlaneAverage qc_ave(&qcloud, geom[level], solverChoice.ave_plane); - qc_ave.compute_averages(ZDir(), qc_ave.field()); + if (q_size >= 2) { + MultiFab qcloud (*(qmoist[level][1]), make_alias, 0, 1); + PlaneAverage qc_ave(&qcloud, geom[level], solverChoice.ave_plane); + qc_ave.compute_averages(ZDir(), qc_ave.field()); - Gpu::HostVector qc_h(ncell); + Gpu::HostVector qc_h(ncell); - qc_ave.line_average(0, qc_h); - Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin()); - } + qc_ave.line_average(0, qc_h); + Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin()); + } - if (q_size >= 3) { - MultiFab qice (*(qmoist[level]), make_alias, 2, 1); - PlaneAverage qi_ave(&qice , geom[level], solverChoice.ave_plane); - qi_ave.compute_averages(ZDir(), qi_ave.field()); + if (q_size >= 3) { + MultiFab qice (*(qmoist[level][2]), make_alias, 0, 1); + PlaneAverage qi_ave(&qice , geom[level], solverChoice.ave_plane); + qi_ave.compute_averages(ZDir(), qi_ave.field()); - Gpu::HostVector qi_h(ncell); + Gpu::HostVector qi_h(ncell); - qi_ave.line_average(0, qi_h); - Gpu::copyAsync(Gpu::hostToDevice, qi_h.begin(), qi_h.end(), qi_d.begin()); - Gpu::streamSynchronize(); + qi_ave.line_average(0, qi_h); + Gpu::copyAsync(Gpu::hostToDevice, qi_h.begin(), qi_h.end(), qi_d.begin()); + Gpu::streamSynchronize(); + } } #include "TI_no_substep_fun.H" diff --git a/Source/TimeIntegration/ERF_advance_microphysics.cpp b/Source/TimeIntegration/ERF_advance_microphysics.cpp index 93373efcf..e726bab88 100644 --- a/Source/TimeIntegration/ERF_advance_microphysics.cpp +++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp @@ -7,11 +7,8 @@ void ERF::advance_microphysics (int lev, const Real& dt_advance) { if (solverChoice.moisture_type != MoistureType::None) { - micro.Init(cons, *(qmoist[lev]), - grids[lev], - Geom(lev), - dt_advance); - micro.Advance(); - micro.Update(cons, *(qmoist[lev])); + micro.Update_Micro_Vars_Lev(lev, cons); + micro.Advance(lev, dt_advance); + micro.Update_State_Vars_Lev(lev, cons); } } diff --git a/Source/TimeIntegration/ERF_advance_radiation.cpp b/Source/TimeIntegration/ERF_advance_radiation.cpp index c42222b60..826c0d934 100644 --- a/Source/TimeIntegration/ERF_advance_radiation.cpp +++ b/Source/TimeIntegration/ERF_advance_radiation.cpp @@ -14,7 +14,8 @@ void ERF::advance_radiation (int lev, bool do_snow_opt {true}; bool is_cmip6_volcano {true}; - rad.initialize(cons, qmoist[lev], + // TODO: Only passing qv from the qmoist vector!!! + rad.initialize(cons, *(qmoist[lev][0]), grids[lev], Geom(lev), dt_advance, diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index 22c540c4c..c6dea21cb 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -36,7 +36,7 @@ using namespace amrex; void make_buoyancy (Vector& S_data, const MultiFab& S_prim, MultiFab& buoyancy, - MultiFab* qmoist, + Vector qmoist, Gpu::DeviceVector qv_d, Gpu::DeviceVector qc_d, Gpu::DeviceVector qi_d, @@ -184,9 +184,9 @@ void make_buoyancy (Vector& S_data, // Base state density const Array4& r0_arr = r0->const_array(mfi); - const Array4 & qv_data = qmoist->array(mfi,0); - const Array4 & qc_data = qmoist->array(mfi,1); - const Array4 & qi_data = qmoist->array(mfi,2); + const Array4 & qv_data = qmoist[0]->array(mfi); + const Array4 & qc_data = qmoist[1]->array(mfi); + const Array4 & qi_data = qmoist[2]->array(mfi); amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -252,9 +252,9 @@ void make_buoyancy (Vector& S_data, const Array4 & cell_data = S_data[IntVar::cons].array(mfi); const Array4 & cell_prim = S_prim.array(mfi); - const Array4 & qv_data = qmoist->const_array(mfi,0); - const Array4 & qc_data = qmoist->const_array(mfi,1); - const Array4 & qi_data = qmoist->const_array(mfi,2); + const Array4 & qv_data = qmoist[0]->const_array(mfi); + const Array4 & qc_data = qmoist[1]->const_array(mfi); + const Array4 & qi_data = qmoist[2]->const_array(mfi); amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -318,9 +318,9 @@ void make_buoyancy (Vector& S_data, const Array4 & cell_data = S_data[IntVar::cons].array(mfi); const Array4 & cell_prim = S_prim.array(mfi); - const Array4 & qv_data = qmoist->const_array(mfi,0); - const Array4 & qc_data = qmoist->const_array(mfi,1); - const Array4 & qi_data = qmoist->const_array(mfi,2); + const Array4 & qv_data = qmoist[0]->const_array(mfi); + const Array4 & qc_data = qmoist[1]->const_array(mfi); + const Array4 & qi_data = qmoist[2]->const_array(mfi); amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 4b673ac08..0a0e042eb 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -80,7 +80,7 @@ void erf_slow_rhs_pre (int level, int finest_level, const MultiFab& xvel, const MultiFab& yvel, const MultiFab& zvel, - const MultiFab* qv, + Vector qmoist, std::unique_ptr& z_t_mf, MultiFab& Omega, const MultiFab& source, @@ -543,7 +543,7 @@ void erf_slow_rhs_pre (int level, int finest_level, FArrayBox pprime; pprime.resize(gbx,1); Elixir pp_eli = pprime.elixir(); const Array4 & pp_arr = pprime.array(); - const Array4 & qv_arr = (use_moisture) ? qv->const_array(mfi) : Array4 {}; + const Array4 & qv_arr = (use_moisture) ? qmoist[0]->const_array(mfi) : Array4 {}; { BL_PROFILE("slow_rhs_pre_pprime"); ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H index 01e72c9b4..6e72b11d4 100644 --- a/Source/TimeIntegration/TI_headers.H +++ b/Source/TimeIntegration/TI_headers.H @@ -22,7 +22,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, const amrex::MultiFab& zvel, - const amrex::MultiFab* qmoist, + amrex::Vector qmoist, std::unique_ptr& z_t, amrex::MultiFab& Omega, const amrex::MultiFab& source, @@ -218,7 +218,7 @@ void make_fast_coeffs (int level, void make_buoyancy (amrex::Vector< amrex::MultiFab>& S_data, const amrex::MultiFab & S_prim, amrex::MultiFab& buoyancy, - amrex::MultiFab* qmoist, + amrex::Vector qmoist, const amrex::Gpu::DeviceVector qv_d, const amrex::Gpu::DeviceVector qc_d, const amrex::Gpu::DeviceVector qi_d, diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 5ac193aa0..f783597e5 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -90,12 +90,12 @@ MultiFab* p0_new = &p_hse_new; make_buoyancy(S_data, S_prim, buoyancy, - qmoist[level].get(), qv_d, qc_d, qi_d, + qmoist[level], qv_d, qc_d, qi_d, 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, - qmoist[level].get(), + qmoist[level], z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss, @@ -183,12 +183,12 @@ // If not moving_terrain make_buoyancy(S_data, S_prim, buoyancy, - qmoist[level].get(), qv_d, qc_d, qi_d, + qmoist[level], qv_d, qc_d, qi_d, 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, - qmoist[level].get(), + qmoist[level], z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss,