diff --git a/src/model/common/non_linear_solver/non_linear_solver_petsc.cc b/src/model/common/non_linear_solver/non_linear_solver_petsc.cc
index 29a63530d..7a956190c 100644
--- a/src/model/common/non_linear_solver/non_linear_solver_petsc.cc
+++ b/src/model/common/non_linear_solver/non_linear_solver_petsc.cc
@@ -1,301 +1,213 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "non_linear_solver_petsc.hh"
#include "dof_manager_petsc.hh"
#include "mpi_communicator_data.hh"
#include "solver_callback.hh"
#include "solver_vector_petsc.hh"
#include "sparse_matrix_petsc.hh"
/* -------------------------------------------------------------------------- */
#include "petscsnes.h"
/* -------------------------------------------------------------------------- */
namespace akantu {
NonLinearSolverPETSc::NonLinearSolverPETSc(
DOFManagerPETSc & dof_manager, const ModelSolverOptions & solver_options,
const ID & id)
: NonLinearSolver(dof_manager, solver_options, id) {
if (solver_options.sparse_solver_type != SparseSolverType::_petsc)
AKANTU_EXCEPTION(
"petsc non linear solver works only with petsc sparse solver");
- // std::unordered_map<NonLinearSolverType, SNESType>
- // petsc_non_linear_solver_types{
- // {NonLinearSolverType::_newton_raphson, SNESNEWTONLS},
- // {NonLinearSolverType::_linear, SNESKSPONLY},
- // {NonLinearSolverType::_gmres, SNESNGMRES},
- // {NonLinearSolverType::_bfgs, SNESQN},
- // {NonLinearSolverType::_cg, SNESNCG}};
-
this->has_internal_set_param = true;
- // for (const auto & pair : petsc_non_linear_solver_types) {
- // supported_type.insert(pair.first);
- // }
supported_type.insert(NonLinearSolverType::_petsc_snes);
this->checkIfTypeIsSupported();
auto && mpi_comm = dof_manager.getMPIComm();
SNESCreate(mpi_comm, &snes);
- // auto it = petsc_non_linear_solver_types.find(non_linear_solver_type);
- // if (it != petsc_non_linear_solver_types.end()) {
- // SNESSetType(snes, it->second);
- // }
SNESSetType(snes, SNESNEWTONLS);
-
SNESSetFromOptions(snes);
this->registerParam("max_iterations", max_iterations, 10, _pat_parsmod,
"Max number of iterations");
this->registerParam("threshold", convergence_criteria, 1e-10, _pat_parsmod,
"Threshold to consider results as converged");
this->registerParam("convergence_type", convergence_criteria_type,
SolveConvergenceCriteria::_solution, _pat_parsmod,
"Type of convergence criteria");
}
/* -------------------------------------------------------------------------- */
NonLinearSolverPETSc::~NonLinearSolverPETSc() { SNESDestroy(&snes); }
/* -------------------------------------------------------------------------- */
void NonLinearSolverPETSc::saveSolution() { AKANTU_TO_IMPLEMENT(); }
/* -------------------------------------------------------------------------- */
void NonLinearSolverPETSc::restoreSolution() { AKANTU_TO_IMPLEMENT(); }
/* -------------------------------------------------------------------------- */
-void NonLinearSolverPETSc::corrector(Vec & x) {
- PetscInt iteration;
- SNESGetIterationNumber(snes, &iteration);
-
- // if (prev_iteration == iteration) {
- // return;
- // }
-
- prev_iteration = iteration;
-
- // SolverVectorPETSc dx(*this->x, "dx");
- // VecWAXPY(dx, -1., *solution_prev, x); // w = alpha x + y.
-
- // std::cout << "x= ";
- // PetscPrint(x);
- // std::cout << std::endl;
- // std::cout << "solution_prev= " << *solution_prev << std::endl;
- // std::cout << "dx= " << dx << std::endl;
+void NonLinearSolverPETSc::corrector(Vec x) {
auto & solution = aka::as_type<SolverVectorPETSc>(dof_manager.getSolution());
-
- // VecCopy(solution, *solution_prev);
- // VecCopy(dx, solution);
if (x != solution.getVec())
VecCopy(x, solution);
dof_manager.splitSolutionPerDOFs();
- auto & model_x = this->dof_manager.getDOFs("displacement");
-
callback->restoreLastConvergedStep();
callback->corrector();
- // VecCopy(*solution_prev, solution);
}
+/* -------------------------------------------------------------------------- */
+
void NonLinearSolverPETSc::assembleResidual(Vec x, Vec f) {
- // saveSolution();
corrector(x);
auto & residual =
dynamic_cast<SolverVectorPETSc &>(dof_manager.getResidual());
- // if (residual.getVec() != f) {
- // VecCopy(residual, *residual_prev);
- // }
callback->assembleResidual();
const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFsIndexes();
-
std::vector<Real> zeros_to_set(blocked_dofs.size());
-
VecSetValuesLocal(residual, blocked_dofs.size(), blocked_dofs.data(),
zeros_to_set.data(), INSERT_VALUES);
// for PETSc F is -akantu::residual
VecScale(residual, -1);
if (residual.getVec() != f) {
VecCopy(residual, f);
- // VecCopy(*residual_prev, residual);
}
- std::cout << "x= ";
- PetscPrint(x);
- std::cout << " f = ";
- PetscPrint(f);
- std::cout << std::endl;
-
- // restoreSolution();
}
-void NonLinearSolverPETSc::assembleJacobian(Vec x) {
- // corrector(x);
+/* -------------------------------------------------------------------------- */
+
+void NonLinearSolverPETSc::assembleJacobian(Vec x, Mat J) {
+ corrector(x);
callback->assembleMatrix("J");
- PetscPrint(aka::as_type<SparseMatrixPETSc>(this->dof_manager.getMatrix("J")));
+ auto & _J = aka::as_type<SparseMatrixPETSc>(dof_manager.getMatrix("J"));
+ if (_J.getMat() != J) {
+ MatCopy(_J, J, SAME_NONZERO_PATTERN);
+ }
}
-// void NonLinearSolverPETSc::reset() { prev_iteration = -1; }
-
-// void NonLinearSolverPETSc::setInitialSolution(SolverVectorPETSc & x) {
-// VecCopy(x, x_prev);
-// }
-
-// void NonLinearSolverPETSc::setCallback(SolverCallback & callback) {
-// this->callback = &callback;
-// }
-
/* -------------------------------------------------------------------------- */
+
PetscErrorCode NonLinearSolverPETSc::FormFunction(SNES /*snes*/, Vec x, Vec f,
void * ctx) {
- auto * _this = reinterpret_cast<NonLinearSolverPETSc *>(ctx);
- _this->assembleResidual(x, f);
-
+ reinterpret_cast<NonLinearSolverPETSc *>(ctx)->assembleResidual(x, f);
return 0;
}
/* -------------------------------------------------------------------------- */
PetscErrorCode NonLinearSolverPETSc::FormJacobian(SNES /*snes*/, Vec x, Mat J,
- Mat P, void * ctx) {
- std::cout << "J= ";
- PetscPrint(J);
- std::cout << std::endl;
- std::cout << "P= ";
- PetscPrint(P);
- std::cout << std::endl;
- auto * _this = reinterpret_cast<NonLinearSolverPETSc *>(ctx);
- _this->assembleJacobian(x);
- std::cout << "Jafter= ";
- PetscPrint(J);
- std::cout << std::endl;
- std::cout << "Pafter= ";
- PetscPrint(P);
- std::cout << std::endl;
+ Mat /*P*/, void * ctx) {
+ reinterpret_cast<NonLinearSolverPETSc *>(ctx)->assembleJacobian(x, J);
return 0;
}
/* -------------------------------------------------------------------------- */
void NonLinearSolverPETSc::solve(SolverCallback & callback) {
callback.beforeSolveStep();
this->dof_manager.updateGlobalBlockedDofs();
callback.assembleMatrix("J");
- auto & global_x =
- dynamic_cast<SolverVectorPETSc &>(dof_manager.getSolution());
- global_x.zero();
-
- if (not x) {
- x = std::make_unique<SolverVectorPETSc>(global_x, "temporary_solution");
- solution_prev =
- std::make_unique<SolverVectorPETSc>(global_x, "previous_solution");
- residual_prev =
- std::make_unique<SolverVectorPETSc>(global_x, "previous_residual");
- }
+ auto & x = dynamic_cast<SolverVectorPETSc &>(dof_manager.getSolution());
+ x.zero();
- *x = global_x;
-
- // this->reset();
- prev_iteration = -1;
-
- // this->setCallback(callback);
this->callback = &callback;
- // this->setInitialSolution(global_x);
- VecCopy(*x, *solution_prev);
-
auto & rhs = aka::as_type<SolverVectorPETSc>(dof_manager.getResidual());
+ rhs.zero();
+
auto & J = aka::as_type<SparseMatrixPETSc>(dof_manager.getMatrix("J"));
SNESSetFunction(snes, rhs, NonLinearSolverPETSc::FormFunction, this);
SNESSetJacobian(snes, J, J, NonLinearSolverPETSc::FormJacobian, this);
- rhs.zero();
-
callback.predictor();
- // callback.assembleResidual();
- SNESView(snes, PETSC_VIEWER_STDOUT_WORLD);
- SNESSolve(snes, nullptr, *x);
+ // SNESView(snes, PETSC_VIEWER_STDOUT_WORLD);
+ SNESSolve(snes, nullptr, x);
SNESGetConvergedReason(snes, &reason);
SNESGetIterationNumber(snes, &n_iter);
- auto & model_x = this->dof_manager.getDOFs("displacement");
+ // access the model solution counter part: only for debug
+ // auto & model_x = this->dof_manager.getDOFs("displacement");
dof_manager.splitSolutionPerDOFs();
callback.restoreLastConvergedStep();
callback.corrector();
bool converged = reason >= 0;
callback.afterSolveStep(converged);
if (not converged) {
PetscReal atol;
PetscReal rtol;
PetscReal stol;
PetscInt maxit;
PetscInt maxf;
SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
AKANTU_CUSTOM_EXCEPTION(debug::SNESNotConvergedException(
this->reason, this->n_iter, stol, atol, rtol, maxit));
}
}
/* -------------------------------------------------------------------------- */
void NonLinearSolverPETSc::updateInternalParameters() {
std::map<ID, ID> akantu_to_petsc_option = {{"max_iterations", "snes_max_it"},
{"threshold", "snes_stol"}};
for (auto && [param, param_akantu] : akantu_to_petsc_option) {
auto & value = this->get(param);
PetscOptionsSetValue(nullptr, ("-" + param_akantu).c_str(),
value.to_string().c_str());
}
SNESSetFromOptions(snes);
PetscOptionsClear(nullptr);
}
/* -------------------------------------------------------------------------- */
void NonLinearSolverPETSc::parseSection(const ParserSection & section) {
auto parameters = section.getParameters();
for (auto && param : range(parameters.first, parameters.second)) {
PetscOptionsSetValue(nullptr, param.getName().c_str(),
param.getValue().c_str());
}
SNESSetFromOptions(snes);
PetscOptionsClear(nullptr);
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/model/common/non_linear_solver/non_linear_solver_petsc.hh b/src/model/common/non_linear_solver/non_linear_solver_petsc.hh
index ec9cddc45..0cf10f655 100644
--- a/src/model/common/non_linear_solver/non_linear_solver_petsc.hh
+++ b/src/model/common/non_linear_solver/non_linear_solver_petsc.hh
@@ -1,112 +1,108 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "non_linear_solver.hh"
#include "solver_vector_petsc.hh"
/* -------------------------------------------------------------------------- */
#include <petscsnes.h>
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_NON_LINEAR_SOLVER_PETSC_HH_
#define AKANTU_NON_LINEAR_SOLVER_PETSC_HH_
namespace akantu {
class DOFManagerPETSc;
class SolverVectorPETSc;
} // namespace akantu
namespace akantu {
class NonLinearSolverPETSc : public NonLinearSolver {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
NonLinearSolverPETSc(DOFManagerPETSc & dof_manager,
const ModelSolverOptions & solver_options,
const ID & id = "non_linear_solver_petsc");
~NonLinearSolverPETSc() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// solve the system described by the jacobian matrix, and rhs contained in
/// the dof manager
void solve(SolverCallback & callback) override;
/// parse the arguments from the input file
void parseSection(const ParserSection & section) override;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
static PetscErrorCode FormFunction(SNES snes, Vec dx, Vec f, void * ctx);
static PetscErrorCode FormJacobian(SNES snes, Vec dx, Mat J, Mat P,
void * ctx);
- void corrector(Vec & x);
+ void corrector(Vec x);
void assembleResidual(Vec x, Vec f);
- void assembleJacobian(Vec x);
+ void assembleJacobian(Vec x, Mat J);
void updateInternalParameters() override;
void saveSolution();
void restoreSolution();
/// PETSc non linear solver
SNES snes;
SNESConvergedReason reason;
SolverCallback * callback{nullptr};
std::unique_ptr<SolverVectorPETSc> x;
- std::unique_ptr<SolverVectorPETSc> solution_prev;
- std::unique_ptr<SolverVectorPETSc> residual_prev;
Int n_iter{0};
Int max_iterations;
/// Type of convergence criteria
SolveConvergenceCriteria convergence_criteria_type;
/// convergence threshold
Real convergence_criteria;
- // save previous iteration
- PetscInt prev_iteration{-1};
};
namespace debug {
class SNESNotConvergedException : public NLSNotConvergedException {
public:
SNESNotConvergedException(SNESConvergedReason reason, Int niter, Real error,
Real absolute_tolerance, Real relative_tolerance,
Int max_iterations)
: NLSNotConvergedException(relative_tolerance, niter, error),
reason(reason), absolute_tolerance(absolute_tolerance),
max_iterations(max_iterations) {}
SNESConvergedReason reason;
Real absolute_tolerance;
Int max_iterations;
};
} // namespace debug
} // namespace akantu
#endif /* AKANTU_NON_LINEAR_SOLVER_PETSC_HH_ */