From 7d394aa60e1ba9b102571c8e220a797ba111d857 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri, 24 Feb 2017 13:53:05 +0100
Subject: [PATCH] [LF] Added an optional input of reference temperature for
 LiquidFlow

---
 .../LIQUID_FLOW/t_reference_temperature.md    |  3 +++
 .../LiquidFlow/CreateLiquidFlowProcess.cpp    | 13 +++++++++++--
 .../LiquidFlowLocalAssembler-impl.h           | 19 +++++++------------
 .../LiquidFlow/LiquidFlowLocalAssembler.h     |  3 +++
 ProcessLib/LiquidFlow/LiquidFlowProcess.cpp   |  5 ++++-
 ProcessLib/LiquidFlow/LiquidFlowProcess.h     |  2 ++
 6 files changed, 30 insertions(+), 15 deletions(-)
 create mode 100644 Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/t_reference_temperature.md

diff --git a/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/t_reference_temperature.md b/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/t_reference_temperature.md
new file mode 100644
index 00000000000..914f3123da7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/t_reference_temperature.md
@@ -0,0 +1,3 @@
+An optional input of a reference temperature for liquid flow process when the
+ flow process is not coupled to heat transport process, and a value other than
+ the default one (room temperature of 273.15+18 K) is selected.
diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
index b1c0981b5ba..e81a4749bb6 100644
--- a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
@@ -15,6 +15,7 @@
 
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 
+#include "MaterialLib/PhysicalConstant.h"
 #include "ProcessLib/Utils/ParseSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
 
@@ -67,6 +68,13 @@ std::unique_ptr<Process> createLiquidFlowProcess(
     assert(g >= 0.);
     const int gravity_axis_id = (g == 0.) ? -1 : gravity_axis_id_input;
 
+    auto const& refT =
+        //! \ogs_file_param{prj__processes__process__LIQUID_FLOW__reference_temperature}
+        config.getConfigParameterOptional<double>("reference_temperature");
+    const double reference_temperature =
+        refT ? *refT
+             : MaterialLib::PhysicalConstant::CelsiusZeroInKelvin + 18.0;
+
     //! \ogs_file_param{prj__processes__process__LIQUID_FLOW__material_property}
     auto const& mat_config = config.getConfigSubtree("material_property");
 
@@ -80,7 +88,7 @@ std::unique_ptr<Process> createLiquidFlowProcess(
             mesh, std::move(jacobian_assembler), parameters, integration_order,
             std::move(process_variables), std::move(secondary_variables),
             std::move(named_function_caller), *mat_ids, has_material_ids,
-            gravity_axis_id, g, mat_config}};
+            gravity_axis_id, g, reference_temperature, mat_config}};
     }
     else
     {
@@ -101,7 +109,8 @@ std::unique_ptr<Process> createLiquidFlowProcess(
             mesh, std::move(jacobian_assembler), parameters, integration_order,
             std::move(process_variables), std::move(secondary_variables),
             std::move(named_function_caller), *dummy_property_vector,
-            has_material_ids, gravity_axis_id, g, mat_config}};
+            has_material_ids, gravity_axis_id, g, reference_temperature,
+            mat_config}};
     }
 }
 
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 50f6db7ead6..a7fb7140b89 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -14,7 +14,6 @@
 
 #include "LiquidFlowLocalAssembler.h"
 
-#include "MaterialLib/PhysicalConstant.h"
 #include "NumLib/Function/Interpolation.h"
 
 #include "ProcessLib/HeatConduction/HeatConductionProcess.h"
@@ -133,10 +132,6 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     //       the integration loop for non-constant porosity and storage models.
     double porosity_variable = 0.;
     double storage_variable = 0.;
-    // Reference temperature for fluid property models because of no coupled
-    // heat transport process. TODO: Add an optional input for this.
-    const double temperature =
-        MaterialLib::PhysicalConstant::CelsiusZeroInKelvin + 18.0;
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         auto const& sm = _shape_matrices[ip];
@@ -151,15 +146,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         // Assemble mass matrix, M
         local_M.noalias() += _material_properties.getMassCoefficient(
                                  material_id, t, pos, porosity_variable,
-                                 storage_variable, p, temperature) *
+                                 storage_variable, p, _reference_temperature) *
                              sm.N.transpose() * sm.N * integration_factor;
 
         // Compute density:
         const double rho_g =
