[NL] pass solution to getKnownSolutions()

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)

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

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

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

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
     {