From 46d6acefd97c7a6cb3e69882ab89a01abe09b128 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Wed, 18 Jan 2017 13:41:22 +0100
Subject: [PATCH] [PL/HT] Use porosity model from MaterialLib.
---
ProcessLib/HT/CreateHTProcess.cpp | 14 ++++++------
ProcessLib/HT/HTFEM.h | 36 +++++++++++++++++--------------
ProcessLib/HT/HTProcessData.h | 9 ++++----
3 files changed, 32 insertions(+), 27 deletions(-)
diff --git a/ProcessLib/HT/CreateHTProcess.cpp b/ProcessLib/HT/CreateHTProcess.cpp
index d2793614eb7..fd73cdcc7eb 100644
--- a/ProcessLib/HT/CreateHTProcess.cpp
+++ b/ProcessLib/HT/CreateHTProcess.cpp
@@ -17,6 +17,7 @@
#include "MaterialLib/Fluid/Density/createFluidDensityModel.h"
#include "MaterialLib/Fluid/Viscosity/createViscosityModel.h"
+#include "MaterialLib/PorousMedium/Porosity/createPorosityModel.h"
namespace ProcessLib
{
@@ -48,12 +49,11 @@ std::unique_ptr<Process> createHTProcess(
//! \ogs_file_param_special{prj__processes__process__HT__process_variables__pressure}
"pressure"});
- // Porosity parameter.
- auto& porosity = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HT__porosity}
- "porosity", parameters, 1);
- DBUG("Use \'%s\' as porosity parameter.", porosity.name.c_str());
+ auto const& porous_medium_config = config.getConfigSubtree("porous_medium");
+ auto const& porosity_conf =
+ porous_medium_config.getConfigSubtree("porosity");
+ auto porosity_model =
+ MaterialLib::PorousMedium::createPorosityModel(porosity_conf);
// Parameter for the intrinsic permeability (only one scalar per element,
// i.e., the isotropic case is handled at the moment)
@@ -157,7 +157,7 @@ std::unique_ptr<Process> createHTProcess(
std::copy_n(b.data(), b.size(), specific_body_force.data());
HTProcessData process_data{
- porosity,
+ std::move(porosity_model),
intrinsic_permeability,
specific_storage,
std::move(viscosity_model),
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index ba546eeca89..266e29947f7 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -127,31 +127,36 @@ public:
auto const intrinsic_permeability =
_process_data.intrinsic_permeability(t, pos)[0];
- auto const porosity = _process_data.porosity(t, pos)[0];
-
- auto const specific_heat_capacity_solid =
- _process_data.specific_heat_capacity_solid(t, pos)[0];
- auto const specific_heat_capacity_fluid =
- _process_data.specific_heat_capacity_fluid(t, pos)[0];
- double const heat_capacity =
- density_solid * specific_heat_capacity_solid * (1 - porosity) +
- fluid_reference_density * specific_heat_capacity_fluid * porosity;
-
auto const thermal_conductivity_solid =
_process_data.thermal_conductivity_solid(t, pos)[0];
auto const thermal_conductivity_fluid =
_process_data.thermal_conductivity_fluid(t, pos)[0];
+ auto const& sm = _shape_matrices[ip];
+ double T_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ // Order matters: First T, then P!
+ NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt, p_int_pt);
+
+ // TODO the first argument has to be changed for non constant
+ // porosity model
+ auto const porosity =
+ _process_data.porosity_model->getValue(0.0, T_int_pt);
+
Eigen::MatrixXd const thermal_conductivity =
(thermal_conductivity_solid * (1 - porosity) +
thermal_conductivity_fluid * porosity) * unit_mat;
+ auto const specific_heat_capacity_solid =
+ _process_data.specific_heat_capacity_solid(t, pos)[0];
+ auto const specific_heat_capacity_fluid =
+ _process_data.specific_heat_capacity_fluid(t, pos)[0];
+
auto const thermal_dispersivity_longitudinal =
_process_data.thermal_dispersivity_longitudinal(t, pos)[0];
auto const thermal_dispersivity_transversal =
_process_data.thermal_dispersivity_transversal(t, pos)[0];
- auto const& sm = _shape_matrices[ip];
auto const& wp = _integration_method.getWeightedPoint(ip);
auto Ktt = local_K.template block<num_nodes, num_nodes>(0, 0);
auto Mtt = local_M.template block<num_nodes, num_nodes>(0, 0);
@@ -161,11 +166,6 @@ public:
num_nodes);
auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
- double T_int_pt = 0.0;
- double p_int_pt = 0.0;
- // Order matters: First T, then P!
- NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt, p_int_pt);
-
// Use the fluid density model to compute the density
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
@@ -194,6 +194,10 @@ public:
auto const hydrodynamic_thermodispersion =
thermal_conductivity + thermal_dispersivity;
+ double const heat_capacity =
+ density_solid * specific_heat_capacity_solid * (1 - porosity) +
+ fluid_reference_density * specific_heat_capacity_fluid * porosity;
+
auto const integral_term =
sm.integralMeasure * sm.detJ * wp.getWeight();
// matrix assembly
diff --git a/ProcessLib/HT/HTProcessData.h b/ProcessLib/HT/HTProcessData.h
index 7ed768d97f3..1a79cad410a 100644
--- a/ProcessLib/HT/HTProcessData.h
+++ b/ProcessLib/HT/HTProcessData.h
@@ -12,6 +12,7 @@
#include <memory>
#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
namespace MeshLib
{
@@ -28,7 +29,7 @@ namespace HT
struct HTProcessData
{
HTProcessData(
- ProcessLib::Parameter<double> const& porosity_,
+ std::unique_ptr<MaterialLib::PorousMedium::Porosity>&& porosity_model_,
ProcessLib::Parameter<double> const& intrinsic_permeability_,
ProcessLib::Parameter<double> const& specific_storage_,
std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& viscosity_model_,
@@ -43,7 +44,7 @@ struct HTProcessData
ProcessLib::Parameter<double> const& thermal_conductivity_fluid_,
Eigen::Vector3d const& specific_body_force_,
bool const has_gravity_)
- : porosity(porosity_),
+ : porosity_model(std::move(porosity_model_)),
intrinsic_permeability(intrinsic_permeability_),
specific_storage(specific_storage_),
viscosity_model(std::move(viscosity_model_)),
@@ -62,7 +63,7 @@ struct HTProcessData
}
HTProcessData(HTProcessData&& other)
- : porosity(other.porosity),
+ : porosity_model(other.porosity_model.release()),
intrinsic_permeability(other.intrinsic_permeability),
specific_storage(other.specific_storage),
viscosity_model(other.viscosity_model.release()),
@@ -91,7 +92,7 @@ struct HTProcessData
//! Assignments are not needed.
void operator=(HTProcessData&&) = delete;
- Parameter<double> const& porosity;
+ std::unique_ptr<MaterialLib::PorousMedium::Porosity> porosity_model;
Parameter<double> const& intrinsic_permeability;
Parameter<double> const& specific_storage;
std::unique_ptr<MaterialLib::Fluid::FluidProperty> viscosity_model;
--
GitLab