From c25fdba3a01c577340d5d32c21ca75a5979e7c32 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 7 Mar 2024 17:43:39 -0800 Subject: [PATCH 1/2] Removing FastEddy completely (#1483) * Removing FastEddy * Updating CMake to reflect FastEddy removal * Removing Plotfile.cpp.convergence --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- CMake/BuildERFExe.cmake | 5 - Exec/Make.ERF | 5 - Exec/RegTests/Bubble/inputs_BF02_moist_bubble | 1 - Source/DataStructs/DataStruct.H | 6 +- Source/ERF.cpp | 3 - Source/IO/Plotfile.cpp.convergence | 1068 ----------------- Source/Microphysics/FastEddy/Diagnose_FE.cpp | 8 - Source/Microphysics/FastEddy/FastEddy.H | 156 --- Source/Microphysics/FastEddy/FastEddy.cpp | 70 -- Source/Microphysics/FastEddy/Init_FE.cpp | 91 -- Source/Microphysics/FastEddy/Make.package | 5 - Source/Microphysics/FastEddy/Update_FE.cpp | 38 - Source/Microphysics/Microphysics.H | 1 - Source/TimeIntegration/ERF_make_buoyancy.cpp | 2 +- 14 files changed, 3 insertions(+), 1456 deletions(-) delete mode 100644 Source/IO/Plotfile.cpp.convergence delete mode 100644 Source/Microphysics/FastEddy/Diagnose_FE.cpp delete mode 100644 Source/Microphysics/FastEddy/FastEddy.H delete mode 100644 Source/Microphysics/FastEddy/FastEddy.cpp delete mode 100644 Source/Microphysics/FastEddy/Init_FE.cpp delete mode 100644 Source/Microphysics/FastEddy/Make.package delete mode 100644 Source/Microphysics/FastEddy/Update_FE.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index d97720169..64f6f947e 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -162,10 +162,6 @@ 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 ${SRC_DIR}/LandSurfaceModel/SLM/SLM.cpp ${SRC_DIR}/LandSurfaceModel/MM5/MM5.cpp ) @@ -212,7 +208,6 @@ 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) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/LandSurfaceModel) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/LandSurfaceModel/Null) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/LandSurfaceModel/SLM) diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 0478133d0..370c29038 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -112,11 +112,6 @@ 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/Exec/RegTests/Bubble/inputs_BF02_moist_bubble b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble index 9288ab3e9..af09852f0 100644 --- a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble +++ b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble @@ -53,7 +53,6 @@ erf.moistscal_vert_adv_type = "Upwind_3rd" # PHYSICS OPTIONS erf.les_type = "None" erf.pbl_type = "None" -#erf.moisture_model = "FastEddy" erf.moisture_model = "Kessler_NoRain" erf.buoyancy_type = 1 erf.use_moist_background = true diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 1ef2e096a..a87161180 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -29,7 +29,7 @@ enum struct TerrainType { }; enum struct MoistureType { - Kessler, SAM, FastEddy, Kessler_NoRain, None + Kessler, SAM, Kessler_NoRain, None }; enum struct WindFarmType { @@ -91,8 +91,6 @@ struct SolverChoice { moisture_type = MoistureType::Kessler; }else if (moisture_model_string == "Kessler_NoRain") { moisture_type = MoistureType::Kessler_NoRain; - } else if (moisture_model_string == "FastEddy") { - moisture_type = MoistureType::FastEddy; } else { moisture_type = MoistureType::None; } @@ -100,7 +98,7 @@ struct SolverChoice { // TODO: should we set default for dry?? // Set a different default for moist vs dry if (moisture_type != MoistureType::None) { - if (moisture_type == MoistureType::FastEddy || + if (moisture_type == MoistureType::Kessler_NoRain || moisture_type == MoistureType::SAM) { buoyancy_type = 1; // asserted in make buoyancy } else { diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 5156692f3..aa6e842f6 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1232,9 +1232,6 @@ ERF::ReadParameters () solverChoice.moisture_type == MoistureType::Kessler_NoRain) { micro.SetModel(); amrex::Print() << "Kessler moisture model!\n"; - } else if (solverChoice.moisture_type == MoistureType::FastEddy) { - micro.SetModel(); - amrex::Print() << "FastEddy moisture model!\n"; } 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/Plotfile.cpp.convergence b/Source/IO/Plotfile.cpp.convergence deleted file mode 100644 index 82af98b5c..000000000 --- a/Source/IO/Plotfile.cpp.convergence +++ /dev/null @@ -1,1068 +0,0 @@ -#include -#include -#include "AMReX_Interp_3D_C.H" -#include "AMReX_PlotFileUtil.H" -#include "TerrainMetrics.H" -#include "ERF_Constants.H" - -using namespace amrex; - -template -bool containerHasElement(const V& iterable, const T& query) { - return std::find(iterable.begin(), iterable.end(), query) != iterable.end(); -} - -void -ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector& plot_var_names) -{ - ParmParse pp(pp_prefix); - - if (pp.contains(pp_plot_var_names.c_str())) - { - std::string nm; - - int nPltVars = pp.countval(pp_plot_var_names.c_str()); - - for (int i = 0; i < nPltVars; i++) - { - pp.get(pp_plot_var_names.c_str(), nm, i); - - // Add the named variable to our list of plot variables - // if it is not already in the list - if (!containerHasElement(plot_var_names, nm)) { - plot_var_names.push_back(nm); - } - } - } else { - // - // The default is to add none of the variables to the list - // - plot_var_names.clear(); - } - - // Get state variables in the same order as we define them, - // since they may be in any order in the input list - Vector tmp_plot_names; - for (int i = 0; i < Cons::NumVars; ++i) { - if ( containerHasElement(plot_var_names, cons_names[i]) ) { - tmp_plot_names.push_back(cons_names[i]); - } - } - // check for velocity since it's not in cons_names - // if we are asked for any velocity component, we will need them all - if (containerHasElement(plot_var_names, "x_velocity") || - containerHasElement(plot_var_names, "y_velocity") || - containerHasElement(plot_var_names, "z_velocity")) { - tmp_plot_names.push_back("x_velocity"); - tmp_plot_names.push_back("y_velocity"); - tmp_plot_names.push_back("z_velocity"); - } - for (int i = 0; i < derived_names.size(); ++i) { - if ( containerHasElement(plot_var_names, derived_names[i]) ) { - if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { - tmp_plot_names.push_back(derived_names[i]); - } - } - } - - // Check to see if we found all the requested variables - for (const auto& plot_name : plot_var_names) { - if (!containerHasElement(tmp_plot_names, plot_name)) { - if (amrex::ParallelDescriptor::IOProcessor()) { - Warning("\nWARNING: Requested to plot variable '" + plot_name + "' but it is not available"); - } - } - } - plot_var_names = tmp_plot_names; -} - -// set plotfile variable names -Vector -ERF::PlotFileVarNames ( Vector plot_var_names ) -{ - Vector names; - - names.insert(names.end(), plot_var_names.begin(), plot_var_names.end()); - - return names; - -} - -// Write plotfile to disk -void -ERF::WritePlotFile (int which, Vector plot_var_names) -{ - const Vector varnames = PlotFileVarNames(plot_var_names); - const int ncomp_mf = varnames.size(); - - // We fillpatch here because some of the derived quantities require derivatives - // which require ghost cells to be filled - for (int lev = 0; lev <= finest_level; ++lev) { - FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], - &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}); - } - - if (ncomp_mf == 0) - return; - - Vector mf(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - mf[lev].define(grids[lev], dmap[lev], ncomp_mf, 0); - } - - Vector mf_nd(finest_level+1); - if (solverChoice.use_terrain) { - for (int lev = 0; lev <= finest_level; ++lev) { - BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes(); - mf_nd[lev].define(nodal_grids, dmap[lev], 3, 0); - mf_nd[lev].setVal(0.); - } - } - - for (int lev = 0; lev <= finest_level; ++lev) { - int mf_comp = 0; - - // First, copy any of the conserved state variables into the output plotfile - AMREX_ALWAYS_ASSERT(cons_names.size() == Cons::NumVars); - for (int i = 0; i < Cons::NumVars; ++i) { - if (containerHasElement(plot_var_names, cons_names[i])) { - MultiFab::Copy(mf[lev],vars_new[lev][Vars::cons],i,mf_comp,1,0); - mf_comp++; - } - } - - // Next, check for velocities and if desired, output them -- note we output none or all, not just some - if (containerHasElement(plot_var_names, "x_velocity") || - containerHasElement(plot_var_names, "y_velocity") || - containerHasElement(plot_var_names, "z_velocity")) { - - average_face_to_cellcenter(mf[lev],mf_comp, - Array{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]}); - mf_comp += AMREX_SPACEDIM; - } - - // Finally, check for any derived quantities and compute them, inserting - // them into our output multifab - auto calculate_derived = [&](const std::string& der_name, - decltype(derived::erf_dernull)& der_function) - { - if (containerHasElement(plot_var_names, der_name)) { - MultiFab dmf(mf[lev], make_alias, mf_comp, 1); -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - auto& dfab = dmf[mfi]; - auto& sfab = vars_new[lev][Vars::cons][mfi]; - der_function(bx, dfab, 0, 1, sfab, Geom(lev), t_new[0], nullptr, lev); - } - - mf_comp++; - } - }; - - // Note: All derived variables must be computed in order of "derived_names" defined in ERF.H - calculate_derived("pressure", derived::erf_derpres); - calculate_derived("soundspeed", derived::erf_dersoundspeed); - calculate_derived("temp", derived::erf_dertemp); - calculate_derived("theta", derived::erf_dertheta); - calculate_derived("KE", derived::erf_derKE); - calculate_derived("QKE", derived::erf_derQKE); - calculate_derived("scalar", derived::erf_derscalar); - - 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 - if (containerHasElement(plot_var_names, "pres_hse")) - { - // p_0 is second component of base_state - MultiFab::Copy(mf[lev],p_hse,0,mf_comp,1,0); - mf_comp += 1; - } - if (containerHasElement(plot_var_names, "dens_hse")) - { - // r_0 is first component of base_state - MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0); - mf_comp += 1; - } - if (containerHasElement(plot_var_names, "pert_pres")) - { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - 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); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - - const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); - derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta) - p0_arr(i,j,k); - }); - } - mf_comp += 1; - } - if (containerHasElement(plot_var_names, "pert_dens")) - { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - 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& r0_arr = r_hse.const_array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - derdat(i, j, k, mf_comp) = S_arr(i,j,k,Rho_comp) - r0_arr(i,j,k); - }); - } - mf_comp ++; - } - - int klo = geom[lev].Domain().smallEnd(2); - int khi = geom[lev].Domain().bigEnd(2); - - if (containerHasElement(plot_var_names, "dpdx")) - { - auto dxInv = geom[lev].InvCellSizeArray(); - MultiFab pres(vars_new[lev][Vars::cons].boxArray(), vars_new[lev][Vars::cons].DistributionMap(), 1, 1); -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // First define pressure on grown box - const Box& gbx = mfi.growntilebox(1); - const Array4 & p_arr = pres.array(mfi); - const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - p_arr(i,j,k) = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp)); - }); - } - pres.FillBoundary(geom[lev].periodicity()); - -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // Now compute pressure gradient on valid box - const Box& bx = mfi.tilebox(); - const Array4& derdat = mf[lev].array(mfi); - const Array4 & p_arr = pres.array(mfi); - - if (solverChoice.use_terrain) { - const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - - // Pgrad at lower I face - Real met_h_xi_lo = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); - Real met_h_zeta_lo = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); - Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k)); - Real gp_zeta_on_iface_lo; - if(k == klo) { - gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( - p_arr(i-1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i-1,j,k ) - p_arr(i,j,k ) ); - } else if (k == khi) { - gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( - p_arr(i-1,j,k ) + p_arr(i,j,k ) - - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); - } else { - gp_zeta_on_iface_lo = 0.25 * dxInv[2] * ( - p_arr(i-1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); - } - amrex::Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo; - - // Pgrad at higher I face - Real met_h_xi_hi = Compute_h_xi_AtIface (i+1, j, k, dxInv, z_nd); - Real met_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j, k, dxInv, z_nd); - Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k)); - Real gp_zeta_on_iface_hi; - if(k == klo) { - gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( - p_arr(i+1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i+1,j,k ) - p_arr(i,j,k ) ); - } else if (k == khi) { - gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( - p_arr(i+1,j,k ) + p_arr(i,j,k ) - - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); - } else { - gp_zeta_on_iface_hi = 0.25 * dxInv[2] * ( - p_arr(i+1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); - } - amrex::Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi; - - // Average P grad to CC - derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi); - }); - } else { - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i+1,j,k) - p_arr(i-1,j,k)) * dxInv[0]; - }); - } - } // mfi - mf_comp ++; - } // dpdx - - if (containerHasElement(plot_var_names, "dpdy")) - { - auto dxInv = geom[lev].InvCellSizeArray(); - - MultiFab pres(vars_new[lev][Vars::cons].boxArray(), vars_new[lev][Vars::cons].DistributionMap(), 1, 1); -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // First define pressure on grown box - const Box& gbx = mfi.growntilebox(1); - const Array4 & p_arr = pres.array(mfi); - const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - p_arr(i,j,k) = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp)); - }); - } - pres.FillBoundary(geom[lev].periodicity()); - -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // Now compute pressure gradient on valid box - const Box& bx = mfi.tilebox(); - const Array4& derdat = mf[lev].array(mfi); - const Array4 & p_arr = pres.array(mfi); - - if (solverChoice.use_terrain) { - const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - - Real met_h_eta_lo = Compute_h_eta_AtJface (i, j, k, dxInv, z_nd); - Real met_h_zeta_lo = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); - Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k)); - Real gp_zeta_on_jface_lo; - if (k == klo) { - gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( - p_arr(i,j,k+1) + p_arr(i,j-1,k+1) - - p_arr(i,j,k ) - p_arr(i,j-1,k ) ); - } else if (k == khi) { - gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( - p_arr(i,j,k ) + p_arr(i,j-1,k ) - - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); - } else { - gp_zeta_on_jface_lo = 0.25 * dxInv[2] * ( - p_arr(i,j,k+1) + p_arr(i,j-1,k+1) - - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); - } - amrex::Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo; - - Real met_h_eta_hi = Compute_h_eta_AtJface (i, j+1, k, dxInv, z_nd); - Real met_h_zeta_hi = Compute_h_zeta_AtJface(i, j+1, k, dxInv, z_nd); - Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k)); - Real gp_zeta_on_jface_hi; - if (k == klo) { - gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( - p_arr(i,j+1,k+1) + p_arr(i,j,k+1) - - p_arr(i,j+1,k ) - p_arr(i,j,k ) ); - } else if (k == khi) { - gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( - p_arr(i,j+1,k ) + p_arr(i,j,k ) - - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); - } else { - gp_zeta_on_jface_hi = 0.25 * dxInv[2] * ( - p_arr(i,j+1,k+1) + p_arr(i,j,k+1) - - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); - } - amrex::Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi; - - derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi); - }); - } else { - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i,j+1,k) - p_arr(i,j-1,k)) * dxInv[1]; - }); - } - } // mf - mf_comp ++; - } // dpdy - - if (containerHasElement(plot_var_names, "pres_hse_x")) - { - auto dxInv = geom[lev].InvCellSizeArray(); -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - const Array4& derdat = mf[lev].array(mfi); - const Array4& p_arr = p_hse.const_array(mfi); - - //USE_TERRAIN POSSIBLE ISSUE HERE - const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real met_h_xi_lo = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); - Real met_h_zeta_lo = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); - Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k)); - Real gp_zeta_on_iface_lo; - if (k == klo) { - gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( - p_arr(i-1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i-1,j,k ) - p_arr(i,j,k ) ); - } else if (k == khi) { - gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( - p_arr(i-1,j,k ) + p_arr(i,j,k ) - - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); - } else { - gp_zeta_on_iface_lo = 0.25 * dxInv[2] * ( - p_arr(i-1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); - } - amrex::Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo; - - Real met_h_xi_hi = Compute_h_xi_AtIface (i+1, j, k, dxInv, z_nd); - Real met_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j, k, dxInv, z_nd); - Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k)); - Real gp_zeta_on_iface_hi; - if (k == klo) { - gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( - p_arr(i+1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i+1,j,k ) - p_arr(i,j,k ) ); - } else if (k == khi) { - gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( - p_arr(i+1,j,k ) + p_arr(i,j,k ) - - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); - } else { - gp_zeta_on_iface_hi = 0.25 * dxInv[2] * ( - p_arr(i+1,j,k+1) + p_arr(i,j,k+1) - - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); - } - amrex::Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi; - - derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi); - }); - } - mf_comp += 1; - } // pres_hse_x - - if (containerHasElement(plot_var_names, "pres_hse_y")) - { - auto dxInv = geom[lev].InvCellSizeArray(); -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - const Array4& derdat = mf[lev].array(mfi); - const Array4& p_arr = p_hse.const_array(mfi); - const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real met_h_eta_lo = Compute_h_eta_AtJface (i, j, k, dxInv, z_nd); - Real met_h_zeta_lo = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); - Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k)); - Real gp_zeta_on_jface_lo; - if (k == klo) { - gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( - p_arr(i,j,k+1) + p_arr(i,j-1,k+1) - - p_arr(i,j,k ) - p_arr(i,j-1,k ) ); - } else if (k == khi) { - gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( - p_arr(i,j,k ) + p_arr(i,j-1,k ) - - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); - } else { - gp_zeta_on_jface_lo = 0.25 * dxInv[2] * ( - p_arr(i,j,k+1) + p_arr(i,j-1,k+1) - - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); - } - amrex::Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo; - - Real met_h_eta_hi = Compute_h_eta_AtJface (i, j+1, k, dxInv, z_nd); - Real met_h_zeta_hi = Compute_h_zeta_AtJface(i, j+1, k, dxInv, z_nd); - Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k)); - Real gp_zeta_on_jface_hi; - if (k == klo) { - gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( - p_arr(i,j+1,k+1) + p_arr(i,j,k+1) - - p_arr(i,j+1,k ) - p_arr(i,j,k ) ); - } else if (k == khi) { - gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( - p_arr(i,j+1,k ) + p_arr(i,j,k ) - - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); - } else { - gp_zeta_on_jface_hi = 0.25 * dxInv[2] * ( - p_arr(i,j+1,k+1) + p_arr(i,j,k+1) - - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); - } - amrex::Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi; - - derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi); - }); - } - mf_comp += 1; - } // pres_hse_y - - if (solverChoice.use_terrain) { - if (containerHasElement(plot_var_names, "z_phys")) - { - MultiFab::Copy(mf[lev],*z_phys_cc[lev],0,mf_comp,1,0); - mf_comp ++; - } - - if (containerHasElement(plot_var_names, "detJ")) - { - MultiFab::Copy(mf[lev],*detJ_cc[lev],0,mf_comp,1,0); - mf_comp ++; - } - } // use_terrain - - if (containerHasElement(plot_var_names, "mapfac")) { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - const Array4& derdat = mf[lev].array(mfi); - const Array4& mf_m = mapfac_m[lev]->array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - derdat(i ,j ,k, mf_comp) = mf_m(i,j,0); - }); - } - mf_comp ++; - } - - if (solverChoice.moisture_type != MoistureType::None) { - calculate_derived("qt", derived::erf_derQt); - calculate_derived("qp", derived::erf_derQp); - - MultiFab qv_mf(qmoist[lev], make_alias, 0, 1); - MultiFab qc_mf(qmoist[lev], make_alias, 1, 1); - MultiFab qi_mf(qmoist[lev], make_alias, 2, 1); - MultiFab qr_mf(qmoist[lev], make_alias, 3, 1); - MultiFab qs_mf(qmoist[lev], make_alias, 4, 1); - MultiFab qg_mf(qmoist[lev], make_alias, 5, 1); - - if (containerHasElement(plot_var_names, "qv")) - { - MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); - mf_comp += 1; - } - - if (containerHasElement(plot_var_names, "qc")) - { - MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); - mf_comp += 1; - } - - if (containerHasElement(plot_var_names, "qi")) - { - MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); - mf_comp += 1; - } - - if (containerHasElement(plot_var_names, "qrain")) - { - MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); - mf_comp += 1; - } - - if (containerHasElement(plot_var_names, "qsnow")) - { - MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0); - mf_comp += 1; - } - - if (containerHasElement(plot_var_names, "qgraup")) - { - MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); - mf_comp += 1; - } - } - -#ifdef ERF_COMPUTE_ERROR - // Next, check for error in velocities and if desired, output them -- note we output none or all, not just some - if (containerHasElement(plot_var_names, "xvel_err") || - containerHasElement(plot_var_names, "yvel_err") || - containerHasElement(plot_var_names, "zvel_err")) - { - // - // Moving terrain ANALYTICAL - // - Real H = geom[lev].ProbHi()[2]; - Real Ampl = 0.16; - Real wavelength = 100.; - Real kp = 2. * PI / wavelength; - Real g = CONST_GRAV; - Real omega = std::sqrt(g * kp); - Real omega_t = omega * t_new[lev]; - - const auto dx = geom[lev].CellSizeArray(); - -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - Box xbx(bx); xbx.surroundingNodes(0); - const Array4 xvel_arr = vars_new[lev][Vars::xvel].array(mfi); - const Array4 zvel_arr = vars_new[lev][Vars::zvel].array(mfi); - - const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); - - ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real x = i * dx[0]; - Real z = 0.25 * (z_nd(i,j,k) + z_nd(i,j+1,k) + z_nd(i,j,k+1) + z_nd(i,j+1,k+1)); - - Real z_base = Ampl * std::sin(kp * x - omega_t); - z -= z_base; - - Real fac = std::cosh( kp * (z - H) ) / std::sinh(kp * H); - - xvel_arr(i,j,k) -= -Ampl * omega * fac * std::sin(kp * x - omega_t); - }); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real x = (i + 0.5) * dx[0]; - Real z = 0.25 * ( z_nd(i,j,k) + z_nd(i+1,j,k) + z_nd(i,j+1,k) + z_nd(i+1,j+1,k)); - - Real z_base = Ampl * std::sin(kp * x - omega_t); - z -= z_base; - - Real fac = std::sinh( kp * (z - H) ) / std::sinh(kp * H); - - zvel_arr(i,j,k) -= Ampl * omega * fac * std::cos(kp * x - omega_t); - }); - } - - MultiFab temp_mf(mf[lev].boxArray(), mf[lev].DistributionMap(), AMREX_SPACEDIM, 0); - average_face_to_cellcenter(temp_mf,0, - Array{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]}); - - if (containerHasElement(plot_var_names, "xvel_err")) { - MultiFab::Copy(mf[lev],temp_mf,0,mf_comp,1,0); - mf_comp += 1; - } - if (containerHasElement(plot_var_names, "yvel_err")) { - MultiFab::Copy(mf[lev],temp_mf,1,mf_comp,1,0); - mf_comp += 1; - } - if (containerHasElement(plot_var_names, "zvel_err")) { - MultiFab::Copy(mf[lev],temp_mf,2,mf_comp,1,0); - mf_comp += 1; - } - - // Now restore the velocities to what they were -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - Box xbx(bx); xbx.surroundingNodes(0); - - const Array4 xvel_arr = vars_new[lev][Vars::xvel].array(mfi); - const Array4 zvel_arr = vars_new[lev][Vars::zvel].array(mfi); - - const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); - - ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real x = i * dx[0]; - Real z = 0.25 * (z_nd(i,j,k) + z_nd(i,j+1,k) + z_nd(i,j,k+1) + z_nd(i,j+1,k+1)); - Real z_base = Ampl * std::sin(kp * x - omega_t); - - z -= z_base; - - Real fac = std::cosh( kp * (z - H) ) / std::sinh(kp * H); - xvel_arr(i,j,k) += -Ampl * omega * fac * std::sin(kp * x - omega_t); - }); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real x = (i + 0.5) * dx[0]; - Real z = 0.25 * ( z_nd(i,j,k) + z_nd(i+1,j,k) + z_nd(i,j+1,k) + z_nd(i+1,j+1,k)); - Real z_base = Ampl * std::sin(kp * x - omega_t); - - z -= z_base; - Real fac = std::sinh( kp * (z - H) ) / std::sinh(kp * H); - - zvel_arr(i,j,k) += Ampl * omega * fac * std::cos(kp * x - omega_t); - }); - } - } // end xvel_err, yvel_err, zvel_err - - if (containerHasElement(plot_var_names, "pp_err")) - { - // Moving terrain ANALYTICAL -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - 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 auto dx = geom[lev].CellSizeArray(); -#if 0 - const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); - const Array4& p0_arr = p_hse.const_array(mfi); - const Array4& r0_arr = r_hse.const_array(mfi); - - Real H = geom[lev].ProbHi()[2]; - Real Ampl = 0.16; - Real wavelength = 100.; - Real kp = 2. * PI / wavelength; - Real g = CONST_GRAV; - Real omega = std::sqrt(g * kp); - Real omega_t = omega * t_new[lev]; - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); - derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta) - p0_arr(i,j,k); - - Real rho_hse = r0_arr(i,j,k); - - Real x = (i + 0.5) * dx[0]; - Real z = 0.125 * ( z_nd(i,j,k ) + z_nd(i+1,j,k ) + z_nd(i,j+1,k ) + z_nd(i+1,j+1,k ) - +z_nd(i,j,k+1) + z_nd(i+1,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k+1) ); - Real z_base = Ampl * std::sin(kp * x - omega_t); - - z -= z_base; - Real fac = std::cosh( kp * (z - H) ) / std::sinh(kp * H); - Real pprime_exact = -(Ampl * omega * omega / kp) * fac * - std::sin(kp * x - omega_t) * r0_arr(i,j,k); - - derdat(i,j,k,mf_comp) -= pprime_exact; - }); -#endif - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - Real x = (i + 0.5) * dx[0]; - Real y = (j + 0.5) * dx[1]; - Real S_exact = cos(PI * x); - - derdat(i,j,k,mf_comp) = S_arr(i,j,k,RhoScalar_comp) - S_exact; - }); - } // mfi -#endif - mf_comp += 1; - - int nx = geom[lev].Domain().length(0); - int nv = std::sqrt(static_cast(geom[lev].Domain().numPts())); - Real dx = 2.0 / static_cast(nx); - Real norm_0 = mf[0].norm0(mf_comp-1,0); - Real norm_2 = mf[0].norm2(mf_comp-1) / nv; - amrex::Print() << " \n ERR: dx = " << dx << " " << norm_2 << std::endl; -// << " order = " << solverChoice.horiz_spatial_order -// << " 2-norm = " << norm_2 -// << " 0-norm = " << norm_0 << "\n" << std::endl; - } // pp_err - } // lev - - std::string plotfilename; - if (which == 1) - plotfilename = Concatenate(plot_file_1, istep[0], 5); - else if (which == 2) - plotfilename = Concatenate(plot_file_2, istep[0], 5); - - if (finest_level == 0) - { - if (plotfile_type == "amrex") { - amrex::Print() << "Writing plotfile " << plotfilename << "\n"; - if (solverChoice.use_terrain) { - // We started with mf_nd holding 0 in every component; here we fill only the offset in z - int lev = 0; - MultiFab::Copy(mf_nd[lev],*z_phys_nd[lev],0,2,1,0); - Real dz = Geom()[lev].CellSizeArray()[2]; - for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); - Array4< Real> mf_arr = mf_nd[lev].array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - mf_arr(i,j,k,2) -= k * dz; - }); - } - WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, - GetVecOfConstPtrs(mf), - GetVecOfConstPtrs(mf_nd), - varnames, - t_new[0], istep); - } else { - WriteMultiLevelPlotfile(plotfilename, finest_level+1, - GetVecOfConstPtrs(mf), - varnames, - Geom(), t_new[0], istep, refRatio()); - } - writeJobInfo(plotfilename); -#ifdef ERF_USE_HDF5 - } else if (plotfile_type == "hdf5" || plotfile_type == "HDF5") { - amrex::Print() << "Writing plotfile " << plotfilename+"d01.h5" << "\n"; - WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1, - GetVecOfConstPtrs(mf), - varnames, - Geom(), t_new[0], istep, refRatio()); -#endif -#ifdef ERF_USE_NETCDF - } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { - int lev = 0; - int l_which = 0; - writeNCPlotFile(lev, l_which, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]); -#endif - } else { - amrex::Print() << "User specified plot_filetype = " << plotfile_type << std::endl; - amrex::Abort("Dont know this plot_filetype"); - } - - } else { // multilevel - - Vector r2(finest_level); - Vector g2(finest_level+1); - Vector mf2(finest_level+1); - - mf2[0].define(grids[0], dmap[0], ncomp_mf, 0); - - // Copy level 0 as is - MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0); - - // Define a new multi-level array of Geometry's so that we pass the new "domain" at lev > 0 - Array periodicity = - {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)}; - g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data()); - - if (plotfile_type == "amrex") { - r2[0] = IntVect(1,1,ref_ratio[0][0]); - for (int lev = 1; lev <= finest_level; ++lev) { - if (lev > 1) { - r2[lev-1][0] = 1; - r2[lev-1][1] = 1; - r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0]; - } - - mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0); - - // Set the new problem domain - Box d2(Geom()[lev].Domain()); - d2.refine(r2[lev-1]); - - g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data()); - } - - // Do piecewise interpolation of mf into mf2 - for (int lev = 1; lev <= finest_level; ++lev) { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(mf2[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); - pcinterp_interp(bx,mf2[lev].array(mfi), 0, mf[lev].nComp(), mf[lev].const_array(mfi),0,r2[lev-1]); - } - } - - // Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile - Vector rr(finest_level); - for (int lev = 0; lev < finest_level; ++lev) { - rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]); - } - - WriteMultiLevelPlotfile(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), varnames, - g2, t_new[0], istep, rr); - writeJobInfo(plotfilename); -#ifdef ERF_USE_NETCDF - } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { - for (int lev = 0; lev <= finest_level; ++lev) { - for (int which_box = 0; which_box < num_boxes_at_level[lev]; which_box++) { - writeNCPlotFile(lev, which_box, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]); - } - } -#endif - } - } // end multi-level -} - -void -ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nlevels, - const Vector& mf, - const Vector& mf_nd, - const Vector& varnames, - Real time, - const Vector& level_steps, - const std::string &versionName, - const std::string &levelPrefix, - const std::string &mfPrefix, - const Vector& extra_dirs) const -{ - BL_PROFILE("WriteMultiLevelPlotfileWithTerrain()"); - - BL_ASSERT(nlevels <= mf.size()); - BL_ASSERT(nlevels <= ref_ratio.size()+1); - BL_ASSERT(nlevels <= level_steps.size()); - BL_ASSERT(mf[0]->nComp() == varnames.size()); - - bool callBarrier(false); - PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier); - if (!extra_dirs.empty()) { - for (const auto& d : extra_dirs) { - const std::string ed = plotfilename+"/"+d; - PreBuildDirectorHierarchy(ed, levelPrefix, nlevels, callBarrier); - } - } - ParallelDescriptor::Barrier(); - - if (ParallelDescriptor::MyProc() == ParallelDescriptor::NProcs()-1) { - Vector boxArrays(nlevels); - for(int level(0); level < boxArrays.size(); ++level) { - boxArrays[level] = mf[level]->boxArray(); - } - - auto f = [=]() { - VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); - std::string HeaderFileName(plotfilename + "/Header"); - std::ofstream HeaderFile; - HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); - HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out | - std::ofstream::trunc | - std::ofstream::binary); - if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName); - WriteGenericPlotfileHeaderWithTerrain(HeaderFile, nlevels, boxArrays, varnames, - time, level_steps, versionName, - levelPrefix, mfPrefix); - }; - - if (AsyncOut::UseAsyncOut()) { - AsyncOut::Submit(std::move(f)); - } else { - f(); - } - } - - std::string mf_nodal_prefix = "Nu_nd"; - for (int level = 0; level <= finest_level; ++level) - { - if (AsyncOut::UseAsyncOut()) { - VisMF::AsyncWrite(*mf[level], - MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix), - true); - VisMF::AsyncWrite(*mf_nd[level], - MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix), - true); - } else { - const MultiFab* data; - std::unique_ptr mf_tmp; - if (mf[level]->nGrowVect() != 0) { - mf_tmp = std::make_unique(mf[level]->boxArray(), - mf[level]->DistributionMap(), - mf[level]->nComp(), 0, MFInfo(), - mf[level]->Factory()); - MultiFab::Copy(*mf_tmp, *mf[level], 0, 0, mf[level]->nComp(), 0); - data = mf_tmp.get(); - } else { - data = mf[level]; - } - VisMF::Write(*data , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix)); - VisMF::Write(*mf_nd[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix)); - } - } -} - -void -ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile, - int nlevels, - const Vector &bArray, - const Vector &varnames, - Real time, - const Vector &level_steps, - const std::string &versionName, - const std::string &levelPrefix, - const std::string &mfPrefix) const -{ - BL_ASSERT(nlevels <= bArray.size()); - BL_ASSERT(nlevels <= ref_ratio.size()+1); - BL_ASSERT(nlevels <= level_steps.size()); - - HeaderFile.precision(17); - - // ---- this is the generic plot file type name - HeaderFile << versionName << '\n'; - - HeaderFile << varnames.size() << '\n'; - - for (int ivar = 0; ivar < varnames.size(); ++ivar) { - HeaderFile << varnames[ivar] << "\n"; - } - HeaderFile << AMREX_SPACEDIM << '\n'; - HeaderFile << time << '\n'; - HeaderFile << finest_level << '\n'; - for (int i = 0; i < AMREX_SPACEDIM; ++i) { - HeaderFile << geom[0].ProbLo(i) << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i < AMREX_SPACEDIM; ++i) { - HeaderFile << geom[0].ProbHi(i) << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i < finest_level; ++i) { - HeaderFile << ref_ratio[i][0] << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i <= finest_level; ++i) { - HeaderFile << geom[i].Domain() << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i <= finest_level; ++i) { - HeaderFile << level_steps[i] << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i <= finest_level; ++i) { - for (int k = 0; k < AMREX_SPACEDIM; ++k) { - HeaderFile << geom[i].CellSize()[k] << ' '; - } - HeaderFile << '\n'; - } - HeaderFile << (int) geom[0].Coord() << '\n'; - HeaderFile << "0\n"; - - for (int level = 0; level <= finest_level; ++level) { - HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n'; - HeaderFile << level_steps[level] << '\n'; - - const IntVect& domain_lo = geom[level].Domain().smallEnd(); - for (int i = 0; i < bArray[level].size(); ++i) - { - // Need to shift because the RealBox ctor we call takes the - // physical location of index (0,0,0). This does not affect - // the usual cases where the domain index starts with 0. - const Box& b = shift(bArray[level][i], -domain_lo); - RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo()); - for (int n = 0; n < AMREX_SPACEDIM; ++n) { - HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n'; - } - } - - HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n'; - } - HeaderFile << "1" << "\n"; - HeaderFile << "3" << "\n"; - HeaderFile << "amrexvec_nu_x" << "\n"; - HeaderFile << "amrexvec_nu_y" << "\n"; - HeaderFile << "amrexvec_nu_z" << "\n"; - std::string mf_nodal_prefix = "Nu_nd"; - for (int level = 0; level <= finest_level; ++level) { - HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n'; - } -} diff --git a/Source/Microphysics/FastEddy/Diagnose_FE.cpp b/Source/Microphysics/FastEddy/Diagnose_FE.cpp deleted file mode 100644 index c4ae826d7..000000000 --- a/Source/Microphysics/FastEddy/Diagnose_FE.cpp +++ /dev/null @@ -1,8 +0,0 @@ -#include "FastEddy.H" - -/** - * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid - * from the existing Microphysics variables. - */ -void FastEddy::Diagnose () { -} diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H deleted file mode 100644 index 80431088f..000000000 --- a/Source/Microphysics/FastEddy/FastEddy.H +++ /dev/null @@ -1,156 +0,0 @@ -#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" -#include "NullMoist.H" - -namespace MicVar_FE { - enum { - // independent variables - rho=0, // density - theta, // liquid/ice water potential temperature - tabs, // temperature - pres, // pressure - // non-precipitating vars - qt, // total cloud - qv, // cloud vapor - qc, // cloud water - // derived vars - 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, - 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; - - 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 (const amrex::Real& dt_advance, - const SolverChoice& /*solverChoice*/) 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; } - - int - Qstate_Size () override { return FastEddy::m_qstate_size; } - -private: - // Number of qmoist variables (qt, qv, qc) - int m_qmoist_size = 3; - - // Number of qstate variables - int m_qstate_size = 2; - - // 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; - - // 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; - - // independent variables - amrex::Array mic_fab_vars; -}; -#endif diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp deleted file mode 100644 index 34003313e..000000000 --- a/Source/Microphysics/FastEddy/FastEddy.cpp +++ /dev/null @@ -1,70 +0,0 @@ -#include -#include -#include "FastEddy.H" - -using namespace amrex; - -/** - * Compute Precipitation-related Microphysics quantities. - */ -void FastEddy::AdvanceFE () -{ - auto tabs = mic_fab_vars[MicVar_FE::tabs]; - - // get the temperature, dentisy, theta, qt and qc from input - for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); - auto qt_array = mic_fab_vars[MicVar_FE::qt]->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 pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi); - - const auto& box3d = mfi.tilebox(); - - // Expose for GPU - Real d_fac_cond = m_fac_cond; - - ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - - qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); - - //------- Autoconversion/accretion - Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; - - Real pressure = pres_array(i,j,k); - 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_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 qc_array(i,j,k) > 0.0){ - dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); - } - - qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor; - qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor; - - theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor); - - qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k)); - qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); - - qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k); - }); - } -} diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp deleted file mode 100644 index fa0e8a914..000000000 --- a/Source/Microphysics/FastEddy/Init_FE.cpp +++ /dev/null @@ -1,91 +0,0 @@ -#include -#include "FastEddy.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, - 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::qt, 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 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_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 - 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); - 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), - qv_array(i,j,k)); - pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k))/100.; - }); - } -} diff --git a/Source/Microphysics/FastEddy/Make.package b/Source/Microphysics/FastEddy/Make.package deleted file mode 100644 index 06893d229..000000000 --- a/Source/Microphysics/FastEddy/Make.package +++ /dev/null @@ -1,5 +0,0 @@ -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 deleted file mode 100644 index 735f15ad0..000000000 --- a/Source/Microphysics/FastEddy/Update_FE.cpp +++ /dev/null @@ -1,38 +0,0 @@ -#include "FastEddy.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::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 - 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/Microphysics.H b/Source/Microphysics/Microphysics.H index dc7e722ef..36ade90c5 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -4,7 +4,6 @@ #include #include #include -#include #include class Microphysics { diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index 122de4e5e..baee1b36e 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -137,7 +137,7 @@ void make_buoyancy (Vector& S_data, // ****************************************************************************************** if (solverChoice.moisture_type != MoistureType::None) { - if (solverChoice.moisture_type == MoistureType::FastEddy) + if (solverChoice.moisture_type == MoistureType::Kessler_NoRain) AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); if (solverChoice.moisture_type == MoistureType::SAM) From fad11b616abdcc12bebb59a4fb0d411b238f20df Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 7 Mar 2024 18:37:29 -0800 Subject: [PATCH 2/2] Plot rain accumulation (#1478) Co-authored-by: Mahesh Natarajan Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> Co-authored-by: Ann Almgren --- Exec/SquallLine_2D/inputs_moisture | 2 +- Exec/SquallLine_2D/inputs_moisture_Gabersek | 2 +- Source/IO/Plotfile.cpp | 7 +++++++ Source/Microphysics/Kessler/Init_Kessler.cpp | 2 +- 4 files changed, 10 insertions(+), 3 deletions(-) diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 967c65a20..27acd6c3f 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -38,7 +38,7 @@ amr.check_int = 1000 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # root name of plotfile erf.plot_int_1 = 60 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qv qc qrain pert_dens +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qv qc qrain rain_accum pert_dens # SOLVER CHOICE erf.use_gravity = true diff --git a/Exec/SquallLine_2D/inputs_moisture_Gabersek b/Exec/SquallLine_2D/inputs_moisture_Gabersek index bba8be6de..13a7e2d73 100644 --- a/Exec/SquallLine_2D/inputs_moisture_Gabersek +++ b/Exec/SquallLine_2D/inputs_moisture_Gabersek @@ -50,7 +50,7 @@ amr.check_int = 1000 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # root name of plotfile erf.plot_int_1 = 100 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qv qc qrain pert_dens +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qv qc qrain rain_accum pert_dens # SOLVER CHOICE erf.use_gravity = true diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 2ab8b59de..1a40d414b 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -766,6 +766,13 @@ ERF::WritePlotFile (int which, Vector plot_var_names) MultiFab::Divide(mf[lev], Sm, Rho_comp , mf_comp, 1, 0); mf_comp += 1; } + + if (containerHasElement(plot_var_names, "rain_accum")) + { + MultiFab rain_accum_mf(*(qmoist[lev][0]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],rain_accum_mf,0,mf_comp,1,0); + mf_comp += 1; + } } #ifdef ERF_USE_PARTICLES diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp index ab9084356..553a3cee9 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -29,7 +29,7 @@ void Kessler::Init (const MultiFab& cons_in, m_gtoe = grids; MicVarMap.resize(m_qmoist_size); - MicVarMap = {MicVar_Kess::qt, MicVar_Kess::qv, MicVar_Kess::qcl, MicVar_Kess::qp, MicVar_Kess::rain_accum}; + MicVarMap = {MicVar_Kess::rain_accum, MicVar_Kess::qt, MicVar_Kess::qv, MicVar_Kess::qcl, MicVar_Kess::qp}; // initialize microphysics variables for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) {