From 12519e6ab234b5847b02951b118a35a47590f168 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 18 Oct 2017 17:16:23 +0200
Subject: [PATCH] [HT] Modified the velocity calculation in order that it can
be used for both the monolithic and staggered schemes.
---
ProcessLib/HT/HTFEM.h | 130 +++++++++++-----------
ProcessLib/HT/HTLocalAssemblerInterface.h | 21 ++++
ProcessLib/HT/HTProcess.cpp | 18 ++-
ProcessLib/HT/HTProcess.h | 8 +-
ProcessLib/HT/MonolithicHTFEM.h | 20 ++++
ProcessLib/HT/StaggeredHTFEM-impl.h | 24 +++-
ProcessLib/HT/StaggeredHTFEM.h | 6 +
ProcessLib/UncoupledProcessesTimeLoop.cpp | 7 +-
8 files changed, 152 insertions(+), 82 deletions(-)
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index da24a85a767..1b0eed69354 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -49,7 +49,8 @@ public:
unsigned const integration_order,
HTMaterialProperties const& material_properties,
const unsigned dof_per_node)
- : _element(element),
+ : HTLocalAssemblerInterface(),
+ _element(element),
_material_properties(material_properties),
_integration_method(integration_order)
{
@@ -86,72 +87,6 @@ public:
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
- std::vector<double> const& getIntPtDarcyVelocity(
- const double t,
- GlobalVector const& current_solution,
- NumLib::LocalToGlobalIndexMap const& dof_table,
- std::vector<double>& cache) const override
- {
- auto const n_integration_points =
- _integration_method.getNumberOfPoints();
-
- auto const indices = NumLib::getIndices(_element.getID(), dof_table);
- assert(!indices.empty());
- auto const local_x = current_solution.get(indices);
-
- cache.clear();
- auto cache_mat = MathLib::createZeroedMatrix<
- Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, GlobalDim, n_integration_points);
-
- SpatialPosition pos;
- pos.setElementID(_element.getID());
-
- MaterialLib::Fluid::FluidProperty::ArrayType vars;
-
- 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;
-
- pos.setIntegrationPoint(ip);
-
- 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 K =
- _material_properties.porous_media_properties.getIntrinsicPermeability(
- t, pos).getValue(t, pos, 0.0, T_int_pt);
-
- auto const mu = _material_properties.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
- GlobalDimMatrixType const K_over_mu = K / mu;
-
- cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
-
- if (_material_properties.has_gravity)
- {
- auto const rho_w =
- _material_properties.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::Density, vars);
- auto const b = _material_properties.specific_body_force;
- // here it is assumed that the vector b is directed 'downwards'
- cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
- }
- }
-
- return cache;
- }
-
protected:
MeshLib::Element const& _element;
HTMaterialProperties const& _material_properties;
@@ -212,6 +147,67 @@ protected:
return thermal_conductivity * I + thermal_dispersivity;
}
+
+ std::vector<double> const& getIntPtDarcyVelocityLocal(
+ const double t, std::vector<double> const& local_p,
+ std::vector<double> const& local_T, std::vector<double>& cache) const
+ {
+ auto const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, GlobalDim, n_integration_points);
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+ &local_p[0], 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;
+
+ pos.setIntegrationPoint(ip);
+
+ double T_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
+ NumLib::shapeFunctionInterpolate(local_T, N, T_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 K = _material_properties.porous_media_properties
+ .getIntrinsicPermeability(t, pos)
+ .getValue(t, pos, 0.0, T_int_pt);
+
+ auto const mu = _material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+ GlobalDimMatrixType const K_over_mu = K / mu;
+
+ cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
+
+ if (_material_properties.has_gravity)
+ {
+ auto const rho_w =
+ _material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ auto const b = _material_properties.specific_body_force;
+ // here it is assumed that the vector b is directed 'downwards'
+ cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
+ }
+ }
+
+ return cache;
+ }
};
} // namespace HT
diff --git a/ProcessLib/HT/HTLocalAssemblerInterface.h b/ProcessLib/HT/HTLocalAssemblerInterface.h
index 82408a4d99f..ecb4d10faf5 100644
--- a/ProcessLib/HT/HTLocalAssemblerInterface.h
+++ b/ProcessLib/HT/HTLocalAssemblerInterface.h
@@ -17,6 +17,8 @@
namespace ProcessLib
{
+struct CoupledSolutionsForStaggeredScheme;
+
namespace HT
{
template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
@@ -42,11 +44,30 @@ class HTLocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
public NumLib::ExtrapolatableElement
{
public:
+ HTLocalAssemblerInterface() : _coupled_solutions(nullptr) {}
+ void setStaggeredCoupledSolutions(
+ std::size_t const /*mesh_item_id*/,
+ CoupledSolutionsForStaggeredScheme* const coupling_term)
+ {
+ _coupled_solutions = coupling_term;
+ }
+
virtual std::vector<double> const& getIntPtDarcyVelocity(
const double /*t*/,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
+
+protected:
+ // TODO: remove _coupled_solutions or move integration point data from local
+ // assembler class to a new class to make local assembler unique for each
+ // process.
+ /** Pointer to CoupledSolutionsForStaggeredScheme that is set in a member of
+ * Process class, setCoupledTermForTheStaggeredSchemeToLocalAssemblers.
+ * It is used for calculate the secondary variables like velocity for
+ * coupled processes.
+ */
+ CoupledSolutionsForStaggeredScheme* _coupled_solutions;
};
} // namespace HT
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index c17570c121a..f76fa2a9fa6 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -104,9 +104,9 @@ void HTProcess::assembleWithJacobianConcreteProcess(
}
void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
- const double /*t*/,
- const double /*delta_t*/,
- const int process_id)
+ const double /*t*/,
+ const double /*delta_t*/,
+ const int process_id)
{
assert(process_id < 2);
@@ -123,11 +123,19 @@ void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
auto& x0 = *_xs_previous_timestep[process_id];
MathLib::LinAlg::copy(x, x0);
}
-
+
auto& x0 = *_xs_previous_timestep[process_id];
MathLib::LinAlg::setLocalAccessibleVector(x0);
}
+void HTProcess::setCoupledTermForTheStaggeredSchemeToLocalAssemblers()
+{
+ DBUG("Set the coupled term for the staggered scheme to local assembers.");
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &HTLocalAssemblerInterface::setStaggeredCoupledSolutions,
+ _local_assemblers, _coupled_solutions);
+}
+
} // namespace HT
} // namespace ProcessLib
-
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 61b399d572a..0d5672d187a 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -67,6 +67,8 @@ public:
return _xs_previous_timestep[process_id].get();
}
+ void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() override;
+
private:
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -82,9 +84,9 @@ private:
const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
- void preTimestepConcreteProcess(
- GlobalVector const& x, double const t, double const dt,
- const int process_id) override;
+ void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+ double const dt,
+ const int process_id) override;
const std::unique_ptr<HTMaterialProperties> _material_properties;
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 6aed470778d..775b56f1370 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -183,6 +183,26 @@ public:
* in buoyancy effects */
}
}
+
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const override
+ {
+ auto const indices =
+ NumLib::getIndices(this->_element.getID(), dof_table);
+ assert(!indices.empty());
+ auto local_x = current_solution.get(indices);
+
+ std::vector<double> local_T(
+ std::make_move_iterator(local_x.begin() + local_x.size() / 2),
+ std::make_move_iterator(local_x.end()));
+ // only p is kept in local_x
+ local_x.erase(local_x.begin() + local_x.size() / 2, local_x.end());
+
+ return this->getIntPtDarcyVelocityLocal(t, local_x, local_T, cache);
+ }
};
} // namespace HT
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 5a2a916180d..7389c7a63bc 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -30,13 +30,13 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
{
if (coupled_term.variable_id == 0)
{
- assembleHydraulicEquation(t, local_M_data, local_K_data, local_b_data,
- coupled_term);
+ assembleHeatTransportEquation(t, local_M_data, local_K_data,
+ local_b_data, coupled_term);
return;
}
- assembleHeatTransportEquation(t, local_M_data, local_K_data, local_b_data,
- coupled_term);
+ assembleHydraulicEquation(t, local_M_data, local_K_data, local_b_data,
+ coupled_term);
}
template <typename ShapeFunction, typename IntegrationMethod,
@@ -267,5 +267,21 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
}
}
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::vector<double> const&
+StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
+ getIntPtDarcyVelocity(const double t,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const
+{
+ auto const indices = NumLib::getIndices(this->_element.getID(), dof_table);
+ assert(!indices.empty());
+ auto const local_xs =
+ getCurrentLocalSolutions(*(this->_coupled_solutions), indices);
+
+ return this->getIntPtDarcyVelocityLocal(t, local_xs[0], local_xs[1], cache);
+}
} // namespace HT
} // namespace ProcessLib
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index ca1217cee25..b2ee534b49b 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -67,6 +67,12 @@ public:
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
LocalCoupledSolutions const& coupled_term) override;
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const override;
+
private:
void assembleHydraulicEquation(double const t,
std::vector<double>& local_M_data,
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index 92f8a612967..00f19e43d7f 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -987,11 +987,12 @@ void UncoupledProcessesTimeLoop::outputSolutions(
spd->timestepper->getTimeStep().dt(), pcs_idx);
if (is_staggered_coupling)
{
- CoupledSolutionsForStaggeredScheme coupled_xs(
+ CoupledSolutionsForStaggeredScheme coupled_solutions(
_solutions_of_coupled_processes, 0.0, pcs_idx);
- spd->process.setCoupledSolutionsForStaggeredScheme(&coupled_xs);
- spd->process.setStaggeredCouplingTermToLocalAssemblers();
+ spd->process.setCoupledSolutionsForStaggeredScheme(
+ &coupled_solutions);
+ spd->process.setCoupledTermForTheStaggeredSchemeToLocalAssemblers();
(output_object.*output_class_member)(pcs, spd->process_output,
timestep, t, x);
}
--
GitLab