diff --git a/ProcessLib/HT/CreateHTProcess.cpp b/ProcessLib/HT/CreateHTProcess.cpp
index 0adb84aba85e7036164a621c3b6c29520a2dd1dd..f8584ed408eeb5227fe09150317e433b2b01fe9f 100644
--- a/ProcessLib/HT/CreateHTProcess.cpp
+++ b/ProcessLib/HT/CreateHTProcess.cpp
@@ -88,6 +88,8 @@ std::unique_ptr<Process> createHTProcess(
auto const dispersion_config =
//! \ogs_file_param{prj__processes__process__HT__thermal_dispersivity}
config.getConfigSubtreeOptional("thermal_dispersivity");
+ bool const has_fluid_thermal_dispersivity =
+ dispersion_config ? true : false;
if (dispersion_config)
{
thermal_dispersivity_longitudinal = &findParameter<double>(
@@ -154,6 +156,7 @@ std::unique_ptr<Process> createHTProcess(
auto const solid_config =
//! \ogs_file_param{prj__processes__process__HT__solid_thermal_expansion}
config.getConfigSubtreeOptional("solid_thermal_expansion");
+ const bool has_fluid_thermal_expansion = solid_config ? true : false;
if (solid_config)
{
solid_thermal_expansion = &findParameter<double>(
@@ -172,13 +175,15 @@ std::unique_ptr<Process> createHTProcess(
std::move(porous_media_properties),
density_solid,
std::move(fluid_properties),
+ has_fluid_thermal_dispersivity,
*thermal_dispersivity_longitudinal,
*thermal_dispersivity_transversal,
specific_heat_capacity_solid,
thermal_conductivity_solid,
thermal_conductivity_fluid,
- default_solid_thermal_expansion,
- default_biot_constant,
+ has_fluid_thermal_expansion,
+ *solid_thermal_expansion,
+ *biot_constant,
specific_body_force,
has_gravity);
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 1b0eed69354c85f9dee209300200a5d916432b3b..72888f29e4b22b8132f2c3a860c0d02378423d3a 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -58,6 +58,7 @@ public:
// matrices.
assert(local_matrix_size == ShapeFunction::NPOINTS * dof_per_node);
(void)local_matrix_size;
+ (void)dof_per_node;
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -128,15 +129,14 @@ protected:
thermal_conductivity_solid * (1 - porosity) +
thermal_conductivity_fluid * porosity;
- auto const thermal_dispersivity_longitudinal =
+ if (!material_properties.has_fluid_thermal_dispersivity)
+ return thermal_conductivity * I;
+
+ double const thermal_dispersivity_longitudinal =
material_properties.thermal_dispersivity_longitudinal(t, pos)[0];
auto const thermal_dispersivity_transversal =
material_properties.thermal_dispersivity_transversal(t, pos)[0];
- if (thermal_dispersivity_longitudinal == 0.0 &&
- thermal_dispersivity_transversal == 0.0)
- return thermal_conductivity * I;
-
double const velocity_magnitude = velocity.norm();
GlobalDimMatrixType const thermal_dispersivity =
fluid_density * specific_heat_capacity_fluid *
diff --git a/ProcessLib/HT/HTMaterialProperties.h b/ProcessLib/HT/HTMaterialProperties.h
index f667c0bad728f186c863f935e24e2aaf2ed6bbf3..995b3c4d4512bb399ab687706bc2e2b464a528fa 100644
--- a/ProcessLib/HT/HTMaterialProperties.h
+++ b/ProcessLib/HT/HTMaterialProperties.h
@@ -30,11 +30,13 @@ struct HTMaterialProperties
ProcessLib::Parameter<double> const& density_solid_,
std::unique_ptr<MaterialLib::Fluid::FluidProperties>&&
fluid_properties_,
+ bool const has_fluid_thermal_dispersivity_,
ProcessLib::Parameter<double> const& thermal_dispersivity_longitudinal_,
ProcessLib::Parameter<double> const& thermal_dispersivity_transversal_,
ProcessLib::Parameter<double> const& specific_heat_capacity_solid_,
ProcessLib::Parameter<double> const& thermal_conductivity_solid_,
ProcessLib::Parameter<double> const& thermal_conductivity_fluid_,
+ bool const has_fluid_thermal_expansion_,
Parameter<double> const& solid_thermal_expansion_,
Parameter<double> const& biot_constant_,
Eigen::VectorXd specific_body_force_,
@@ -43,10 +45,12 @@ struct HTMaterialProperties
density_solid(density_solid_),
fluid_properties(std::move(fluid_properties_)),
specific_heat_capacity_solid(specific_heat_capacity_solid_),
+ has_fluid_thermal_dispersivity(has_fluid_thermal_dispersivity_),
thermal_dispersivity_longitudinal(thermal_dispersivity_longitudinal_),
thermal_dispersivity_transversal(thermal_dispersivity_transversal_),
thermal_conductivity_solid(thermal_conductivity_solid_),
thermal_conductivity_fluid(thermal_conductivity_fluid_),
+ has_fluid_thermal_expansion(has_fluid_thermal_expansion_),
solid_thermal_expansion(solid_thermal_expansion_),
biot_constant(biot_constant_),
specific_body_force(std::move(specific_body_force_)),
@@ -58,11 +62,13 @@ struct HTMaterialProperties
Parameter<double> const& density_solid;
std::unique_ptr<MaterialLib::Fluid::FluidProperties> fluid_properties;
Parameter<double> const& specific_heat_capacity_solid;
+ bool const has_fluid_thermal_dispersivity;
Parameter<double> const& thermal_dispersivity_longitudinal;
Parameter<double> const& thermal_dispersivity_transversal;
Parameter<double> const& thermal_conductivity_solid;
Parameter<double> const& thermal_conductivity_fluid;
+ bool const has_fluid_thermal_expansion;
Parameter<double> const& solid_thermal_expansion;
Parameter<double> const& biot_constant;
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index f76fa2a9fa628e6c81ef783e22583852d8f9a405..98e305be1df611797368ba8e2e78adc0d379955e 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -75,10 +75,26 @@ void HTProcess::assembleConcreteProcess(const double t,
GlobalMatrix& K,
GlobalVector& b)
{
- DBUG("Assemble HTProcess.");
-
- if (!_is_monolithic_scheme)
+ if (_is_monolithic_scheme)
+ {
+ DBUG("Assemble HTProcess.");
+ }
+ else
+ {
+ if (_coupled_solutions->process_id == 0)
+ {
+ DBUG(
+ "Assemble the equations of heat transport process within "
+ "HTProcess.");
+ }
+ else
+ {
+ DBUG(
+ "Assemble the equations of single phase fully saturated "
+ "fluid flow process within HTProcess.");
+ }
setCoupledSolutionsOfPreviousTimeStep();
+ }
// Call global assembler for each local assembly item.
GlobalExecutor::executeMemberDereferenced(
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 7389c7a63bcc88bd1bfced93813fbdd4a3acbba8..a5add61c701ced6d0f97d51b55e9fb727d6d50f8 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -26,17 +26,17 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& local_b_data,
- LocalCoupledSolutions const& coupled_term)
+ LocalCoupledSolutions const& coupled_xs)
{
- if (coupled_term.variable_id == 0)
+ if (coupled_xs.process_id == 0)
{
assembleHeatTransportEquation(t, local_M_data, local_K_data,
- local_b_data, coupled_term);
+ local_b_data, coupled_xs);
return;
}
assembleHydraulicEquation(t, local_M_data, local_K_data, local_b_data,
- coupled_term);
+ coupled_xs);
}
template <typename ShapeFunction, typename IntegrationMethod,
@@ -45,17 +45,17 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
assembleHydraulicEquation(double const t, std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& local_b_data,
- LocalCoupledSolutions const& coupled_term)
+ LocalCoupledSolutions const& coupled_xs)
{
- auto const& local_p = coupled_term.local_coupled_xs[0];
+ auto const& local_p = coupled_xs.local_coupled_xs[1];
auto const local_matrix_size = local_p.size();
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
assert(local_matrix_size == ShapeFunction::NPOINTS);
- auto const& local_T1 = coupled_term.local_coupled_xs[1];
- auto const& local_T0 = coupled_term.local_coupled_xs0[1];
- const double dt = coupled_term.dt;
+ auto const& local_T1 = coupled_xs.local_coupled_xs[0];
+ auto const& local_T0 = coupled_xs.local_coupled_xs0[0];
+ const double dt = coupled_xs.dt;
auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
local_M_data, local_matrix_size, local_matrix_size);
@@ -71,9 +71,6 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
auto const& b = material_properties.specific_body_force;
- GlobalDimMatrixType const& I(
- GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
-
MaterialLib::Fluid::FluidProperty::ArrayType vars;
unsigned const n_integration_points =
@@ -110,10 +107,6 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
material_properties.fluid_properties->getdValue(
MaterialLib::Fluid::FluidPropertyType::Density, vars,
MaterialLib::Fluid::PropertyVariableType::p);
- const double dfluid_density_dT =
- material_properties.fluid_properties->getdValue(
- MaterialLib::Fluid::FluidPropertyType::Density, vars,
- MaterialLib::Fluid::PropertyVariableType::T);
// Use the viscosity model to compute the viscosity
auto const viscosity = material_properties.fluid_properties->getValue(
@@ -139,11 +132,23 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
local_K.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
+ if (material_properties.has_gravity)
+ {
+ local_b.noalias() +=
+ w * fluid_density * dNdx.transpose() * K_over_mu * b;
+ }
+
+ if (!material_properties.has_fluid_thermal_expansion)
+ return;
+
// Add the thermal expansion term
- auto const solid_thermal_expansion =
- material_properties.solid_thermal_expansion(t, pos)[0];
- if (solid_thermal_expansion > 0.0)
{
+ auto const solid_thermal_expansion =
+ material_properties.solid_thermal_expansion(t, pos)[0];
+ const double dfluid_density_dT =
+ material_properties.fluid_properties->getdValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars,
+ MaterialLib::Fluid::PropertyVariableType::T);
double T0_at_xi = 0.;
NumLib::shapeFunctionInterpolate(local_T0, N, T0_at_xi);
auto const biot_constant =
@@ -154,12 +159,6 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
local_b.noalias() +=
eff_thermal_expansion * (T1_at_xi - T0_at_xi) * w * N / dt;
}
-
- if (material_properties.has_gravity)
- {
- local_b.noalias() +=
- w * fluid_density * dNdx.transpose() * K_over_mu * b;
- }
}
}
@@ -170,9 +169,9 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& /*local_b_data*/,
- LocalCoupledSolutions const& coupled_term)
+ LocalCoupledSolutions const& coupled_xs)
{
- auto const& local_p = coupled_term.local_coupled_xs[0];
+ auto const& local_p = coupled_xs.local_coupled_xs[1];
auto const local_matrix_size = local_p.size();
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
@@ -181,7 +180,7 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
auto local_p_Eigen_type =
Eigen::Map<const NodalVectorType>(&local_p[0], local_matrix_size);
- auto const& local_T1 = coupled_term.local_coupled_xs[1];
+ auto const& local_T1 = coupled_xs.local_coupled_xs[0];
auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
local_M_data, local_matrix_size, local_matrix_size);
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index b2ee534b49ba19417d93d8bc37f22ee001da216f..a8b37f16bde70ea8051d991f0208124e9f5f0a0c 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -65,7 +65,7 @@ public:
void assembleWithCoupledTerm(
double const t, std::vector<double>& local_M_data,
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
- LocalCoupledSolutions const& coupled_term) override;
+ LocalCoupledSolutions const& coupled_xs) override;
std::vector<double> const& getIntPtDarcyVelocity(
const double t,
@@ -78,12 +78,12 @@ private:
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& local_b_data,
- LocalCoupledSolutions const& coupled_term);
+ LocalCoupledSolutions const& coupled_xs);
void assembleHeatTransportEquation(
double const t, std::vector<double>& local_M_data,
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
- LocalCoupledSolutions const& coupled_term);
+ LocalCoupledSolutions const& coupled_xs);
};
} // namespace HT