Commit 764c3045 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch '142-add-parameterization-for-the-effective-hydrodynamic-dispersion-tensor'

Resolve "Add parameterization for the effective hydrodynamic dispersion tensor"

Closes #142

See merge request !141
parents 114eb2fd ee58535b
......@@ -40,6 +40,7 @@
* Coupling between transient models for water flow and solute transport !96
* Initial conditions generated from H5 input data !130
* Generic Python VTK file reader !143
* Parameterizations for hydrodynamic dispersion in solute transport !141
### Changed
* `Simulation` is renamed `RichardsSimulation` and moved to
......
......@@ -197,7 +197,8 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
to the `CMAKE_FLAGS` in the above command, replacing `<hdf5-path>` with the
path chosen as installation prefix when configuring HDF5.
#### Experimental Features
### Experimental Features
The local operator implementing Richards equation's discretization supports
multiple scheme settings. Setting these via the config file is disabled by
default. You can enable this feature by reconfiguring DORiE with the CMake flag
......
......@@ -119,7 +119,7 @@ endfunction()
# copy BC & parameter files
file(COPY 2d_infiltr.bcdat 3d_infiltr.bcdat
2d_solute.bcdat 3d_solute.bcdat
param.yml
richards_param.yml transport_param.yml
DESTINATION .)
# Random field generator
......
......@@ -32,10 +32,15 @@ adding an empty line, make text **bold** or ``monospaced``.
<dorie>
<category name="simulation">
<parameter name="mode">
<definition> Mode of simulation in DORiE. </definition>
<definition> Sets the simulation mode of DORiE.
(``richards``) mode calculates the matric head with a
DG scheme and produce water fluxes with Raviart Thomas reconstruction.
(``richards+transport``) mode calculates (``richards``) mode and use
the generated water flux and saturation at each step to calculate the
solute transport model for unsaturated media with a DG scheme.
</definition>
<suggestion> richards </suggestion>
<values> richards, richards+transport </values>
<comment> richards, richards+transport </comment>
</parameter>
</category>
......
volumes:
sand:
index: 0
type: MvG
parameters:
alpha: -2.3
n: 4.17
k0: 2.2E-5
theta_r: 0.03
theta_s: 0.31
tau: -1.1
silt:
index: 1
type: MvG
parameters:
alpha: -0.7
n: 1.3
k0: 1.0E-5
theta_r: 0.01
theta_s: 0.41
tau: 0.0
scaling:
type: None
......@@ -227,6 +227,7 @@ adding an empty line, make text **bold** or ``monospaced``.
<definition> YAML file containing the parameter definitions.
</definition>
<values> path </values>
<suggestion> richards_param.yml </suggestion>
</parameter>
</category>
......
volumes:
sand:
index: 0
richards:
type: MvG
parameters:
alpha: -2.3
n: 4.17
k0: 2.2E-5
theta_r: 0.03
theta_s: 0.31
tau: -1.1
silt:
index: 1
richards:
type: MvG
parameters:
alpha: -0.7
n: 1.3
k0: 1.0E-5
theta_r: 0.01
theta_s: 0.41
tau: 0.0
scaling:
type: None
......@@ -67,6 +67,16 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> endOfTransportStep, endOfRichardsStep, none </values>
</parameter>
<parameter name="writeDispersionTensor">
<definition> Defines whether VTK files should include the hydrodynamic
dispersion tensor. Tensors are written in 3D and have 9 componentents
independently of the world dimension. This can be easily be visualizated
in Paraview with the ``Tensor Glyph`` filter.
</definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="subsamplingLevel">
<definition> Plot VTK files with virtually refined grids. VTK only
supports bilinear triangulations and displays higher-order solutions
......@@ -109,6 +119,13 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="dirichletMode">
<definition> Type of the input value for the dirichlet condition.
</definition>
<values> soluteConcentration, totalSolute </values>
<suggestion> soluteConcentration </suggestion>
</parameter>
</category>
<category name="initial">
......@@ -222,11 +239,12 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
</category>
<category name="parameters">
<parameter name="effHydrDips">
<definition> Effective hydrodynamic dispersion </definition>
<values> float &gt; 0 </values>
<suggestion> 2E-9 </suggestion>
<category name="parameters">
<parameter name="file">
<definition> YAML file containing the parameter definitions.
</definition>
<values> path </values>
<suggestion> transport_param.yml </suggestion>
</parameter>
</category>
......@@ -253,8 +271,7 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
<parameter name="FEorder">
<definition> Order of the finite element method used. Values &gt; 1 are not
thoroughly tested. </definition>
<definition> Order of the finite element method used. </definition>
<suggestion> 0 </suggestion>
<values> 0 </values>
</parameter>
......
volumes:
sand:
index: 0
transport:
type: Deff_MQ1 + Dhm_iso
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
silt:
index: 1
transport:
type: Deff_MQ1 + Dhm_iso
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
\ No newline at end of file
......@@ -64,11 +64,16 @@ Output Files
.. object:: Transport
- ``solute``: Solute concentration in water phase
:math:`c_w \, [\mathrm{kg}/\mathrm{m}^3]`
:math:`c_w \, [\mathrm{kg}\mathrm{m}^{-3}]`
- ``solute_total``: Total solute concentration
:math:`c_t \, [\mathrm{kg}/\mathrm{m}^3]`
:math:`c_t \, [\mathrm{kg}\mathrm{m}^{-3}]`
- ``micro_peclet``: Microscopic peclet number
:math:`\mathsf{pe_\mu} \, [-]`
- ``flux_RT<k-1>``: Reconstructed solute flux
:math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}/\mathrm{m}^{-2}\mathrm{s}^{-1}]`
:math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}\,\mathrm{m}^{-2}\mathrm{s}^{-1}]`
(if enabled!)
- ``eff_hd_dispersion``: Hydrodynamic dispersion tensor
:math:`\mathsf{D}_\mathsf{hd} \, [\mathrm{m}^2\mathrm{s}^{-1}]`
(if enabled!)
* VTK Parallel Collection Files (``.pvtu``): Merging multiple VTK files in case
......
......@@ -7,7 +7,7 @@ Cook Book
Introduction
============
This part of the documentation is intended for first-time DORiE users. It explaines the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement.
This part of the documentation is intended for first-time DORiE users. It explains the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement.
Prerequisites
-------------
......@@ -57,7 +57,7 @@ Let's begin with the ``output``. We want the program to give us at least some ou
.. code-block:: none
[output]
verbose = 1
logLevel = info
outputPath = ./sand
fileName = sand
......
......@@ -4,11 +4,16 @@ Soil Parameterization
Specifying the domain properties requires input of the
:doc:`soil architecture <grid>`
(in relation to the grid) and of the soil parameters. The latter incorporate
the *subscale physics*, the representation of the soil hydraulic properties
below the REV scale. DORiE requires several input files for retrieving this
information, depending on the type of grid used for the computation.
the *subscale physics*, the representation of the soil hydraulic and solute
properties below the REV scale. DORiE requires several input files for
retrieving this information, depending on the type of grid used for the
computation.
A YAML_ parameter file is always required. It specifies a set of soil
.. contents::
:depth: 3
:local:
A YAML_ parameter file is always required. It specifies a set of
parameterizations, along with a medium index. This index is used to reference
the medium specified in the architecture files, or with the "global" indices
given in the config file.
......@@ -24,13 +29,18 @@ YAML Parameter File
-------------------
This file is used to specify the parameterization for each medium inside the
simulated domain. It follows a simple hierarchical syntax. The file path and
name must be specified via the ``parameters.file`` key of the
name must be specified via the ``<model>.parameters.file`` key of the
:doc:`config file <config-file>`.
The top-level mapping must contain the key ``volumes``. This key contains a
sequence of arbitrary parameterization names. These, in turn, contain the
parameterization ``type``, the medium ``index``, and the actual ``parameters``.
The medium index must be an **integer**.
medium ``index``, and the model type (``richards`` or ``transport``).
The medium index must be an **integer**. Each model key contains the
parameterization ``type``, and the actual ``parameters``.
.. note:: Parameterization data for model ``transport`` is only required if
the model is actually enabled in the :doc:`config file settings
<config-file>`, via ``simulation.mode = richards+transport``.
Heterogeneities throughout the entire domain are set via the top level key
``scaling``. It must contain a supported scaling ``type``, which may default
......@@ -38,6 +48,10 @@ to ``none``. In case an actual scaling is to be applied, the section must
contain the key ``data``, which specifies the datasets required for this type
of scaling.
.. warning::
Transport parameters like `porosity` or `characteristic length` currently
are not affected by small scale heterogoeneities.
Every ``scaling_section`` has the same keys: The H5 filepath is given by
``file``, and the internal dataset path by ``dataset``. You can then choose an
``interpolation`` method for the dataset, which requires you to insert values
......@@ -51,11 +65,19 @@ shapes.
volumes:
<name-0>:
type: <string>
index: <int>
parameters:
<param-name-0>: <float>
# more parameters ...
richards:
type: <string>
parameters:
<param-name-0>: <float>
# more parameters ...
# if transport mode is enabled ...
transport:
type: <string>
parameters:
<param-name-0>: <float or sequence>
# more parameters ...
# more volumes ...
......@@ -74,9 +96,9 @@ shapes.
You find a documentation of the parameterizations implemented in DORiE
along with the parameters they require below.
The parameterization ``type`` must match the static ``type`` member of one of
them. Likewise, the parameters must match the names of parameters
these objects use.
The parameterization ``type`` must match a static parameterization name or a
valid combination of them. Likewise, the parameters must match the names of
parameters these objects use.
.. _man-parameter_scaling:
......@@ -101,8 +123,8 @@ This section lists the available types for parameterizations, scalings, and
interpolators in a understandable format. You can also dig through the code
documentation below, or the code itself, to retrieve this information.
Parameterizations
^^^^^^^^^^^^^^^^^
Richards Parameterizations
~~~~~~~~~~~~~~~~~~~~~~~~~~
As the *soil porosity* :math:`\phi \, [-]` and the *residual water content*
:math:`\theta_s \, [-]` vary between soil types, it is convenient to define the
......@@ -111,13 +133,15 @@ As the *soil porosity* :math:`\phi \, [-]` and the *residual water content*
.. math::
\Theta = \frac{\theta_w - \theta_r}{\phi - \theta_r},
\quad \Theta \in \left[ 0, 1 \right]
\quad \Theta \in \left[ 0, 1 \right],
where :math:`\theta_w \, [-]` is the volumetric soil water content.
One typically assumes that the entire pore space can be filled up with water,
hence we set the *saturated water content* :math:`\theta_s \, [-]` equal to the
porosity, :math:`\theta_s = \phi`.
.. object:: Mualem–van Genuchten parameterization
Mualem–van Genuchten
""""""""""""""""""""
Implements the following parameterization functions:
......@@ -136,8 +160,8 @@ porosity, :math:`\theta_s = \phi`.
Parameters:
* ``theta_r``: Residual water content :math:`\theta_r`
* ``theta_s``: Saturated water content / Porosity :math:`\theta_s`
* ``k0``: Saturated conductivity :math:`K_0 \, [\mathrm{ms}^{-1}]`
* ``alpha``: Air-entry value :math:`\alpha \, [\mathrm{m}^{-1}]`
* ``k0``: Saturated conductivity :math:`K_0 \, [\text{ms}^{-1}]`
* ``alpha``: Air-entry value :math:`\alpha \, [\text{m}^{-1}]`
* ``n``: Retention curve shape factor :math:`n`
* ``tau``: Skew factor :math:`\tau`
......@@ -145,6 +169,7 @@ porosity, :math:`\theta_s = \phi`.
.. code-block:: yaml
richards:
type: MvG
parameters:
theta_r:
......@@ -154,8 +179,323 @@ porosity, :math:`\theta_s = \phi`.
n:
tau:
Transport Parameterizations
~~~~~~~~~~~~~~~~~~~~~~~~~~~
Regardless of the parameterization, the transport simulation always computes
the microscopic péclet number, for which it requires the characteristic pore
length, the molecular diffusion, and the fluid velocity. The latter is directly
provided by the richards simulation while the other two have to be specified
for each volume:
Permanent parameters:
* ``mol_diff``: Molecular diffusion
:math:`D_m \, [\text{m}^2\text{s}^{-1}]`
* ``char_length``: Characteristic pore length :math:`\ell \, [\text{m}]`
.. note::
We support two options for specifying tensors through the YAML syntax.
You may either specify every entry of the tensor with a dedicated key, like
.. code-block:: yaml
<param>_xx: <value_xx>
<param>_xy: <value_xy>
# ...
or give an entire sequence that will be mapped to the respective entries,
.. code-block:: yaml
<param>: [<value_xx>, <value_xy> # , ...
]
The sequence is interpreted as a flattened tensor with row-major layout.
.. warning::
You must omit any component containing the spatial direction ``z`` in a
2D setup.
The hydrodynamic dispersion tensor :math:`\mathsf{D}_\text{hd} \,
[\text{m}^2\text{s}^{-1}]` is the main parameter to provide in the
transport simulation. Below you will find several parameterization for this.
Constant
""""""""
In this case, the hydrodynamic dispersion tensor is inserted directly
component-by-compoment.
.. note::
From a physical point of view, the hydrodynamic tensor must be
symmetric, but the user input is not verified by DORiE in this regard.
.. math::
\mathsf{D}_\text{hd} = \text{const}.
* ``type: Dhd_const``
Parameters:
* ``hydrodynamic_disp_<ij>``: (i, j)-th component of the
hydrodynamic dispersion tensor,
:math:`\left( \mathsf{D}_\text{hd} \right)_{ij} \, [\text{m}^2\text{s}^{-1}]`,
**or**
* ``hydrodynamic_disp``: Flattened hydrodynamic dispersion tensor
:math:`\mathsf{D}_\text{hd} \, [\text{m}^2\text{s}^{-1}]`.
YAML template:
.. code-block:: yaml
transport:
type: Dhd_const
parameters:
mol_diff:
char_length:
# sequence notation, or...
hydrodynamic_disp: [ ]
# component notation
hydrodynamic_disp_xx:
hydrodynamic_disp_xy:
hydrodynamic_disp_xz:
hydrodynamic_disp_yx:
hydrodynamic_disp_yy:
hydrodynamic_disp_yz:
hydrodynamic_disp_zx:
hydrodynamic_disp_zy:
hydrodynamic_disp_zz:
Power Law
"""""""""
Implements the following parameterization function:
.. math::
D_\text{hd} = \gamma D_m \operatorname{pe}^\alpha.
* ``type: Dhd_pl``
Parameters:
* ``gamma``: Scale the power law :math:`\gamma \, [-]`
* ``alpha``: Exponent of the power law :math:`\alpha \, [-]`
* ``mol_diff``: Molecular diffusion
:math:`D_m \, [\text{m}^2\text{s}^{-1}]`
The Peclét number :math:`\operatorname{pe}` is specified in the
:doc:`config file <config-file>`.
YAML template:
.. code-block:: yaml
transport:
type: Dhd_pl
parameters:
mol_diff:
char_length:
alpha:
gamma:
Superposition
"""""""""""""
The hydrodynamic dispersion tensor is said to be the superposition of
several diffusion/dispersion processes. In DORiE this possible by summing
several valid parameterizations types. Currently we provide
parameterizations for the *effective diffusion*
:math:`D_\text{eff}` and for the *effective hydromechanic tensor*
:math:`\mathsf{D}_\text{hm}` identified by the key prefixes ``Deff`` and
``Dhm`` respectively.
.. math::
\mathsf{D}_\text{hd} = \mathsf{D}_\text{hm} + D_\text{eff}.
* ``type: <Dhm_type> + <Deff_type>``
Each of the types are further parameterized and can be found below.
Effective Diffusion
^^^^^^^^^^^^^^^^^^^
Constant Effective Diffusion
############################
In this case, the effective diffusion is inserted directly,
.. math::
D_\text{eff} = \text{const}.
* ``Deff_type: Deff_const``
Parameters:
* ``eff_diff``: Effective diffusion
:math:`D_\text{eff} \, [\text{m}^2\text{s}^{-1}]`
YAML template:
.. code-block:: yaml
transport:
type: <Dhm_type> + Deff_const
parameters:
mol_diff:
char_length:
eff_diff:
# <Dhm_type> parameters ...
Milligton-Quirk I Effective Diffusion
#####################################
Implements the following parameterization function:
.. math::
D_\text{eff} = D_m \frac{\theta_w^{7/3}}{\phi^{2/3}}.
where the volumetric water content :math:`\theta_w \, [-]`
is automatically obtained from the Richards simulation.
* ``Deff_type: Deff_MQ1``
Parameters:
* ``mol_diff``: Molecular diffusion
:math:`D_m \, [\text{m}^2\text{s}^{-1}]`
* ``phi``: Soil porosity :math:`\phi \, [-]`
YAML template:
.. code-block:: yaml
transport:
type: <Dhm_type> + Deff_MQ1
parameters:
mol_diff:
char_length:
phi:
# <Dhm_type> parameters ...
Milligton-Quirk II Effective Diffusion
######################################
Implements the following parameterization function:
.. math::
D_\text{eff} = D_m \frac{\theta_w}{\phi^{2/3}}.
where the volumetric water content :math:`\theta_w \, [-]`
is automatically obtained from the Richards simulation.
* ``Deff_type: Deff_MQ2``
Parameters:
* ``mol_diff``: Molecular diffusion
:math:`D_m \, [\text{m}^2\text{s}^{-1}]`
* ``phi``: Soil porosity :math:`\phi \, [-]`
YAML template:
.. code-block:: yaml
transport:
type: <Dhm_type> + Deff_MQ1
parameters:
mol_diff:
char_length:
phi:
# <Dhm_type> parameters ...
Effective Hydromechanic Dispersion
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Constant Effective Hydromechanic Dispersion Tensor
##################################################
In this case, the effective hydromechanic dispersion tensor is inserted
directly.
.. math::
\mathsf{D}_\text{hm} = \text{const}.
* ``Dhm_type: Dhm_const``
Parameters:
* ``eff_hydromechanic_disp_<ij>``: (i, j)-th component of the
hydromechanic dispersion tensor,
:math:`\left(\mathsf{D}_\text{hm}\right)_{ij} \, [\text{m}^2\text{s}^{-1}]`,
**or**
* ``eff_hydromechanic_disp``: Flattened hydromechanic dispersion tensor
:math:`\mathsf{D}_\text{hm} \, [\text{m}^2\text{s}^{-1}]`.
YAML template:
.. code-block:: yaml
transport:
type: Dhm_const + <Deff_type>
parameters:
mol_diff:
char_length:
# sequence notation, or...
eff_hydromechanic_disp: [ ]
# component notation
eff_hydromechanic_disp_xx:
eff_hydromechanic_disp_xy:
eff_hydromechanic_disp_xz:
eff_hydromechanic_disp_yx:
eff_hydromechanic_disp_yy:
eff_hydromechanic_disp_yz:
eff_hydromechanic_disp_zx:
eff_hydromechanic_disp_zy:
eff_hydromechanic_disp_zz:
# <Deff_type> parameters ...
Effective Hydromechanic Dispersion Tensor for Isotropic Media
#############################################################
Implements the following parameterization function:
.. math::
\left( \mathsf{D}_\text{hm} \right)_{ij}
= (\lambda_l-\lambda_t)\frac{v_i v_j}{\lvert \mathbf{v} \rvert}
+ \delta_{ij}\lambda_t \lvert \mathbf{v} \rvert,
where :math:`\mathbf{v} \, [\text{ms}^{-1}]` is the local fluid velocity
and :math:`\delta_{ij}` is the Kronecker delta.
* ``Dhm_type: Dhm_iso``
Parameters:
* ``lambda_l``: Longitudinal dispersivity :math:`\lambda_l \, [\text{m}^2\text{s}^{-1}]`
* ``lambda_t``: Transverse dispersivity :math:`\lambda_t \, [\text{m}^2\text{s}^{-1}]`
YAML template:
.. code-block:: yaml
transport:
type: Dhm_iso + <Deff_type>
parameters:
mol_diff:
char_length:
lambda_l:
lambda_t:
# <Deff_type> parameters ...
Scalings
^^^^^^^^
~~~~~~~~
Every ``scaling_section`` has the following layout:
.. code-block:: yaml
......@@ -178,35 +518,36 @@ coordinates in the respective spatial dimensions.
omitted. Only omitting the respective *values* of both keys will lead to
a parser error.
.. object:: Miller scaling
Miller Scaling
""""""""""""""
Assumes geometric similarity between scaled regions. Applies the following
scaling:
Assumes geometric similarity between scaled regions. Applies the following
scaling:
.. math ::
.. math ::
\Theta &= \Theta (h_m \cdot \xi_M)\\
K &= K (\Theta) \cdot \xi_M^2
\Theta &= \Theta (h_m \cdot \xi_M)\\
K &= K (\Theta) \cdot \xi_M^2
* ``type: Miller``
* ``type: Miller``
Scaling sections:
* ``scale_miller``: Scaling factor :math:`\xi_M` to be applied onto
matric head and conductivity simultaneously.
Scaling sections:
* ``scale_miller``: Scaling factor :math:`\xi_M` to be applied onto
matric head and conductivity simultaneously.
YAML template:
YAML template:
.. code-block:: yaml
.. code-block:: yaml
type: Miller
data: