From 1f81164117ecad41ecefc081579aa0e406f6f974 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Thu, 8 Feb 2018 12:17:29 +0100
Subject: [PATCH] [NL] ODE; Add PETSc SNES nonlinear solver.
---
NumLib/ODESolver/PETScNonlinearSolver.cpp | 183 ++++++++++++++++++++++
NumLib/ODESolver/PETScNonlinearSolver.h | 75 +++++++++
2 files changed, 258 insertions(+)
create mode 100644 NumLib/ODESolver/PETScNonlinearSolver.cpp
create mode 100644 NumLib/ODESolver/PETScNonlinearSolver.h
diff --git a/NumLib/ODESolver/PETScNonlinearSolver.cpp b/NumLib/ODESolver/PETScNonlinearSolver.cpp
new file mode 100644
index 00000000000..6ea95db32bc
--- /dev/null
+++ b/NumLib/ODESolver/PETScNonlinearSolver.cpp
@@ -0,0 +1,183 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifdef USE_PETSC
+
+#include "PETScNonlinearSolver.h"
+
+#include <petscmat.h>
+#include <petscvec.h>
+
+#include "BaseLib/RunTime.h"
+
+namespace
+{
+struct PetscContext
+{
+ using System = NumLib::NonlinearSystem<NumLib::NonlinearSolverTag::Newton>;
+ System* system;
+ std::vector<GlobalVector*>& x;
+ std::vector<GlobalVector*> const& x_prev;
+ GlobalVector* r;
+ GlobalMatrix* J;
+ int const process_id;
+};
+
+PetscErrorCode updateResidual(SNES /*snes*/, Vec x, Vec petsc_r,
+ void* petsc_context)
+{
+ auto context = static_cast<PetscContext*>(petsc_context);
+
+ DBUG("PETScNonlinearSolver: residual callback called.");
+
+ VecCopy(x, context->x[context->process_id]->getRawVector());
+
+ // Assemble in ogs context.
+ BaseLib::RunTime time_assembly;
+ time_assembly.start();
+ context->system->assemble(context->x, context->x_prev, context->process_id);
+
+ INFO("[time] Assembly took {} s.", time_assembly.elapsed());
+ context->system->getResidual(*context->x[context->process_id],
+ *context->x_prev[context->process_id],
+ *context->r);
+ context->J->finalizeAssembly();
+
+ context->system->getJacobian(*context->J);
+ context->system->applyKnownSolutionsNewton(
+ *context->J, *context->r, *context->x[context->process_id]);
+
+ VecCopy(context->r->getRawVector(), petsc_r);
+
+ return 0;
+}
+
+PetscErrorCode updateJacobian(SNES /*snes*/, Vec /*x*/, Mat J,
+ Mat /*same as J*/, void* petsc_context)
+{
+ DBUG("PETScNonlinearSolver: jacobian callback called.");
+ // Assume the system is already assembled.
+ // Copy the results into petsc vectors.
+
+ auto context = static_cast<PetscContext*>(petsc_context);
+ MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
+ MatCopy(context->J->getRawMatrix(), J, DIFFERENT_NONZERO_PATTERN);
+ return 0;
+}
+} // namespace
+
+namespace NumLib
+{
+PETScNonlinearSolver::PETScNonlinearSolver(
+ GlobalLinearSolver& /*linear_solver*/)
+{
+ SNESCreate(PETSC_COMM_WORLD, &_snes_solver);
+ SNESSetFromOptions(_snes_solver);
+#ifndef NDEBUG
+ PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_WORLD);
+#endif // NDEBUG
+}
+
+void PETScNonlinearSolver::setEquationSystem(System& eq,
+ ConvergenceCriterion& conv_crit)
+{
+ _equation_system = &eq;
+ _convergence_criterion = &conv_crit;
+}
+
+void PETScNonlinearSolver::compensateNonEquilibriumInitialResiduum(
+ bool const value)
+{
+ _compensate_non_equilibrium_initial_residuum = value;
+}
+
+void PETScNonlinearSolver::calculateNonEquilibriumInitialResiduum(
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<GlobalVector*> const& /*x_prev*/, int const /*process_id*/)
+{
+ if (!_compensate_non_equilibrium_initial_residuum)
+ {
+ return;
+ }
+
+ OGS_FATAL(
+ "Non-equilibrium initial residuum is not implemented for the "
+ "PETScNonlinearSolver.");
+}
+
+NonlinearSolverStatus PETScNonlinearSolver::solve(
+ std::vector<GlobalVector*>& x,
+ std::vector<GlobalVector*> const& x_prev,
+ std::function<void(
+ int,
+ std::vector<GlobalVector*> const&)> const& /*postIterationCallback*/,
+ int const process_id)
+{
+ DBUG("PETScNonlinearSolver: solve()");
+ using TimeDiscretizedSystem =
+ TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Newton>;
+
+ auto* system = static_cast<TimeDiscretizedSystem*>(_equation_system);
+
+ DBUG("PETScNonlinearSolver: create vectors");
+ // r and J on which the ogs assembly operates.
+ auto& J = NumLib::GlobalMatrixProvider::provider.getMatrix(
+ system->getMatrixSpecifications(process_id), _jacobian_id);
+ auto& r = NumLib::GlobalVectorProvider::provider.getVector(
+ system->getMatrixSpecifications(process_id), _residual_id);
+
+ // Vectors and matrices used by SNES for solutions. These will be locked
+ // during the SNES' solve call.
+ auto& J_snes = NumLib::GlobalMatrixProvider::provider.getMatrix(
+ system->getMatrixSpecifications(process_id), _petsc_jacobian_id);
+ auto& r_snes = NumLib::GlobalVectorProvider::provider.getVector(
+ system->getMatrixSpecifications(process_id), _petsc_residual_id);
+ auto& x_snes = NumLib::GlobalVectorProvider::provider.getVector(
+ system->getMatrixSpecifications(process_id), _petsc_x_id);
+ MathLib::LinAlg::copy(*x[process_id], x_snes); // Initial guess.
+
+ PetscContext petsc_context{_equation_system, x, x_prev, &r, &J, process_id};
+
+ DBUG("PETScNonlinearSolver: set function");
+ SNESSetFunction(_snes_solver, r_snes.getRawVector(), updateResidual,
+ &petsc_context);
+
+ DBUG("PETScNonlinearSolver: set jacobian");
+ // The jacobian and the preconditioner matrices point to the same location.
+ SNESSetJacobian(_snes_solver, J_snes.getRawMatrix(), J_snes.getRawMatrix(),
+ updateJacobian, &petsc_context);
+
+ DBUG("PETScNonlinearSolver: call SNESSolve");
+ SNESSolve(_snes_solver, nullptr, x_snes.getRawVector());
+
+ SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
+ SNESGetConvergedReason(_snes_solver, &reason);
+ INFO("PETSsSNES convergence reason {}.", reason);
+
+ PetscInt iterations;
+ SNESGetIterationNumber(_snes_solver, &iterations);
+ INFO("PETScSNES used {} iterations.", iterations);
+
+ // Copy back the solution.
+ MathLib::LinAlg::copy(x_snes, *x[process_id]);
+
+ NumLib::GlobalVectorProvider::provider.releaseVector(x_snes);
+ NumLib::GlobalVectorProvider::provider.releaseVector(r_snes);
+ NumLib::GlobalMatrixProvider::provider.releaseMatrix(J_snes);
+ NumLib::GlobalVectorProvider::provider.releaseVector(r);
+ NumLib::GlobalMatrixProvider::provider.releaseMatrix(J);
+
+ return {reason >= 0, iterations};
+}
+
+} // namespace NumLib
+#endif // USE_PETSC
diff --git a/NumLib/ODESolver/PETScNonlinearSolver.h b/NumLib/ODESolver/PETScNonlinearSolver.h
new file mode 100644
index 00000000000..6074a553071
--- /dev/null
+++ b/NumLib/ODESolver/PETScNonlinearSolver.h
@@ -0,0 +1,75 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+#ifdef USE_PETSC
+
+#include <petscsnes.h>
+
+#include "ConvergenceCriterion.h"
+#include "NonlinearSolver.h"
+#include "TimeDiscretizedODESystem.h"
+
+namespace NumLib
+{
+//! \addtogroup ODESolver
+//! @{
+class PETScNonlinearSolver final : public NonlinearSolverBase
+{
+public:
+ //! Type of the nonlinear equation system to be solved.
+ using System = NonlinearSystem<NonlinearSolverTag::Newton>;
+
+ /*! Constructs a new instance.
+ *
+ * \param linear_solver the linear solver used by this nonlinear solver.
+ */
+ explicit PETScNonlinearSolver(GlobalLinearSolver& linear_solver);
+
+ //! Set the nonlinear equation system that will be solved.
+ void setEquationSystem(System& eq, ConvergenceCriterion& conv_crit);
+
+ void compensateNonEquilibriumInitialResiduum(bool const value);
+
+ void calculateNonEquilibriumInitialResiduum(
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev,
+ int const process_id) override;
+
+ NonlinearSolverStatus solve(
+ std::vector<GlobalVector*>& x,
+ std::vector<GlobalVector*> const& x_prev,
+ std::function<void(int, std::vector<GlobalVector*> const&)> const&
+ postIterationCallback,
+ int const process_id) override;
+
+private:
+ SNES _snes_solver;
+
+ System* _equation_system = nullptr;
+ ConvergenceCriterion* _convergence_criterion = nullptr;
+
+ std::size_t _residual_id = 0u; //!< ID of the residual vector.
+ std::size_t _jacobian_id = 0u; //!< ID of the Jacobian matrix.
+
+ std::size_t _petsc_x_id = 0u;
+ std::size_t _petsc_jacobian_id = 0u;
+ std::size_t _petsc_residual_id = 0u;
+
+ // clang-format off
+ /// \copydoc NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_non_equilibrium_initial_residuum
+ bool _compensate_non_equilibrium_initial_residuum = false;
+ // clang-format on
+};
+
+//! @}
+
+} // namespace NumLib
+#endif // USE_PETSC
--
GitLab