diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index cf1d3202ceba566957b2f74dd1635d15ba75c251..c3216b1cc56db1714c0b3ab62148e3a73840d122 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -22,13 +22,13 @@
namespace NumLib
{
void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
- GlobalVector const& x, int const process_id) const
+ std::vector<GlobalVector*> const& x, int const process_id) const
{
_equation_system->assemble(x, process_id);
}
NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
- GlobalVector& x,
+ std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const& postIterationCallback,
int const process_id)
{
@@ -39,11 +39,13 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(
_rhs_id);
- auto& x_new = NumLib::GlobalVectorProvider::provider.getVector(_x_new_id);
- bool error_norms_met = false;
+ std::vector<GlobalVector*> x_new{x};
+ x_new[process_id] =
+ &NumLib::GlobalVectorProvider::provider.getVector(_x_new_id);
+ LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
- LinAlg::copy(x, x_new); // set initial guess
+ bool error_norms_met = false;
_convergence_criterion->preFirstIteration();
@@ -58,11 +60,11 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
time_iteration.start();
timer_dirichlet.start();
- sys.computeKnownSolutions(x_new, process_id);
- sys.applyKnownSolutions(x_new);
+ sys.computeKnownSolutions(*x_new[process_id], process_id);
+ sys.applyKnownSolutions(*x_new[process_id]);
time_dirichlet += timer_dirichlet.elapsed();
- sys.preIteration(iteration, x_new);
+ sys.preIteration(iteration, *x_new[process_id]);
BaseLib::RunTime time_assembly;
time_assembly.start();
@@ -72,20 +74,21 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
timer_dirichlet.start();
- sys.applyKnownSolutionsPicard(A, rhs, x_new);
+ sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
time_dirichlet += timer_dirichlet.elapsed();
INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
GlobalVector res;
- LinAlg::matMult(A, x_new, res); // res = A * x_new
+ LinAlg::matMult(A, *x_new[process_id], res); // res = A * x_new
LinAlg::axpy(res, -1.0, rhs); // res -= rhs
_convergence_criterion->checkResidual(res);
}
BaseLib::RunTime time_linear_solver;
time_linear_solver.start();
- bool iteration_succeeded = _linear_solver.solve(A, rhs, x_new);
+ bool iteration_succeeded =
+ _linear_solver.solve(A, rhs, *x_new[process_id]);
INFO("[time] Linear solver took %g s.", time_linear_solver.elapsed());
if (!iteration_succeeded)
@@ -96,10 +99,10 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
{
if (postIterationCallback)
{
- postIterationCallback(iteration, x_new);
+ postIterationCallback(iteration, *x_new[process_id]);
}
- switch (sys.postIteration(x_new))
+ switch (sys.postIteration(*x_new[process_id]))
{
case IterationResult::SUCCESS:
// Don't copy here. The old x might still be used further
@@ -112,13 +115,15 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
// Copy new solution to x.
// Thereby the failed solution can be used by the caller for
// debugging purposes.
- LinAlg::copy(x_new, x);
+ LinAlg::copy(*x_new[process_id], *x[process_id]);
break;
case IterationResult::REPEAT_ITERATION:
INFO(
"Picard: The postIteration() hook decided that this "
"iteration has to be repeated.");
- LinAlg::copy(x, x_new); // throw the iteration result away
+ LinAlg::copy(
+ *x[process_id],
+ *x_new[process_id]); // throw the iteration result away
continue;
}
}
@@ -134,17 +139,18 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
error_norms_met = true;
} else {
if (_convergence_criterion->hasDeltaXCheck()) {
- GlobalVector minus_delta_x(x);
+ GlobalVector minus_delta_x(*x[process_id]);
LinAlg::axpy(minus_delta_x, -1.0,
- x_new); // minus_delta_x = x - x_new
- _convergence_criterion->checkDeltaX(minus_delta_x, x_new);
+ *x_new[process_id]); // minus_delta_x = x - x_new
+ _convergence_criterion->checkDeltaX(minus_delta_x,
+ *x_new[process_id]);
}
error_norms_met = _convergence_criterion->isSatisfied();
}
// Update x s.t. in the next iteration we will compute the right delta x
- LinAlg::copy(x_new, x);
+ LinAlg::copy(*x_new[process_id], *x[process_id]);
INFO("[time] Iteration #%u took %g s.", iteration,
time_iteration.elapsed());
@@ -171,13 +177,13 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
NumLib::GlobalMatrixProvider::provider.releaseMatrix(A);
NumLib::GlobalVectorProvider::provider.releaseVector(rhs);
- NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
+ NumLib::GlobalVectorProvider::provider.releaseVector(*x_new[process_id]);
return {error_norms_met, iteration};
}
void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
- GlobalVector const& x, int const process_id) const
+ std::vector<GlobalVector*> const& x, int const process_id) const
{
_equation_system->assemble(x, process_id);
// TODO if the equation system would be reset to nullptr after each
@@ -186,7 +192,7 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
}
NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
- GlobalVector& x,
+ std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const& postIterationCallback,
int const process_id)
{
@@ -205,7 +211,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
// TODO be more efficient
// init minus_delta_x to the right size
- LinAlg::copy(x, minus_delta_x);
+ LinAlg::copy(*x[process_id], minus_delta_x);
_convergence_criterion->preFirstIteration();
@@ -220,16 +226,16 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
time_iteration.start();
timer_dirichlet.start();
- sys.computeKnownSolutions(x, process_id);
- sys.applyKnownSolutions(x);
+ sys.computeKnownSolutions(*x[process_id], process_id);
+ sys.applyKnownSolutions(*x[process_id]);
time_dirichlet += timer_dirichlet.elapsed();
- sys.preIteration(iteration, x);
+ sys.preIteration(iteration, *x[process_id]);
BaseLib::RunTime time_assembly;
time_assembly.start();
sys.assemble(x, process_id);
- sys.getResidual(x, res);
+ sys.getResidual(*x[process_id], res);
sys.getJacobian(J);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
@@ -259,8 +265,8 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
// TODO could be solved in a better way
// 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);
+ auto& x_new = NumLib::GlobalVectorProvider::provider.getVector(
+ *x[process_id], _x_new_id);
LinAlg::axpy(x_new, -_damping, minus_delta_x);
if (postIterationCallback)
@@ -287,7 +293,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
continue; // That throws the iteration result away.
}
- LinAlg::copy(x_new, x); // copy new solution to x
+ LinAlg::copy(x_new, *x[process_id]); // copy new solution to x
NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
}
@@ -303,7 +309,8 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
} else {
if (_convergence_criterion->hasDeltaXCheck()) {
// Note: x contains the new solution!
- _convergence_criterion->checkDeltaX(minus_delta_x, x);
+ _convergence_criterion->checkDeltaX(minus_delta_x,
+ *x[process_id]);
}
error_norms_met = _convergence_criterion->isSatisfied();
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 4ec739bbcc01e09c87a57691d271114e2b83638d..bcffe0a4a476e2cbd6a2b108d461396ec2847c57 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -46,7 +46,7 @@ public:
* assembled.
* \param process_id the process' id which will be assembled.
*/
- virtual void assemble(GlobalVector const& x,
+ virtual void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const = 0;
/*! Assemble and solve the equation system.
@@ -59,7 +59,7 @@ public:
* \retval false otherwise
*/
virtual NonlinearSolverStatus solve(
- GlobalVector& x,
+ std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const&
postIterationCallback,
int const process_id) = 0;
@@ -110,10 +110,11 @@ public:
_convergence_criterion = &conv_crit;
}
- void assemble(GlobalVector const& x, int const process_id) const override;
+ void assemble(std::vector<GlobalVector*> const& x,
+ int const process_id) const override;
NonlinearSolverStatus solve(
- GlobalVector& x,
+ std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const&
postIterationCallback,
int const process_id) override;
@@ -171,10 +172,11 @@ public:
_convergence_criterion = &conv_crit;
}
- void assemble(GlobalVector const& x, int const process_id) const override;
+ void assemble(std::vector<GlobalVector*> const& x,
+ int const process_id) const override;
NonlinearSolverStatus solve(
- GlobalVector& x,
+ std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const&
postIterationCallback,
int const process_id) override;
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index d5367eae3cd81075e24b247a6af0031983de2cc7..196b00b7def03abb79e32ea6d9e3d81f67e35828 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -37,7 +37,8 @@ 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, int const process_id) = 0;
+ virtual void assemble(std::vector<GlobalVector*> const& x,
+ int const process_id) = 0;
/*! Writes the residual at point \c x to \c res.
*
@@ -83,7 +84,8 @@ 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, int const process_id) = 0;
+ virtual void assemble(std::vector<GlobalVector*> const& x,
+ int const process_id) = 0;
//! Writes the linearized equation system matrix to \c A.
//! \pre assemble() must have been called before.
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index d57c682f41219d6809c587f72a682ec3fb388506..dd45a353b0bd2af70967df9387d2abfcb560f919 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -69,21 +69,21 @@ TimeDiscretizedODESystem<
NumLib::GlobalVectorProvider::provider.releaseVector(*_b);
}
-void TimeDiscretizedODESystem<
- ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Newton>::assemble(const GlobalVector& x_new_timestep,
- int const process_id)
+void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Newton>::
+ assemble(std::vector<GlobalVector*> const& x_new_timestep,
+ int const process_id)
{
namespace LinAlg = MathLib::LinAlg;
auto const t = _time_disc.getCurrentTime();
auto const dt = _time_disc.getCurrentTimeIncrement();
- auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
+ auto const& x_curr = _time_disc.getCurrentX(*x_new_timestep[process_id]);
auto const dxdot_dx = _time_disc.getNewXWeight();
auto const dx_dx = _time_disc.getDxDx();
auto& xdot = NumLib::GlobalVectorProvider::provider.getVector(_xdot_id);
- _time_disc.getXdot(x_new_timestep, xdot);
+ _time_disc.getXdot(*x_new_timestep[process_id], xdot);
_M->setZero();
_K->setZero();
@@ -188,16 +188,16 @@ TimeDiscretizedODESystem<
NumLib::GlobalVectorProvider::provider.releaseVector(*_b);
}
-void TimeDiscretizedODESystem<
- ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Picard>::assemble(const GlobalVector& x_new_timestep,
- int const process_id)
+void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Picard>::
+ assemble(std::vector<GlobalVector*> const& x_new_timestep,
+ int const process_id)
{
namespace LinAlg = MathLib::LinAlg;
auto const t = _time_disc.getCurrentTime();
auto const dt = _time_disc.getCurrentTimeIncrement();
- auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
+ auto const& x_curr = _time_disc.getCurrentX(*x_new_timestep[process_id]);
_M->setZero();
_K->setZero();
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 9fa14e095fe7000f4ed7cb3e990b2de79317c573..3c533b04adca4c623ff0c858a59141babe2dccd1 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -83,7 +83,7 @@ public:
~TimeDiscretizedODESystem() override;
- void assemble(const GlobalVector& x_new_timestep,
+ void assemble(std::vector<GlobalVector*> const& x_new_timestep,
int const process_id) override;
void getResidual(GlobalVector const& x_new_timestep,
@@ -181,7 +181,7 @@ public:
~TimeDiscretizedODESystem() override;
- void assemble(const GlobalVector& x_new_timestep,
+ void assemble(std::vector<GlobalVector*> const& x_new_timestep,
int const process_id) override;
void getA(GlobalMatrix& A) const override
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 96e35a35b795db881f8bef16e7d8fb78aca353e9..f0753cf8a5f6cd4a7c950e5a7cd1c0d213881802 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -157,7 +157,7 @@ std::vector<GlobalVector*> setInitialConditions(
auto const nl_tag = process_data->nonlinear_solver_tag;
setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
- nonlinear_solver.assemble(x0, process_id);
+ nonlinear_solver.assemble(process_solutions, process_id);
time_disc.pushState(
t0, x0,
mat_strg); // TODO: that might do duplicate work
@@ -168,7 +168,7 @@ std::vector<GlobalVector*> setInitialConditions(
}
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
- GlobalVector& x, std::size_t const timestep, double const t,
+ std::vector<GlobalVector*>& x, std::size_t const timestep, double const t,
double const delta_t, ProcessData const& process_data,
Output& output_control)
{
@@ -200,7 +200,7 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
if (nonlinear_solver_status.error_norms_met)
{
- process.postNonLinearSolver(x, t, delta_t, process_id);
+ process.postNonLinearSolver(*x[process_id], t, delta_t, process_id);
}
return nonlinear_solver_status;
@@ -535,7 +535,8 @@ void preTimestepForAllProcesses(
static NumLib::NonlinearSolverStatus solveMonolithicProcess(
const double t, const double dt, const std::size_t timestep_id,
- ProcessData const& process_data, GlobalVector& x, Output& output)
+ ProcessData const& process_data, std::vector<GlobalVector*>& x,
+ Output& output)
{
BaseLib::RunTime time_timestep_process;
time_timestep_process.start();
@@ -588,9 +589,8 @@ NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
for (auto& process_data : _per_process_data)
{
auto const process_id = process_data->process_id;
- nonlinear_solver_status =
- solveMonolithicProcess(t, dt, timestep_id, *process_data,
- *_process_solutions[process_id], *_output);
+ nonlinear_solver_status = solveMonolithicProcess(
+ t, dt, timestep_id, *process_data, _process_solutions, *_output);
process_data->nonlinear_solver_status = nonlinear_solver_status;
if (!nonlinear_solver_status.error_norms_met)
@@ -662,8 +662,9 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
process_data->process.setCoupledSolutionsForStaggeredScheme(
&coupled_solutions);
- nonlinear_solver_status = solveOneTimeStepOneProcess(
- x, timestep_id, t, dt, *process_data, *_output);
+ nonlinear_solver_status =
+ solveOneTimeStepOneProcess(_process_solutions, timestep_id, t,
+ dt, *process_data, *_output);
process_data->nonlinear_solver_status = nonlinear_solver_status;
INFO(