From 697294ed3f2869d0dde93b53f62b6dd36f0035ef Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 16 May 2022 16:24:22 +0200
Subject: [PATCH 1/6] [NumLib] Added a file of NumericalType.h, which firstly
 contains

enum class UpwindingScheme
---
 NumLib/NumericalType.h | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)
 create mode 100644 NumLib/NumericalType.h

diff --git a/NumLib/NumericalType.h b/NumLib/NumericalType.h
new file mode 100644
index 00000000000..65b9009a889
--- /dev/null
+++ b/NumLib/NumericalType.h
@@ -0,0 +1,21 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on May 16, 2022, 10:53 AM
+ */
+
+#pragma once
+
+namespace NumLib
+{
+enum class UpwindingScheme
+{
+    none = 0,
+    full_upwinding = 1
+};
+}
-- 
GitLab


From be5d15c62f3db3b71430c62bb1c9c7aad3576a83 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 16 May 2022 16:25:49 +0200
Subject: [PATCH 2/6] [Process/HT] Added upwinding type data and their input.

---
 ProcessLib/HT/CreateHTProcess.cpp | 51 ++++++++++++++++++++++++++++++-
 ProcessLib/HT/HTProcessData.h     |  9 ++++--
 2 files changed, 57 insertions(+), 3 deletions(-)

diff --git a/ProcessLib/HT/CreateHTProcess.cpp b/ProcessLib/HT/CreateHTProcess.cpp
index ed8cc003b17..4376c3c8f4b 100644
--- a/ProcessLib/HT/CreateHTProcess.cpp
+++ b/ProcessLib/HT/CreateHTProcess.cpp
@@ -16,6 +16,7 @@
 #include "MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MeshLib/IO/readMeshFromFile.h"
+#include "NumLib/NumericalType.h"
 #include "ParameterLib/ConstantParameter.h"
 #include "ParameterLib/Utils.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
@@ -80,6 +81,53 @@ std::unique_ptr<Process> createHTProcess(
     const bool use_monolithic_scheme =
         !(coupling_scheme && (*coupling_scheme == "staggered"));
 
+    NumLib::UpwindingScheme upwinding_scheme_T = NumLib::UpwindingScheme::none;
+    NumLib::UpwindingScheme upwinding_scheme_H = NumLib::UpwindingScheme::none;
+    auto const upwinding_config =
+        //! \ogs_file_param{prj__processes__process__HT__upwinding_scheme}
+        config.getConfigSubtreeOptional("upwinding_scheme");
+    //
+    if (upwinding_config)
+    {
+        auto getUpwindingScheme =
+            [](std::string const& scheme_name, std::string const& equation_name)
+        {
+            if (scheme_name == "FullUpwinding")
+            {
+                INFO("Full upwinding scheme is applied for the {:s}.",
+                     equation_name);
+                return NumLib::UpwindingScheme::full_upwinding;
+            }
+            else
+            {
+                OGS_FATAL(
+                    "The specified upwinding scheme of {:s} for the {:s} "
+                    "is not available so far.",
+                    scheme_name, equation_name);
+            }
+        };
+
+        auto const upwinding_scheme_T_ptr =
+            //! \ogs_file_param{prj__processes__process__HT__upwinding_scheme__energy_balance_equation}
+            upwinding_config->getConfigParameterOptional<std::string>(
+                "energy_balance_equation");
+        if (upwinding_scheme_T_ptr)
+        {
+            upwinding_scheme_T = getUpwindingScheme(*upwinding_scheme_T_ptr,
+                                                    "energy balance equation");
+        }
+
+        auto const upwinding_scheme_H_ptr =
+            //! \ogs_file_param{prj__processes__process__HT__upwinding_scheme__mass_balance_equation}
+            upwinding_config->getConfigParameterOptional<std::string>(
+                "mass_balance_equation");
+        if (upwinding_scheme_H_ptr)
+        {
+            upwinding_scheme_H = getUpwindingScheme(*upwinding_scheme_H_ptr,
+                                                    "mass balance equation");
+        }
+    }
+
     /// \section processvariablesht Process Variables
 
     //! \ogs_file_param{prj__processes__process__HT__process_variables}
@@ -183,7 +231,8 @@ std::unique_ptr<Process> createHTProcess(
         std::move(media_map),      has_fluid_thermal_expansion,
         *solid_thermal_expansion,  *biot_constant,
         specific_body_force,       has_gravity,
-        heat_transport_process_id, hydraulic_process_id};
+        heat_transport_process_id, hydraulic_process_id,
+        upwinding_scheme_T,        upwinding_scheme_H};
 
     SecondaryVariableCollection secondary_variables;
 
