Skip to content

Commit

Permalink
Forest canopy (#1980)
Browse files Browse the repository at this point in the history
* 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 <[email protected]>
Co-authored-by: Aaron M. Lattanzi <[email protected]>
  • Loading branch information
3 people authored Nov 25, 2024
1 parent d63c6ae commit dade055
Show file tree
Hide file tree
Showing 17 changed files with 386 additions and 16 deletions.
1 change: 1 addition & 0 deletions CMake/BuildERFExe.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 3 additions & 0 deletions Docs/sphinx_doc/Inputs.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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) | | |
+-------------------------------------+------------------------+-------------------+---------------------+
Expand Down
1 change: 1 addition & 0 deletions Docs/sphinx_doc/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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::
Expand Down
31 changes: 31 additions & 0 deletions Docs/sphinx_doc/theory/Forest.rst
Original file line number Diff line number Diff line change
@@ -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)
<https://doi.org/10.1175/1520-0450(2004)043<0641:AERDLD>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.
4 changes: 4 additions & 0 deletions Exec/ABL/erf_forest_def
Original file line number Diff line number Diff line change
@@ -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
64 changes: 64 additions & 0 deletions Exec/ABL/inputs_canopy
Original file line number Diff line number Diff line change
@@ -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
6 changes: 3 additions & 3 deletions Source/BoundaryConditions/ERF_ABLMost.H
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
3 changes: 3 additions & 0 deletions Source/DataStructs/ERF_DataStruct.H
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 2 additions & 0 deletions Source/ERF.H
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@
#include <ERF_PhysBCFunct.H>
#include <ERF_FillPatcher.H>
#include <ERF_SampleData.H>
#include <ERF_ForestDrag.H>

#ifdef ERF_USE_PARTICLES
#include "ERF_ParticleData.H"
Expand Down Expand Up @@ -1123,6 +1124,7 @@ private:
std::unique_ptr<WriteBndryPlanes> m_w2d = nullptr;
std::unique_ptr<ReadBndryPlanes> m_r2d = nullptr;
std::unique_ptr<ABLMost> m_most = nullptr;
amrex::Vector<std::unique_ptr<ForestDrag>> m_forest;

//
// Holds info for dynamically generated tagging criteria
Expand Down
16 changes: 15 additions & 1 deletion Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -300,6 +304,7 @@ ERF::ERF_shared ()
Hwave_onegrid[lev] = nullptr;
Lwave_onegrid[lev] = nullptr;
}

// Theta prim for MOST
Theta_prim.resize(nlevs_max);

Expand Down Expand Up @@ -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<ForestDrag>(forestfile);
}
}

// AmrMesh iterate on grids?
bool iterate(true);
pp_amr.query("iterate_grids",iterate);
Expand Down
15 changes: 15 additions & 0 deletions Source/ERF_make_new_level.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
// *******************************************************************************************
Expand Down Expand Up @@ -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
// *******************************************************************************************
Expand Down Expand Up @@ -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
// *****************************************************************************************************
Expand Down
40 changes: 40 additions & 0 deletions Source/SourceTerms/ERF_ForestDrag.H
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
#ifndef ERF_FORESTDRAG_H_
#define ERF_FORESTDRAG_H_

#include <memory>
#include <AMReX_MultiFab.H>

/*
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<amrex::Real> m_type_forest;
amrex::Vector<amrex::Real> m_x_forest;
amrex::Vector<amrex::Real> m_y_forest;
amrex::Vector<amrex::Real> m_height_forest;
amrex::Vector<amrex::Real> m_diameter_forest;
amrex::Vector<amrex::Real> m_cd_forest;
amrex::Vector<amrex::Real> m_lai_forest;
amrex::Vector<amrex::Real> m_laimax_forest;
std::unique_ptr<amrex::MultiFab> m_forest_drag;
};
#endif
Loading

0 comments on commit dade055

Please sign in to comment.