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/LinearSolvers.rst b/Docs/sphinx_doc/LinearSolvers.rst index 6b6a01f2a..1d34c2bb8 100644 --- a/Docs/sphinx_doc/LinearSolvers.rst +++ b/Docs/sphinx_doc/LinearSolvers.rst @@ -17,8 +17,9 @@ Multigrid can also be used when the union of grids at a level is not in itself r For simulations using grid stretching in the vertical but flat terrain, we must use the hybrid FFT solver in which we perform 2D transforms only in the lateral directions and couple the solution in the vertical direction with a tridiagonal solve. In both these cases we use a 7-point stencil. + To solve the Poisson equation on terrain-fitted coordinates with general terrain, we rely on the FFT-preconditioned GMRES solver since the stencil effectively has variable coefficients and requires 19 points. - .. note:: - **Currently only doubly periodic lateral boundary conditions are supported by the hybrid FFT, and therefore by the GMRES solver. More general boundary conditions are a work in progress.** +- .. note:: +- **The FFT solver / preconditioner can only be used when the union of grids at a level is itself rectangular. 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/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index b353fa291..4d4b28ea7 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,11 +7,24 @@ .. _sec:Buoyancy: +ERF has several options for how to define the buoyancy force. + Buoyancy ========= -ERF has three options for how to define the buoyancy force. Even in the absence of moisture these -expressions are not equivalent. +When using the anelastic formulation, the expression for buoyancy used by ERF is + +.. math:: + \mathbf{B} = -\rho^\prime \mathbf{g} \approx -\rho_0 \mathbf{g} ( \frac{\theta_v^\prime}{\overline{\theta_0}} + +where we define + +.. math:: + \theta_v = \theta_d (1 - (1 - \frac{R_v}{R_d}) (q_v + q_c) - \frac{R_v}{R_d} q_c) + + +When using the compressible formulation, the following choices are available. + Density of the mixture ----------------------- @@ -52,56 +65,10 @@ for the background state. For eg., a usual scenario is that of a background stat the total background density :math:`\rho_0 = \rho_{d_0}(1 + q_{v_0})`, where :math:`\rho_{d_0}` and :math:`q_{v_0}` are the background dry density and vapor mixing ratio respectively. As a check, we observe that :math:`\rho^\prime_0 = 0`, which means that the background state is not buoyant. -Type 2 ------- - -The second option for the buoyancy force is - -.. math:: - \mathbf{B} = -\rho_0 \mathbf{g} ( 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime - + \frac{T^\prime}{\bar{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ) - -To derive this expression, we define :math:`T_v = T (1 + 0.61 q_v − q_c − q_i - q_p)`, then we can write - -.. math:: - p = \rho (R_d q_d + R_v q_v) T = \rho R_d T (1 + 0.61 q_v − q_c − q_i - q_p ) = \rho R_d T_v - - -Starting from :math:`p = \rho R_d T_v` and neglecting :math:`\frac{p^\prime}{\bar{p}}`, we now write - -.. math:: - \frac{\rho^\prime}{\overline{\rho}} = -\frac{T_v^\prime}{\overline{T_v}} - -and define - -.. math:: - - T_v^\prime = T_v - \overline{T_v} \approx \overline{T} [ 0.61 q_v^\prime - (q_c^\prime + q_i^\prime + q_p^\prime)] + - (T - \overline{T}) [1+ 0.61 \bar{q_v} - \bar{q_c} - \bar{q_i} - \bar{q_p} ] . - -where we have retained only first order terms in perturbational quantities. - -Then - -.. math:: - - \mathbf{B} = \rho^\prime \mathbf{g} = -\overline{\rho} \frac{\overline{T}}{\overline{T_v}} \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime ) + \frac{T^\prime}{\overline{T_v}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ] - -where the overbar represents a horizontal average of the current state and the perturbation is defined relative to that average. - -Again keeping only the first order terms in the mass mixing ratios, we can simplify this to - -.. math:: - \mathbf{B} = \rho^\prime \mathbf{g} = -\rho_0 \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime + q_i^\prime + q_p^\prime - + \frac{T^\prime}{\overline{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ] - -We note that this reduces to Type 3 if the horizontal averages of the moisture terms are all zero. - Type 3 ------ -The third formulation of the buoyancy term assumes that the horizontal averages of the moisture quantities are negligible, -which removes the need to compute horizontal averages of these quantities. This reduces the Type 2 expression to the following: +This formulation of the buoyancy term assumes that the horizontal averages of the moisture quantities are negligible. .. math:: \mathbf{B} = \rho^\prime \mathbf{g} \approx -\rho_0 \mathbf{g} ( \frac{T^\prime}{\overline{T}} @@ -120,7 +87,7 @@ This expression for buoyancy is from `khairoutdinov2003cloud`_ and `bryan2002ben .. math:: \begin{equation} - \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg)\mathbf{g}, + \mathbf{B} = \rho'\mathbf{g} \approx -\rho \mathbf{g} \Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg) \end{equation} The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as 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..41707c247 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -126,14 +126,15 @@ struct SolverChoice { if (moisture_type == MoistureType::Kessler_NoRain || moisture_type == MoistureType::SAM || moisture_type == MoistureType::SAM_NoIce || - moisture_type == MoistureType::SAM_NoPrecip_NoIce) { - buoyancy_type = 1; // asserted in make buoyancy + moisture_type == MoistureType::SAM_NoPrecip_NoIce) + { + buoyancy_type = 1; // uses Rhoprime } else { buoyancy_type = 2; // uses Tprime } } - // Which expression (1,2 or 3) to use for buoyancy + // Which expression (1,2/3 or 4) to use for buoyancy pp.query("buoyancy_type", buoyancy_type); // What type of land surface model to use @@ -196,7 +197,7 @@ struct SolverChoice { if (any_anelastic == 1) { constant_density = true; project_initial_velocity = true; - buoyancy_type = 2; // uses Tprime + buoyancy_type = 3; // (This isn't actually used when anelastic is set) } else { pp.query("project_initial_velocity", project_initial_velocity); @@ -644,5 +645,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 8a89c39bc..07055ed5a 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" @@ -1124,6 +1125,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 11449624e..7b427e1ad 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); @@ -1590,6 +1595,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_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H index fd18827f3..967a954fc 100644 --- a/Source/SourceTerms/ERF_buoyancy_utils.H +++ b/Source/SourceTerms/ERF_buoyancy_utils.H @@ -17,16 +17,11 @@ buoyancy_dry_anelastic (int& i, int& j, int& k, amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp); amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp); - amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi); + amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi); + amrex::Real theta_d0_wface = amrex::Real(0.5) * (th0_arr(i,j,k) + th0_arr(i,j,k-1)); + amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - amrex::Real theta_d_0_lo = th0_arr(i,j,k-1); - amrex::Real theta_d_0_hi = th0_arr(i,j,k ); - - amrex::Real theta_d_0_wface = amrex::Real(0.5) * (theta_d_0_lo + theta_d_0_hi); - - amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - - return (-rho0_wface * grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface); + return (-rho0_wface * grav_gpu * (theta_d_wface - theta_d0_wface) / theta_d0_wface); } AMREX_GPU_DEVICE @@ -51,39 +46,11 @@ buoyancy_moist_anelastic (int& i, int& j, int& k, amrex::Real qt_hi = qv_hi + qc_hi; amrex::Real theta_v_hi = theta_d_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi); - amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi); - - amrex::Real theta_v_0_lo = th0_arr(i,j,k-1) * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo); - amrex::Real theta_v_0_hi = th0_arr(i,j,k ) * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi); - - amrex::Real theta_v_0_wface = amrex::Real(0.5) * (theta_v_0_lo + theta_v_0_hi); - - amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - - return (-rho0_wface * grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface); -} - -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -amrex::Real -buoyancy_dry_default (int& i, int& j, int& k, - amrex::Real const& grav_gpu, - amrex::Real const& rd_over_cp, - const amrex::Array4& r0_arr, - const amrex::Array4& p0_arr, - const amrex::Array4& th0_arr, - const amrex::Array4& cell_data) -{ - amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); - amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); - amrex::Real qplus = (t_hi-t0_hi)/t0_hi; - - amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp); - amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - amrex::Real qminus = (t_lo-t0_lo)/t0_lo; + amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi); + amrex::Real theta_v0_wface = amrex::Real(0.5) * (th0_arr(i,j,k) + th0_arr(i,j,k-1)); + amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus); - return (-r0_q_avg * grav_gpu); + return (-rho0_wface * grav_gpu * (theta_v_wface - theta_v0_wface) / theta_v0_wface); } AMREX_GPU_DEVICE @@ -104,165 +71,138 @@ buoyancy_rhopert (int& i, int& j, int& k, return( grav_gpu * amrex::Real(0.5) * ( rhop_hi + rhop_lo ) ); } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_type2 (int& i, int& j, int& k, - const int& n_qstate, - amrex::Real const& grav_gpu, - amrex::Real* rho_d_ptr, - amrex::Real* theta_d_ptr, - amrex::Real* qv_d_ptr, - amrex::Real* qc_d_ptr, - amrex::Real* qp_d_ptr, - const amrex::Array4& cell_prim, - const amrex::Array4& cell_data) +buoyancy_dry_Tpert (int& i, int& j, int& k, + amrex::Real const& grav_gpu, + amrex::Real const& rd_over_cp, + const amrex::Array4& r0_arr, + const amrex::Array4& p0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& cell_data) { - amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]); - amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]); - - amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), - cell_data(i,j,k ,RhoTheta_comp), - cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp)); - amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), - cell_data(i,j,k-1,RhoTheta_comp), - cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp)); + amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); + amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); - amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; - amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; + amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp); + amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; - amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; + amrex::Real tprime_hi = (t_hi-t0_hi)/t0_hi; + amrex::Real tprime_lo = (t_lo-t0_lo)/t0_lo; - amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; - amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; + amrex::Real tp_avg = amrex::Real(0.5) * (tprime_hi + tprime_lo); + amrex::Real r0_avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - amrex::Real qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) - - ( qc_plus - qc_d_ptr[k] + - qp_plus - qp_d_ptr[k] ) - + (tempp3d-tempp1d)/tempp1d*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qp_d_ptr[k]); + return ( -r0_avg * grav_gpu * tp_avg); +} - amrex::Real qminus = 0.61 * ( qv_minus - qv_d_ptr[k-1] ) - - ( qc_minus - qc_d_ptr[k-1] + - qp_minus - qp_d_ptr[k-1] ) - + (tempm3d-tempm1d)/tempm1d*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]); +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +buoyancy_dry_Thpert (int& i, int& j, int& k, + amrex::Real const& grav_gpu, + const amrex::Array4& r0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& cell_prim) +{ + // + // TODO: we are currently using Theta_prime/Theta_0 to replace T_prime/T_0 - P_prime/P_0 + // but I don't think that is quite right. + // + amrex::Real thetaprime_hi = (cell_prim(i,j,k ,PrimTheta_comp) - th0_arr(i,j,k )) / th0_arr(i,j,k ); + amrex::Real thetaprime_lo = (cell_prim(i,j,k-1,PrimTheta_comp) - th0_arr(i,j,k-1)) / th0_arr(i,j,k-1); - amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus); - amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); + amrex::Real thp_avg = amrex::Real(0.5) * (thetaprime_hi + thetaprime_lo); + amrex::Real r0_avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - return (-qavg * r0avg * grav_gpu); + return ( -r0_avg * grav_gpu * thp_avg); } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_type3 (int& i, int& j, int& k, - const int& n_qstate, - amrex::Real const& grav_gpu, - amrex::Real* rho_d_ptr, - amrex::Real* theta_d_ptr, - amrex::Real* qv_d_ptr, - const amrex::Array4& cell_prim, - const amrex::Array4& cell_data) +buoyancy_moist_Tpert (int& i, int& j, int& k, + const int& n_qstate, + amrex::Real const& grav_gpu, + amrex::Real const& rd_over_cp, + const amrex::Array4& r0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& p0_arr, + const amrex::Array4& cell_prim, + const amrex::Array4& cell_data) { - amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]); - amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]); + // + // Note: this currently assumes the base state qv0 is identically zero + // TODO: generalize this to allow for moist base state + // + amrex::Real qv_hi = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; + amrex::Real qv_lo = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; - amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), - cell_data(i,j,k ,RhoTheta_comp), - cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp)); - amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), - cell_data(i,j,k-1,RhoTheta_comp), - cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp)); + amrex::Real qc_hi = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; + amrex::Real qc_lo = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; - amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; - amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; + amrex::Real qp_hi = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; + amrex::Real qp_lo = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; - amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; - amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; + amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); + amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); - amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; - amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; + amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp), qv_hi); + amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp), qv_lo); - amrex::Real qplus = 0.61 * qv_plus - (qc_plus + qp_plus) + (tempp3d-tempp1d)/tempp1d; + amrex::Real q_hi = 0.61 * qv_hi - (qc_hi + qp_hi) + (t_hi-t0_hi)/t0_hi; + amrex::Real q_lo = 0.61 * qv_lo - (qc_lo + qp_lo) + (t_lo-t0_lo)/t0_lo; - amrex::Real qminus = 0.61 * qv_minus - (qc_minus + qp_minus) + (tempm3d-tempm1d)/tempm1d; + amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo); + amrex::Real r0avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus); - amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); - - return ( -qavg * r0avg * grav_gpu ); + return ( -r0avg * grav_gpu * qavg); } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_type4 (int& i, int& j, int& k, - const int& n_qstate, - amrex::Real const& grav_gpu, - amrex::Real* rho_d_ptr, - amrex::Real* theta_d_ptr, - amrex::Real* qv_d_ptr, - amrex::Real* qc_d_ptr, - amrex::Real* qp_d_ptr, - const amrex::Array4& cell_prim, - const amrex::Array4& cell_data) +buoyancy_moist_Thpert (int& i, int& j, int& k, + const int& n_qstate, + amrex::Real const& grav_gpu, + const amrex::Array4& r0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& cell_prim) { - amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; - amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; + // + // Note: this currently assumes the base state qv0 is identically zero + // TODO: generalize this to allow for moist base state + // + amrex::Real qv_hi = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; + amrex::Real qv_lo = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; - amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; - amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; + amrex::Real qc_hi = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; + amrex::Real qc_lo = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; - amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; - amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; + amrex::Real qp_hi = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; + amrex::Real qp_lo = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; - amrex::Real qplus = amrex::Real(0.61) * ( qv_plus - qv_d_ptr[k] ) - - ( qc_plus - qc_d_ptr[k] + - qp_plus - qp_d_ptr[k] ) - + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ]; + // + // TODO: we are currently using Theta_prime/Theta_0 to replace T_prime/T_0 - P_prime/P_0 + // but I don't think that is quite right. + // + amrex::Real q_hi = amrex::Real(0.61) * qv_hi - (qc_hi + qp_hi) + + (cell_prim(i,j,k ,PrimTheta_comp) - th0_arr(i,j,k )) / th0_arr(i,j,k ); - amrex::Real qminus = amrex::Real(0.61) * ( qv_minus - qv_d_ptr[k-1] ) - - ( qc_minus - qc_d_ptr[k-1] + - qp_minus - qp_d_ptr[k-1] ) - + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1]; + amrex::Real q_lo = amrex::Real(0.61) * qv_lo - (qc_lo + qp_lo) + + (cell_prim(i,j,k-1,PrimTheta_comp) - th0_arr(i,j,k-1)) / th0_arr(i,j,k-1); - amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus); - amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); + amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo); + amrex::Real r0avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - return (-qavg * r0avg * grav_gpu); + return ( -r0avg * grav_gpu * qavg); } // ************************************************************************************** // Routines below this line are not currently used // ************************************************************************************** -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -amrex::Real -buoyancy_dry_Tpert (int& i, int& j, int& k, - amrex::Real const& grav_gpu, - amrex::Real const& rd_over_cp, - const amrex::Array4& r0_arr, - const amrex::Array4& p0_arr, - const amrex::Array4& th0_arr, - const amrex::Array4& cell_data) -{ - amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); - amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); - amrex::Real qplus = (t_hi-t0_hi)/t0_hi; - - amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp); - amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - amrex::Real qminus = (t_lo-t0_lo)/t0_lo; - - amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus); - return (-r0_q_avg * grav_gpu); -} - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -276,14 +216,14 @@ buoyancy_dry_anelastic_T (int& i, int& j, int& k, amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k)); amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp); amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp); - amrex::Real qplus = (t_hi-t0_hi)/t0_hi; + amrex::Real q_hi = (t_hi-t0_hi)/t0_hi; amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)); amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp); amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp); - amrex::Real qminus = (t_lo-t0_lo)/t0_lo; + amrex::Real q_lo = (t_lo-t0_lo)/t0_lo; - amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus); + amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * q_hi + r0_arr(i,j,k-1) * q_lo); return (-r0_q_avg * grav_gpu); } diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index 416c79a71..c42882319 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -53,267 +53,133 @@ void make_buoyancy (Vector& S_data, MultiFab p0 (base_state, make_alias, BaseState::p0_comp , 1); MultiFab th0(base_state, make_alias, BaseState::th0_comp, 1); - if (anelastic == 1) { #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); - // We don't compute a source term for z-momentum on the bottom or top boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); + // We don't compute a source term for z-momentum on the bottom or top boundary + if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); + if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + const Array4 & cell_data = S_data[IntVars::cons].array(mfi); + const Array4 & cell_prim = S_prim.array(mfi); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - // Base state density and pressure - const Array4& r0_arr = r0.const_array(mfi); - const Array4& th0_arr = th0.const_array(mfi); + // Base state density and pressure + const Array4& r0_arr = r0.const_array(mfi); + const Array4& p0_arr = p0.const_array(mfi); + const Array4& th0_arr = th0.const_array(mfi); + + if ( anelastic && (solverChoice.moisture_type == MoistureType::None) ) + { + // ****************************************************************************************** + // Dry anelastic + // ****************************************************************************************** + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // + // Return -rho0 g (thetaprime / theta0) + // + buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,grav_gpu[2], + r0_arr,th0_arr,cell_data); + }); + } + else if ( anelastic && (solverChoice.moisture_type != MoistureType::None) ) + { + // ****************************************************************************************** + // Moist anelastic + // ****************************************************************************************** + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // + // Return -rho0 g (thetaprime / theta0) + // + buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k,grav_gpu[2],rv_over_rd, + r0_arr,th0_arr,cell_data); + }); + } + else if ( !anelastic && (solverChoice.moisture_type == MoistureType::None) ) + { + // ****************************************************************************************** + // Dry compressible + // ****************************************************************************************** + int n_q_dry = 0; + if (solverChoice.buoyancy_type == 1) { - if (solverChoice.moisture_type == MoistureType::None) { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // // Return -rho0 g (thetaprime / theta0) // - buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k, - grav_gpu[2], - r0_arr,th0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2], + r0_arr,cell_data); }); - } else { - // NOTE: For decomposition in the vertical direction, klo may not - // reside in the valid box and this call will yield an out - // of bounds error since it depends upon the surface theta_l + } + else if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 3) + { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // - // Return -rho0 g (thetaprime / theta0) + // Return -rho0 g (Tprime / T0) // - buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k, - grav_gpu[2],rv_over_rd, - r0_arr,th0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_dry_Tpert(i,j,k,grav_gpu[2],rd_over_cp, + r0_arr,p0_arr,th0_arr,cell_data); }); } - } // mfi - } - else - { - // ****************************************************************************************** - // Dry versions of buoyancy expressions (type 1 and type 2/3 -- types 2 and 3 are equivalent) - // ****************************************************************************************** - if (solverChoice.moisture_type == MoistureType::None) - { - int n_q_dry = 0; - if (solverChoice.buoyancy_type == 1) { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - // Base state density - const Array4& r0_arr = r0.const_array(mfi); - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data); - }); - } // mfi - - } - else // (buoyancy_type != 1) + else if (solverChoice.buoyancy_type == 4) { - // We now use the base state rather than planar average because - // 1) we don't want to average over the limited region of the fine level if doing multilevel. - // 2) it's cheaper to use the base state than to compute the horizontal averages - // 3) when running in a smallish domain, the horizontal average may evolve over time, - // which is not necessarily the intended behavior - // -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - // Base state density and pressure - const Array4& r0_arr = r0.const_array(mfi); - const Array4& p0_arr = p0.const_array(mfi); - const Array4& th0_arr = th0.const_array(mfi); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_dry_default(i,j,k, - grav_gpu[2],rd_over_cp, - r0_arr,p0_arr,th0_arr,cell_data); - }); - } // mfi - } // buoyancy_type - } // moisture type - else + // + // Return -rho0 g (Theta_prime / Theta_0) + // + buoyancy_fab(i, j, k) = buoyancy_dry_Thpert(i,j,k,grav_gpu[2], + r0_arr,th0_arr,cell_data); + }); + } // buoyancy_type for dry compressible + } + else // if ( !anelastic && (solverChoice.moisture_type != MoistureType::None) ) { - // ****************************************************************************************** - // Moist versions of buoyancy expressions - // ****************************************************************************************** + // ****************************************************************************************** + // Moist compressible + // ****************************************************************************************** - if ( (solverChoice.moisture_type == MoistureType::Kessler_NoRain) || - (solverChoice.moisture_type == MoistureType::SAM) || - (solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) ) - { - AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); - } - - if (solverChoice.buoyancy_type == 1) { - -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + if ( (solverChoice.moisture_type == MoistureType::Kessler_NoRain) || + (solverChoice.moisture_type == MoistureType::SAM) || + (solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) ) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - // Base state density - const Array4& r0_arr = r0.const_array(mfi); + AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); + } + if (solverChoice.buoyancy_type == 1) + { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate, - grav_gpu[2],r0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate,grav_gpu[2], + r0_arr,cell_data); }); - } // mfi - - } else { - - PlaneAverage state_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane); - PlaneAverage prim_ave(&S_prim , geom, solverChoice.ave_plane); - - // Compute horizontal averages of all components of each field - state_ave.compute_averages(ZDir(), state_ave.field()); - prim_ave.compute_averages(ZDir(), prim_ave.field()); - - int ncell = state_ave.ncell_line(); - - Gpu::HostVector rho_h(ncell), theta_h(ncell); - Gpu::DeviceVector rho_d(ncell), theta_d(ncell); - - state_ave.line_average(Rho_comp, rho_h); - Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - - prim_ave.line_average(PrimTheta_comp, theta_h); - Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); - - Real* rho_d_ptr = rho_d.data(); - Real* theta_d_ptr = theta_d.data(); - - // Average valid moisture vars - Gpu::HostVector qv_h(ncell) , qc_h(ncell) , qp_h(ncell); - Gpu::DeviceVector qv_d(ncell,0.0), qc_d(ncell,0.0), qp_d(ncell,0.0); - if (n_qstate >=1) { - prim_ave.line_average(PrimQ1_comp, qv_h); - Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); - } - if (n_qstate >=2) { - prim_ave.line_average(PrimQ2_comp, qc_h); - Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin()); - } - if (n_qstate >=3) { - prim_ave.line_average(PrimQ3_comp, qp_h); - Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin()); } - Real* qv_d_ptr = qv_d.data(); - Real* qc_d_ptr = qc_d.data(); - Real* qp_d_ptr = qp_d.data(); - - if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { + else if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 3) + { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4 & cell_prim = S_prim.array(mfi); - - // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp) - if (solverChoice.buoyancy_type == 2) { - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_type2(i,j,k,n_qstate,grav_gpu[2], - rho_d_ptr,theta_d_ptr, - qv_d_ptr,qc_d_ptr,qp_d_ptr, - cell_prim,cell_data); - }); - } else { - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_type4(i,j,k,n_qstate,grav_gpu[2], - rho_d_ptr,theta_d_ptr, - qv_d_ptr,qc_d_ptr,qp_d_ptr, - cell_prim,cell_data); - }); - } - } // mfi - - } else if (solverChoice.buoyancy_type == 3) { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + buoyancy_fab(i, j, k) = buoyancy_moist_Tpert(i,j,k,n_qstate,grav_gpu[2],rd_over_cp, + r0_arr,th0_arr,p0_arr, + cell_prim,cell_data); + }); + } + else if (solverChoice.buoyancy_type == 4) + { + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4 & cell_prim = S_prim.array(mfi); - - // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp) - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_type3(i,j,k,n_qstate,grav_gpu[2], - rho_d_ptr,theta_d_ptr,qv_d_ptr, - cell_prim,cell_data); + buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu[2], + r0_arr,th0_arr,cell_prim); }); - } // mfi - } // buoyancy_type - } // not buoyancy_type == 1 - } // has moisture - } // anelastic? + } + } // moist compressible + } // mfi } 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..00951a969 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -15,6 +15,8 @@ // // Define primitive variables for all later RK stages // (We have already done this for the first RK step) + // Note that it is essential this happen before the call to make_mom_sources + // because some of the buoyancy routines use the primitive variables // if (nrk > 0) { int ng_cons = S_data[IntVars::cons].nGrow(); @@ -58,6 +60,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 +129,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 +238,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 +430,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 +447,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, diff --git a/Submodules/AMReX b/Submodules/AMReX index f7f3ee9f5..1c8b54530 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit f7f3ee9f57170e227af1fee1c6b95e5b41a45af8 +Subproject commit 1c8b5453023acb66bf65580b56e2f2f1dfd40c8e