From dade0554988c56d49997927cab6344077b4d8372 Mon Sep 17 00:00:00 2001 From: Harish Date: Mon, 25 Nov 2024 14:05:57 -0700 Subject: [PATCH] Forest canopy (#1980) * compiles but allocation issue, needs docs, needs testing. * This ran. * Test Case Works * Added documentation * Fix when we don't use canopy model. --------- Co-authored-by: AMLattanzi Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- CMake/BuildERFExe.cmake | 1 + Docs/sphinx_doc/Inputs.rst | 3 + Docs/sphinx_doc/index.rst | 1 + Docs/sphinx_doc/theory/Forest.rst | 31 +++++ Exec/ABL/erf_forest_def | 4 + Exec/ABL/inputs_canopy | 64 +++++++++ Source/BoundaryConditions/ERF_ABLMost.H | 6 +- Source/DataStructs/ERF_DataStruct.H | 3 + Source/ERF.H | 2 + Source/ERF.cpp | 16 ++- Source/ERF_make_new_level.cpp | 15 +++ Source/SourceTerms/ERF_ForestDrag.H | 40 ++++++ Source/SourceTerms/ERF_ForestDrag.cpp | 133 +++++++++++++++++++ Source/SourceTerms/ERF_Src_headers.H | 2 + Source/SourceTerms/ERF_make_mom_sources.cpp | 59 +++++++- Source/SourceTerms/Make.package | 2 + Source/TimeIntegration/ERF_TI_slow_rhs_fun.H | 20 ++- 17 files changed, 386 insertions(+), 16 deletions(-) create mode 100644 Docs/sphinx_doc/theory/Forest.rst 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 bcb15eff8..c7850f48d 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -169,6 +169,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/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/Docs/sphinx_doc/index.rst b/Docs/sphinx_doc/index.rst index aefbf3798..cea2b5e37 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..1fbbe8b21 --- /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 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..01bf340c8 --- /dev/null +++ b/Exec/ABL/inputs_canopy @@ -0,0 +1,64 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 50 + +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 = 2.0 # fixed time step depending on grid resolution +erf.cfl = 0.5 + +# 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 = 50 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +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 +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/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.H b/Source/ERF.H index 777625041..e8fcc4a2a 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" @@ -1123,6 +1124,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 573d8c7bf..21fbb1ff9 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); @@ -1588,6 +1593,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; + 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); + } + } + // 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..c7c529e3c 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) + // ******************************************************************************************** + if (solverChoice.do_forest) { 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) + // ******************************************************************************************** + if (solverChoice.do_forest) { 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) + // ******************************************************************************************** + 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_ForestDrag.H b/Source/SourceTerms/ERF_ForestDrag.H new file mode 100644 index 000000000..a95b24a33 --- /dev/null +++ b/Source/SourceTerms/ERF_ForestDrag.H @@ -0,0 +1,40 @@ +#ifndef ERF_FORESTDRAG_H_ +#define ERF_FORESTDRAG_H_ + +#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: + + explicit ForestDrag (std::string forestfile); + + ~ForestDrag () = default; + + void + define_drag_field (const amrex::BoxArray& ba, + const 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; +}; +#endif diff --git a/Source/SourceTerms/ERF_ForestDrag.cpp b/Source/SourceTerms/ERF_ForestDrag.cpp new file mode 100644 index 000000000..fb3b35e32 --- /dev/null +++ b/Source/SourceTerms/ERF_ForestDrag.cpp @@ -0,0 +1,133 @@ +#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 (const BoxArray& ba, + const 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); + 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..50f33f93a 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*/, @@ -87,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 @@ -210,8 +213,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 +226,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 +488,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 (l_do_forest) { + 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..499120896 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -58,6 +58,10 @@ dptr_wbar_sub, d_rayleigh_ptrs_at_lev, input_sounding_data, turbPert); + // Canopy data for mom sources + MultiFab* forest_drag = nullptr; + if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); } + // Moving terrain if ( solverChoice.terrain_type == TerrainType::Moving ) { @@ -123,9 +127,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 +236,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 +428,10 @@ 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 = 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) qheating_rates[level], @@ -437,9 +445,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,