[NL] SNES Add inequality constraints callback.

This allows a process to set lower or upper bounds
to the solution vector for SNESVI solvers.

Only two types of SNES solver support this.

NumLib/ODESolver/NonlinearSystem.h

@@ -71,6 +71,10 @@ public:
     virtual void applyKnownSolutionsNewton(
         GlobalMatrix& Jac, GlobalVector& res,
         GlobalVector& minus_delta_x) const = 0;
+
+    virtual void updateConstraints(GlobalVector& /*lower*/,
+                                   GlobalVector& /*upper*/,
+                                   int /*process_id*/) = 0;
 };
 
 /*! A System of nonlinear equations to be solved with the Picard fixpoint

NumLib/ODESolver/ODESystem.h

@@ -66,6 +66,10 @@ public:
     {
         return nullptr;  // by default there are no known solutions
     }
+
+    virtual void updateConstraints(GlobalVector& /*lower*/,
+                                   GlobalVector& /*upper*/,
+                                   int const /*process_id*/){};
 };
 
 /*! Interface for a first-order implicit quasi-linear ODE.

NumLib/ODESolver/PETScNonlinearSolver.cpp

@@ -156,6 +156,29 @@ NonlinearSolverStatus PETScNonlinearSolver::solve(
     SNESSetJacobian(_snes_solver, J_snes.getRawMatrix(), J_snes.getRawMatrix(),
                     updateJacobian, &petsc_context);
 
+    std::unique_ptr<GlobalVector> xl = nullptr;
+    std::unique_ptr<GlobalVector> xu = nullptr;
+
+    SNESType snes_type;
+    SNESGetType(_snes_solver, &snes_type);
+    if ((std::strcmp(snes_type, SNESVINEWTONRSLS) == 0) ||
+        (std::strcmp(snes_type, SNESVINEWTONSSLS) == 0))
+    {
+        // Set optional constraints via callback.
+        DBUG("PETScNonlinearSolver: set constraints");
+        xl = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+            system->getMatrixSpecifications(process_id));
+        xu = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+            system->getMatrixSpecifications(process_id));
+
+        system->updateConstraints(*xl, *xu, process_id);
+        MathLib::finalizeVectorAssembly(*xl);
+        MathLib::finalizeVectorAssembly(*xu);
+
+        SNESVISetVariableBounds(_snes_solver, xl->getRawVector(),
+                                xu->getRawVector());
+    }
+
     DBUG("PETScNonlinearSolver: call SNESSolve");
     SNESSolve(_snes_solver, nullptr, x_snes.getRawVector());
 

NumLib/ODESolver/TimeDiscretizedODESystem.h

@@ -100,6 +100,12 @@ public:
     void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
                                    GlobalVector& minus_delta_x) const override;
 
+    void updateConstraints(GlobalVector& lower, GlobalVector& upper,
+                           int const process_id) override
+    {
+        _ode.updateConstraints(lower, upper, process_id);
+    }
+
     bool isLinear() const override { return _ode.isLinear(); }
 
     void preIteration(const unsigned iter, GlobalVector const& x) override