diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index f1e4c472c962cfe950d9d76d941f3ae5da8517a7..ba32dd836d7c5b19c54abe5f20dd98f2686609af 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -209,7 +209,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
BaseLib::RunTime time_dirichlet;
time_dirichlet.start();
- sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
+ sys.applyKnownSolutionsNewton(J, res, minus_delta_x, x);
INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
@@ -230,8 +230,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
// cf.
// http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
auto& x_new =
- NumLib::GlobalVectorProvider::provider.getVector(
- x, _x_new_id);
+ NumLib::GlobalVectorProvider::provider.getVector(x, _x_new_id);
LinAlg::axpy(x_new, -_damping, minus_delta_x);
if (postIterationCallback)
@@ -252,16 +251,14 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
"iteration"
" has to be repeated.");
// TODO introduce some onDestroy hook.
- NumLib::GlobalVectorProvider::provider
- .releaseVector(x_new);
+ NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
continue; // That throws the iteration result away.
}
// TODO could be done via swap. Note: that also requires swapping
// the ids. Same for the Picard scheme.
LinAlg::copy(x_new, x); // copy new solution to x
- NumLib::GlobalVectorProvider::provider.releaseVector(
- x_new);
+ NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
}
if (!iteration_succeeded)
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index 91269cb7151948a6d1006ce9919977e880651880..b3bf35c59a6a99be5ab17ef1d0f69226b62f9efc 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -59,7 +59,8 @@ public:
//! Apply known solutions to the linearized equation system
//! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
virtual void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x) = 0;
+ GlobalVector& minus_delta_x,
+ GlobalVector& x) = 0;
};
/*! A System of nonlinear equations to be solved with the Picard fixpoint
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index 8d16dd2b235dc9848eb252feddbef614267eaada..64005e7206fbdf952cf05678ab4030cd3974b0d7 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -48,18 +48,19 @@ public:
virtual void preAssemble(const double t, GlobalVector const& x) = 0;
//! 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;
+ virtual void assemble(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) = 0;
using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
//! Provides known solutions (Dirichlet boundary conditions) vector for
//! the ode system at the given time \c t.
virtual std::vector<NumLib::IndexValueVector<Index>> const*
- getKnownSolutions(double const t) const
+ getKnownSolutions(double const t, GlobalVector const& x) const
{
(void)t;
+ (void)x;
return nullptr; // by default there are no known solutions
}
};
@@ -76,7 +77,6 @@ class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Picard>
{
public:
-
//! Calls process' pre-assembly with the provided state (\c t, \c x).
virtual void preAssemble(const double t, GlobalVector const& x) = 0;
@@ -124,11 +124,12 @@ public:
* \mathtt{Jac} \f$.
* \endparblock
*/
- virtual 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) = 0;
+ virtual 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) = 0;
};
//! @}
-}
+} // namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 5c597a88e5f55e557feb2cf37291a4d7d642ae02..617c722a65fe4f41124d7b7712fcf457b35bb1e7 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -34,14 +34,14 @@ void applyKnownSolutions(std::vector<Solutions> const* const known_solutions,
}
MathLib::LinAlg::finalizeAssembly(x);
}
-}
+} // namespace detail
namespace NumLib
{
TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
- TimeDiscretizedODESystem(const int process_id,
- ODE& ode, TimeDisc& time_discretization)
+ TimeDiscretizedODESystem(const int process_id, ODE& ode,
+ TimeDisc& time_discretization)
: _ode(ode),
_time_disc(time_discretization),
_mat_trans(createMatrixTranslator<ODETag>(time_discretization))
@@ -66,9 +66,9 @@ TimeDiscretizedODESystem<
NumLib::GlobalVectorProvider::provider.releaseVector(*_b);
}
-void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Newton>::
- assemble(const GlobalVector& x_new_timestep)
+void TimeDiscretizedODESystem<
+ ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Newton>::assemble(const GlobalVector& x_new_timestep)
{
namespace LinAlg = MathLib::LinAlg;
@@ -125,16 +125,16 @@ void TimeDiscretizedODESystem<
NonlinearSolverTag::Newton>::applyKnownSolutions(GlobalVector& x) const
{
::detail::applyKnownSolutions(
- _ode.getKnownSolutions(_time_disc.getCurrentTime()), x);
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
}
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x)
+ GlobalVector& minus_delta_x, GlobalVector& x)
{
auto const* known_solutions =
- _ode.getKnownSolutions(_time_disc.getCurrentTime());
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
if (!known_solutions || known_solutions->empty())
return;
@@ -149,6 +149,8 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
// For the Newton method the values must be zero
std::vector<double> values(ids.size(), 0);
MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
+
+ ::detail::applyKnownSolutions(known_solutions, x);
}
TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
@@ -176,9 +178,9 @@ TimeDiscretizedODESystem<
NumLib::GlobalVectorProvider::provider.releaseVector(*_b);
}
-void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Picard>::
- assemble(const GlobalVector& x_new_timestep)
+void TimeDiscretizedODESystem<
+ ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Picard>::assemble(const GlobalVector& x_new_timestep)
{
namespace LinAlg = MathLib::LinAlg;
@@ -202,7 +204,7 @@ void TimeDiscretizedODESystem<
NonlinearSolverTag::Picard>::applyKnownSolutions(GlobalVector& x) const
{
::detail::applyKnownSolutions(
- _ode.getKnownSolutions(_time_disc.getCurrentTime()), x);
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
}
void TimeDiscretizedODESystem<
@@ -212,7 +214,7 @@ void TimeDiscretizedODESystem<
GlobalVector& x)
{
auto const* known_solutions =
- _ode.getKnownSolutions(_time_disc.getCurrentTime());
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
if (known_solutions)
{
@@ -229,4 +231,4 @@ void TimeDiscretizedODESystem<
}
}
-} // NumLib
+} // namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 30e8ab1948be14bf453a3963ead8c43a36a792f1..59ead507782483e9247ee37c7fd2d91f13e20287 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -92,7 +92,8 @@ public:
void applyKnownSolutions(GlobalVector& x) const override;
void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x) override;
+ GlobalVector& minus_delta_x,
+ GlobalVector& x) override;
bool isLinear() const override
{