From 7b71f23d2c6e577de66eddf08a486c574f91b127 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 13 Nov 2017 17:49:01 +0100
Subject: [PATCH] [PCS] Separate the process variables, BCs and STs by process
in class Process
---
ProcessLib/Process.cpp | 139 ++++++++++++++++++++++-------------------
ProcessLib/Process.h | 49 +++++++++------
2 files changed, 105 insertions(+), 83 deletions(-)
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index ad810458bbe..5eabe1688e9 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -24,18 +24,28 @@ Process::Process(
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
unsigned const integration_order,
- std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
SecondaryVariableCollection&& secondary_variables,
- NumLib::NamedFunctionCaller&& named_function_caller)
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ const bool use_monolithic_scheme)
: _mesh(mesh),
_secondary_variables(std::move(secondary_variables)),
_named_function_caller(std::move(named_function_caller)),
_global_assembler(std::move(jacobian_assembler)),
- _is_monolithic_scheme(true),
+ _use_monolithic_scheme(use_monolithic_scheme),
_coupled_solutions(nullptr),
_integration_order(integration_order),
_process_variables(std::move(process_variables)),
- _boundary_conditions(parameters)
+ _boundary_conditions([&]() -> std::vector<BoundaryConditionCollection> {
+ std::vector<BoundaryConditionCollection> pcs_BCs;
+ pcs_BCs.reserve(process_variables.size());
+ for (std::size_t i = 0; i < process_variables.size(); i++)
+ {
+ pcs_BCs.emplace_back(BoundaryConditionCollection(parameters));
+ }
+ return pcs_BCs;
+ }())
{
}
@@ -58,56 +68,69 @@ void Process::initialize()
finishNamedFunctionsInitialization();
DBUG("Initialize boundary conditions.");
- _boundary_conditions.addBCsForProcessVariables(
- _process_variables, *_local_to_global_index_map, _integration_order);
-
- for (int variable_id = 0;
- variable_id < static_cast<int>(_process_variables.size());
- ++variable_id)
+ for (std::size_t pcs_id = 0; pcs_id < _process_variables.size(); pcs_id++)
{
- ProcessVariable& pv = _process_variables[variable_id];
- auto sts =
- pv.createSourceTerms(*_local_to_global_index_map, variable_id,
- _integration_order);
+ auto const& per_process_variables = _process_variables[pcs_id];
+ auto& per_process_BCs = _boundary_conditions[pcs_id];
+ std::vector<std::unique_ptr<NodalSourceTerm>> per_process_source_terms;
+ for (std::size_t variable_id = 0;
+ variable_id < per_process_variables.size();
+ variable_id++)
+ {
+ per_process_BCs.addBCsForProcessVariables(
+ per_process_variables, *_local_to_global_index_map,
+ _integration_order);
+
+ ProcessVariable& pv = per_process_variables[variable_id];
+ auto sts = pv.createSourceTerms(*_local_to_global_index_map, 0,
+ _integration_order);
- std::move(sts.begin(), sts.end(),
- std::back_inserter(_source_terms));
+ std::move(sts.begin(), sts.end(),
+ std::back_inserter(per_process_source_terms));
+ }
+ _source_terms.push_back(std::move(per_process_source_terms));
}
}
-void Process::setVariableInitialCondition(const int variable_id, double const t,
- GlobalVector& x)
+void Process::setInitialConditions(const unsigned pcs_id, double const t,
+ GlobalVector& x)
{
- SpatialPosition pos;
+ auto const& per_process_variables = _process_variables[pcs_id];
+ for (std::size_t variable_id = 0;
+ variable_id < per_process_variables.size();
+ variable_id++)
+ {
+ SpatialPosition pos;
- auto const& pv = _process_variables[variable_id];
- DBUG("Set the initial condition of variable %s.",
- pv.get().getName().data());
+ auto const& pv = per_process_variables[variable_id];
+ DBUG("Set the initial condition of variable %s of process %d.",
+ pv.get().getName().data(), pcs_id);
- auto const& ic = pv.get().getInitialCondition();
+ auto const& ic = pv.get().getInitialCondition();
- auto const num_comp = pv.get().getNumberOfComponents();
+ auto const num_comp = pv.get().getNumberOfComponents();
- const int mesh_subset_id = _is_monolithic_scheme ? variable_id : 0;
+ const int mesh_subset_id = _use_monolithic_scheme ? variable_id : 0;
- for (int component_id = 0; component_id < num_comp; ++component_id)
- {
- auto const& mesh_subsets = _local_to_global_index_map->getMeshSubsets(
- mesh_subset_id, component_id);
- for (auto const& mesh_subset : mesh_subsets)
+ for (int component_id = 0; component_id < num_comp; ++component_id)
{
- auto const mesh_id = mesh_subset->getMeshID();
- for (auto const* node : mesh_subset->getNodes())
+ auto const& mesh_subsets =
+ _local_to_global_index_map->getMeshSubsets(mesh_subset_id,
+ component_id);
+ for (auto const& mesh_subset : mesh_subsets)
{
- MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
- node->getID());
-
- pos.setNodeID(node->getID());
- auto const& ic_value = ic(t, pos);
-
- auto global_index =
- std::abs(_local_to_global_index_map->getGlobalIndex(
- l, mesh_subset_id, component_id));
+ auto const mesh_id = mesh_subset->getMeshID();
+ for (auto const* node : mesh_subset->getNodes())
+ {
+ MeshLib::Location const l(
+ mesh_id, MeshLib::MeshItemType::Node, node->getID());
+
+ pos.setNodeID(node->getID());
+ auto const& ic_value = ic(t, pos);
+
+ auto global_index =
+ std::abs(_local_to_global_index_map->getGlobalIndex(
+ l, mesh_subset_id, component_id));
#ifdef USE_PETSC
// The global indices of the ghost entries of the global
// matrix or the global vectors need to be set as negative
@@ -120,28 +143,13 @@ void Process::setVariableInitialCondition(const int variable_id, double const t,
if (global_index == x.size())
global_index = 0;
#endif
- x.set(global_index, ic_value[component_id]);
+ x.set(global_index, ic_value[component_id]);
+ }
}
}
}
}
-void Process::setInitialConditions(const unsigned pcs_id, double const t,
- GlobalVector& x)
-{
- if (!_is_monolithic_scheme)
- {
- setVariableInitialCondition(pcs_id, t, x);
- return;
- }
-
- for (std::size_t variable_id = 0; variable_id < _process_variables.size();
- variable_id++)
- {
- setVariableInitialCondition(variable_id, t, x);
- }
-}
-
MathLib::MatrixSpecifications Process::getMatrixSpecifications() const
{
auto const& l = *_local_to_global_index_map;
@@ -161,9 +169,12 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
assembleConcreteProcess(t, x, M, K, b);
- _boundary_conditions.applyNaturalBC(t, x, K, b);
+ const auto pcs_id =
+ (_coupled_solutions) ? _coupled_solutions->process_id : 0;
+ _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b);
- for (auto const& st : _source_terms)
+ const auto _source_terms_per_pcs = _source_terms[pcs_id];
+ for (auto const& st : _source_terms_per_pcs)
{
st->integrateNodalSourceTerm(t, b);
}
@@ -182,7 +193,9 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x,
Jac);
// TODO apply BCs to Jacobian.
- _boundary_conditions.applyNaturalBC(t, x, K, b);
+ const auto pcs_id =
+ (_coupled_solutions) ? _coupled_solutions->process_id : 0;
+ _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b);
}
void Process::constructDofTable()
@@ -196,7 +209,7 @@ void Process::constructDofTable()
// Vector of the number of variable components
std::vector<int> vec_var_n_components;
- if (_is_monolithic_scheme)
+ if (_use_monolithic_scheme)
{
// Collect the mesh subsets in a vector.
for (ProcessVariable const& pv : _process_variables)
@@ -243,7 +256,7 @@ void Process::initializeExtrapolator()
bool manage_storage;
if (_local_to_global_index_map->getNumberOfComponents() == 1 ||
- !_is_monolithic_scheme)
+ !_use_monolithic_scheme)
{
// For single-variable-single-component processes reuse the existing DOF
// table.
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 39b2cc50b3b..9bdd143265e 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -46,10 +46,11 @@ public:
std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
unsigned const integration_order,
- std::vector<std::reference_wrapper<ProcessVariable>>&&
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
process_variables,
SecondaryVariableCollection&& secondary_variables,
- NumLib::NamedFunctionCaller&& named_function_caller);
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ const bool use_monolithic_scheme = true);
/// Preprocessing before starting assembly for new timestep.
void preTimestep(GlobalVector const& x, const double t,
@@ -78,14 +79,9 @@ public:
_coupled_solutions = coupled_solutions;
}
- void setDecouplingSchemeType(const bool is_monolithic_scheme)
- {
- _is_monolithic_scheme = is_monolithic_scheme;
- }
- bool useMonolithicScheme() const { return _is_monolithic_scheme; }
-
- virtual void setCoupledSolutionsForStaggeredSchemeToLocalAssemblers() {}
+ bool isMonolithicSchemeUsed() const { return _use_monolithic_scheme; }
+ virtual void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() {}
void preAssemble(const double t, GlobalVector const& x) override final;
void assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b) final;
@@ -99,7 +95,9 @@ public:
std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
getKnownSolutions(double const t) const final
{
- return _boundary_conditions.getKnownSolutions(t);
+ const auto pcs_id =
+ (_coupled_solutions) ? _coupled_solutions->process_id : 0;
+ return _boundary_conditions[pcs_id].getKnownSolutions(t);
}
NumLib::LocalToGlobalIndexMap const& getDOFTable() const
@@ -111,7 +109,9 @@ public:
std::vector<std::reference_wrapper<ProcessVariable>> const&
getProcessVariables() const
{
- return _process_variables;
+ const auto pcs_id =
+ (_coupled_solutions) ? _coupled_solutions->process_id : 0;
+ return _process_variables[pcs_id];
}
SecondaryVariableCollection const& getSecondaryVariables() const
@@ -120,6 +120,7 @@ public:
}
// Get the solution of the previous time step.
+
virtual GlobalVector* getPreviousTimeStepSolution(
const int /*process_id*/) const
{
@@ -127,6 +128,7 @@ public:
}
// Used as a call back for CalculateSurfaceFlux process.
+
virtual std::vector<double> getFlux(std::size_t /*element_id*/,
MathLib::Point3d const& /*p*/,
GlobalVector const& /*x*/) const
@@ -204,9 +206,6 @@ private:
/// DOF-table.
void computeSparsityPattern();
- void setVariableInitialCondition(const int variable_id, double const t,
- GlobalVector& x);
-
protected:
MeshLib::Mesh& _mesh;
std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
@@ -222,7 +221,7 @@ protected:
VectorMatrixAssembler _global_assembler;
- mutable bool _is_monolithic_scheme;
+ const bool _use_monolithic_scheme;
/// Pointer to CoupledSolutionsForStaggeredScheme, which contains the
/// references to the solutions of the coupled processes.
@@ -240,11 +239,21 @@ protected:
private:
GlobalSparsityPattern _sparsity_pattern;
- /// Variables used by this process.
- std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
-
- BoundaryConditionCollection _boundary_conditions;
- std::vector<std::unique_ptr<NodalSourceTerm>> _source_terms;
+ /// Variables used by this process. For the monolithic scheme or a
+ /// single process, the size of the outer vector is one. For the
+ /// staggered scheme, the size of the outer vector is the number of the
+ /// coupled processes.
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+ _process_variables;
+
+ /// Vector for boundary conditions. For the monolithic scheme or a
+ /// single process, the size of the vector is one. For the staggered
+ /// scheme, the size of vector is the number of the coupled processes.
+ std::vector<BoundaryConditionCollection> _boundary_conditions;
+
+ /// Vector for nodal source terms. The outer vector is for processes,
+ /// which has the same size as that for boundary conditions.
+ std::vector<std::vector<std::unique_ptr<NodalSourceTerm>>> _source_terms;
ExtrapolatorData _extrapolator_data;
};
--
GitLab