From 252d7e2a6a93c7c9cc089f39d7edf22f4fb7e380 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 19 Dec 2023 14:31:54 -0800 Subject: [PATCH] Fix index error in last commit and fix plotfile. --- Source/IO/Plotfile.cpp | 54 ++-- Source/Microphysics/FastEddy/FastEddy.H | 8 +- Source/Microphysics/FastEddy/FastEddy.cpp | 2 +- Source/Microphysics/FastEddy/Init_FE.cpp | 98 +------ .../Microphysics/Kessler/Diagnose_Kessler.cpp | 2 +- Source/Microphysics/Kessler/Init_Kessler.cpp | 257 +----------------- Source/Microphysics/Kessler/Kessler.H | 39 +-- Source/Microphysics/SAM/Init_SAM.cpp | 243 +---------------- Source/Microphysics/SAM/SAM.H | 9 +- 9 files changed, 47 insertions(+), 665 deletions(-) diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 2efe8ad1e..8fde84b7f 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -605,78 +605,66 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } - // TODO: The introduction of Q3 changes these plots - // Should we make the micro model own qt and qp - // since it knows whether qt = qv + qc or qt = qv + qc + qi? - //----------------------------------------------------------- // NOTE: Protect against accessing non-existent data if (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)) + if (containerHasElement(plot_var_names, "qt") && (q_size >= 1)) { - 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); + MultiFab qt_mf(*(qmoist[lev][0]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qt_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qp") && (q_size >= 6)) + if (containerHasElement(plot_var_names, "qv") && (q_size >= 2)) { - 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); + MultiFab qv_mf(*(qmoist[lev][1]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qv")) + if (containerHasElement(plot_var_names, "qc") && (q_size >= 3)) { - MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); + MultiFab qc_mf(*(qmoist[lev][2]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qc") && (q_size >= 2)) + if (containerHasElement(plot_var_names, "qi") && (q_size >= 4)) { - MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1); - MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); + MultiFab qi_mf(*(qmoist[lev][3]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qi") && (q_size >= 3)) + if (containerHasElement(plot_var_names, "qp") && (q_size >= 5)) { - MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1); - MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); + MultiFab qp_mf(*(qmoist[lev][4]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qp_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qrain") && (q_size >= 4)) + if (containerHasElement(plot_var_names, "qrain") && (q_size >= 6)) { - MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1); + MultiFab qr_mf(*(qmoist[lev][5]), 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)) + if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 7)) { - MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1); + MultiFab qs_mf(*(qmoist[lev][6]), 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)) + if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 8)) { - MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1); + MultiFab qg_mf(*(qmoist[lev][7]), make_alias, 0, 1); MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); mf_comp += 1; } } - //----------------------------------------------------------- #ifdef ERF_USE_PARTICLES if (containerHasElement(plot_var_names, "tracer_particle_count")) diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H index 6dc7d6613..9b8907f23 100644 --- a/Source/Microphysics/FastEddy/FastEddy.H +++ b/Source/Microphysics/FastEddy/FastEddy.H @@ -22,9 +22,11 @@ namespace MicVar_FE { // independent variables qv = 0, qc, + qt, + rho, // density theta, // liquid/ice water potential temperature tabs, // temperature - rho, // density + pres, // pressure NumVars }; } @@ -114,8 +116,8 @@ public: Qstate_Size () { return FastEddy::m_qstate_size; } private: - // Number of qmoist variables - int m_qmoist_size = 2; + // Number of qmoist variables (qt, qv, qc) + int m_qmoist_size = 3; // Number of qstate variables int m_qstate_size = 2; diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp index fc1aef273..27210ba84 100644 --- a/Source/Microphysics/FastEddy/FastEddy.cpp +++ b/Source/Microphysics/FastEddy/FastEddy.cpp @@ -21,7 +21,7 @@ void FastEddy::AdvanceFE () 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); + auto pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi); const auto& box3d = mfi.tilebox(); diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp index f16a4c65b..7ff34df97 100644 --- a/Source/Microphysics/FastEddy/Init_FE.cpp +++ b/Source/Microphysics/FastEddy/Init_FE.cpp @@ -1,4 +1,3 @@ - #include #include "Microphysics.H" #include "IndexDefines.H" @@ -29,7 +28,7 @@ void FastEddy::Init (const MultiFab& cons_in, m_gtoe = grids; MicVarMap.resize(m_qmoist_size); - MicVarMap = {MicVar_FE::qv, MicVar_FE::qc}; + MicVarMap = {MicVar_FE::qt, MicVar_FE::qv, MicVar_FE::qc}; // initialize microphysics variables for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) { @@ -65,13 +64,14 @@ void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in) auto states_array = cons_in.array(mfi); - auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); - auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); + auto qt_array = mic_fab_vars[MicVar_FE::qt]->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_Kess::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); - auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi); - auto pres_array = mic_fab_vars[MicVar::pres]->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 tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi); // Get pressure, theta, temperature, density, and qt, qp amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -80,6 +80,7 @@ void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in) 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); + qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k); tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp), states_array(i,j,k,RhoTheta_comp), @@ -88,84 +89,3 @@ void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in) }); } } - - - -#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. - * - * @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, MultiFab& qmoist, - const BoxArray& grids, - const Geometry& geom, - const Real& dt_advance) -{ - m_geom = geom; - m_gtoe = grids; - - 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.); - } - - // 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; - } - - // 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); - - 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)); - }); - } -} -#endif diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp index e9908effd..e63a4a84d 100644 --- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -35,7 +35,7 @@ void Kessler::Diagnose () 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;; + 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 40767e7f5..6f626550f 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -1,4 +1,3 @@ - #include #include "Microphysics.H" #include "IndexDefines.H" @@ -30,8 +29,8 @@ void Kessler::Init (const MultiFab& cons_in, 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}; + MicVarMap = {MicVar_Kess::qt, MicVar_Kess::qv, MicVar_Kess::qcl, MicVar_Kess::qci, + MicVar_Kess::qp, MicVar_Kess::qpl, MicVar_Kess::qpi}; // initialize microphysics variables for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) { @@ -95,255 +94,3 @@ void Kessler::Copy_State_to_Micro (const MultiFab& cons_in) } } -#if 0 -/** - * 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, - 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.; - }); - } - - - // 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(); - -#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; - }); -#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); - }); - -} -#endif diff --git a/Source/Microphysics/Kessler/Kessler.H b/Source/Microphysics/Kessler/Kessler.H index 16f5cc5f1..26c7c147d 100644 --- a/Source/Microphysics/Kessler/Kessler.H +++ b/Source/Microphysics/Kessler/Kessler.H @@ -134,8 +134,8 @@ public: Qstate_Size () { return Kessler::m_qstate_size; } private: - // Number of qmoist variables - int m_qmoist_size = 6; + // Number of qmoist variables (qt, qv, qcl, qci, qp, qpl, qpi) + int m_qmoist_size = 7; // Number of qstate variables int m_qstate_size = 3; @@ -169,40 +169,5 @@ private: // independent variables amrex::Array mic_fab_vars; - -#if 0 - // 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 - }; #endif diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp index 00a847e73..7db0099a9 100644 --- a/Source/Microphysics/SAM/Init_SAM.cpp +++ b/Source/Microphysics/SAM/Init_SAM.cpp @@ -1,4 +1,3 @@ - #include #include "Microphysics.H" #include "IndexDefines.H" @@ -29,8 +28,8 @@ void SAM::Init (const MultiFab& cons_in, m_gtoe = grids; MicVarMap.resize(m_qmoist_size); - MicVarMap = {MicVar::qv, MicVar::qcl, MicVar::qci, - MicVar::qr, MicVar::qpi, MicVar::qg}; + MicVarMap = {MicVar::qt, MicVar::qv, MicVar::qcl, MicVar::qci, + MicVar::qp, MicVar::qpl, MicVar::qpi, MicVar::qg}; // initialize microphysics variables for (auto ivar = 0; ivar < MicVar::NumVars; ++ivar) { @@ -251,241 +250,3 @@ void SAM::Compute_Coefficients () (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); }); } - -#if 0 -/** - * 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, 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::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 qc_array = qc.array(mfi); - auto qi_array = qi.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 rho_array = mic_fab_vars[MicVar::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar::tabs]->array(mfi); - auto pres_array = mic_fab_vars[MicVar::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); - temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp)); - 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); - }); - - // 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; - }); -#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); - }); -} -#endif diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 6590ff7b3..41a416262 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -33,7 +33,6 @@ namespace MicVar { rho, // density pres, // pressure // derived variables - qr, // rain water qv, // water vapor qn, // cloud condensate (liquid+ice), initial to zero qci, // cloud ice @@ -142,17 +141,17 @@ public: } void - Compute_Coefficients ( ); + Compute_Coefficients (); int - Qmoist_Size ( ) override { return SAM::m_qmoist_size; } + Qmoist_Size () override { return SAM::m_qmoist_size; } int Qstate_Size () { return SAM::m_qstate_size; } private: - // Number of qmoist variables - int m_qmoist_size = 6; + // Number of qmoist variables (qt, qv, qcl, qci, qp, qpl, qpi, qpg) + int m_qmoist_size = 8; // Number of qmoist variables int m_qstate_size = 3;