From 61beb8cbceaf879fb1b2ea8049fb5d110db77072 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 10 Dec 2023 18:52:55 -0800 Subject: [PATCH 1/3] More steps towards making moisture a run-time setting --- CMake/BuildERFExe.cmake | 1 + Exec/Make.ERF | 5 + Source/DataStructs/DataStruct.H | 4 +- Source/ERF.cpp | 4 +- Source/IO/Checkpoint.cpp | 90 ++++--- Source/IO/ERF_Write1DProfiles.cpp | 56 ++-- Source/IO/Plotfile.cpp.convergence | 5 +- Source/IO/ReadFromWRFInput.cpp | 19 +- .../Initialization/ERF_init_from_wrfinput.cpp | 50 ++-- Source/Microphysics/FastEddy/Diagnose_FE.cpp | 41 +++ Source/Microphysics/FastEddy/FastEddy.H | 153 +++++++++++ Source/Microphysics/FastEddy/FastEddy.cpp | 214 +++++++++++++++ Source/Microphysics/FastEddy/Init_FE.cpp | 249 ++++++++++++++++++ Source/Microphysics/FastEddy/Make.package | 5 + Source/Microphysics/FastEddy/Update_FE.cpp | 58 ++++ Source/Microphysics/Kessler/Kessler.cpp | 2 +- Source/Microphysics/Microphysics.H | 1 + 17 files changed, 856 insertions(+), 101 deletions(-) create mode 100644 Source/Microphysics/FastEddy/Diagnose_FE.cpp create mode 100644 Source/Microphysics/FastEddy/FastEddy.H create mode 100644 Source/Microphysics/FastEddy/FastEddy.cpp create mode 100644 Source/Microphysics/FastEddy/Init_FE.cpp create mode 100644 Source/Microphysics/FastEddy/Make.package create mode 100644 Source/Microphysics/FastEddy/Update_FE.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 358149741..ad50a774d 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -203,6 +203,7 @@ function(build_erf_lib erf_lib_name) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/FastEddy) if(ERF_ENABLE_RRTMGP) target_link_libraries(${erf_lib_name} PUBLIC yakl) diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 320ae7927..787166d3c 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -112,6 +112,11 @@ include $(ERF_MOISTURE_NULL_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) INCLUDE_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) +ERF_MOISTURE_FE_DIR = $(ERF_SOURCE_DIR)/Microphysics/FastEddy +include $(ERF_MOISTURE_FE_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_FE_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_FE_DIR) + ERF_MOISTURE_SAM_DIR = $(ERF_SOURCE_DIR)/Microphysics/SAM include $(ERF_MOISTURE_SAM_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index d24d138b7..80540afbc 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -29,7 +29,7 @@ enum struct TerrainType { }; enum struct MoistureType { - Kessler, SAM, None + Kessler, SAM, FastEddy, None }; enum struct Coord { @@ -90,6 +90,8 @@ struct SolverChoice { moisture_type = MoistureType::SAM; } else if (moisture_model_string == "Kessler") { moisture_type = MoistureType::Kessler; + } else if (moisture_model_string == "FastEddy") { + moisture_type = MoistureType::FastEddy; } else { moisture_type = MoistureType::None; } diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 55f7229fb..1efc763a9 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -905,8 +905,6 @@ ERF::init_only (int lev, Real time) lev_new[Vars::xvel].OverrideSync(geom[lev].periodicity()); lev_new[Vars::yvel].OverrideSync(geom[lev].periodicity()); lev_new[Vars::zvel].OverrideSync(geom[lev].periodicity()); - - print_state(vars_new[0][Vars::cons], IntVect(30,30,2)); } // read in some parameters from inputs file @@ -965,6 +963,8 @@ ERF::ReadParameters () micro.SetModel(); } else if (solverChoice.moisture_type == MoistureType::Kessler) { micro.SetModel(); + } else if (solverChoice.moisture_type == MoistureType::FastEddy) { + micro.SetModel(); } else if (solverChoice.moisture_type == MoistureType::None) { micro.SetModel(); amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n"; diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index 8c75dd695..0fd1c7db5 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -106,34 +106,34 @@ ERF::WriteCheckpointFile () const } } - // write the MultiFab data to, e.g., chk00010/Level_0/ - // Here we make copies of the MultiFab with no ghost cells - for (int lev = 0; lev <= finest_level; ++lev) - { - MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); - MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,NVAR,0); - VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell")); - - MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); - MultiFab::Copy(xvel,vars_new[lev][Vars::xvel],0,0,1,0); - VisMF::Write(xvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "XFace")); - - MultiFab yvel(convert(grids[lev],IntVect(0,1,0)),dmap[lev],1,0); - MultiFab::Copy(yvel,vars_new[lev][Vars::yvel],0,0,1,0); - VisMF::Write(yvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "YFace")); - - MultiFab zvel(convert(grids[lev],IntVect(0,0,1)),dmap[lev],1,0); - MultiFab::Copy(zvel,vars_new[lev][Vars::zvel],0,0,1,0); - VisMF::Write(zvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace")); - - IntVect ng; -#ifdef ERF_USE_MOISTURE - // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars - ng = qmoist[lev].nGrowVect(); - MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng); - MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),ng); - VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars")); -#endif + // write the MultiFab data to, e.g., chk00010/Level_0/ + // Here we make copies of the MultiFab with no ghost cells + for (int lev = 0; lev <= finest_level; ++lev) + { + MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); + MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,NVAR,0); + VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell")); + + MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); + MultiFab::Copy(xvel,vars_new[lev][Vars::xvel],0,0,1,0); + VisMF::Write(xvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "XFace")); + + MultiFab yvel(convert(grids[lev],IntVect(0,1,0)),dmap[lev],1,0); + MultiFab::Copy(yvel,vars_new[lev][Vars::yvel],0,0,1,0); + VisMF::Write(yvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "YFace")); + + MultiFab zvel(convert(grids[lev],IntVect(0,0,1)),dmap[lev],1,0); + MultiFab::Copy(zvel,vars_new[lev][Vars::zvel],0,0,1,0); + VisMF::Write(zvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace")); + + IntVect ng; + if (solverChoice.moisture_type != MoistureType::None) { + // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars + ng = qmoist[lev].nGrowVect(); + MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng); + MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),ng); + VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars")); + } // Note that we write the ghost cells of the base state (unlike above) ng = base_state[lev].nGrowVect(); @@ -249,7 +249,14 @@ ERF::ReadCheckpointFile () // conservative, cell-centered vars is >> chk_ncomp; GotoNextLine(is); - AMREX_ASSERT(chk_ncomp == Cons::NumVars); +#if 0 + if (solverChoice.moisture_type != MoistureType::None) { + AMREX_ASSERT(chk_ncomp == Cons::NumVars); + } else { + // This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them + AMREX_ASSERT(chk_ncomp == Cons::NumVars-2); + } +#endif // x-velocity on faces is >> chk_ncomp; @@ -314,7 +321,16 @@ ERF::ReadCheckpointFile () { MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); VisMF::Read(cons, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Cell")); - MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars,0); + + if (solverChoice.moisture_type != MoistureType::None) { + MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars,0); + } else { + // This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them + MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars-2,0); + + // We zero them here so they won't cause problems while we're still debugging + vars_new[lev][Vars::cons].setVal(0.0, RhoQ1_comp, 2); + } MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); VisMF::Read(xvel, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "XFace")); @@ -329,12 +345,14 @@ ERF::ReadCheckpointFile () MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0); IntVect ng; -#ifdef ERF_USE_MOISTURE - 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); -#endif + 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); + } // Note that we read the ghost cells of the base state (unlike above) ng = base_state[lev].nGrowVect(); diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index 37c5c401d..cacc3d0b8 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -194,9 +194,7 @@ ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVector& h_avg_u , Gpu::HostVector& v_cc_arr = v_cc.array(mfi); const Array4& w_cc_arr = w_cc.array(mfi); const Array4& cons_arr = mf_cons.array(mfi); -#if defined(ERF_USE_MOISTURE) - const Array4& qv_arr = qv.array(mfi); -#endif const Array4& p0_arr = p_hse.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -237,18 +232,45 @@ ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVector& 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); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k)); + + p -= p0_arr(i,j,k); + fab_arr(i, j, k,17) = p; // p' + fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k); // p'u + fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k); // p'v + fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k); // p'w + }); + } // mfi + } // use_moisture h_avg_rho = sumToLine(mf_out, 0,1,domain,zdir); h_avg_th = sumToLine(mf_out, 1,1,domain,zdir); diff --git a/Source/IO/Plotfile.cpp.convergence b/Source/IO/Plotfile.cpp.convergence index 8396e370d..f5dc08436 100644 --- a/Source/IO/Plotfile.cpp.convergence +++ b/Source/IO/Plotfile.cpp.convergence @@ -544,7 +544,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } -#if defined(ERF_USE_MOISTURE) + if (solverChoice.moisture_type != MoistureType::None) { calculate_derived("qt", derived::erf_derQt); calculate_derived("qp", derived::erf_derQp); @@ -590,7 +590,8 @@ ERF::WritePlotFile (int which, Vector plot_var_names) MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); mf_comp += 1; } -#elif defined(ERF_USE_WARM_NO_PRECIP) + } +#if defined(ERF_USE_WARM_NO_PRECIP) calculate_derived("qv", derived::erf_derQv); calculate_derived("qc", derived::erf_derQc); #endif diff --git a/Source/IO/ReadFromWRFInput.cpp b/Source/IO/ReadFromWRFInput.cpp index f930af6bb..4af0cfd15 100644 --- a/Source/IO/ReadFromWRFInput.cpp +++ b/Source/IO/ReadFromWRFInput.cpp @@ -1,5 +1,6 @@ #include "NCWpsFile.H" #include "AMReX_FArrayBox.H" +#include "DataStruct.H" using namespace amrex; @@ -16,14 +17,12 @@ read_from_wrfinput (int lev, FArrayBox& NC_MSFM_fab, FArrayBox& NC_SST_fab, FArrayBox& NC_C1H_fab , FArrayBox& NC_C2H_fab, FArrayBox& NC_RDNW_fab, -#if defined(ERF_USE_MOISTURE) FArrayBox& NC_QVAPOR_fab, FArrayBox& NC_QCLOUD_fab, FArrayBox& NC_QRAIN_fab, -#elif defined(ERF_USE_WARM_NO_PRECIP) -#endif FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab, - FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab) + FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab, + MoistureType moisture_type) { amrex::Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; @@ -51,12 +50,12 @@ read_from_wrfinput (int lev, NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 15 NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 16 NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17 -#if defined(ERF_USE_MOISTURE) - NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 18 - NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 19 - NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20 -# elif defined(ERF_USE_WARM_NO_PRECIP) -#endif + + if (moisture_type != MoistureType::None) { + NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 18 + NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 19 + NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20 + } // Read the netcdf file and fill these FABs amrex::Print() << "Building initial FABS from file " << fname << std::endl; diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index b2af0162b..556ed7de8 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -7,6 +7,7 @@ #include #include #include +#include using namespace amrex; @@ -22,14 +23,12 @@ read_from_wrfinput (int lev, const Box& domain, const std::string& fname, FArrayBox& NC_MSFM_fab, FArrayBox& NC_SST_fab, FArrayBox& NC_C1H_fab , FArrayBox& NC_C2H_fab, FArrayBox& NC_RDNW_fab, -#if defined(ERF_USE_MOISTURE) FArrayBox& NC_QVAPOR_fab, FArrayBox& NC_QCLOUD_fab, FArrayBox& NC_QRAIN_fab, -#elif defined(ERF_USE_WARM_NO_PRECIP) -#endif FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab, - FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab); + FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab, + MoistureType moisture_type); Real read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, @@ -60,17 +59,15 @@ void init_state_from_wrfinput (int lev, FArrayBox& state_fab, FArrayBox& x_vel_fab, FArrayBox& y_vel_fab, FArrayBox& z_vel_fab, -#if defined(ERF_USE_MOISTURE) const Vector& NC_QVAPOR_fab, const Vector& NC_QCLOUD_fab, const Vector& NC_QRAIN_fab, -#elif defined(ERF_USE_WARM_NO_PRECIP) -#endif const Vector& NC_xvel_fab, const Vector& NC_yvel_fab, const Vector& NC_zvel_fab, const Vector& NC_rho_fab, - const Vector& NC_rhotheta_fab); + const Vector& NC_rhotheta_fab, + MoistureType moisture_type); void init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab, @@ -117,12 +114,9 @@ ERF::init_from_wrfinput (int lev) Vector NC_PHB_fab ; NC_PHB_fab.resize(num_boxes_at_level[lev]); Vector NC_ALB_fab ; NC_ALB_fab.resize(num_boxes_at_level[lev]); Vector NC_PB_fab ; NC_PB_fab.resize(num_boxes_at_level[lev]); -#if defined(ERF_USE_MOISTURE) Vector NC_QVAPOR_fab; NC_QVAPOR_fab.resize(num_boxes_at_level[lev]); Vector NC_QCLOUD_fab; NC_QCLOUD_fab.resize(num_boxes_at_level[lev]); Vector NC_QRAIN_fab ; NC_QRAIN_fab.resize(num_boxes_at_level[lev]); -#elif defined(ERF_USE_WARM_NO_PRECIP) -#endif // amrex::Print() << "Building initial FABS from file " << nc_init_file[lev][idx] << std::endl; if (nc_init_file.empty()) @@ -135,11 +129,9 @@ ERF::init_from_wrfinput (int lev) NC_rhop_fab[idx], NC_rhoth_fab[idx], NC_MUB_fab[idx], NC_MSFU_fab[idx], NC_MSFV_fab[idx], NC_MSFM_fab[idx], NC_SST_fab[idx], NC_C1H_fab[idx], NC_C2H_fab[idx], NC_RDNW_fab[idx], -#if defined(ERF_USE_MOISTURE) NC_QVAPOR_fab[idx], NC_QCLOUD_fab[idx], NC_QRAIN_fab[idx], -#elif defined(ERF_USE_WARM_NO_PRECIP) -#endif - NC_PH_fab[idx],NC_PHB_fab[idx],NC_ALB_fab[idx],NC_PB_fab[idx]); + NC_PH_fab[idx],NC_PHB_fab[idx],NC_ALB_fab[idx],NC_PB_fab[idx], + solverChoice.moisture_type); } auto& lev_new = vars_new[lev]; @@ -158,12 +150,9 @@ ERF::init_from_wrfinput (int lev) FArrayBox &zvel_fab = lev_new[Vars::zvel][mfi]; init_state_from_wrfinput(lev, cons_fab, xvel_fab, yvel_fab, zvel_fab, -#if defined(ERF_USE_MOISTURE) NC_QVAPOR_fab, NC_QCLOUD_fab, NC_QRAIN_fab, -#elif defined(ERF_USE_WARM_NO_PRECIP) -#endif NC_xvel_fab, NC_yvel_fab, NC_zvel_fab, - NC_rho_fab, NC_rhoth_fab); + NC_rho_fab, NC_rhoth_fab, solverChoice.moisture_type); } // mf #ifdef _OPENMP @@ -280,17 +269,15 @@ init_state_from_wrfinput (int lev, FArrayBox& x_vel_fab, FArrayBox& y_vel_fab, FArrayBox& z_vel_fab, -#if defined(ERF_USE_MOISTURE) const Vector& NC_QVAPOR_fab, const Vector& NC_QCLOUD_fab, const Vector& NC_QRAIN_fab, -#elif defined(ERF_USE_WARM_NO_PRECIP) -#endif const Vector& NC_xvel_fab, const Vector& NC_yvel_fab, const Vector& NC_zvel_fab, const Vector& NC_rho_fab, - const Vector& NC_rhotheta_fab) + const Vector& NC_rhotheta_fab, + MoistureType moisture_type) { int nboxes = NC_xvel_fab.size(); for (int idx = 0; idx < nboxes; idx++) @@ -313,19 +300,18 @@ init_state_from_wrfinput (int lev, // This copies the density state_fab.template copy(NC_rho_fab[idx], 0, Rho_comp, 1); - // This copies (rho*theta) state_fab.template copy(NC_rhotheta_fab[idx], 0, RhoTheta_comp, 1); -#if defined(ERF_USE_MOISTURE) - state_fab.template copy(NC_QVAPOR_fab[idx], 0, RhoQ1_comp, 1); - state_fab.template plus(NC_QCLOUD_fab[idx], 0, RhoQ1_comp, 1); - state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ1_comp, 1); + if (moisture_type != MoistureType::None) + { + state_fab.template copy(NC_QVAPOR_fab[idx], 0, RhoQ1_comp, 1); + state_fab.template plus(NC_QCLOUD_fab[idx], 0, RhoQ1_comp, 1); + state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ1_comp, 1); - state_fab.template copy(NC_QRAIN_fab[idx], 0, RhoQ2_comp, 1); - state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ2_comp, 1); -# elif defined(ERF_USE_WARM_NO_PRECIP) -#endif + state_fab.template copy(NC_QRAIN_fab[idx], 0, RhoQ2_comp, 1); + state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ2_comp, 1); + } } // idx } diff --git a/Source/Microphysics/FastEddy/Diagnose_FE.cpp b/Source/Microphysics/FastEddy/Diagnose_FE.cpp new file mode 100644 index 000000000..b3c1307c1 --- /dev/null +++ b/Source/Microphysics/FastEddy/Diagnose_FE.cpp @@ -0,0 +1,41 @@ +#include "Microphysics.H" + +/** + * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid + * from the existing Microphysics variables. + */ +void FastEddy::Diagnose () { + + auto qt = mic_fab_vars[MicVar_FE::qt]; + auto qp = mic_fab_vars[MicVar_FE::qp]; + auto qv = mic_fab_vars[MicVar_FE::qv]; + auto qn = mic_fab_vars[MicVar_FE::qn]; + auto qcl = mic_fab_vars[MicVar_FE::qcl]; + auto qci = mic_fab_vars[MicVar_FE::qci]; + auto qpl = mic_fab_vars[MicVar_FE::qpl]; + auto qpi = mic_fab_vars[MicVar_FE::qpi]; + auto tabs = mic_fab_vars[MicVar_FE::tabs]; + + for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qt_array = qt->array(mfi); + auto qp_array = qp->array(mfi); + auto qv_array = qv->array(mfi); + auto qn_array = qn->array(mfi); + auto qcl_array = qcl->array(mfi); + auto qci_array = qci->array(mfi); + auto qpl_array = qpl->array(mfi); + auto qpi_array = qpi->array(mfi); + + const auto& box3d = mfi.tilebox(); + + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + qv_array(i,j,k) = std::max(0.0,qt_array(i,j,k) - qn_array(i,j,k)); + amrex::Real omn = 1.0; + qcl_array(i,j,k) = qn_array(i,j,k)*omn; + qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn); + amrex::Real omp = 1.0;; + qpl_array(i,j,k) = qp_array(i,j,k)*omp; + qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp); + }); + } +} diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H new file mode 100644 index 000000000..df54a1736 --- /dev/null +++ b/Source/Microphysics/FastEddy/FastEddy.H @@ -0,0 +1,153 @@ +/* + */ +#ifndef FastEddy_H +#define FastEddy_H + +#include +#include +#include + +#include +#include +#include +#include + +#include "ERF_Constants.H" +#include "Microphysics_Utils.H" +#include "IndexDefines.H" +#include "DataStruct.H" + +namespace MicVar_FE { + enum { + // independent variables + qt = 0, + qp, + theta, // liquid/ice water potential temperature + tabs, // temperature + rho, // density + pres, // pressure + // derived variables + qr, // rain water + qv, // water vapor + qn, // cloud condensate (liquid+ice), initial to zero + qci, // cloud ice + qcl, // cloud water + qpl, // precip rain + qpi, // precip ice + // temporary variable + omega, + NumVars + }; +} + +// +// use MultiFab for 3D data, but table for 1D data +// +class FastEddy : public NullMoist { + + using FabPtr = std::shared_ptr; + + public: + // constructor + FastEddy () {} + + // destructor + virtual ~FastEddy () = default; + + // cloud physics + void AdvanceFE (); + + // diagnose + void Diagnose () override; + + // Set up for first time + void + define (SolverChoice& sc) override + { + docloud = sc.do_cloud; + doprecip = sc.do_precip; + m_fac_cond = lcond / sc.c_p; + m_fac_fus = lfus / sc.c_p; + m_fac_sub = lsub / sc.c_p; + m_gOcp = CONST_GRAV / sc.c_p; + m_axis = sc.ave_plane; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) override; + + // update ERF variables + void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist) override; + + // wrapper to do all the updating + void + Advance ( ) override + { + this->AdvanceFE(); + this->Diagnose(); + } + + private: + // geometry + amrex::Geometry m_geom; + // valid boxes on which to evolve the solution + amrex::BoxArray m_gtoe; + + // timestep + amrex::Real dt; + + // number of vertical levels + int nlev, zlo, zhi; + + // plane average axis + int m_axis; + + // model options + bool docloud, doprecip; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // microphysics parameters/coefficients + amrex::TableData accrrc; + amrex::TableData accrsi; + amrex::TableData accrsc; + amrex::TableData coefice; + amrex::TableData evaps1; + amrex::TableData evaps2; + amrex::TableData accrgi; + amrex::TableData accrgc; + amrex::TableData evapg1; + amrex::TableData evapg2; + amrex::TableData evapr1; + amrex::TableData evapr2; + + // vertical plane average data + amrex::TableData rho1d; + amrex::TableData pres1d; + amrex::TableData tabs1d; + amrex::TableData qt1d; + amrex::TableData qv1d; + amrex::TableData qn1d; + + // independent variables + amrex::Array mic_fab_vars; + + amrex::TableData gamaz; + amrex::TableData zmid; // mid value of vertical coordinate in physical domain + + // data (output) + amrex::TableData qifall; + amrex::TableData tlatqi; +}; +#endif diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp new file mode 100644 index 000000000..e3575ed99 --- /dev/null +++ b/Source/Microphysics/FastEddy/FastEddy.cpp @@ -0,0 +1,214 @@ +/* + * this file is modified from precip_proc from samxx + */ +#include "Microphysics.H" + +#include +#include + +using namespace amrex; + +/** + * Compute Precipitation-related Microphysics quantities. + */ +void FastEddy::AdvanceFE() { + + auto qt = mic_fab_vars[MicVar_FE::qt]; + auto qp = mic_fab_vars[MicVar_FE::qp]; + auto qn = mic_fab_vars[MicVar_FE::qn]; + auto tabs = mic_fab_vars[MicVar_FE::tabs]; + + auto qcl = mic_fab_vars[MicVar_FE::qcl]; + auto theta = mic_fab_vars[MicVar_FE::theta]; + auto qv = mic_fab_vars[MicVar_FE::qv]; + auto rho = mic_fab_vars[MicVar_FE::rho]; + + auto dz = m_geom.CellSize(2); + auto domain = m_geom.Domain(); + int k_lo = domain.smallEnd(2); + int k_hi = domain.bigEnd(2); + + MultiFab fz; + auto ba = tabs->boxArray(); + auto dm = tabs->DistributionMap(); + fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells + + for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_FE::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_FE::qn]->array(mfi); + auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_FE::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; + } + }); + } + + + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + + const auto& box3d = mfi.tilebox(); + auto fz_array = fz.array(mfi); + + 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; + }); + }*/ + + + + + // 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_FE::tabs]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi); + auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); + + auto theta_array = theta->array(mfi); + auto qv_array = qv->array(mfi); + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + auto fz_array = fz.array(mfi); + + // Expose for GPU + Real d_fac_cond = m_fac_cond; + + 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; + } + + //------- Autoconversion/accretion + Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; + + Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; + erf_qsatw(tabs_array(i,j,k), pressure, qsat); + + // If there is precipitating water (i.e. rain), and the cell is not saturated + // then the rain water can evaporate leading to extraction of latent heat, hence + // reducing temperature and creating negative buoyancy + + dq_clwater_to_rain = 0.0; + dq_rain_to_vapor = 0.0; + dq_vapor_to_clwater = 0.0; + dq_clwater_to_vapor = 0.0; + + + //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); + Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); + + // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature + if(qv_array(i,j,k) > qsat){ + dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); + } + + + // If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and + // reducing temperature + if(qv_array(i,j,k) < qsat and qn_array(i,j,k) > 0.0){ + dq_clwater_to_vapor = std::min(qn_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); + } + + if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) { + Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046); + dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/ + (5.4e5 + 2.55e6/(pressure*qsat))*dtn; + // The negative sign is to make this variable (vapor formed from evaporation) + // a poistive quantity (as qv/qs < 1) + dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor}); + + // Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature + } + + // If there is cloud water present then do accretion and autoconversion to rain + + if (qn_array(i,j,k) > 0.0) { + qcc = qn_array(i,j,k); + + autor = 0.0; + if (qcc > qcw0) { + autor = 0.001; + } + + accrr = 0.0; + accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875); + dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001)); + + // If the amount of change is more than the amount of qc present, then dq = qc + dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k)); + } + + + Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; + + qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; + qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; + qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; + + theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor); + + qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); + qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); + + }); + } +} diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp new file mode 100644 index 000000000..033d56136 --- /dev/null +++ b/Source/Microphysics/FastEddy/Init_FE.cpp @@ -0,0 +1,249 @@ + +#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; + + auto dz = m_geom.CellSize(2); + auto lowz = m_geom.ProbLo(2); + + 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; + + // 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_FE::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_FE::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); + }); +} diff --git a/Source/Microphysics/FastEddy/Make.package b/Source/Microphysics/FastEddy/Make.package new file mode 100644 index 000000000..06893d229 --- /dev/null +++ b/Source/Microphysics/FastEddy/Make.package @@ -0,0 +1,5 @@ +CEXE_sources += Init_FE.cpp +CEXE_sources += Diagnose_FE.cpp +CEXE_sources += Update_FE.cpp +CEXE_sources += FastEddy.cpp +CEXE_headers += FastEddy.H diff --git a/Source/Microphysics/FastEddy/Update_FE.cpp b/Source/Microphysics/FastEddy/Update_FE.cpp new file mode 100644 index 000000000..a157279f7 --- /dev/null +++ b/Source/Microphysics/FastEddy/Update_FE.cpp @@ -0,0 +1,58 @@ + +#include "Microphysics.H" +#include "IndexDefines.H" +#include "TileNoZ.H" + +/** + * 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::Update (amrex::MultiFab& cons, + amrex::MultiFab& qmoist) +{ + // copy multifab data to qc, qv, and qi + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qv], 0, 0, 1, mic_fab_vars[MicVar_FE::qv]->nGrowVect()); // vapor + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qcl], 0, 1, 1, mic_fab_vars[MicVar_FE::qcl]->nGrowVect()); // cloud water + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qci], 0, 2, 1, mic_fab_vars[MicVar_FE::qci]->nGrowVect()); // cloud ice + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpl], 0, 3, 1, mic_fab_vars[MicVar_FE::qpl]->nGrowVect()); // rain + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpi], 0, 4, 1, mic_fab_vars[MicVar_FE::qpi]->nGrowVect()); // snow + + // Don't need to copy this since it is filled below + // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpi], 0, 5, 1, mic_fab_vars[MicVar_FE::qci]->nGrowVect()); // graupel + + amrex::MultiFab qgraup_mf(qmoist, amrex::make_alias, 5, 1); + + // Get the temperature, density, theta, qt and qp from input + for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto qt_arr = mic_fab_vars[MicVar_FE::qt]->array(mfi); + auto qp_arr = mic_fab_vars[MicVar_FE::qp]->array(mfi); + + auto qgraup_arr= qgraup_mf.array(mfi); + + const auto& box3d = mfi.tilebox(); + + // get potential total density, temperature, qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + + // Graupel == precip total - rain - snow (but must be >= 0) + qgraup_arr(i,j,k) = 0.0;// + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); + qmoist.FillBoundary(m_geom.periodicity()); +} + + diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 6c32a0ff2..a8624af65 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -202,7 +202,7 @@ void Kessler::AdvanceKessler() { Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; //dq_sed = 0.0; - + qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index 1338f7674..c384ad8a0 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -4,6 +4,7 @@ #include #include #include +#include class Microphysics { From e5d00c597722dae9b12094edac83fbfb2bf6515a Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 10 Dec 2023 18:54:04 -0800 Subject: [PATCH 2/3] fix tabs --- Source/Microphysics/Kessler/Kessler.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index a8624af65..64fe35639 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -197,11 +197,11 @@ void Kessler::AdvanceKessler() { } - if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; - if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; + if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; - if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; - //dq_sed = 0.0; + if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; + //dq_sed = 0.0; qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; From b6ead7d352cb2db8979f1b5f31a8005033e027bd Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 10 Dec 2023 19:03:00 -0800 Subject: [PATCH 3/3] fix cmake build --- CMake/BuildERFExe.cmake | 4 ++++ Docs/sphinx_doc/Inputs.rst | 3 ++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index ad50a774d..b64704f24 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -159,6 +159,10 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Microphysics/Kessler/Kessler.cpp ${SRC_DIR}/Microphysics/Kessler/Diagnose_Kessler.cpp ${SRC_DIR}/Microphysics/Kessler/Update_Kessler.cpp + ${SRC_DIR}/Microphysics/FastEddy/Init_FE.cpp + ${SRC_DIR}/Microphysics/FastEddy/FastEddy.cpp + ${SRC_DIR}/Microphysics/FastEddy/Diagnose_FE.cpp + ${SRC_DIR}/Microphysics/FastEddy/Update_FE.cpp ) if(NOT "${erf_exe_name}" STREQUAL "erf_unit_tests") diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index ebf539b16..dbc767156 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -893,7 +893,8 @@ List of Parameters | Parameter | Definition | Acceptable | Default | | | | Values | | +=============================+==========================+====================+============+ -| **erf.moisture_model** | Name of moisture model | "SAM", "Kessler" | "Null" | +| **erf.moisture_model** | Name of moisture model | "SAM", "Kessler", | "Null" | +| | | "FastEddy" | | +-----------------------------+--------------------------+--------------------+------------+ | **erf.do_cloud** | use basic moisture model | true / false | true | +-----------------------------+--------------------------+--------------------+------------+