From a98cbd778aaec3e2e805c9e3882fdadfd2ed9cee Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Wed, 7 Jun 2017 12:48:00 +0200
Subject: [PATCH] [PL/HT] Compute secondary vars. using actual vals.
---
ProcessLib/HT/HTFEM.h | 48 ++++++++++++++++++++++++++++++++++++-
ProcessLib/HT/HTProcess.cpp | 10 ++++++++
ProcessLib/HT/HTProcess.h | 3 +++
3 files changed, 60 insertions(+), 1 deletion(-)
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 8d9a61ae3b5..72f0d5a7580 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -261,8 +261,54 @@ public:
Bp += w * density_water_T * dNdx.transpose() * perm_over_visc * b;
/* with Oberbeck-Boussing assumption density difference only exists
* in buoyancy effects */
+ }
+ }
+
+ void computeSecondaryVariableConcrete(
+ double const t, std::vector<double> const& local_x) override
+ {
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+
+ auto const K =
+ _process_data.porous_media_properties.getIntrinsicPermeability(t,
+ pos);
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+ &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ auto const& ip_data = _ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+
+ double T_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+
+ auto const mu = _process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+ GlobalDimMatrixType const K_over_mu = K / mu;
+
+ GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values;
+ if (_process_data.has_gravity)
+ {
+ auto const rho_w = _process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ auto const b = _process_data.specific_body_force;
+ // here it is assumed that the vector b is directed 'downwards'
+ velocity += K_over_mu * rho_w * b;
+ }
- // velocity computed for output.
for (unsigned d = 0; d < GlobalDim; ++d)
{
_darcy_velocities[d][ip] = velocity[d];
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 3cebe46de02..7a70c5ad767 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -98,6 +98,16 @@ void HTProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac, coupling_term);
}
+void HTProcess::computeSecondaryVariableConcrete(
+ double const t, GlobalVector const& x,
+ StaggeredCouplingTerm const& coupling_term)
+{
+ DBUG("Compute the Darcy velocity for HTProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &HTLocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
+ *_local_to_global_index_map, t, x, coupling_term);
+}
+
} // namespace HT
} // namespace ProcessLib
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index a5048e295cd..6dcdaaeaf26 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -56,6 +56,9 @@ public:
bool isLinear() const override { return false; }
//! @}
+ void computeSecondaryVariableConcrete(
+ double const t, GlobalVector const& x,
+ StaggeredCouplingTerm const& coupling_term) override;
private:
void initializeConcreteProcess(
--
GitLab