nox_solver.h

// This file is part of HermesCommon
//
// Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
// Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
//
// Hermes2D is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation; either version 2 of the License,
// or (at your option) any later version.
//
// Hermes2D 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Hermes2D; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#ifndef __H2D_NEWTON_SOLVER_NOX_H_
#define __H2D_NEWTON_SOLVER_NOX_H_

#include "solvers/linear_matrix_solver.h"
#include "solvers/interfaces/epetra.h"
#if(defined HAVE_NOX && defined HAVE_EPETRA && defined HAVE_TEUCHOS)
#include <NOX.H>
#ifdef _POSIX_C_SOURCE
# undef _POSIX_C_SOURCE  // pyconfig.h included by NOX_Epetra defines it
#endif
#ifdef _XOPEN_SOURCE
# undef _XOPEN_SOURCE  // pyconfig.h included by NOX_Epetra defines it
#endif
#undef interface
#include <NOX_Epetra.H>
#include "exceptions.h"
#include "discrete_problem.h"

namespace Hermes
{
  namespace Hermes2D
  {
    template <typename Scalar> class NewtonSolverNOX;

    template <typename Scalar>
    class HERMES_API DiscreteProblemNOX :
      public Hermes::Hermes2D::DiscreteProblem<Scalar>,
      public NOX::Epetra::Interface::Required,
      public NOX::Epetra::Interface::Jacobian,
      public NOX::Epetra::Interface::Preconditioner
    {
    public:
      DiscreteProblemNOX(WeakFormSharedPtr<Scalar> wf, std::vector<SpaceSharedPtr<Scalar> > spaces);
      DiscreteProblemNOX(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space);
      DiscreteProblemNOX();
      ~DiscreteProblemNOX();

      void set_precond(Teuchos::RCP<Hermes::Preconditioners::EpetraPrecond<Scalar> > &pc);

      Teuchos::RCP<Hermes::Preconditioners::EpetraPrecond<Scalar> > get_precond();

      // that is generated by the Interface class.
      EpetraMatrix<Scalar> *get_jacobian();

      virtual bool computeF(const Epetra_Vector &x, Epetra_Vector &f, FillType flag = Residual);

      virtual bool computeJacobian(const Epetra_Vector &x, Epetra_Operator &op);

      virtual bool computePreconditioner(const Epetra_Vector &x, Epetra_Operator &m,
        Teuchos::ParameterList *precParams = 0);

    private:
      EpetraMatrix<Scalar>* jacobian;
      Teuchos::RCP<Hermes::Preconditioners::EpetraPrecond<Scalar> > precond;

      friend class NewtonSolverNOX < Scalar > ;
    };

    template <typename Scalar>
    class HERMES_API NewtonSolverNOX : public Hermes::Mixins::Loggable
    {
    private:
      Teuchos::RCP<Teuchos::ParameterList> nl_pars;
    public:
      NewtonSolverNOX(DiscreteProblemNOX<Scalar> *problem);
      NewtonSolverNOX(WeakFormSharedPtr<Scalar> wf, std::vector<SpaceSharedPtr<Scalar> > spaces);
      NewtonSolverNOX(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space);

      void init();

      virtual ~NewtonSolverNOX();

      virtual void set_time(double time);
      virtual void set_time_step(double time_step);

      virtual void solve(Scalar* coeff_vec = nullptr);
      virtual void solve(std::vector<MeshFunctionSharedPtr<Scalar> > initial_guess);
      virtual void solve(MeshFunctionSharedPtr<Scalar> initial_guess);

      Scalar* get_sln_vector();

      virtual int get_num_iters();
      virtual double get_residual();
      int get_num_lin_iters();
      double get_achieved_tol();

      void set_output_flags(int flags);

      void set_ls_type(const char *type);
      void set_ls_max_iters(int iters);
      void set_ls_tolerance(double tolerance);
      void set_ls_sizeof_krylov_subspace(int size);

      void set_norm_type(NOX::Abstract::Vector::NormType type);
      void set_scale_type(NOX::StatusTest::NormF::ScaleType type);
      void set_conv_iters(int iters);
      void set_conv_abs_resid(double resid);
      void set_conv_rel_resid(double resid);
      void disable_abs_resid();
      void disable_rel_resid();
      void set_conv_update(double update);
      void set_conv_wrms(double rtol, double atol);

      void set_precond_reuse(const char * pc_reuse);
      void set_precond_max_age(int max_age);
      virtual void set_precond(Hermes::Preconditioners::EpetraPrecond<Scalar> &pc);
      virtual void set_precond(const char *pc);

      Scalar* sln_vector;
      DiscreteProblemNOX<Scalar> *dp;

      bool own_dp;

    protected:
      int num_iters;
      double residual;
      int num_lin_iters;
      double achieved_tol;

      // convergence params
      struct conv_t
      {
        int max_iters;
        double abs_resid;
        double rel_resid;
        NOX::Abstract::Vector::NormType norm_type;
        NOX::StatusTest::NormF::ScaleType stype;
        double update;
        double wrms_rtol;
        double wrms_atol;
      } conv;

      struct conv_flag_t
      {
        unsigned absresid : 1;
        unsigned relresid : 1;
        unsigned wrms : 1;
        unsigned update : 1;
      } conv_flag;
    };
  }
}
#endif
#endif