diff --git a/ProcessLib/HT/HTProcessData.h b/ProcessLib/HT/HTProcessData.h
index 8378bac3ce1..0395713244f 100644
--- a/ProcessLib/HT/HTProcessData.h
+++ b/ProcessLib/HT/HTProcessData.h
@@ -10,12 +10,12 @@
 
 #pragma once
 
+#include <Eigen/Eigen>
 #include <memory>
 #include <utility>
 
-#include <Eigen/Eigen>
-
 #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "NumLib/NumericalType.h"
 #include "ParameterLib/Parameter.h"
 
 namespace ProcessLib
@@ -33,6 +33,11 @@ struct HTProcessData final
     bool const has_gravity;
     int const heat_transport_process_id;
     int const hydraulic_process_id;
+
+    /// Upwinding scheme for the energy balance equation (T).
+    NumLib::UpwindingScheme upwinding_scheme_T;
+    /// Upwinding scheme for the mass balance equation (H).
+    NumLib::UpwindingScheme upwinding_scheme_H;
 };
 
 }  // namespace HT
-- 
GitLab


From d0e0c087a816fec61ee539edaebd6cde59948f9d Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 17 May 2022 16:50:49 +0200
Subject: [PATCH 3/6] [bib] Added a reference about the fully upwinding scheme

---
 Documentation/bibliography/other.bib | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/Documentation/bibliography/other.bib b/Documentation/bibliography/other.bib
index 974274cd5d2..6eec22e7952 100644
--- a/Documentation/bibliography/other.bib
+++ b/Documentation/bibliography/other.bib
@@ -421,3 +421,15 @@ URL = {https://doi.org/10.1680/geot.2008.58.3.157}
   year={2017},
   publisher={Elsevier}
 }
+
+@article{dalen1979,
+    author = "Dalen, V.",
+    abstract = "This paper summarizes some research that was conducted to construct finite-element models for reservoir flow problems. The models are based on Galerkin's method, but the method is applied in an unorthodox manner to simplify calculation of coefficients and to improve stability. Specifically, techniques of compatibility relaxation, capacity lumping, and upstream mobility weighting are used, and schemes are obtained that seem to combine the simplicity and high stability of conventional finite-difference models with the generality and modeling flexibility of finite-element methods. The development of a model for single-phase gas flow and a two-phase oil/water model is described. Numerical examples are included to illustrate the usefulness of finite elements. In particular, the triangular element with linear interpolation is shown to be an attractive alternative to the standard five-point, finite-difference approximation for two-dimensional analysis.",
+    doi = "10.2118/7196-PA",
+    journal = "Society of Petroleum Engineers Journal",
+    pages = "333-342",
+    title = "{Simplified Finite-Element Models for Reservoir Flow Problems}",
+    url= "https://doi.org/10.2118/7196-PA"
+    volume = "19",
+    year = "1979"
+}
-- 
GitLab


From ff0237805807541e0e854f8b259c1afdcb942e98 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 17 May 2022 16:52:12 +0200
Subject: [PATCH 4/6] [NumLib] Added a description about the full upwinding
 scheme

---
 NumLib/NumericalType.h | 128 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 128 insertions(+)

