Skip to content

Commit

Permalink
update docs for boundary conditions (#1306)
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren authored Nov 18, 2023
1 parent 3b5d301 commit 98734ce
Show file tree
Hide file tree
Showing 5 changed files with 267 additions and 266 deletions.
261 changes: 39 additions & 222 deletions Docs/sphinx_doc/BoundaryConditions.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,40 +2,9 @@
.. role:: cpp(code)
:language: c++

.. _sec:BoundaryConditions:
.. _sec:LateralBoundaryConditions:

Filling Ghost Values
--------------------------
ERF uses an operation called ``FillPatch`` to fill the ghost cells/faces for each grid of data.
The data is filled outside the valid region with a combination of three operations: interpolation
from coarser level, copy from same level, and enforcement of physical boundary conditions.

Interpolation from Coarser level
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Interpolation is controlled by which interpolater we choose to use. The default is
conservative interpolation for cell-centered quantities, and analogous for faces.
The paradigm is that fine faces on a coarse-fine boundary are filled as Dirichlet
boundary conditions from the coarser level; all faces outside the valid region are
similarly filled, while fine faces inside the valid region are not over-written.

Copy from other grids at same level (includes periodic boundaries)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This is part of the ``FillPatch`` operation, but can also be applied independently,
e.g. by the call

::

mf.FillBoundary(geom[lev].periodicity());

would fill all the ghost cells/faces of the grids in MultiFab ``mf``, including those
that occur at periodic boundaries.

In the ``FillPatch`` operation, ``FillBoundary`` always overrides any interpolated values, i.e. if
there is fine data available (except at coarse-fine boundary) we always use it.

Imposition of physical/domain boundary conditions
Ideal Domain Boundary Conditions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

There are two primary types of physical/domain boundary conditions: those which rely only on the
Expand Down Expand Up @@ -71,9 +40,6 @@ for each type; this is summarized in the table below.

.. _sec:dirichlet:

Boundary Conditions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

ERF provides the ability to specify a variety of boundary conditions (BCs) in the inputs file.
We use the following options preceded by ``xlo``, ``xhi``, ``ylo``, ``yhi``, ``zlo``, and ``zhi``:

Expand Down Expand Up @@ -164,211 +130,62 @@ Couette regression test example, in which we specify
We also note that in the case of a "slipwall" boundary condition in a simulation with non-zero
viscosity specified, the "foextrap" boundary condition enforces zero strain at the wall.

The keywork "MOST" is an ERF-specific boundary type and the mapping is described below.

The keywork "MOST" is an ERF-specific boundary type; see :ref:`sec:MOST` for more information.

It is important to note that external Dirichlet boundary data should be specified
as the value on the face of the cell bounding the domain, even for cell-centered
state data.

More general boundary types are a WIP; one type that will be supported soon is the ability
to read in a time sequence of data at a domain boundary and impose this data as "ext_dir"
boundary values using ``FillPatch``.

.. _MostBoundary:

MOST Boundaries
-------------------
Monin-Obukhov similarity theory (MOST) is used to describe the atmospheric surface layer (ASL), the lowest part of the atmospheric boundary layer. The implementation of MOST in ERF follows that in `AMR-Wind <https://github.com/Exawind/amr-wind/>`_, which is based on the surface layer profiles presented in
`P. van der Laan, et al., Wind Energy, 2017 <https://onlinelibrary.wiley.com/doi/10.1002/we.2017>`_ and
`D. Etling, "Modeling the vertical ABL structure", 1999 <https://www.worldscientific.com/doi/abs/10.1142/9789814447164_0003>`_.
MOST theory assumes that the ASL is in a steady state and horizontally homogenous, and kinematic fluxes due to turbulent transport (:math:`\overline{u^{'}w^{'}}`, :math:`\overline{v^{'}w^{'}}`, and :math:`\overline{\theta^{'}w^{'}}`) are constant with height.
:math:`\Phi_m` and :math:`\Phi_h` are the nondimensional wind shear and temperature gradient, respectively, which are assumed to follow universal similarity laws based on dimensional arguments.
With these assumptions, the MOST theory can be written as:

.. math::
\overline{u^{'}} \overline{w^{'}} = const = -u^{2}_{\star},
\overline{w^{'}} \overline{\theta^{'}} = const = -u_{\star}\theta_{\star},
\Phi_{m}(\zeta) = \frac{\kappa z}{u_{\star}} \frac{\partial \overline{u}(z)}{\partial z},
\Phi_{h}(\zeta) = \frac{\kappa z}{u_{\star}} \frac{\partial \overline{\theta}(z)}{\partial z}
where the nondimensional gradients are expressed in terms of the MOST stability parameter, :math:`\zeta = \frac{z}{L} = -\frac{\kappa z}{u_{\star}^{3}} \frac{g}{\overline{\theta}} \overline{w^{'}\theta^{'}}`, which serves as a surface layer scaling parameter.
Here, :math:`L` is the Monin-Obukhov length,
:math:`u_{\star}` is the friction velocity (defined for :math:`u` aligned with the wind direction),
:math:`\theta_{\star}` is the surface layer temperature scale,
:math:`\overline{\theta}` is the reference virtual potential temperature for the ASL,
and :math:`\kappa` is the von Karman constant (taken to be :math:`0.41`).

Integration of the MOST assumption equations give the classical MOST profiles of mean velocity and potential temperature

.. math::
\overline{u}(z) &= \frac{u_{\star}}{\kappa} \left[ \mathrm{ln} \left(\frac{z}{z_0}\right) - \Psi_m(\zeta)\right],
\overline{\theta}(z) - \theta_0 &= \frac{\theta_{\star}}{\kappa} \left[ \mathrm{ln}\left(\frac{z}{z_0}\right) - \Psi_{h}(\zeta) \right]
where :math:`\theta_0` is the surface potential temperature and :math:`z_0` is a characteristic roughness height. The integrated similarity functions,

.. math::
\Psi_{m}(\zeta) &= \int_{0} ^{\frac{z}{L}} [1-\Phi_{m}(\zeta)]\frac{d\zeta}{\zeta},
\Psi_{h}(\zeta) &= \int_{0} ^{\frac{z}{L}} [1-\Phi_{h}(\zeta)]\frac{d\zeta}{\zeta}
are calculated analytically from empirical gradient functions :math:`\Phi_m` and :math:`\Phi_h`, which are
defined piecewise for stable and unstable values of the stability parameter.

Unstable: :math:`(\zeta < 0)`

.. math::
\Phi_{m} &= (1-\gamma_{1}\zeta)^{-\frac{1}{4}}, \quad
\Psi_{m} = \mathrm{ln}\left[\frac{1}{8}(1+\Phi_{m}^{-2})(1+\Phi_{m}^{-1})^{2}\right]-2\arctan(\Phi_{m}^{-1})+\frac{\pi}{2},
\Phi_{h} &= \sigma_{\theta}(1-\gamma_{2}\zeta)^{-\frac{1}{2}}, \quad
\Psi_{h} = (1+\sigma_{\theta}) \mathrm{ln} \left[\frac{1}{2}(1+\Phi_{h}^{-1}) \right]+(1-\sigma_{\theta}) {\mathrm{ln}} \left[\frac{1}{2}(-1+\Phi_{h}^{-1})\right]
Stable: :math:`(\zeta > 0)`

.. math::
\Phi_{m} &= 1+\beta \zeta, \quad \Psi_{m}=-\beta \zeta,
\Phi_{h} &= \sigma_{\theta}+\beta \zeta, \quad \Psi_{h}=(1-\sigma_{\theta})\mathrm{ln}(\zeta)-\beta \zeta,
where the constants take the values proposed in `Dyer, Boundary Layer Meteorology, 1974
<https://link.springer.com/article/10.1007/BF00240838>`_:

.. math::
\sigma_{\theta}=1, \quad \beta = 5, \quad \gamma_{1}=16, \quad \gamma_{2}=16
Inverting the equations above, the MOST stability parameter,

.. math::
\zeta=\frac{z}{L} = -\kappa z \frac{g}{\bar{\theta}} \frac{\theta_{\star}}{u^{2}_{\star}}
is determined by the friction velocity

.. math::
u_{\star} = \kappa \overline{u}/[\mathrm{ln}(z/z_0)-\Psi_{m}({z}/{L})]
and the characteristic surface layer temperature

.. math::
\theta_{\star} = \kappa (\overline{\theta}-\theta_0)/[\mathrm{ln}(z / z_0)-\Psi_{h}(z/L)]
MOST Implementation
~~~~~~~~~~~~~~~~~~~

In ERF, when the MOST boundary condition is applied, velocity and temperature in the ghost cells are set to give stresses that are consistent with the MOST equations laid out above. The code is structured to allow either the surface temperature (:math:`\theta_0`) or surface temperature flux (:math:`\overline{w^{'}\theta^{'}}`) to be enforced. To apply the MOST boundary, the following algorithm is applied:

#. Horizontal (planar) averages :math:`\bar{u}`, :math:`\bar{v}` and :math:`\overline{\theta}` are computed at a reference height :math:`z_{ref}` assumed to be within the surface layer.

#. Initially, neutral conditions (:math:`L=\infty, \zeta=0`) are assumed and used to compute a provisional :math:`u_{\star}` using the equation given above. If :math:`\theta_0` is specified, the above equation for :math:`\theta_{\star}` is applied and then the surface flux is computed :math:`\overline{w^{'}\theta^{'}} = -u_{\star} \theta_{\star}`. If :math:`\overline{w^{'}\theta^{'}}` is specified, :math:`\theta_{\star}` is computed as :math:`-\overline{w^{'}\theta^{'}}/u_{\star}` and the previous equation is inverted to compute :math:`\theta_0`.

#. The stability parameter :math:`\zeta` is recomputed using the equation given above based on the provisional values of :math:`u_{\star}` and :math:`\theta_{\star}`.

#. The previous two steps are repeated iteratively, sequentially updating the values of :math:`u_{\star}` and :math:`\zeta`, until the change in the value of :math:`u_{\star}` on each iteration falls below a specified tolerance.

#. Once the MOST iterations have converged, and the planar average surface flux values are known, the approach from `Moeng, Journal of the Atmospheric Sciences, 1984 <https://journals.ametsoc.org/view/journals/atsc/41/13/1520-0469_1984_041_2052_alesmf_2_0_co_2.xml>`_ is applied to consistently compute local surface-normal stress/flux values (e.g., :math:`\tau_{xz} = - \rho \overline{u^{'}w^{'}}`):

.. math::
\left. \frac{\tau_{xz}}{\rho} \right|_0 &= u_{\star}^{2} \frac{(u - \bar{u})|\mathbf{\bar{u}}| + \bar{u}\sqrt{u^2 + v^2} }{|\mathbf{\bar{u}}|^2},
\left. \frac{\tau_{yz}}{\rho} \right|_0 &= u_{\star}^{2} \frac{(v - \bar{v})|\mathbf{\bar{u}}| + \bar{v}\sqrt{u^2 + v^2} }{|\mathbf{\bar{u}}|^2},
\left. \frac{\tau_{\theta z}}{\rho} \right|_0 &= \theta_\star u_{\star} \frac{|\mathbf{\bar{u}}| ({\theta} - \overline{\theta}) +
\sqrt{u^2+v^2} (\overline{\theta} - \theta_0) }{ |\mathbf{\bar{u}}| (\overline{\theta} -\theta_0) } =
u_{\star} \kappa \frac{|\mathbf{\bar{u}}| ({\theta} - \overline{\theta}) +
\sqrt{u^2+v^2} (\overline{\theta} - \theta_0) }{ |\mathbf{\bar{u}}| [ \mathrm{ln}(z_{ref} / z_0)-\Psi_{h}(z_{ref}/L)] }
where :math:`\bar{u}`, :math:`\bar{v}` and :math:`\overline{\theta}` are the plane averaged values (at :math:`z_{ref}`) of the
two horizontal velocity components and the potential temperature, respectively, and
:math:`|\mathbf{\bar{u}}|` is the plane averaged magnitude of horizontal velocity (plane averaged wind speed). We note a slight variation in the denominator
of the velocity terms from the form of the
equations presented in Moeng to match the form implemented in AMR-Wind.

#. These local flux values are used to populate values in the ghost cells that will lead to appropiate fluxes, assuming the fluxes are computed from the turbulent transport coefficients (in the vertical direction, if applicable) :math:`K_{m,v}` and :math:`K_{\theta,v}` as follows:

.. math::
\tau_{xz} = K_{m,v} \frac{\partial u}{\partial z}
\tau_{yz} = K_{m,v} \frac{\partial v}{\partial z}
\tau_{\theta z} = K_{\theta,v} \frac{\partial \theta}{\partial z}.
This implies that, for example, the value set for the conserved :math:`\rho\theta` variable in the :math:`-n\mathrm{th}` ghost cell is

.. math::
(\rho \theta)_{i,j,-n} = \rho_{i,j,-n} \left[ \frac{(\rho\theta)_{i,j,0}}{\rho_{i,j,0}} - \left. \frac{\tau_{\theta z}}{\rho} \right|_{i,j,0} \frac{\rho_{i,j,0}}{K_{\theta,v,(i,j,0)}} n \Delta z \right]
MOST Inputs
~~~~~~~~~~~~~~~~~~~
To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports three methods for parameterizing the surface roughness: ``constant``, ``charnock``, and ``modified_charnock``. The latter two methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 <https://doi.org/10.1175/JAMC-D-17-0137.1>`_. The rougness calculation method may be specified with

::

erf.most.roughness_type = STRING #Z_0 type (constant, charnock, modified_charnock)

If the ``charnock`` method is employed, the :math:`a` constant may be specified with ``erf.most.charnock_constant`` (defaults to 0.0185). If the ``modified_charnock`` method is employed, the depth :math:`d` may be specified with ``erf.most.modified_charnock_depth`` (defaults to 30 m).

When computing an average :math:`\overline{\phi}` for the MOST boundary, where :math:`\phi` denotes a generic variable, ERF supports a variety of approaches. Specifically, ``planar averages`` and ``local region averages`` may be computed with or without ``time averaging``. With each averaging methodology, the query point :math:`z` may be determined from the following procedures: specified vertical distance :math:`z_{ref}` from the bottom surface, specified :math:`k_{index}`, or (when employing terrain-fit coordinates) specified normal vector length :math:`z_{ref}`. The available inputs to the MOST boundary and their associated data types are

::
Specified Domain Boundary Conditions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

erf.most.average_policy = INT #POLICY FOR AVERAGING
erf.most.use_normal_vector = BOOL #USE NORMAL VECTOR W/ TERRAIN?
erf.most.use_interpolation = BOOL #INTERPOLATE QUERY POINT W/ TERRAIN?
erf.most.time_average = BOOL #USE TIME AVERAGING?
erf.most.z0 = FLOAT #SURFACE ROUGHNESS
erf.most.zref = FLOAT #QUERY DISTANCE (HEIGHT OR NORM LENGTH)
erf.most.surf_temp = FLOAT #SPECIFIED SURFACE TEMP
erf.most.surf_temp_flux = FLOAT #SPECIFIED SURFACE FLUX
erf.most.k_arr_in = INT #SPECIFIED K INDEX ARRAY (MAXLEV)
erf.most.radius = INT #SPECIFIED REGION RADIUS
erf.most.time_window = FLOAT #WINDOW FOR TIME AVG
When we use specified lateral boundary conditions, we read time-dependent Dirichlet data
from a file. The user may specify (in the inputs file)
the total width of the interior Dirichlet and relaxation region with
``erf.wrfbdy_width = <Int>`` (yellow + blue)
and analogously the width of the interior Dirichlet region may be specified with
``erf.wrfbdy_set_width = <Int>`` (yellow).
We note that all dycore variables are set and relaxed, but moisture and other scalars
are only set in the yellow region if present in the boundary file.

We now consider two concrete examples. To employ an instantaneous ``planar average`` at a specified vertical height above the bottom surface, one would specify:
.. |wrfbdy| image:: figures/wrfbdy_BCs.png
:width: 600

::
.. _fig:Lateral BCs
erf.most.average_policy = 0
erf.most.use_normal_vector = false
erf.most.time_average = false
erf.most.z0 = 0.1
erf.most.zref = 1.0
.. table:: Lateral boundaries

By contrast, ``local region averaging`` would be employed in conjunction with ``time averaging`` for the following inputs:
+-----------------------------------------------------+
| |wrfbdy| |
+-----------------------------------------------------+
| Image taken from `Skamarock et al. (2021)`_ |
+-----------------------------------------------------+

::
.. _`Skamarock et al. (2021)`: http://dx.doi.org/10.5065/1dfh-6p97

erf.most.average_policy = 1
erf.most.use_normal_vector = true
erf.most.use_interpolation = true
erf.most.time_average = true
erf.most.z0 = 0.1
erf.most.zref = 1.0
erf.most.surf_temp_flux = 0.0
erf.most.radius = 1
erf.most.time_window = 10.0
Here we describe the relaxation algorithm.

In the above case, ``use_normal_vector`` utilizes the a local surface-normal vector with length :math:`z_{ref}` to construct the positions of the query points. Each query point, and surrounding points that are within ``erf.most.radius`` from the query point, are interpolated to and averaged; for a radius of 1, 27 points are averaged. The ``time average`` is completed by way of an exponential filter function whose peak coincides with the current time step and tail extends backwards in time
Within the interior Dirichlet region (yellow), the RHS is exactly 0.
However, within the relaxation region (blue), the RHS (:math:`F`) is given by the following:

.. math::
\frac{1}{\tau} \int_{-\infty}^{0} \exp{\left(t/\tau\right)} \, f(t) \; \rm{d}t.
\begin{align}
F &= G + R, \\
\psi^{\prime} &= \psi^{n} + \Delta t \; G, \\
R &= H_{1} \left( \psi^{FP} - \psi^{\prime} \right) - H_{2} \Delta^2 \left( \psi^{FP} - \psi^{\prime} \right), \\
H_{1} &= \frac{1}{10 \Delta t} \frac{{\rm SpecWidth} + {\rm RelaxWidth} - n}{{\rm RelaxWidth} - 1}, \\
H_{2} &= \frac{1}{50 \Delta t} \frac{{\rm SpecWidth} + {\rm RelaxWidth} - n}{{\rm RelaxWidth} - 1},
\end{align}
Due to the form of the above integral, it is advantageous to consider :math:`\tau` as a multiple of the simulation time step :math:`\Delta t`, which is specified by ``erf.most.time_window``. As ``erf.most.time_window`` is reduced to 0, the exponential filter function tends to a Dirac delta function (prior averages are irrelevant). Increasing ``erf.most.time_window`` extends the tail of the exponential and more heavily weights prior averages.
where :math:`G` is the RHS of the evolution equations, :math:`\psi^{\prime}` is the predicted update without
relaxation, :math:`\psi^{FP}` is the fine data obtained from spatial and temporal interpolation of the
coarse data, and :math:`n` is the minimum number of grid points from a lateral boundary. The specified and
relaxation regions are applied to all dycore variables :math:`\left[\rho \; \rho\Theta \; U\; V\; W \right]`
on the fine mesh.

Sponge zone boundary conditions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

ERF provides the capability to apply sponge zones at the boundaries to prevent spurious reflections that otherwise occur at the domain boundaries if standard extrapolation boundary condition is used. The sponge zone is implemented as a source term in the governing equations, which are active in a volumteric region at the boundaries that is specified by the user in the inputs file. Currently the target condition to which the sponge zones should be forced towards is to be specifed by the user in the inputs file.

Expand Down
Loading

0 comments on commit 98734ce

Please sign in to comment.