-            _material_properties.getLiquidDensity(p, temperature) *
+            _material_properties.getLiquidDensity(p, _reference_temperature) *
             _gravitational_acceleration;
         // Compute viscosity:
-        const double mu = _material_properties.getViscosity(p, temperature);
+        const double mu =
+            _material_properties.getViscosity(p, _reference_temperature);
 
         // Assemble Laplacian, K, and RHS by the gravitational term
         LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
@@ -286,8 +282,6 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
 
-    const double temperature =
-        MaterialLib::PhysicalConstant::CelsiusZeroInKelvin + 18.0;
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         auto const& sm = _shape_matrices[ip];
@@ -295,10 +289,11 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
 
         const double rho_g =
-            _material_properties.getLiquidDensity(p, temperature) *
+            _material_properties.getLiquidDensity(p, _reference_temperature) *
             _gravitational_acceleration;
         // Compute viscosity:
-        const double mu = _material_properties.getViscosity(p, temperature);
+        const double mu =
+            _material_properties.getViscosity(p, _reference_temperature);
 
         LaplacianGravityVelocityCalculator::calculateVelocity(
             _darcy_velocities, local_p_vec, sm, permeability, ip, mu, rho_g,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 17850159952..d45f13017ee 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -71,6 +71,7 @@ public:
         unsigned const integration_order,
         int const gravitational_axis_id,
         double const gravitational_acceleration,
+        double const reference_temperature,
         LiquidFlowMaterialProperties const& material_propertries)
         : _element(element),
           _integration_method(integration_order),
@@ -79,6 +80,7 @@ public:
               element, is_axially_symmetric, _integration_method)),
           _gravitational_axis_id(gravitational_axis_id),
           _gravitational_acceleration(gravitational_acceleration),
+          _reference_temperature(reference_temperature),
           _material_properties(material_propertries)
     {
     }
@@ -219,6 +221,7 @@ private:
 
     const int _gravitational_axis_id;
     const double _gravitational_acceleration;
+    const double _reference_temperature;
     const LiquidFlowMaterialProperties& _material_properties;
 };
 
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 74c9bde7697..3eaaa0b739b 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -41,12 +41,14 @@ LiquidFlowProcess::LiquidFlowProcess(
     bool const has_material_ids,
     int const gravitational_axis_id,
     double const gravitational_acceleration,
+    double const reference_temperature,
     BaseLib::ConfigTree const& config)
     : Process(mesh, std::move(jacobian_assembler), parameters,
               integration_order, std::move(process_variables),
               std::move(secondary_variables), std::move(named_function_caller)),
       _gravitational_axis_id(gravitational_axis_id),
       _gravitational_acceleration(gravitational_acceleration),
+      _reference_temperature(reference_temperature),
       _material_properties(createLiquidFlowMaterialProperties(
           config, parameters, has_material_ids, material_ids))
 {
@@ -63,7 +65,8 @@ void LiquidFlowProcess::initializeConcreteProcess(
         mesh.getDimension(), mesh.getElements(), dof_table,
         pv.getShapeFunctionOrder(), _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _gravitational_axis_id,
-        _gravitational_acceleration, *_material_properties);
+        _gravitational_acceleration, _reference_temperature,
+        *_material_properties);
 
     _secondary_variables.addSecondaryVariable(
         "darcy_velocity_x", 1,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index dd5836da3ac..fca99047492 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -70,6 +70,7 @@ public:
         bool const has_material_ids,
         int const gravitational_axis_id,
         double const gravitational_acceleration,
+        double const reference_temperature,
         BaseLib::ConfigTree const& config);
 
     void computeSecondaryVariableConcrete(
@@ -106,6 +107,7 @@ private:
 
     const int _gravitational_axis_id;
     const double _gravitational_acceleration;
+    const double _reference_temperature;
     const std::unique_ptr<LiquidFlowMaterialProperties> _material_properties;
 
     std::vector<std::unique_ptr<LiquidFlowLocalAssemblerInterface>>
-- 
GitLab