From 7c4d508e496ae0f8d4ea3ccb05314e20046bdb76 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 19 Dec 2016 09:40:49 +0100
Subject: [PATCH] [Coupling] Added staggered coupling computation.
---
NumLib/ODESolver/NonlinearSolver.cpp | 18 +-
NumLib/ODESolver/NonlinearSolver.h | 17 +-
NumLib/ODESolver/NonlinearSystem.h | 10 +-
NumLib/ODESolver/ODESystem.h | 10 +-
NumLib/ODESolver/TimeDiscretizedODESystem.cpp | 10 +-
NumLib/ODESolver/TimeDiscretizedODESystem.h | 10 +-
ProcessLib/Process.cpp | 6 +-
ProcessLib/Process.h | 9 +-
ProcessLib/StaggeredCouplingTerm.cpp | 26 ++
ProcessLib/StaggeredCouplingTerm.h | 42 +++
ProcessLib/UncoupledProcessesTimeLoop.cpp | 289 ++++++++++++++----
ProcessLib/UncoupledProcessesTimeLoop.h | 20 +-
Tests/NumLib/ODEs.h | 31 +-
Tests/NumLib/TimeLoopSingleODE.h | 8 +-
14 files changed, 416 insertions(+), 90 deletions(-)
create mode 100644 ProcessLib/StaggeredCouplingTerm.cpp
create mode 100644 ProcessLib/StaggeredCouplingTerm.h
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index bf55b1a5d8d..2353746d8cb 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -24,13 +24,14 @@
namespace NumLib
{
void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
- GlobalVector const& x) const
+ GlobalVector const& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term) const
{
- _equation_system->assemble(x);
+ _equation_system->assemble(x, coupled_term);
}
bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
- GlobalVector& x,
+ GlobalVector& x, ProcessLib::StaggeredCouplingTerm const& coupled_term,
std::function<void(unsigned, GlobalVector const&)> const& postIterationCallback)
{
namespace LinAlg = MathLib::LinAlg;
@@ -61,7 +62,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
BaseLib::RunTime time_assembly;
time_assembly.start();
- sys.assemble(x);
+ sys.assemble(x, coupled_term);
sys.getA(A);
sys.getRhs(rhs);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
@@ -161,16 +162,17 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
}
void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
- GlobalVector const& x) const
+ GlobalVector const& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term) const
{
- _equation_system->assemble(x);
+ _equation_system->assemble(x, coupled_term);
// TODO if the equation system would be reset to nullptr after each
// assemble() or solve() call, the user would be forced to set the
// equation every time and could not forget it.
}
bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
- GlobalVector& x,
+ GlobalVector& x, ProcessLib::StaggeredCouplingTerm const& coupled_term,
std::function<void(unsigned, GlobalVector const&)> const& postIterationCallback)
{
namespace LinAlg = MathLib::LinAlg;
@@ -204,7 +206,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
BaseLib::RunTime time_assembly;
time_assembly.start();
- sys.assemble(x);
+ sys.assemble(x, coupled_term);
sys.getResidual(x, res);
sys.getJacobian(J);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index c4a9be12b97..1eeb874bdf7 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -13,6 +13,8 @@
#include <utility>
#include <logog/include/logog.hpp>
+#include "ProcessLib/StaggeredCouplingTerm.h"
+
#include "ConvergenceCriterion.h"
#include "NonlinearSystem.h"
#include "Types.h"
@@ -42,7 +44,9 @@ public:
*
* \param x the state at which the equation system will be assembled.
*/
- virtual void assemble(GlobalVector const& x) const = 0;
+ virtual void assemble(GlobalVector const& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term
+ ) const = 0;
/*! Assemble and solve the equation system.
*
@@ -53,6 +57,7 @@ public:
* \retval false otherwise
*/
virtual bool solve(GlobalVector& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term,
std::function<void(unsigned, GlobalVector const&)> const&
postIterationCallback) = 0;
@@ -101,9 +106,12 @@ public:
_equation_system = &eq;
_convergence_criterion = &conv_crit;
}
- void assemble(GlobalVector const& x) const override;
+ void assemble(GlobalVector const& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term
+ ) const override;
bool solve(GlobalVector& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term,
std::function<void(unsigned, GlobalVector const&)> const&
postIterationCallback) override;
@@ -156,9 +164,12 @@ public:
_equation_system = &eq;
_convergence_criterion = &conv_crit;
}
- void assemble(GlobalVector const& x) const override;
+ void assemble(GlobalVector const& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term
+ ) const override;
bool solve(GlobalVector& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term,
std::function<void(unsigned, GlobalVector const&)> const&
postIterationCallback) override;
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index edfa53160ee..69b1c116977 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -9,6 +9,8 @@
#pragma once
+#include "ProcessLib/StaggeredCouplingTerm.h"
+
#include "EquationSystem.h"
#include "Types.h"
@@ -36,7 +38,9 @@ public:
//! Assembles the linearized equation system at the point \c x.
//! The linearized system is \f$A(x) \cdot x = b(x)\f$. Here the matrix
//! \f$A(x)\f$ and the vector \f$b(x)\f$ are assembled.
- virtual void assemble(GlobalVector const& x) = 0;
+ virtual void assemble(GlobalVector const& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term
+ ) = 0;
/*! Writes the residual at point \c x to \c res.
*
@@ -75,7 +79,9 @@ public:
//! Assembles the linearized equation system at the point \c x.
//! The linearized system is \f$J(x) \cdot \Delta x = (x)\f$. Here the
//! residual vector \f$r(x)\f$ and its Jacobian \f$J(x)\f$ are assembled.
- virtual void assemble(GlobalVector const& x) = 0;
+ virtual void assemble(GlobalVector const& x,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term
+ ) = 0;
//! Writes the linearized equation system matrix to \c A.
//! \pre assemble() must have been called before.
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index aafe142c7a5..f1a40fa7bb9 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -12,6 +12,8 @@
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "NumLib/IndexValueVector.h"
+#include "ProcessLib/StaggeredCouplingTerm.h"
+
#include "EquationSystem.h"
#include "Types.h"
@@ -47,7 +49,9 @@ public:
//! Assemble \c M, \c K and \c b at the provided state (\c t, \c x).
virtual void assemble(const double t, GlobalVector const& x,
GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b) = 0;
+ GlobalVector& b,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term)
+ = 0;
using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
@@ -115,7 +119,9 @@ public:
GlobalVector const& xdot,
const double dxdot_dx, const double dx_dx,
GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b, GlobalMatrix& Jac) = 0;
+ GlobalVector& b, GlobalMatrix& Jac,
+ ProcessLib::StaggeredCouplingTerm
+ const& coupled_term) = 0;
};
//! @}
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index a1b6f29685c..2298e5c0f99 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -67,7 +67,8 @@ TimeDiscretizedODESystem<
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
- assemble(const GlobalVector& x_new_timestep)
+ assemble(const GlobalVector& x_new_timestep,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term)
{
namespace LinAlg = MathLib::LinAlg;
@@ -85,7 +86,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
_Jac->setZero();
_ode.assembleWithJacobian(t, x_curr, xdot, dxdot_dx, dx_dx, *_M, *_K, *_b,
- *_Jac);
+ *_Jac, coupled_term);
LinAlg::finalizeAssembly(*_M);
LinAlg::finalizeAssembly(*_K);
@@ -175,7 +176,8 @@ TimeDiscretizedODESystem<
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Picard>::
- assemble(const GlobalVector& x_new_timestep)
+ assemble(const GlobalVector& x_new_timestep,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term)
{
namespace LinAlg = MathLib::LinAlg;
@@ -186,7 +188,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
_K->setZero();
_b->setZero();
- _ode.assemble(t, x_curr, *_M, *_K, *_b);
+ _ode.assemble(t, x_curr, *_M, *_K, *_b, coupled_term);
LinAlg::finalizeAssembly(*_M);
LinAlg::finalizeAssembly(*_K);
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 4a9856d0847..2db7ddc4edb 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -11,6 +11,8 @@
#include <memory>
+#include "ProcessLib/StaggeredCouplingTerm.h"
+
#include "MatrixTranslator.h"
#include "NonlinearSystem.h"
#include "ODESystem.h"
@@ -80,7 +82,9 @@ public:
~TimeDiscretizedODESystem();
- void assemble(const GlobalVector& x_new_timestep) override;
+ void assemble(const GlobalVector& x_new_timestep,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term)
+ override;
void getResidual(GlobalVector const& x_new_timestep,
GlobalVector& res) const override;
@@ -172,7 +176,9 @@ public:
~TimeDiscretizedODESystem();
- void assemble(const GlobalVector& x_new_timestep) override;
+ void assemble(const GlobalVector& x_new_timestep,
+ ProcessLib::StaggeredCouplingTerm const& coupled_term)
+ override;
void getA(GlobalMatrix& A) const override
{
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 82e75ec037b..4e69b98a00c 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -119,7 +119,8 @@ MathLib::MatrixSpecifications Process::getMatrixSpecifications() const
}
void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
- GlobalMatrix& K, GlobalVector& b)
+ GlobalMatrix& K, GlobalVector& b,
+ StaggeredCouplingTerm const& coupled_term)
{
assembleConcreteProcess(t, x, M, K, b);
_boundary_conditions.applyNaturalBC(t, x, K, b);
@@ -129,7 +130,8 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x,
GlobalVector const& xdot,
const double dxdot_dx, const double dx_dx,
GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b, GlobalMatrix& Jac)
+ GlobalVector& b, GlobalMatrix& Jac,
+ StaggeredCouplingTerm const& coupled_term)
{
assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac);
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 874be468a4c..44df54f18ca 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -22,6 +22,7 @@
#include "CachedSecondaryVariable.h"
#include "AbstractJacobianAssembler.h"
#include "ProcessType.h"
+#include "StaggeredCouplingTerm.h"
#include "VectorMatrixAssembler.h"
namespace MeshLib
@@ -74,13 +75,17 @@ public:
const override final;
void assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
- GlobalMatrix& K, GlobalVector& b) override final;
+ GlobalMatrix& K, GlobalVector& b,
+ StaggeredCouplingTerm const& coupled_term)
+ override final;
void assembleWithJacobian(const double t, GlobalVector const& x,
GlobalVector const& xdot, const double dxdot_dx,
const double dx_dx, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b,
- GlobalMatrix& Jac) override final;
+ GlobalMatrix& Jac,
+ StaggeredCouplingTerm const& coupled_term)
+ override final;
std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
getKnownSolutions(double const t) const override final
diff --git a/ProcessLib/StaggeredCouplingTerm.cpp b/ProcessLib/StaggeredCouplingTerm.cpp
new file mode 100644
index 00000000000..b588710308b
--- /dev/null
+++ b/ProcessLib/StaggeredCouplingTerm.cpp
@@ -0,0 +1,26 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file StaggeredCouplingTerm.cpp
+ *
+ * Created on November 7, 2016, 12:14 PM
+ */
+
+#include "StaggeredCouplingTerm.h"
+#include "Process.h"
+
+namespace ProcessLib
+{
+const StaggeredCouplingTerm createVoidStaggeredCouplingTerm()
+{
+ std::map<ProcessType, Process const&> coupled_pcsss;
+ std::map<ProcessType, GlobalVector const*> coupled_xs;
+ return StaggeredCouplingTerm(coupled_pcsss, coupled_xs);
+}
+
+} // end of ProcessLib
+
diff --git a/ProcessLib/StaggeredCouplingTerm.h b/ProcessLib/StaggeredCouplingTerm.h
new file mode 100644
index 00000000000..c7a36166d8e
--- /dev/null
+++ b/ProcessLib/StaggeredCouplingTerm.h
@@ -0,0 +1,42 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file StaggeredCouplingTerm.h
+ *
+ * Created on November 7, 2016, 12:14 PM
+ */
+
+#ifndef OGS_STAGGERED_COUPLING_TERM_H
+#define OGS_STAGGERED_COUPLING_TERM_H
+
+#include <map>
+
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
+
+#include "ProcessType.h"
+
+namespace ProcessLib
+{
+class Process;
+struct StaggeredCouplingTerm
+{
+ StaggeredCouplingTerm(
+ std::map<ProcessType, Process const&> const& coupled_processes_,
+ std::map<ProcessType, GlobalVector const*> const& coupled_xs_)
+ : coupled_processes(coupled_processes_), coupled_xs(coupled_xs_)
+ {
+ }
+
+ std::map<ProcessType, Process const&> const& coupled_processes;
+ std::map<ProcessType, GlobalVector const*> const& coupled_xs;
+};
+
+const StaggeredCouplingTerm createVoidStaggeredCouplingTerm();
+
+} // end of ProcessLib
+#endif /* OGS_STAGGERED_COUPLING_TERM_H */
+
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index 4640820555b..c7489c21488 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -15,6 +15,8 @@
#include "NumLib/ODESolver/ConvergenceCriterionPerComponent.h"
#include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
+#include "MathLib/LinAlg/LinAlg.h"
+
std::unique_ptr<NumLib::ITimeStepAlgorithm> createTimeStepper(
BaseLib::ConfigTree const& config)
{
@@ -51,7 +53,6 @@ std::unique_ptr<ProcessLib::Output> createOutput(
return ProcessLib::Output::newInstance(config, output_directory);
}
-
//! Sets the EquationSystem for the given nonlinear solver,
//! which is Picard or Newton depending on the NLTag.
template <NumLib::NonlinearSolverTag NLTag>
@@ -82,12 +83,10 @@ static void setEquationSystem(NumLib::NonlinearSolverBase& nonlinear_solver,
switch (nl_tag)
{
case Tag::Picard:
- setEquationSystem<Tag::Picard>(nonlinear_solver, eq_sys,
- conv_crit);
+ setEquationSystem<Tag::Picard>(nonlinear_solver, eq_sys, conv_crit);
break;
case Tag::Newton:
- setEquationSystem<Tag::Newton>(nonlinear_solver, eq_sys,
- conv_crit);
+ setEquationSystem<Tag::Newton>(nonlinear_solver, eq_sys, conv_crit);
break;
}
}
@@ -133,6 +132,7 @@ struct SingleProcessData
std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
Process& process_,
+ std::map<ProcessType, Process const&>&& coupled_processes_,
ProcessOutput&& process_output_);
SingleProcessData(SingleProcessData&& spd);
@@ -150,6 +150,7 @@ struct SingleProcessData
NumLib::InternalMatrixStorage* mat_strg = nullptr;
Process& process;
+ std::map<ProcessType, Process const&> const coupled_processes;
ProcessOutput process_output;
};
@@ -159,12 +160,14 @@ SingleProcessData::SingleProcessData(
std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
Process& process_,
+ std::map<ProcessType, Process const&>&& coupled_processes_,
ProcessOutput&& process_output_)
: nonlinear_solver_tag(NLTag),
nonlinear_solver(nonlinear_solver),
conv_crit(std::move(conv_crit_)),
time_disc(std::move(time_disc_)),
process(process_),
+ coupled_processes(coupled_processes_),
process_output(std::move(process_output_))
{
}
@@ -177,6 +180,7 @@ SingleProcessData::SingleProcessData(SingleProcessData&& spd)
tdisc_ode_sys(std::move(spd.tdisc_ode_sys)),
mat_strg(spd.mat_strg),
process(spd.process),
+ coupled_processes(spd.coupled_processes),
process_output(std::move(spd.process_output))
{
spd.mat_strg = nullptr;
@@ -240,6 +244,7 @@ std::unique_ptr<SingleProcessData> makeSingleProcessData(
Process& process,
std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit,
+ std::map<ProcessType, Process const&>&& coupled_processes,
ProcessOutput&& process_output)
{
using Tag = NumLib::NonlinearSolverTag;
@@ -250,7 +255,8 @@ std::unique_ptr<SingleProcessData> makeSingleProcessData(
{
return std::unique_ptr<SingleProcessData>{new SingleProcessData{
*nonlinear_solver_picard, std::move(conv_crit),
- std::move(time_disc), process, std::move(process_output)}};
+ std::move(time_disc), process, std::move(coupled_processes),
+ std::move(process_output)}};
}
else if (auto* nonlinear_solver_newton =
dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
@@ -258,16 +264,18 @@ std::unique_ptr<SingleProcessData> makeSingleProcessData(
{
return std::unique_ptr<SingleProcessData>{new SingleProcessData{
*nonlinear_solver_newton, std::move(conv_crit),
- std::move(time_disc), process, std::move(process_output)}};
- } else {
+ std::move(time_disc), process, std::move(coupled_processes),
+ std::move(process_output)}};
+ }
+ else
+ {
OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
}
}
std::vector<std::unique_ptr<SingleProcessData>> createPerProcessData(
BaseLib::ConfigTree const& config,
- const std::map<std::string, std::unique_ptr<Process>>&
- processes,
+ const std::map<std::string, std::unique_ptr<Process>>& processes,
std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>> const&
nonlinear_solvers)
{
@@ -297,12 +305,39 @@ std::vector<std::unique_ptr<SingleProcessData>> createPerProcessData(
//! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion}
pcs_config.getConfigSubtree("convergence_criterion"));
+ auto const& cpl_pcses
+ //! \ogs_file_param{prj__time_loop__processes__process__coupled_processes}
+ = pcs_config.getConfigSubtreeOptional("coupled_processes");
+ std::map<ProcessType, Process const&> coupled_processes;
+ if (cpl_pcses)
+ {
+ for (
+ auto const cpl_pcs_name :
+ //! \ogs_file_param{prj__time_loop__processes__process__coupled_processes__coupled_process}
+ cpl_pcses->getConfigParameterList<std::string>(
+ "coupled_process"))
+ {
+ auto& cpl_pcs_ptr = BaseLib::getOrError(
+ processes, cpl_pcs_name,
+ "A process with the given name has not been defined.");
+
+ auto const cpl_pcs = cpl_pcs_ptr.get();
+ auto const inserted = coupled_processes.emplace(
+ cpl_pcs_ptr->getProcessType(), *cpl_pcs);
+ if (!inserted.second)
+ { // insertion failed, i.e., key already exists
+ OGS_FATAL("Coupled process `%s' already exists.",
+ cpl_pcs_name.data());
+ }
+ }
+ }
+
//! \ogs_file_param{prj__time_loop__processes__process__output}
ProcessOutput process_output{pcs_config.getConfigSubtree("output")};
per_process_data.emplace_back(makeSingleProcessData(
nl_slv, pcs, std::move(time_disc), std::move(conv_crit),
- std::move(process_output)));
+ std::move(coupled_processes), std::move(process_output)));
}
if (per_process_data.size() != processes.size())
@@ -315,11 +350,26 @@ std::vector<std::unique_ptr<SingleProcessData>> createPerProcessData(
std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
BaseLib::ConfigTree const& config, std::string const& output_directory,
- const std::map<std::string, std::unique_ptr<Process>>&
- processes,
+ const std::map<std::string, std::unique_ptr<Process>>& processes,
const std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>>&
nonlinear_solvers)
{
+ auto const& coupling_config
+ //! \ogs_file_param{prj__time_loop__global_process_coupling}
+ = config.getConfigSubtreeOptional("global_process_coupling");
+
+ std::unique_ptr<NumLib::ConvergenceCriterion> coupling_conv_crit = nullptr;
+ unsigned max_coupling_iterations = 1;
+ if (coupling_config)
+ {
+ max_coupling_iterations
+ //! \ogs_file_param{prj__time_loop__global_process_coupling__max_iteration}
+ = coupling_config->getConfigParameter<unsigned>("max_iteration");
+ coupling_conv_crit = NumLib::createConvergenceCriterion(
+ //! \ogs_file_param{prj__time_loop__global_process_coupling__max_iteration__convergence_criterion}
+ coupling_config->getConfigSubtree("convergence_criterion"));
+ }
+
auto timestepper =
//! \ogs_file_param{prj__time_loop__time_stepping}
createTimeStepper(config.getConfigSubtree("time_stepping"));
@@ -333,9 +383,10 @@ std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
config.getConfigSubtree("processes"), processes, nonlinear_solvers);
return std::unique_ptr<UncoupledProcessesTimeLoop>{
- new UncoupledProcessesTimeLoop{std::move(timestepper),
- std::move(output),
- std::move(per_process_data)}};
+ new UncoupledProcessesTimeLoop{
+ std::move(timestepper), std::move(output),
+ std::move(per_process_data), max_coupling_iterations,
+ std::move(coupling_conv_crit)}};
}
std::vector<GlobalVector*> setInitialConditions(
@@ -371,7 +422,8 @@ std::vector<GlobalVector*> setInitialConditions(
auto& conv_crit = *spd->conv_crit;
setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
- nonlinear_solver.assemble(x0);
+ nonlinear_solver.assemble(
+ x0, ProcessLib::createVoidStaggeredCouplingTerm());
time_disc.pushState(
t0, x0, mat_strg); // TODO: that might do duplicate work
}
@@ -385,6 +437,7 @@ std::vector<GlobalVector*> setInitialConditions(
bool solveOneTimeStepOneProcess(GlobalVector& x, std::size_t const timestep,
double const t, double const delta_t,
SingleProcessData& process_data,
+ StaggeredCouplingTerm const& coupled_term,
Output const& output_control)
{
auto& process = process_data.process;
@@ -405,35 +458,90 @@ bool solveOneTimeStepOneProcess(GlobalVector& x, std::size_t const timestep,
applyKnownSolutions(ode_sys, nl_tag, x);
- process.preTimestep(x, t, delta_t);
-
- auto const post_iteration_callback = [&](
- unsigned iteration, GlobalVector const& x) {
+ auto const post_iteration_callback = [&](unsigned iteration,
+ GlobalVector const& x) {
output_control.doOutputNonlinearIteration(
process, process_data.process_output, timestep, t, x, iteration);
};
bool nonlinear_solver_succeeded =
- nonlinear_solver.solve(x, post_iteration_callback);
+ nonlinear_solver.solve(x, coupled_term, post_iteration_callback);
auto& mat_strg = *process_data.mat_strg;
time_disc.pushState(t, x, mat_strg);
- process.postTimestep(x);
-
return nonlinear_solver_succeeded;
}
UncoupledProcessesTimeLoop::UncoupledProcessesTimeLoop(
std::unique_ptr<NumLib::ITimeStepAlgorithm>&& timestepper,
std::unique_ptr<Output>&& output,
- std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data)
+ std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data,
+ const unsigned global_coupling_max_iterations,
+ std::unique_ptr<NumLib::ConvergenceCriterion>&& glb_coupling_conv_crit)
: _timestepper{std::move(timestepper)},
_output(std::move(output)),
- _per_process_data(std::move(per_process_data))
+ _per_process_data(std::move(per_process_data)),
+ _global_coupling_max_iterations(global_coupling_max_iterations),
+ _global_coupling_conv_crit(std::move(glb_coupling_conv_crit))
{
}
+bool UncoupledProcessesTimeLoop::setCoupledSolutions()
+{
+ if (!_global_coupling_conv_crit && _global_coupling_max_iterations == 1)
+ return false;
+
+ unsigned pcs_idx = 0;
+ _solutions_of_coupled_processes.resize(_per_process_data.size());
+ for (auto& spd : _per_process_data)
+ {
+ BaseLib::RunTime time_timestep_process;
+ time_timestep_process.start();
+
+ auto const& coupled_processes = spd->coupled_processes;
+ std::map<ProcessType, GlobalVector const*> coupled_xs;
+ auto it = coupled_processes.begin();
+ while (it != coupled_processes.end())
+ {
+ ProcessLib::Process const& coupled_process = it->second;
+ auto const found_item = std::find_if(
+ _per_process_data.begin(),
+ _per_process_data.end(),
+ [&coupled_process](
+ std::unique_ptr<SingleProcessData> const& item) {
+ return coupled_process.getProcessType() ==
+ item->process.getProcessType();
+ });
+
+ if (found_item != _per_process_data.end())
+ {
+ // Id of the coupled process:
+ const std::size_t c_id =
+ found_item - _per_process_data.begin();
+
+ BaseLib::insertMapIfKeyUniqueElseError(coupled_xs, it->first,
+ _process_solutions[c_id],
+ "global_coupled_x");
+ }
+ it++;
+ }
+ _solutions_of_coupled_processes[pcs_idx] = coupled_xs;
+
+ const auto x = _process_solutions[pcs_idx];
+
+ auto x1 = &NumLib::GlobalVectorProvider::provider.getVector(*x);
+ MathLib::LinAlg::copy(*x, *x1);
+
+ // append a solution vector of suitable size
+ _solutions_of_last_cpl_iteration.emplace_back(x1);
+
+ ++pcs_idx;
+ }
+
+ return true;
+}
+
bool UncoupledProcessesTimeLoop::loop()
{
// initialize output, convergence criterion, etc.
@@ -448,7 +556,8 @@ bool UncoupledProcessesTimeLoop::loop()
if (auto* conv_crit =
dynamic_cast<NumLib::ConvergenceCriterionPerComponent*>(
- spd->conv_crit.get())) {
+ spd->conv_crit.get()))
+ {
conv_crit->setDOFTable(pcs.getDOFTable(), pcs.getMesh());
}
@@ -475,6 +584,8 @@ bool UncoupledProcessesTimeLoop::loop()
}
}
+ const bool is_staggered_coupling = setCoupledSolutions();
+
double t = t0;
std::size_t timestep = 1; // the first timestep really is number one
bool nonlinear_solver_succeeded = true;
@@ -492,41 +603,109 @@ bool UncoupledProcessesTimeLoop::loop()
INFO("=== timestep #%u (t=%gs, dt=%gs) ==============================",
timestep, t, delta_t);
- // TODO use process name
- unsigned pcs_idx = 0;
- for (auto& spd : _per_process_data)
+ bool coupling_iteration_converged = true;
+ for (unsigned i = 0; i < _global_coupling_max_iterations; i++)
{
- auto& pcs = spd->process;
- BaseLib::RunTime time_timestep_process;
- time_timestep_process.start();
-
- auto& x = *_process_solutions[pcs_idx];
-
- nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
- x, timestep, t, delta_t, *spd, *_output);
-
- pcs.computeSecondaryVariable(t, x);
+ // TODO use process name
+ unsigned pcs_idx = 0;
+ for (auto& spd : _per_process_data)
+ {
+ auto& pcs = spd->process;
+ BaseLib::RunTime time_timestep_process;
+ time_timestep_process.start();
+
+ auto& x = *_process_solutions[pcs_idx];
+ if (i==0)
+ pcs.preTimestep(x, t, delta_t);
+
+ if (is_staggered_coupling)
+ {
+ if (i == 0)
+ _global_coupling_conv_crit->preFirstIteration();
+ else
+ _global_coupling_conv_crit->setNoFirstIteration();
+ StaggeredCouplingTerm coupled_term(
+ spd->coupled_processes,
+ _solutions_of_coupled_processes[pcs_idx]);
+
+ nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
+ x, timestep, t, delta_t, *spd, coupled_term, *_output);
+ }
+ else
+ {
+ nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
+ x, timestep, t, delta_t, *spd,
+ ProcessLib::createVoidStaggeredCouplingTerm(),
+ *_output);
+ pcs.postTimestep(x);
+ }
+
+ pcs.computeSecondaryVariable(t, x);
+
+ INFO(
+ "[time] Solving process #%u took %g s in timestep #%u "
+ "uncouling iteration #%u",
+ pcs_idx, time_timestep.elapsed(), timestep, i);
+
+ if (!nonlinear_solver_succeeded)
+ {
+ ERR("The nonlinear solver failed in timestep #%u at t = %g "
+ "s"
+ " for process #%u.",
+ timestep, t, pcs_idx);
+
+ // save unsuccessful solution
+ _output->doOutputAlways(pcs, spd->process_output, timestep,
+ t, x);
+
+ break;
+ }
+ else
+ {
+ if (!is_staggered_coupling)
+ _output->doOutput(pcs, spd->process_output, timestep, t,
+ x);
+ }
+
+ // Check the convergence of the coupling iteration
+ if (is_staggered_coupling)
+ {
+ auto& x_old = *_solutions_of_last_cpl_iteration[pcs_idx];
+ MathLib::LinAlg::axpy(x_old, -1.0, x);
+ _global_coupling_conv_crit->checkResidual(x_old);
+ coupling_iteration_converged =
+ coupling_iteration_converged &&
+ _global_coupling_conv_crit->isSatisfied();
+
+ if (coupling_iteration_converged)
+ break;
+ MathLib::LinAlg::copy(x, x_old);
+ }
+
+ ++pcs_idx;
+ }
- INFO("[time] Solving process #%u took %g s in timestep #%u.",
- pcs_idx, time_timestep.elapsed(), timestep);
+ if (is_staggered_coupling && coupling_iteration_converged)
+ break;
if (!nonlinear_solver_succeeded)
{
- ERR("The nonlinear solver failed in timestep #%u at t = %g s"
- " for process #%u.",
- timestep, t, pcs_idx);
-
- // save unsuccessful solution
- _output->doOutputAlways(pcs, spd->process_output, timestep, t, x);
-
break;
}
- else
+ }
+
+ if (is_staggered_coupling && coupling_iteration_converged)
+ {
+ unsigned pcs_idx = 0;
+ for (auto& spd : _per_process_data)
{
+ auto& pcs = spd->process;
+ auto& x = *_process_solutions[pcs_idx];
+ pcs.postTimestep(x);
+
_output->doOutput(pcs, spd->process_output, timestep, t, x);
+ ++pcs_idx;
}
-
- ++pcs_idx;
}
INFO("[time] Timestep #%u took %g s.", timestep,
@@ -544,7 +723,8 @@ bool UncoupledProcessesTimeLoop::loop()
{
auto& pcs = spd->process;
auto const& x = *_process_solutions[pcs_idx];
- _output->doOutputLastTimestep(pcs, spd->process_output, timestep, t, x);
+ _output->doOutputLastTimestep(pcs, spd->process_output, timestep, t,
+ x);
++pcs_idx;
}
@@ -557,6 +737,9 @@ UncoupledProcessesTimeLoop::~UncoupledProcessesTimeLoop()
{
for (auto* x : _process_solutions)
NumLib::GlobalVectorProvider::provider.releaseVector(*x);
+
+ for (auto* x : _solutions_of_last_cpl_iteration)
+ NumLib::GlobalVectorProvider::provider.releaseVector(*x);
}
} // namespace ProcessLib
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.h b/ProcessLib/UncoupledProcessesTimeLoop.h
index 3d6f74dde46..ef829722aeb 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.h
+++ b/ProcessLib/UncoupledProcessesTimeLoop.h
@@ -18,6 +18,11 @@
#include "Output.h"
#include "Process.h"
+namespace NumLib
+{
+class ConvergenceCriterion;
+}
+
namespace ProcessLib
{
struct SingleProcessData;
@@ -29,17 +34,30 @@ public:
explicit UncoupledProcessesTimeLoop(
std::unique_ptr<NumLib::ITimeStepAlgorithm>&& timestepper,
std::unique_ptr<Output>&& output,
- std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data);
+ std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data,
+ const unsigned global_coupling_max_iterations,
+ std::unique_ptr<NumLib::ConvergenceCriterion>&& glb_coupling_conv_crit);
bool loop();
~UncoupledProcessesTimeLoop();
+ bool setCoupledSolutions();
+
private:
std::vector<GlobalVector*> _process_solutions;
std::unique_ptr<NumLib::ITimeStepAlgorithm> _timestepper;
std::unique_ptr<Output> _output;
std::vector<std::unique_ptr<SingleProcessData>> _per_process_data;
+
+ /// Maximum iterations of the global coupling.
+ const unsigned _global_coupling_max_iterations;
+ /// Convergence criteria of the global coupling iterations.
+ std::unique_ptr<NumLib::ConvergenceCriterion> _global_coupling_conv_crit;
+ std::vector<std::map<ProcessType, GlobalVector const*>>
+ _solutions_of_coupled_processes;
+ /// Solutions of the previous coupling iteration.
+ std::vector<GlobalVector*> _solutions_of_last_cpl_iteration;
};
//! Builds an UncoupledProcessesTimeLoop from the given configuration.
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index 848b62cda93..da01e8c8c06 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -13,6 +13,7 @@
#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
#include "NumLib/ODESolver/ODESystem.h"
+#include "ProcessLib/StaggeredCouplingTerm.h"
// debug
//#include <iostream>
@@ -27,7 +28,9 @@ class ODE1 final : public NumLib::ODESystem<
{
public:
void assemble(const double /*t*/, GlobalVector const& /*x*/,
- GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+ ProcessLib::StaggeredCouplingTerm const& /*coupled_term*/
+ ) override
{
MathLib::setMatrix(M, { 1.0, 0.0, 0.0, 1.0 });
MathLib::setMatrix(K, { 0.0, 1.0, -1.0, 0.0 });
@@ -39,11 +42,13 @@ public:
GlobalVector const& /*xdot*/, const double dxdot_dx,
const double dx_dx, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b,
- GlobalMatrix& Jac) override
+ GlobalMatrix& Jac,
+ ProcessLib::StaggeredCouplingTerm
+ const& coupled_term) override
{
namespace LinAlg = MathLib::LinAlg;
- assemble(t, x_curr, M, K, b);
+ assemble(t, x_curr, M, K, b, coupled_term);
// compute Jac = M*dxdot_dx + dx_dx*K
LinAlg::finalizeAssembly(M);
@@ -103,7 +108,9 @@ class ODE2 final : public NumLib::ODESystem<
{
public:
void assemble(const double /*t*/, GlobalVector const& x, GlobalMatrix& M,
- GlobalMatrix& K, GlobalVector& b) override
+ GlobalMatrix& K, GlobalVector& b,
+ ProcessLib::StaggeredCouplingTerm const& /*coupled_term*/
+ ) override
{
MathLib::setMatrix(M, {1.0});
MathLib::setMatrix(K, {x[0]});
@@ -114,9 +121,11 @@ public:
GlobalVector const& /*xdot*/,
const double dxdot_dx, const double dx_dx,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
- GlobalMatrix& Jac) override
+ GlobalMatrix& Jac,
+ ProcessLib::StaggeredCouplingTerm const&
+ coupled_term) override
{
- assemble(t, x, M, K, b);
+ assemble(t, x, M, K, b, coupled_term);
namespace LinAlg = MathLib::LinAlg;
@@ -187,7 +196,9 @@ class ODE3 final : public NumLib::ODESystem<
{
public:
void assemble(const double /*t*/, GlobalVector const& x_curr, GlobalMatrix& M,
- GlobalMatrix& K, GlobalVector& b) override
+ GlobalMatrix& K, GlobalVector& b,
+ ProcessLib::StaggeredCouplingTerm const& /*coupled_term*/
+ ) override
{
auto const u = x_curr[0];
auto const v = x_curr[1];
@@ -202,9 +213,11 @@ public:
GlobalVector const& xdot, const double dxdot_dx,
const double dx_dx, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b,
- GlobalMatrix& Jac) override
+ GlobalMatrix& Jac,
+ ProcessLib::StaggeredCouplingTerm const&
+ coupled_term) override
{
- assemble(t, x_curr, M, K, b);
+ assemble(t, x_curr, M, K, b, coupled_term);
auto const u = x_curr[0];
auto const v = x_curr[1];
diff --git a/Tests/NumLib/TimeLoopSingleODE.h b/Tests/NumLib/TimeLoopSingleODE.h
index 1c67bfe6edc..531ef5432da 100644
--- a/Tests/NumLib/TimeLoopSingleODE.h
+++ b/Tests/NumLib/TimeLoopSingleODE.h
@@ -13,6 +13,7 @@
#include "NumLib/ODESolver/TimeDiscretizedODESystem.h"
#include "NumLib/ODESolver/NonlinearSolver.h"
+#include "ProcessLib/StaggeredCouplingTerm.h"
namespace NumLib
{
@@ -96,7 +97,8 @@ bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
if (time_disc.needsPreload())
{
- _nonlinear_solver->assemble(x);
+ _nonlinear_solver->assemble(x,
+ ProcessLib::createVoidStaggeredCouplingTerm());
time_disc.pushState(t0, x0,
_ode_sys); // TODO: that might do duplicate work
}
@@ -112,7 +114,9 @@ bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
// INFO("time: %e, delta_t: %e", t, delta_t);
time_disc.nextTimestep(t, delta_t);
- nl_slv_succeeded = _nonlinear_solver->solve(x, nullptr);
+ nl_slv_succeeded = _nonlinear_solver->solve(x,
+ ProcessLib::createVoidStaggeredCouplingTerm(),
+ nullptr);
if (!nl_slv_succeeded)
break;
--
GitLab