From 3890775670d9d203dd3cea599683e50c3cb53191 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Fri, 22 Nov 2024 15:50:44 -0800 Subject: [PATCH 1/5] compiles but allocation issue, needs docs, needs testing. --- CMake/BuildERFExe.cmake | 1 + Exec/ABL/erf_forest_def | 4 + Exec/ABL/inputs_canopy | 63 +++++++++ Source/BoundaryConditions/ERF_ABLMost.H | 6 +- Source/ERF.H | 2 + Source/ERF.cpp | 16 ++- Source/ERF_make_new_level.cpp | 15 ++ Source/SourceTerms/ERF_ForestDrag.H | 44 ++++++ Source/SourceTerms/ERF_ForestDrag.cpp | 139 +++++++++++++++++++ Source/SourceTerms/ERF_Src_headers.H | 2 + Source/SourceTerms/ERF_make_mom_sources.cpp | 58 +++++++- Source/SourceTerms/Make.package | 2 + Source/TimeIntegration/ERF_TI_slow_rhs_fun.H | 18 ++- 13 files changed, 354 insertions(+), 16 deletions(-) create mode 100644 Exec/ABL/erf_forest_def create mode 100644 Exec/ABL/inputs_canopy create mode 100644 Source/SourceTerms/ERF_ForestDrag.H create mode 100644 Source/SourceTerms/ERF_ForestDrag.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 884e4da91..45254c449 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -172,6 +172,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/SourceTerms/ERF_make_sources.cpp ${SRC_DIR}/SourceTerms/ERF_moist_set_rhs.cpp ${SRC_DIR}/SourceTerms/ERF_NumericalDiffusion.cpp + ${SRC_DIR}/SourceTerms/ERF_ForestDrag.cpp ${SRC_DIR}/TimeIntegration/ERF_ComputeTimestep.cpp ${SRC_DIR}/TimeIntegration/ERF_Advance.cpp ${SRC_DIR}/TimeIntegration/ERF_TimeStep.cpp diff --git a/Exec/ABL/erf_forest_def b/Exec/ABL/erf_forest_def new file mode 100644 index 000000000..5013cadb6 --- /dev/null +++ b/Exec/ABL/erf_forest_def @@ -0,0 +1,4 @@ +1 1024 512 45 200 0.2 6 0.8 +1 1024 1024 35 200 0.2 6 0.8 +1 1024 1224 75 200 0.2 6 0.8 +2 1024 1524 120 200 0.2 10 0.8 \ No newline at end of file diff --git a/Exec/ABL/inputs_canopy b/Exec/ABL/inputs_canopy new file mode 100644 index 000000000..a3b0ab8f4 --- /dev/null +++ b/Exec/ABL/inputs_canopy @@ -0,0 +1,63 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 4000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 2048 2048 1024 +amr.n_cell = 64 64 32 + +geometry.is_periodic = 1 1 0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.1 # fixed time step depending on grid resolution + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 10 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 + +erf.init_type = "uniform" + +erf.forest_file = erf_forest_def + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + +# Higher values of perturbations lead to instability +# Instability seems to be coming from BC +prob.U_0_Pert_Mag = 0.08 +prob.V_0_Pert_Mag = 0.08 # +prob.W_0_Pert_Mag = 0.0 diff --git a/Source/BoundaryConditions/ERF_ABLMost.H b/Source/BoundaryConditions/ERF_ABLMost.H index e821f2e46..177d90c5a 100644 --- a/Source/BoundaryConditions/ERF_ABLMost.H +++ b/Source/BoundaryConditions/ERF_ABLMost.H @@ -431,15 +431,15 @@ public: get_t_surf (const int& lev) { return t_surf[lev].get(); } amrex::Real - get_zref () {return m_ma.get_zref();} + get_zref () { return m_ma.get_zref(); } const amrex::FArrayBox* - get_z0 (const int& lev) {return &z_0[lev];} + get_z0 (const int& lev) { return &z_0[lev]; } bool have_variable_sea_roughness() { return m_var_z0; } const amrex::iMultiFab* - get_lmask(const int& lev) {return m_lmask_lev[lev][0];} + get_lmask(const int& lev) { return m_lmask_lev[lev][0]; } int lmask_min_reduce (amrex::iMultiFab& lmask, const int& nghost) diff --git a/Source/ERF.H b/Source/ERF.H index aa8c7a8fa..80b45a260 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -45,6 +45,7 @@ #include #include #include +#include #ifdef ERF_USE_PARTICLES #include "ERF_ParticleData.H" @@ -1119,6 +1120,7 @@ private: std::unique_ptr m_w2d = nullptr; std::unique_ptr m_r2d = nullptr; std::unique_ptr m_most = nullptr; + amrex::Vector> m_forest; // // Holds info for dynamically generated tagging criteria diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 2d056fbbd..ebcc8850a 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -30,7 +30,7 @@ Real ERF::change_max = 1.1; int ERF::fixed_mri_dt_ratio = 0; // Dictate verbosity in screen output -int ERF::verbose = 0; +int ERF::verbose = 0; int ERF::mg_verbose = 0; bool ERF::use_fft = false; @@ -129,6 +129,10 @@ ERF::ERF_shared () lsm_data.resize(nlevs_max); lsm_flux.resize(nlevs_max); + // NOTE: size canopy model before readparams (if file exists, we construct) + m_forest.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr; } + ReadParameters(); initializeMicrophysics(nlevs_max); @@ -300,6 +304,7 @@ ERF::ERF_shared () Hwave_onegrid[lev] = nullptr; Lwave_onegrid[lev] = nullptr; } + // Theta prim for MOST Theta_prim.resize(nlevs_max); @@ -1581,6 +1586,15 @@ ERF::ReadParameters () pp.query("cf_width", cf_width); pp.query("cf_set_width", cf_set_width); + // Query the canopy model file name + std::string forestfile; + auto do_forest = pp.query("forest_file", forestfile); + if (do_forest) { + for (int lev = 0; lev < max_level; ++lev) { + m_forest[lev] = std::make_unique(forestfile); + } + } + // AmrMesh iterate on grids? bool iterate(true); pp_amr.query("iterate_grids",iterate); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index d2a69f3bd..254262810 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -131,6 +131,11 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, } } + // ******************************************************************************************** + // Build the data structures for canopy model (depends upon z_phys) + // ******************************************************************************************** + m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); + //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -205,6 +210,11 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, } } + // ******************************************************************************************** + // Build the data structures for canopy model (depends upon z_phys) + // ******************************************************************************************** + m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); + //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -329,6 +339,11 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp } } + // ******************************************************************************************** + // Build the data structures for canopy model (depends upon z_phys) + // ******************************************************************************************** + m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); + // ***************************************************************************************************** // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch // ***************************************************************************************************** diff --git a/Source/SourceTerms/ERF_ForestDrag.H b/Source/SourceTerms/ERF_ForestDrag.H new file mode 100644 index 000000000..009cafe0b --- /dev/null +++ b/Source/SourceTerms/ERF_ForestDrag.H @@ -0,0 +1,44 @@ +#ifndef ERF_FORESTDRAG_H_ +#define ERF_FORESTDRAG_H_ + +#include +#include +#include +#include +#include + +/* + ForestDrag flow physics adapted from: + Lalic & Mihailovic (2004) + https://doi.org/10.1175/1520-0450(2004)043<0641:AERDLD>2.0.CO;2 + */ +class ForestDrag +{ +public: + + ForestDrag (std::string forestfile); + + ~ForestDrag () = default; + + void + define_drag_field (amrex::BoxArray ba, + amrex::DistributionMapping dm, + amrex::Geometry& geom, + amrex::MultiFab* z_phys_nd); + + amrex::MultiFab* + get_drag_field () { return m_forest_drag.get(); } + +private: + amrex::Vector m_type_forest; + amrex::Vector m_x_forest; + amrex::Vector m_y_forest; + amrex::Vector m_height_forest; + amrex::Vector m_diameter_forest; + amrex::Vector m_cd_forest; + amrex::Vector m_lai_forest; + amrex::Vector m_laimax_forest; + std::unique_ptr m_forest_drag; + std::unique_ptr iptr; +}; +#endif diff --git a/Source/SourceTerms/ERF_ForestDrag.cpp b/Source/SourceTerms/ERF_ForestDrag.cpp new file mode 100644 index 000000000..bf5fa510c --- /dev/null +++ b/Source/SourceTerms/ERF_ForestDrag.cpp @@ -0,0 +1,139 @@ +#include + +using namespace amrex; + +/* + Constructor to get the forest parameters: + TreeType xc, yc, height, diameter, cd, lai, laimax +*/ +ForestDrag::ForestDrag (std::string forestfile) +{ + std::ifstream file(forestfile, std::ios::in); + if (!file.good()) { + Abort("Cannot find forest file: " + forestfile); + } + // TreeType xc yc height diameter cd lai laimax + Real value1, value2, value3, value4, value5, value6, value7, value8; + while (file >> value1 >> value2 >> value3 >> value4 >> value5 >> value6 >> + value7 >> value8) { + m_type_forest.push_back(value1); + m_x_forest.push_back(value2); + m_y_forest.push_back(value3); + m_height_forest.push_back(value4); + m_diameter_forest.push_back(value5); + m_cd_forest.push_back(value6); + m_lai_forest.push_back(value7); + m_laimax_forest.push_back(value8); + } + file.close(); +} + +void +ForestDrag::define_drag_field (BoxArray ba, + DistributionMapping dm, + Geometry& geom, + MultiFab* z_phys_nd) +{ + // Geometry params + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + + // Allocate the forest drag MF + // NOTE: 1 ghost cell for averaging to faces + //m_forest_drag.reset(); + //m_forest_drag = std::make_unique(ba,dm,1,1); + Print() << "CHECK: " << ba << ' ' + << dm << "\n"; + MultiFab Test; + Test.define(ba,dm,1,1); + iptr = std::make_unique(2); + exit(0); + + m_forest_drag->setVal(0.); + + // Loop over forest types and pre-compute factors + for (unsigned ii = 0; ii < m_x_forest.size(); ++ii) { + + // Expose CPU data for GPU capture + Real af; // Depends upon the type of forest (tf) + Real treeZm = 0.0; // Only for forest type 2 + int tf = int(m_type_forest[ii]); + Real hf = m_height_forest[ii]; + Real xf = m_x_forest[ii]; + Real yf = m_y_forest[ii]; + Real df = m_diameter_forest[ii]; + Real cdf = m_cd_forest[ii]; + Real laif = m_lai_forest[ii]; + Real laimaxf = m_laimax_forest[ii]; + if (tf == 1) { + // Constant factor + af = laif / hf; + } else { + // Discretize integral with 100 points and pre-compute + int nk = 100; + Real ztree = 0; + Real expFun = 0; + Real ratio = 0; + const Real dz = hf / Real(nk); + treeZm = laimaxf * hf; + for (int k(0); k& levelDrag = m_forest_drag->array(mfi); + const Array4& z_nd = (z_phys_nd) ? z_phys_nd->const_array(mfi) : + Array4{}; + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Physical positions of cell-centers + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + Real z = prob_lo[2] + (k + 0.5) * dx[2]; + if (z_nd) { + 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) ); + } + z = std::max(z,0.0); + + // Proximity to the forest + const Real radius = std::sqrt((x - xf) * (x - xf) + + (y - yf) * (y - yf)); + + // Hit for canopy region + Real factor = 1; + if ((z <= hf) && (radius <= (0.5 * df))) { + if (tf == 2) { + Real ratio = (hf - treeZm) / (hf - z); + if (z < treeZm) { + factor = std::pow(ratio, 6.0) * + std::exp(6.0 * (1.0 - ratio)); + } else if (z <= hf) { + factor = std::pow(ratio, 0.5) * + std::exp(0.5 * (1.0 - ratio)); + } + } + levelDrag(i, j, k) = cdf * af * factor; + } + }); + } // mfi + } // ii (forest type) + + // Fillboundary for periodic ghost cell copy + m_forest_drag->FillBoundary(geom.periodicity()); + +} // init_drag_field diff --git a/Source/SourceTerms/ERF_Src_headers.H b/Source/SourceTerms/ERF_Src_headers.H index 667a9dfe7..5ea93913d 100644 --- a/Source/SourceTerms/ERF_Src_headers.H +++ b/Source/SourceTerms/ERF_Src_headers.H @@ -53,10 +53,12 @@ void make_mom_sources (int level, int nrk, std::unique_ptr& z_phys_cc, const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, + const amrex::MultiFab& wvel, amrex::MultiFab& xmom_source, amrex::MultiFab& ymom_source, amrex::MultiFab& zmom_source, const amrex::MultiFab& base_state, + amrex::MultiFab* forest_drag, const amrex::Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& mapfac_m, diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index 4f25d2cf3..b19c258bc 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -45,12 +45,14 @@ void make_mom_sources (int level, const MultiFab & S_prim, std::unique_ptr& z_phys_nd, std::unique_ptr& z_phys_cc, - const MultiFab & /*xvel*/, - const MultiFab & /*yvel*/, + const MultiFab & xvel, + const MultiFab & yvel, + const MultiFab & wvel, MultiFab & xmom_src, MultiFab & ymom_src, MultiFab & zmom_src, const MultiFab& base_state, + MultiFab* forest_drag, const Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& /*mapfac_m*/, @@ -210,8 +212,9 @@ void make_mom_sources (int level, const Array4& rho_v = S_data[IntVars::ymom].array(mfi); const Array4& rho_w = S_data[IntVars::zmom].array(mfi); - //const Array4& u = xvel.array(mfi); - //const Array4& v = yvel.array(mfi); + const Array4& u = xvel.array(mfi); + const Array4& v = yvel.array(mfi); + const Array4& w = wvel.array(mfi); const Array4< Real>& xmom_src_arr = xmom_src.array(mfi); const Array4< Real>& ymom_src_arr = ymom_src.array(mfi); @@ -222,8 +225,13 @@ void make_mom_sources (int level, //const Array4& mf_u = mapfac_u->const_array(mfi); //const Array4& mf_v = mapfac_v->const_array(mfi); - const Array4& z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& z_cc_arr = (use_terrain) ? z_phys_cc->const_array(mfi) : Array4{}; + const Array4& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) : + Array4{}; + + const Array4& z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : + Array4{}; + const Array4& z_cc_arr = (use_terrain) ? z_phys_cc->const_array(mfi) : + Array4{}; // ***************************************************************************** // 2. Add CORIOLIS forcing (this assumes east is +x, north is +y) @@ -479,5 +487,43 @@ void make_mom_sources (int level, xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w); } + // ***************************************************************************** + // 9. Add CANOPY source terms + // ***************************************************************************** + if (forest_drag) { + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real ux = u(i, j, k); + const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k ) + + v(i, j+1, k ) + v(i-1, j+1, k ) ); + const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k ) + + w(i, j , k+1) + w(i-1, j , k+1) ); + const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k)); + xmom_src_arr(i, j, k) -= f_drag * ux * windspeed; + }); + ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k ) + + u(i+1, j , k ) + u(i+1, j-1, k ) ); + const Real uy = v(i, j, k); + const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k ) + + w(i , j , k+1) + w(i , j-1, k+1) ); + const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k)); + ymom_src_arr(i, j, k) -= f_drag * uy * windspeed; + }); + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k ) + + u(i , j , k-1) + u(i+1, j , k-1) ); + const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k ) + + v(i , j , k-1) + v(i , j+1, k-1) ); + const amrex::Real uz = w(i, j, k); + const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1)); + zmom_src_arr(i, j, k) -= f_drag * uz * windspeed; + }); + } } // mfi } diff --git a/Source/SourceTerms/Make.package b/Source/SourceTerms/Make.package index 2ed03c5a6..f5d607302 100644 --- a/Source/SourceTerms/Make.package +++ b/Source/SourceTerms/Make.package @@ -5,6 +5,7 @@ CEXE_sources += ERF_make_sources.cpp CEXE_sources += ERF_ApplySpongeZoneBCs.cpp CEXE_sources += ERF_ApplySpongeZoneBCs_ReadFromFile.cpp CEXE_sources += ERF_NumericalDiffusion.cpp +CEXE_sources += ERF_ForestDrag.cpp ifeq ($(USE_NETCDF),TRUE) CEXE_sources += ERF_moist_set_rhs.cpp @@ -13,3 +14,4 @@ endif CEXE_headers += ERF_NumericalDiffusion.H CEXE_headers += ERF_Src_headers.H CEXE_headers += ERF_buoyancy_utils.H +CEXE_headers += ERF_ForestDrag.H diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index 2a982d95a..53663fcc9 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -58,6 +58,9 @@ dptr_wbar_sub, d_rayleigh_ptrs_at_lev, input_sounding_data, turbPert); + // Canopy data for mom sources + MultiFab* forest_drag = m_forest[level]->get_drag_field(); + // Moving terrain if ( solverChoice.terrain_type == TerrainType::Moving ) { @@ -123,9 +126,9 @@ make_mom_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, z_phys_nd[level], z_phys_cc[level], - xvel_new, yvel_new, + xvel_new, yvel_new, zvel_new, xmom_src, ymom_src, zmom_src, - base_state_new[level], fine_geom, solverChoice, + base_state_new[level], forest_drag, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, @@ -232,9 +235,9 @@ make_mom_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, z_phys_nd[level], z_phys_cc[level], - xvel_new, yvel_new, + xvel_new, yvel_new, zvel_new, xmom_src, ymom_src, zmom_src, - base_state[level], fine_geom, solverChoice, + base_state[level], forest_drag, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, @@ -424,6 +427,9 @@ Real* dptr_u_geos = solverChoice.have_geo_wind_profile ? d_u_geos[level].data(): nullptr; Real* dptr_v_geos = solverChoice.have_geo_wind_profile ? d_v_geos[level].data(): nullptr; + // Canopy data for mom sources + MultiFab* forest_drag = m_forest[level]->get_drag_field(); + make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level], #if defined(ERF_USE_RRTMGP) qheating_rates[level], @@ -437,9 +443,9 @@ int n_qstate = micro->Get_Qstate_Size(); make_mom_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, z_phys_nd[level], z_phys_cc[level], - xvel_new, yvel_new, + xvel_new, yvel_new, zvel_new, xmom_src, ymom_src, zmom_src, - base_state[level], fine_geom, solverChoice, + base_state[level], forest_drag, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, From d6bbb9a23dc128371ec1c000babc89636fa0ec43 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Fri, 22 Nov 2024 16:48:24 -0800 Subject: [PATCH 2/5] This ran. --- Docs/sphinx_doc/Inputs.rst | 3 +++ Source/ERF.cpp | 2 +- Source/SourceTerms/ERF_ForestDrag.H | 10 +++------- Source/SourceTerms/ERF_ForestDrag.cpp | 16 +++++----------- 4 files changed, 12 insertions(+), 19 deletions(-) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 0148228c8..29136d252 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -1096,6 +1096,9 @@ List of Parameters | **erf.input_sounding_file** | Name(s) of the | String(s) | input_sounding_file | | | input sounding file(s) | | | +-------------------------------------+------------------------+-------------------+---------------------+ +| **erf.forest_file** | Name(s) of the | String | None | +| | canopy forest file | | | ++-------------------------------------+------------------------+-------------------+---------------------+ | **erf.input_sounding_time** | Time(s) of the | Real(s) | false | | | input sounding file(s) | | | +-------------------------------------+------------------------+-------------------+---------------------+ diff --git a/Source/ERF.cpp b/Source/ERF.cpp index ebcc8850a..6f88f2f62 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1590,7 +1590,7 @@ ERF::ReadParameters () std::string forestfile; auto do_forest = pp.query("forest_file", forestfile); if (do_forest) { - for (int lev = 0; lev < max_level; ++lev) { + for (int lev = 0; lev <= max_level; ++lev) { m_forest[lev] = std::make_unique(forestfile); } } diff --git a/Source/SourceTerms/ERF_ForestDrag.H b/Source/SourceTerms/ERF_ForestDrag.H index 009cafe0b..a95b24a33 100644 --- a/Source/SourceTerms/ERF_ForestDrag.H +++ b/Source/SourceTerms/ERF_ForestDrag.H @@ -2,9 +2,6 @@ #define ERF_FORESTDRAG_H_ #include -#include -#include -#include #include /* @@ -16,13 +13,13 @@ class ForestDrag { public: - ForestDrag (std::string forestfile); + explicit ForestDrag (std::string forestfile); ~ForestDrag () = default; void - define_drag_field (amrex::BoxArray ba, - amrex::DistributionMapping dm, + define_drag_field (const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, amrex::Geometry& geom, amrex::MultiFab* z_phys_nd); @@ -39,6 +36,5 @@ private: amrex::Vector m_lai_forest; amrex::Vector m_laimax_forest; std::unique_ptr m_forest_drag; - std::unique_ptr iptr; }; #endif diff --git a/Source/SourceTerms/ERF_ForestDrag.cpp b/Source/SourceTerms/ERF_ForestDrag.cpp index bf5fa510c..fb3b35e32 100644 --- a/Source/SourceTerms/ERF_ForestDrag.cpp +++ b/Source/SourceTerms/ERF_ForestDrag.cpp @@ -29,8 +29,8 @@ ForestDrag::ForestDrag (std::string forestfile) } void -ForestDrag::define_drag_field (BoxArray ba, - DistributionMapping dm, +ForestDrag::define_drag_field (const BoxArray& ba, + const DistributionMapping& dm, Geometry& geom, MultiFab* z_phys_nd) { @@ -40,15 +40,8 @@ ForestDrag::define_drag_field (BoxArray ba, // Allocate the forest drag MF // NOTE: 1 ghost cell for averaging to faces - //m_forest_drag.reset(); - //m_forest_drag = std::make_unique(ba,dm,1,1); - Print() << "CHECK: " << ba << ' ' - << dm << "\n"; - MultiFab Test; - Test.define(ba,dm,1,1); - iptr = std::make_unique(2); - exit(0); - + m_forest_drag.reset(); + m_forest_drag = std::make_unique(ba,dm,1,1); m_forest_drag->setVal(0.); // Loop over forest types and pre-compute factors @@ -137,3 +130,4 @@ ForestDrag::define_drag_field (BoxArray ba, m_forest_drag->FillBoundary(geom.periodicity()); } // init_drag_field + From 4f6c76f1a9308a0090f009f94eb789e2ae194388 Mon Sep 17 00:00:00 2001 From: Harish Date: Fri, 22 Nov 2024 22:57:43 -0700 Subject: [PATCH 3/5] Test Case Works --- Exec/ABL/inputs_canopy | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Exec/ABL/inputs_canopy b/Exec/ABL/inputs_canopy index a3b0ab8f4..293b620c4 100644 --- a/Exec/ABL/inputs_canopy +++ b/Exec/ABL/inputs_canopy @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 4000 +max_step = 50 amrex.fpe_trap_invalid = 1 @@ -15,7 +15,8 @@ zlo.type = "SlipWall" zhi.type = "SlipWall" # TIME STEP CONTROL -erf.fixed_dt = 0.1 # fixed time step depending on grid resolution +erf.fixed_dt = 2.0 # fixed time step depending on grid resolution +erf.cfl = 0.5 # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -27,11 +28,11 @@ amr.max_level = 0 # maximum level number allowed # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file -erf.check_int = 100 # number of timesteps between checkpoints +erf.check_int = 50 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 10 # number of timesteps between plotfiles +erf.plot_int_1 = 50 # number of timesteps between plotfiles erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta # SOLVER CHOICE From e62a057588bc5e73577fb9d9513ab797eb8f81fd Mon Sep 17 00:00:00 2001 From: Harish Date: Sat, 23 Nov 2024 11:27:53 -0700 Subject: [PATCH 4/5] Added documentation --- Docs/sphinx_doc/index.rst | 1 + Docs/sphinx_doc/theory/Forest.rst | 31 +++++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+) create mode 100644 Docs/sphinx_doc/theory/Forest.rst diff --git a/Docs/sphinx_doc/index.rst b/Docs/sphinx_doc/index.rst index aefbf3798..7633ba5cf 100644 --- a/Docs/sphinx_doc/index.rst +++ b/Docs/sphinx_doc/index.rst @@ -57,6 +57,7 @@ In addition to this documentation, there is API documentation for ERF generated theory/Forcings.rst Particles.rst theory/WindFarmModels.rst + theory/Forest.rst theory/UnitsAndConstants.rst .. toctree:: diff --git a/Docs/sphinx_doc/theory/Forest.rst b/Docs/sphinx_doc/theory/Forest.rst new file mode 100644 index 000000000..6cece77f2 --- /dev/null +++ b/Docs/sphinx_doc/theory/Forest.rst @@ -0,0 +1,31 @@ +Forest Model +-------------- + +The forest model provides an option to include the drag from forested regions to be included in the momentum equation. The +drag force is calculated as follows: + +.. math:: + + F_i= - C_d L(x,y,z) U_i | U_i | + + +Here :math:`C_d` is the coefficient of drag for the forested region and :math:`L(x,y,z)` is the leaf area density (LAD) for the +forested region. A three-dimensional model for the LAD is usually unavailable and is also cumbersome to use if there are thousands +of trees. Two different models are available as an alternative: + +.. math:: + L=\frac{LAI}{h} + +.. math:: + L(z)=L_m \left(\frac{h - z_m}{h - z}\right)^n exp\left[n \left(1 -\frac{h - z_m}{h - z}\right )\right] + +Here :math:`LAI` is the leaf area index and is available from measurements, :math:`h` is the height of the tree, :math:`z_m` is the location +of the maximum LAD, :math:`L_m` is the maximum value of LAD at :math:`z_m` and :math:`n` is a model constant with values 6 (below :math:`z_m`) and 0.5 +(above :math:`z_m`), respectively. :math:`L_m` is computed by integrating the following equation (see `Lalic and Mihailovic (2004) +2.0.CO;2>`_): + +.. math:: + LAI = \int_{0}^{h} L(z) dz + +The simplified model with uniform LAD is recommended for forested regions with no knowledge of the individual trees. LAI values can be used from +climate model look-up tables for different regions around the world if no local remote sensing data is available. \ No newline at end of file From d32f92f26d17ba8b33d7eacc4068c483a9bc4bb2 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Mon, 25 Nov 2024 10:20:07 -0800 Subject: [PATCH 5/5] Fix when we don't use canopy model. --- Docs/sphinx_doc/index.rst | 2 +- Docs/sphinx_doc/theory/Forest.rst | 22 ++++++++++---------- Exec/ABL/inputs_canopy | 2 +- Source/DataStructs/ERF_DataStruct.H | 3 +++ Source/ERF.cpp | 4 ++-- Source/ERF_make_new_level.cpp | 6 +++--- Source/SourceTerms/ERF_make_mom_sources.cpp | 3 ++- Source/TimeIntegration/ERF_TI_slow_rhs_fun.H | 6 ++++-- 8 files changed, 27 insertions(+), 21 deletions(-) diff --git a/Docs/sphinx_doc/index.rst b/Docs/sphinx_doc/index.rst index 7633ba5cf..cea2b5e37 100644 --- a/Docs/sphinx_doc/index.rst +++ b/Docs/sphinx_doc/index.rst @@ -57,7 +57,7 @@ In addition to this documentation, there is API documentation for ERF generated theory/Forcings.rst Particles.rst theory/WindFarmModels.rst - theory/Forest.rst + theory/Forest.rst theory/UnitsAndConstants.rst .. toctree:: diff --git a/Docs/sphinx_doc/theory/Forest.rst b/Docs/sphinx_doc/theory/Forest.rst index 6cece77f2..1fbbe8b21 100644 --- a/Docs/sphinx_doc/theory/Forest.rst +++ b/Docs/sphinx_doc/theory/Forest.rst @@ -1,31 +1,31 @@ Forest Model -------------- -The forest model provides an option to include the drag from forested regions to be included in the momentum equation. The -drag force is calculated as follows: +The forest model provides an option to include the drag from forested regions to be included in the momentum equation. The +drag force is calculated as follows: .. math:: F_i= - C_d L(x,y,z) U_i | U_i | -Here :math:`C_d` is the coefficient of drag for the forested region and :math:`L(x,y,z)` is the leaf area density (LAD) for the +Here :math:`C_d` is the coefficient of drag for the forested region and :math:`L(x,y,z)` is the leaf area density (LAD) for the forested region. A three-dimensional model for the LAD is usually unavailable and is also cumbersome to use if there are thousands -of trees. Two different models are available as an alternative: +of trees. Two different models are available as an alternative: .. math:: L=\frac{LAI}{h} -.. math:: +.. math:: L(z)=L_m \left(\frac{h - z_m}{h - z}\right)^n exp\left[n \left(1 -\frac{h - z_m}{h - z}\right )\right] -Here :math:`LAI` is the leaf area index and is available from measurements, :math:`h` is the height of the tree, :math:`z_m` is the location -of the maximum LAD, :math:`L_m` is the maximum value of LAD at :math:`z_m` and :math:`n` is a model constant with values 6 (below :math:`z_m`) and 0.5 +Here :math:`LAI` is the leaf area index and is available from measurements, :math:`h` is the height of the tree, :math:`z_m` is the location +of the maximum LAD, :math:`L_m` is the maximum value of LAD at :math:`z_m` and :math:`n` is a model constant with values 6 (below :math:`z_m`) and 0.5 (above :math:`z_m`), respectively. :math:`L_m` is computed by integrating the following equation (see `Lalic and Mihailovic (2004) -2.0.CO;2>`_): +2.0.CO;2>`_): .. math:: - LAI = \int_{0}^{h} L(z) dz + LAI = \int_{0}^{h} L(z) dz -The simplified model with uniform LAD is recommended for forested regions with no knowledge of the individual trees. LAI values can be used from -climate model look-up tables for different regions around the world if no local remote sensing data is available. \ No newline at end of file +The simplified model with uniform LAD is recommended for forested regions with no knowledge of the individual trees. LAI values can be used from +climate model look-up tables for different regions around the world if no local remote sensing data is available. \ No newline at end of file diff --git a/Exec/ABL/inputs_canopy b/Exec/ABL/inputs_canopy index 293b620c4..01bf340c8 100644 --- a/Exec/ABL/inputs_canopy +++ b/Exec/ABL/inputs_canopy @@ -46,7 +46,7 @@ erf.Cs = 0.1 erf.init_type = "uniform" -erf.forest_file = erf_forest_def +#erf.forest_file = erf_forest_def # PROBLEM PARAMETERS prob.rho_0 = 1.0 diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index a2056cbf3..d8f05797b 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -644,5 +644,8 @@ struct SolverChoice { amrex::Real turb_disk_angle = -1.0; amrex::Real windfarm_x_shift = -1.0; amrex::Real windfarm_y_shift = -1.0; + + // Flag for valid canopy model + bool do_forest {false}; }; #endif diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 6f88f2f62..d8fabb58f 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1588,8 +1588,8 @@ ERF::ReadParameters () // Query the canopy model file name std::string forestfile; - auto do_forest = pp.query("forest_file", forestfile); - if (do_forest) { + solverChoice.do_forest = pp.query("forest_file", forestfile); + if (solverChoice.do_forest) { for (int lev = 0; lev <= max_level; ++lev) { m_forest[lev] = std::make_unique(forestfile); } diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 254262810..c7c529e3c 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -134,7 +134,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); + if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics @@ -213,7 +213,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); + if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics @@ -342,7 +342,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); + if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } // ***************************************************************************************************** // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index b19c258bc..50f33f93a 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -89,6 +89,7 @@ void make_mom_sources (int level, // ***************************************************************************** const bool l_use_ndiff = solverChoice.use_NumDiff; const bool use_terrain = solverChoice.terrain_type != TerrainType::None; + const bool l_do_forest = solverChoice.do_forest; // ***************************************************************************** // Data for Coriolis forcing @@ -490,7 +491,7 @@ void make_mom_sources (int level, // ***************************************************************************** // 9. Add CANOPY source terms // ***************************************************************************** - if (forest_drag) { + if (l_do_forest) { ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const Real ux = u(i, j, k); diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index 53663fcc9..499120896 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -59,7 +59,8 @@ input_sounding_data, turbPert); // Canopy data for mom sources - MultiFab* forest_drag = m_forest[level]->get_drag_field(); + MultiFab* forest_drag = nullptr; + if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); } // Moving terrain if ( solverChoice.terrain_type == TerrainType::Moving ) @@ -428,7 +429,8 @@ Real* dptr_v_geos = solverChoice.have_geo_wind_profile ? d_v_geos[level].data(): nullptr; // Canopy data for mom sources - MultiFab* forest_drag = m_forest[level]->get_drag_field(); + MultiFab* forest_drag = nullptr; + if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); } make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level], #if defined(ERF_USE_RRTMGP)