diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
index fe49dfeb5a6fcea9687f833de54c049d23f70a26..b8375887d2d80a72ee763554607da7fe0789d619 100644
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
+++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
@@ -197,6 +197,33 @@ private:
break;
}
}
+
+ //! Applies known solutions to the solution vector \c x, transparently
+ //! for equation systems linearized with either the Picard or Newton method.
+ template<NumLib::NonlinearSolverTag NLTag>
+ static void applyKnownSolutions(
+ EquationSystem const& eq_sys, Vector& x)
+ {
+ using EqSys = NumLib::NonlinearSystem<Matrix, Vector, NLTag>;
+ assert(dynamic_cast<EqSys const*> (&eq_sys) != nullptr);
+ auto& eq_sys_ = static_cast<EqSys const&> (eq_sys);
+
+ eq_sys_.applyKnownSolutions(x);
+ }
+
+ //! Applies known solutions to the solution vector \c x, transparently
+ //! for equation systems linearized with either the Picard or Newton method.
+ static void applyKnownSolutions(
+ EquationSystem const& eq_sys,
+ NumLib::NonlinearSolverTag const nl_tag, Vector& x)
+ {
+ using Tag = NumLib::NonlinearSolverTag;
+ switch (nl_tag)
+ {
+ case Tag::Picard: applyKnownSolutions<Tag::Picard>(eq_sys, x); break;
+ case Tag::Newton: applyKnownSolutions<Tag::Newton>(eq_sys, x); break;
+ }
+ }
};
//! Builds an UncoupledProcessesTimeLoop from the given configuration.
@@ -310,11 +337,17 @@ solveOneTimeStepOneProcess(
setEquationSystem(nonlinear_solver, ode_sys, nl_tag);
- process.preTimestep(x, t, delta_t);
+ // Note: Order matters!
+ // First advance to the next timestep, then set known solutions at that
+ // time, afterwards pass the right solution vector and time to the
+ // preTimestep() hook.
- // INFO("time: %e, delta_t: %e", t, delta_t);
time_disc.nextTimestep(t, delta_t);
+ applyKnownSolutions(ode_sys, nl_tag, x);
+
+ process.preTimestep(x, t, delta_t);
+
bool nonlinear_solver_succeeded = nonlinear_solver.solve(x);
auto& mat_strg = process_data.mat_strg;
diff --git a/MathLib/LinAlg/MatrixVectorTraits.cpp b/MathLib/LinAlg/MatrixVectorTraits.cpp
index d1036ed4347d240bd323dd387d49741e5f10e2c3..dceca5ee11c47bc09b40634d18e577a643ac1238 100644
--- a/MathLib/LinAlg/MatrixVectorTraits.cpp
+++ b/MathLib/LinAlg/MatrixVectorTraits.cpp
@@ -8,6 +8,7 @@
*/
#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "MatrixProviderUser.h"
#include "MatrixVectorTraits.h"
#ifdef OGS_USE_EIGEN
diff --git a/MathLib/LinAlg/MatrixVectorTraits.h b/MathLib/LinAlg/MatrixVectorTraits.h
index 793e6c2bf823c2ba16c6f2e2c8fc257bd481759e..44f9f4892fa2335b4f005077bbc6ce344eb5978a 100644
--- a/MathLib/LinAlg/MatrixVectorTraits.h
+++ b/MathLib/LinAlg/MatrixVectorTraits.h
@@ -12,12 +12,12 @@
#include<memory>
-#include "MatrixProviderUser.h"
-
namespace MathLib
{
template<typename Matrix>
struct MatrixVectorTraits;
+
+struct MatrixSpecifications;
}
#define SPECIALIZE_MATRIX_VECTOR_TRAITS(MATVEC, IDX) \
diff --git a/MathLib/LinAlg/UnifiedMatrixSetters.cpp b/MathLib/LinAlg/UnifiedMatrixSetters.cpp
index 3b762f0a985ff111dca78236e68d45f2e0e2ea2f..00dbc7ce0ab9cd56fcb7613ac58b5b9cbfb8fa1b 100644
--- a/MathLib/LinAlg/UnifiedMatrixSetters.cpp
+++ b/MathLib/LinAlg/UnifiedMatrixSetters.cpp
@@ -67,6 +67,12 @@ void setVector(Eigen::VectorXd& v, std::initializer_list<double> values)
for (std::size_t i=0; i<values.size(); ++i) v[i] = *(it++);
}
+void setVector(Eigen::VectorXd& v, MatrixVectorTraits<Eigen::VectorXd>::Index const index,
+ double const value)
+{
+ v[index] = value;
+}
+
} // namespace MathLib
#endif // OGS_USE_EIGEN
@@ -97,6 +103,12 @@ void setVector(PETScVector& v,
v.set(idcs, vals);
}
+void setVector(PETScVector& v, MatrixVectorTraits<PETScVector>::Index const index,
+ double const value)
+{
+ v.set(index, value); // TODO handle negative indices
+}
+
void setMatrix(PETScMatrix& m,
std::initializer_list<double> values)
{
@@ -173,6 +185,13 @@ void setVector(EigenVector& v,
setVector(v.getRawVector(), values);
}
+void setVector(EigenVector& v, MatrixVectorTraits<EigenVector>::Index const index,
+ double const value)
+{
+ v.getRawVector()[index] = value;
+}
+
+
void setMatrix(EigenMatrix& m,
std::initializer_list<double> values)
{
diff --git a/MathLib/LinAlg/UnifiedMatrixSetters.h b/MathLib/LinAlg/UnifiedMatrixSetters.h
index 4ece2dd9400627a442b4c44f95c585baee9ce50b..d7402b25139b7ade1185489fba30248b763d683e 100644
--- a/MathLib/LinAlg/UnifiedMatrixSetters.h
+++ b/MathLib/LinAlg/UnifiedMatrixSetters.h
@@ -13,6 +13,7 @@
#define MATHLIB_UNIFIED_MATRIX_SETTERS_H
#include <initializer_list>
+#include "MatrixVectorTraits.h"
#ifdef OGS_USE_EIGEN
@@ -35,6 +36,9 @@ double norm(Eigen::VectorXd const& x);
void setVector(Eigen::VectorXd& v, std::initializer_list<double> values);
+void setVector(Eigen::VectorXd& v, MatrixVectorTraits<Eigen::VectorXd>::Index const index,
+ double const value);
+
} // namespace MathLib
#endif // OGS_USE_EIGEN
@@ -55,6 +59,9 @@ double norm(PETScVector const& x);
void setVector(PETScVector& v,
std::initializer_list<double> values);
+void setVector(PETScVector& v, MatrixVectorTraits<PETScVector>::Index const index,
+ double const value);
+
void setMatrix(PETScMatrix& m, Eigen::MatrixXd const& tmp);
void addToMatrix(PETScMatrix& m,
@@ -79,6 +86,9 @@ class EigenMatrix;
void setVector(EigenVector& v,
std::initializer_list<double> values);
+void setVector(EigenVector& v, MatrixVectorTraits<EigenVector>::Index const index,
+ double const value);
+
void setMatrix(EigenMatrix& m,
std::initializer_list<double> values);
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index d466da77e121c99a781dc86007a710bef5ea2c13..17c9f169cd18e6dbd7544697271777468f699624 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -64,6 +64,9 @@ public:
*/
virtual void getJacobian(Matrix& Jac) const = 0;
+ //! Apply known solutions to the solution vector \c x.
+ virtual void applyKnownSolutions(Vector& x) const = 0;
+
//! Apply known solutions to the linearized equation system
//! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
virtual void applyKnownSolutionsNewton(
@@ -93,6 +96,9 @@ public:
//! Writes the linearized equation system right-hand side to \c rhs.
virtual void getRhs(Vector& rhs) const = 0;
+ //! Apply known solutions to the solution vector \c x.
+ virtual void applyKnownSolutions(Vector& x) const = 0;
+
//! Apply known solutions to the linearized equation system
//! \f$ A \cdot x = \mathit{rhs} \f$.
virtual void applyKnownSolutionsPicard(
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 5acc2d6e26eca63f59e4053eacdd23ce716841a7..22d1d1414737fe100e8114ffe547a12eef410edd 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -13,6 +13,7 @@
#include <memory>
#include "MathLib/LinAlg/ApplyKnownSolution.h"
+#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
#include "ProcessLib/DirichletBc.h"
#include "ODESystem.h"
@@ -20,6 +21,23 @@
#include "TimeDiscretization.h"
#include "MatrixTranslator.h"
+namespace detail
+{
+//! Applies known solutions to the solution vector \c x.
+template<typename Solutions, typename Vector>
+void applyKnownSolutions(std::vector<Solutions> const*const known_solutions,
+ Vector& x)
+{
+ if (!known_solutions) return;
+
+ for (auto const& bc : *known_solutions) {
+ for (std::size_t i=0; i<bc.global_ids.size(); ++i) {
+ // TODO that might have bad performance for some Vector types, e.g., PETSc.
+ MathLib::setVector(x, bc.global_ids[i], bc.values[i]);
+ }
+ }
+}
+}
namespace NumLib
{
@@ -174,11 +192,29 @@ public:
_mat_trans->computeJacobian(*_Jac, Jac);
}
+ void applyKnownSolutions(Vector& x) const override
+ {
+ ::detail::applyKnownSolutions(
+ _ode.getKnownSolutions(_time_disc.getCurrentTime()), x);
+ }
+
void applyKnownSolutionsNewton(Matrix& Jac, Vector& res,
Vector& minus_delta_x) override
{
- (void) Jac; (void) res; (void) minus_delta_x;
- INFO("Method applyKnownSolutionsNewton() not implemented."); // TODO implement
+ auto const* known_solutions =
+ _ode.getKnownSolutions(_time_disc.getCurrentTime());
+
+ if (known_solutions) {
+ std::vector<double> values;
+
+ for (auto const& bc : *known_solutions) {
+ // TODO this is the quick and dirty and bad performance solution.
+ values.resize(bc.values.size(), 0.0);
+
+ // TODO maybe it would be faster to apply all at once
+ MathLib::applyKnownSolution(Jac, res, minus_delta_x, bc.global_ids, values);
+ }
+ }
}
bool isLinear() const override
@@ -308,6 +344,12 @@ public:
_mat_trans->computeRhs(*_M, *_K, *_b, rhs);
}
+ void applyKnownSolutions(Vector& x) const override
+ {
+ ::detail::applyKnownSolutions(
+ _ode.getKnownSolutions(_time_disc.getCurrentTime()), x);
+ }
+
void applyKnownSolutionsPicard(Matrix& A, Vector& rhs, Vector& x) override
{
auto const* known_solutions =