diff --git a/NumLib/NumericalType.h b/NumLib/NumericalType.h
index 65b9009a889..7fd20fb9afb 100644
--- a/NumLib/NumericalType.h
+++ b/NumLib/NumericalType.h
@@ -13,6 +13,134 @@
 
 namespace NumLib
 {
+/**
+ *  So far, there is only full upwinding scheme is implemented. <br><br>
+ *   <li>
+ *    <b>Full upwinding scheme</b> <br>
+ *      The implemented full upwinding scheme is based on the scheme that was
+ *      introduced by Dalen \cite dalen1979.
+ *      The concept the full upwinding scheme is to evaluate the fluid velocity
+ *      related weak from terms in the flow upwinding mannar to achieve
+ numerical
+ *      stability.
+ *
+ *      <ol>
+ *      <li>  <b>Mass balance equation</b> <br>
+ *      For the discretised weak form of the mass balance equation, we consider
+ *      the Laplace term:
+ *      \f[
+ *           q_i = \int_{\Omega_e} \nabla \phi_i \rho \mathbf{K}(p)(\nabla p
+ *               + \rho \mathbf{g}) \mathrm{d}\Omega
+ *      \f]
+ *     where \f$\phi_i\f$ is the test function at node \f$i\f$,
+ \f$\mathbf{K}(p)\f$
+ *     is the hydraulic conductivity, \f$p\f$ is the pore pressure, \f$\rho\f$
+ is
+ *      the fluid density, \f$\mathbf g\f$ is the gravitational acceleration
+ vector,
+ *      and \f$\Omega_e\f$ is the element domain.
+ *     Let's consider an approximation of \f$L_i\f$ by choosing a constant
+ *     value of \f$\rho\mathbf{K}(p)\f$ such as
+ *      \f[
+ *           q_i \approx \overline{(\rho \mathbf{ K})} \int_{\Omega_e} \nabla
+ \phi_i  (\nabla p
+ *               + \rho \mathbf{g}) \mathrm{d}\Omega = \overline{(\rho \mathbf{
+ K})} {\hat q}_i
+ *      \f]
+ *      Since \f$q_i\f$ can be treated as the mass flux flowing throw node
+ \f$i\f$,
+ *      \f$\hat{q}_i\f$ is therefore interpreted as a measure of fluid flux
+ throw node \f$i\f$.
+ *      More over the summation of \f$\hat{q}_i\f$ is zero as
+ *      \f[
+ *          \sum_{i=1}^n{\hat q}_i = \sum_{i=1}^n\int_{\Omega_e} \nabla \phi_i
+ (\nabla p
+ *               + \rho \mathbf{g}) \mathrm{d}\Omega
+ *                = \int_{\Omega_e} \nabla (\sum_{i=1}^n \phi_i)  (\nabla p
+ *               + \rho \mathbf{g}) \mathrm{d}\Omega
+ *                =0
+ *      \f]
+ *      because of \f$\sum_{i=1}^n \phi_i=1\f$ with \f$n\f$ the number of
+ element nodes.
+ *      \f$ \sum_{i=1}^n{\hat q}_i=0\f$ guarantees the mass balance within an
+ element. The nodes with \f${\hat q}_i \geq 0\f$ are called upwind nodes because
+ that fluid
+ *       flow from them to the downwind nodes, at where  \f${\hat q}_i < 0 \f$.
+ *       The concept of the full upwinding scheme is to keep the values of
+ *       the nodal flux at the upwind nodes, while to compute the values of the
+ *       nodal flux at the downwind nodes with weight under the condition of
+ *       the mass convervation. That means
+ *        \f{eqnarray*}{
+ *          q_i \approx
+ *           \begin{cases}
+ *             \, \overline{(\rho \mathbf{ K})} {\hat q}_i,\quad \forall  {\hat
+ q}_i \geq 0\\
+ *              w \overline{(\rho \mathbf{ K})} {\hat q}_i,\quad \forall  {\hat
+ q}_i < 0
+ *           \end{cases}
+ *        \f}
+ *         where the \f$w\f$  is the upwinding weight.
+ *         Let
+ *         \f[
+ *              {\hat q}_{up}=\sum_{{\hat q}_i \geq 0}{\hat q}_i
+ *         \f]
+ *         be the summation of the nodal flux at the upwind nodes, and
+ *         \f[
+ *              {\hat q}_{down}=\sum_{{\hat q}_i < 0}{\hat q}_i
+ *         \f]
+ *         be the summation of the nodal flux at the donw nodes. The weight is
+ *         then given by
+ *         \f[
+ *             w = -\dfrac{{\hat q}_{up}}{{\hat q}_{down}}
+ *         \f]
+ *
+ *       Since \f${\hat q}_{up}+{\hat q}_{down}=0\f$, the element mass
+ *       conservation is guaranteed in the fully upwinging scheme.
+ *
+ *       <li>  <b>Heat/mass transport equation</b> <br>
+ *    For the heat/mass transport equation, we can apply the full upwinding
+ scheme
+ *   to  the advection term of
+ *    \f$
+ *          \nabla \cdot ( c u \mathbf{v}_f)
+ *    \f$
+ *     where \f$c\f$ is the coefficent, \f$u\f$ can be temperature \f$ T\f$ or
+ mass
+ *     component concentration \f$C\f$, \f$\mathbf{v}_f\f$ is the fluid velocity
+ *     defined as
+ *    \f$
+ *        \mathbf{v}_f = - \mathbf{K}(p)(\nabla p + \rho \mathbf{g}).
+ *     \f$
+ *    The discretised weak form of that advection term takes the form of
+ *    \f{eqnarray*}{
+ *           q_i &=& \int_{\Omega_e}  \nabla \cdot (c u \mathbf{v}_f)
+ *                \phi_i \mathrm{d}_\Omega\\
+ *               &=&  \int_{\Omega_e}  \nabla \cdot (c u \mathbf{v}_f \phi_i)
+ *                 \mathrm{d}_\Omega - \int_{\Omega_e} \nabla \phi_i
+ *                  \cdot (c u \mathbf{v}_f )  \mathrm{d}_\Omega
+ *    \f}
+ *    where the first term is converted to boundary integration according to
+ *     the Green's theorem  as
+ *          \f[
+ *               \int_{\Omega_e}  \nabla \cdot (c u \mathbf{v}_f \phi_i)
+ *                 \mathrm{d}_\Omega
+ *               = \int_{\partial\Omega_e}  (c u \mathbf{v}_f \phi_i) \mathbf{n}
+ *                 \mathrm{d}_\Gamma
+ *          \f]
+ *        with \f$\mathbf{n}\f$ the normal to the boundary surface. That
+ *        boundary integration is a part of Neumann boundary condition.
+ *     Therefore what left for the element integration is
+ *       \f[
+             - \int_{\Omega_e} \nabla \phi_i
+ *                  \cdot (c u \mathbf{v}_f )  \mathrm{d}_\Omega,
+ *        \f]
+ *      which takes the same form as that of the Laplace term of the weak form
+ *      of the mass balance equation, and for which the same full upwinding
+ *      approach as that is addressed for the mass balance equation can be
+ *    applied.
+ *   </ol>
+ *   </li>
+ */
 enum class UpwindingScheme
 {
     none = 0,
-- 
GitLab


From 84fb5eed5970d4063de40941fc0787608cd4f752 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 17 May 2022 17:15:25 +0200
Subject: [PATCH 5/6] [HT] Documented the input parameters of the upwinding
 scheme

---
 .../process/HT/upwinding_scheme/i_upwinding_scheme.md         | 4 ++++
 .../process/HT/upwinding_scheme/t_energy_balance_equation.md  | 1 +
 .../process/HT/upwinding_scheme/t_mass_balance_equation.md    | 1 +
 3 files changed, 6 insertions(+)
 create mode 100644 Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/i_upwinding_scheme.md
 create mode 100644 Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_energy_balance_equation.md
 create mode 100644 Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_mass_balance_equation.md

diff --git a/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/i_upwinding_scheme.md b/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/i_upwinding_scheme.md
new file mode 100644
index 00000000000..9a2ff8291ad
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/i_upwinding_scheme.md
@@ -0,0 +1,4 @@
+An optional tag to invoke the use of upwinding scheme for the mass balance
+ equation and the energy balance equation, individually.
+
+\copydoc NumLib::UpwindingScheme
diff --git a/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_energy_balance_equation.md b/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_energy_balance_equation.md
new file mode 100644
index 00000000000..848d9083b3c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_energy_balance_equation.md
@@ -0,0 +1 @@
+An optional tag to select an upwinding scheme for the energy balance equation.
diff --git a/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_mass_balance_equation.md b/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_mass_balance_equation.md
new file mode 100644
index 00000000000..a37d10ea7f5
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/HT/upwinding_scheme/t_mass_balance_equation.md
@@ -0,0 +1 @@
+An optional tag to select an upwinding scheme for the mass balance equation.
-- 
GitLab


From 175da5a3901e3cdb73aae1153d1123a9514d2806 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 19 May 2022 17:44:30 +0200
Subject: [PATCH 6/6] temp

---
 ProcessLib/HT/MonolithicHTFEM.h | 88 ++++++++++++++++++++++++++++++---
 1 file changed, 81 insertions(+), 7 deletions(-)

diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index fa8cd7d01be..f085f6764eb 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -12,6 +12,7 @@
 #pragma once
 
 #include <Eigen/Dense>
+#include <numeric>
 #include <vector>
 
 #include "HTFEM.h"
@@ -23,6 +24,7 @@
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/NumericalType.h"
 #include "ParameterLib/Parameter.h"
 
 namespace ProcessLib
@@ -37,6 +39,7 @@ class MonolithicHTFEM
 {
     using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
     using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
 
     using LocalMatrixType = typename ShapeMatricesType::template MatrixType<
         NUM_NODAL_DOF * ShapeFunction::NPOINTS,
@@ -92,6 +95,8 @@ public:
             pressure_index, pressure_index);
         auto Bp = local_b.template block<pressure_size, 1>(pressure_index, 0);
 
+        NodalMatrixType Kpp_average;
+
         ParameterLib::SpatialPosition pos;
         pos.setElementID(this->_element.getID());
 
@@ -114,6 +119,48 @@ public:
         unsigned const n_integration_points =
             this->_integration_method.getNumberOfPoints();
 
+        double averge_density = 0.0;
+        GlobalDimMatrixType K_over_mu_average =
+            GlobalDimMatrixType::Zero(GlobalDim, GlobalDim);
+
+        bool const using_full_upwinding =
+            (process_data.upwinding_scheme_T ==
+                 NumLib::UpwindingScheme::full_upwinding ||
+             process_data.upwinding_scheme_H ==
+                 NumLib::UpwindingScheme::full_upwinding);
+
+        if (using_full_upwinding)
+        {
+            Kpp_average = NodalMatrixType::Zero(ShapeFunction::NPOINTS,
+                                                ShapeFunction::NPOINTS);
+
+            double const T_avg =
+                std::accumulate(local_x.begin() + temperature_index,
+                                local_x.begin() + temperature_size, 0) /
+                temperature_size;
+            vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
+                T_avg;
+
+            double const p_avg =
+                std::accumulate(local_x.begin() + pressure_index,
+                                local_x.begin() + pressure_size, 0) /
+                pressure_size;
+            vars[static_cast<int>(
+                MaterialPropertyLib::Variable::phase_pressure)] = p_avg;
+
+            auto const intrinsic_permeability =
+                MaterialPropertyLib::formEigenTensor<GlobalDim>(
+                    medium
+                        .property(
+                            MaterialPropertyLib::PropertyType::permeability)
+                        .value(vars, pos, t, dt));
+            auto const viscosity =
+                liquid_phase
+                    .property(MaterialPropertyLib::PropertyType::viscosity)
+                    .template value<double>(vars, pos, t, dt);
+            K_over_mu_average = intrinsic_permeability / viscosity;
+        }
+
         for (unsigned ip(0); ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
@@ -183,21 +230,48 @@ public:
                 this->getThermalConductivityDispersivity(
                     vars, fluid_density, specific_heat_capacity_fluid, velocity,
                     I, pos, t, dt);
-            KTT.noalias() +=
-                (dNdx.transpose() * thermal_conductivity_dispersivity * dNdx +
-                 N.transpose() * velocity.transpose() * dNdx * fluid_density *
-                     specific_heat_capacity_fluid) *
-                w;
-            Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
+
+            if (using_full_upwinding)
+            {
+                Kpp_average.noalias() +=
+                    dNdx.transpose() * K_over_mu_average * dNdx * w;
+            }
+
             MTT.noalias() += w *
                              this->getHeatEnergyCoefficient(
                                  vars, porosity, fluid_density,
                                  specific_heat_capacity_fluid, pos, t, dt) *
                              N.transpose() * N;
+            if (process_data.upwinding_scheme_T ==
+                NumLib::UpwindingScheme::full_upwinding)
+            {
+                KTT.noalias() += dNdx.transpose() *
+                                 thermal_conductivity_dispersivity * dNdx * w;
+            }
+            else
+            {
+                KTT.noalias() +=
+                    (dNdx.transpose() * thermal_conductivity_dispersivity *
+                         dNdx +
+                     N.transpose() * velocity.transpose() * dNdx *
+                         fluid_density * specific_heat_capacity_fluid) *
+                    w;
+            }
+
+            if (process_data.upwinding_scheme_H !=
+                NumLib::UpwindingScheme::full_upwinding)
+            {
+                Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
+            }
+
             Mpp.noalias() += w * N.transpose() * specific_storage * N;
             if (process_data.has_gravity)
             {
-                Bp += w * fluid_density * dNdx.transpose() * K_over_mu * b;
+                Bp +=
+                    using_full_upwinding
+                        ? w * fluid_density * dNdx.transpose() *
+                              K_over_mu_average * b
+                        : w * fluid_density * dNdx.transpose() * K_over_mu * b;
             }
             /* with Oberbeck-Boussing assumption density difference only exists
              * in buoyancy effects */
-- 
GitLab