JFNKSolver.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
 *      _______              __
 *     / ___/ /  ___  __ _  / /  ___
 *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
 *    \___/_//_/\___/_/_/_/_.__/\___/
 *    Please refer to Copyright.txt, in Chombo's root directory.
 */
#endif

#ifndef _JFNKSOLVER_H_
#define _JFNKSOLVER_H_


#include "IceVelocitySolver.H"
#include "NonlinearViscousTensor.H"
#include "AMRMultiGrid.H"
#include "MultilevelIceVelOp.H"
#include "ViscousTensorOp.H"
#ifdef CH_USE_PETSC
#include "PetscAMRSolver.H"
#endif

#include "NamespaceHeader.H"


enum LinearizationMode {PICARD_LINEARIZATION_MODE, JFNK_LINEARIZATION_MODE, MAX_LINEARIZATION_MODE};


class JFNKSolver;

class LinearizedVTOp :  public LinearOp<Vector<LevelData<FArrayBox>*> >
{  
  friend class JFNKSolver;
#ifdef CH_USE_PETSC
  friend class PetscAMRSolver;
#endif
  
  NonlinearViscousTensor* m_u; // current state
  NonlinearViscousTensor* m_uPerturbed; // perturbed state

  Vector<LevelData<FArrayBox>*> m_fu; // current value of f(u)
  Vector<LevelData<FArrayBox>*> m_uplushv; // current value of u + v

  //constant perturbation coefficient h
  Real m_h;
  //parameters that appear in the adaptive calculation of h
  Real m_err, m_umin; bool m_hAdaptive;
  //compute a perturbation coefficient h given a direction a_v
  Real finiteh(const Vector<LevelData<FArrayBox>*>& a_v);
  

  // m_mlOp provides the multigrid precondioner
  // and is used to evaluate f(u)
  MultilevelIceVelOp m_mlOp;
  
  //m_perturbedMlOp is used to evaluate f(u + h*v)
  MultilevelIceVelOp m_perturbedMlOp;
  
  static int m_bottom_solver_type;

  static int m_MG_solver_depth;

  Vector<DisjointBoxLayout> m_grids;
  Vector<int> m_refRatio;
  Vector<ProblemDomain> m_domains;
  Vector<RealVect> m_dxs;
  int m_lBase;

  LinearizationMode m_mode;

  //options related to writing out the residual
  bool m_writeResiduals;
  static int m_residualID;

  void writeResidual
  (const Vector<LevelData<FArrayBox> *>& a_u,
   const Vector<LevelData<FArrayBox> *>& a_residual);

public:

  NonlinearViscousTensor& current() const {return *m_u;}

  virtual ~LinearizedVTOp()
  {

    if (m_uPerturbed != NULL)
      {
        delete m_uPerturbed;
        m_uPerturbed = NULL;
      }

    for (int lev = 0; lev < m_fu.size(); ++lev)
      {
        if (m_fu[lev] != NULL)
          delete m_fu[lev];
      }
    for (int lev = 0; lev < m_uplushv.size(); ++lev)
      {
        if (m_uplushv[lev] != NULL)
          delete m_uplushv[lev];        
      }
  }
  LinearizedVTOp(NonlinearViscousTensor* a_currentState , 
         Vector<LevelData<FArrayBox>*>& a_u,
         Real a_h,  Real m_err, Real m_umin, bool m_hAdaptive,
         Vector<DisjointBoxLayout>& a_grids,
         Vector<int>& a_refRatio,
         Vector<ProblemDomain>& a_domains,
         Vector<RealVect>& a_dxs,
         int a_lBase, 
         int a_numMGSmooth, 
         int a_numMGIter,
         LinearizationMode a_mode
         );

  // set the current solution u (so that J * v \approx (f(u+hv)-f(u))/h
  virtual void setU(Vector<LevelData<FArrayBox>*>& a_u);
  
  // a_lhs = (f(u+ h * (a_v) ) -f(u))/ h (in JFNK mode)
  // or L[u](a_v) (in Picard mode)
  virtual void applyOp(Vector<LevelData<FArrayBox>*>& a_lhs, 
                       const Vector<LevelData<FArrayBox>*>& a_v, 
                       bool a_homogeneous = false);

  // a_lhs = (f(u+ h * (a_v) )-f(u))/h - a_rhs (in JFNK mode)
  // or L[u](a_v) (in Picard mode)
  virtual void residual(Vector<LevelData<FArrayBox>*>& a_lhs, 
                        const Vector<LevelData<FArrayBox>*>& a_v, 
                        const Vector<LevelData<FArrayBox>*>& a_rhs, 
                        bool a_homogeneous = false);
   
  virtual void outerResidual(Vector<LevelData<FArrayBox>*>& a_lhs, 
                             const Vector<LevelData<FArrayBox>*>& a_u, 
                             const Vector<LevelData<FArrayBox>*>& a_rhs, 
                             bool a_homogeneous = false);

  
  virtual void outerLHS(Vector<LevelData<FArrayBox>*>& a_lhs, 
                             const Vector<LevelData<FArrayBox>*>& a_u,  
                             bool a_homogeneous = false)
  {
    
    m_mlOp.applyOp(a_lhs, a_u, a_homogeneous);
  }
  


  //the rest of LinearOp can be palmed off on m_mlOp.
  virtual void preCond(Vector<LevelData<FArrayBox>*>& a_cor,
                       const Vector<LevelData<FArrayBox>*>& a_residual);

    
  virtual void create(Vector<LevelData<FArrayBox>*>& a_lhs,
                      const Vector<LevelData<FArrayBox>*>& a_rhs)
  { 
    m_mlOp.create(a_lhs , a_rhs);
  }

  virtual void clear(Vector<LevelData<FArrayBox>*>& a_lhs)
  {
    m_mlOp.clear(a_lhs);
  }

  // helpful function which clears memory allocated by create
  virtual void clearStorage(Vector<LevelData<FArrayBox> * >& a_lhs)
  {
    for (int lev=0; lev<a_lhs.size(); lev++)
      {
        if (a_lhs[lev] != NULL)
          {
            delete a_lhs[lev];
            a_lhs[lev] = NULL;
          }
      }
  }

  virtual void assign(Vector<LevelData<FArrayBox>*>& a_lhs,
                      const Vector<LevelData<FArrayBox>*>& a_rhs)
  { 
    m_mlOp.assign(a_lhs, a_rhs);
  }
    
    
  virtual Real dotProduct(const Vector<LevelData<FArrayBox>*>& a_1,
                          const Vector<LevelData<FArrayBox>*>& a_2)
  {
    return m_mlOp.dotProduct(a_1,a_2);
  }

  virtual void incr(Vector<LevelData<FArrayBox>*>& a_lhs,
                    const Vector<LevelData<FArrayBox>*>& a_x,
                    Real a_scale)
  {
    m_mlOp.incr(a_lhs, a_x, a_scale);
  }

    
  virtual void axby(Vector<LevelData<FArrayBox>*>& a_lhs,
                    const Vector<LevelData<FArrayBox>*>& a_x,
                    const Vector<LevelData<FArrayBox>*>& a_y,
                    Real a_a,
                    Real a_b)
  {
    m_mlOp.axby(a_lhs, a_x, a_y, a_a, a_b);
  }

  virtual void scale(Vector<LevelData<FArrayBox>*>& a_lhs,
                     const Real& a_scale)
  {
    m_mlOp.scale(a_lhs , a_scale);
  }

  virtual Real norm(const Vector<LevelData<FArrayBox>*>& a_rhs,
                    int a_ord)
  {
    return m_mlOp.norm(a_rhs, a_ord);
  }

  virtual void setToZero(Vector<LevelData<FArrayBox>*>& a_lhs)
  {
    m_mlOp.setToZero(a_lhs);
  }
  

};



class JFNKSolver : public IceVelocitySolver
{

public:
  class Configuration
  {

  public:

    enum LinearSolverType {Relax,BiCGStab,GMRES,CG,petsc};
    
    ~Configuration();
    
    Configuration();
    
    void parse(const char* a_prefix = "JFNKSolver");

    LinearSolverType m_linearSolverType;
    int m_maxIter;
    Real m_absTol; 
    Real m_relTol;
    Real m_BiCGStabRelTol;
    int m_maxBiCGStabIter;
    Real m_CGRelTol;
    int m_maxCGIter;
    Real m_GMRESRelTol;
    int m_maxGMRESIter;
    Real m_RelaxRelTol;
    int m_maxRelaxIter;
    Real m_RelaxHang;
    int m_normType; 
    int m_verbosity;
    Real m_vtopSafety;
    Real m_vtopRelaxTol;
    int m_vtopRelaxMinIter;
    int m_numMGSmooth;
    int m_numMGIter;
    Real m_h;
    Real m_err;
    Real m_umin; 
    bool m_hAdaptive;
    Real m_switchRate;

    int m_minPicardIter;
    Real m_minStepFactor;
    bool m_writeResiduals;
    bool m_eliminateFastIce;
    Real m_eliminateFastIceSpeed;
    bool m_eliminateFastIceEdgeOnly;
    Real m_eliminateRemoteIceTol;
    Real m_eliminateRemoteIceMaxIter;
    Real m_muMin;
    Real  m_muMax;
    //dont solve, just evaluate the residual
    bool m_residualOnly;
    // scale factor, f: solve for f*u internally - still return u 
    Real m_scale;
    // artificial drag that applies everywhere
    /* if used, need to ensure that 
       (m_artificial_drag * |u|)^m_artifical_drag_power << rho * g * h * grad(s) ~ 10^4
       when |u| ~ 10^4
    */
    Real m_artificialDragCoef, m_artificialDragPower;
    
  };

private:

  Configuration m_config;

public:
  void parseConfiguration(const char* a_prefix)
  {
    m_config.parse(a_prefix);
  }

  const Configuration& configuration()
  {
    return m_config;
  }

  Configuration& get_evil_configuration()
  {
    return m_config;
  }
  
private:

  Vector<ProblemDomain> m_domains;
  
  Vector<DisjointBoxLayout> m_grids;
  
  Vector<RealVect> m_dxs;
  
  Vector<int> m_refRatios;

  ConstitutiveRelation* m_constRelPtr;
  
  BasalFrictionRelation* m_basalFrictionRelPtr;
  
  IceThicknessIBC* m_bcPtr;

  
  
#ifdef CH_USE_PETSC
  // petsc solvers are more expensive to set up, so keep one around
  PetscAMRSolver *m_petscSolver;
#endif

  LinearSolver<Vector<LevelData<FArrayBox>* > >* newLinearSolver(LinearizedVTOp& a_op, 
                                                                 LinearizationMode a_mode);

  Real lineSearch(Vector<LevelData<FArrayBox>* >& a_u, 
                  Vector<LevelData<FArrayBox>* >& a_residual,
                  const Vector<LevelData<FArrayBox>* >& a_rhs,
                  const Vector<LevelData<FArrayBox>* >& a_du,
                  LinearizedVTOp& a_op,
                  NonlinearViscousTensor& a_nvt,
                  Real a_minW, 
                  bool a_resetOnFail);

  void linearSolve(LinearizedVTOp& a_op,
                   Vector<LevelData<FArrayBox>* >& a_u,
                   const Vector<LevelData<FArrayBox>* >& a_rhs, 
                   LinearizationMode a_mode);


  void eliminateFastIce(Vector<LevelData<FArrayBox>* >& a_velocity,
                        Vector<LevelData<FArrayBox>* >& a_calvedIce,
                        Vector<LevelData<FArrayBox>* >& a_addedIce,
                        Vector<LevelData<FArrayBox>* >& a_removedIce,
                        Vector<LevelData<FArrayBox>* >& a_rhs,
                        Vector<RefCountedPtr<LevelSigmaCS > >& a_coordSys,
                        IceNonlinearViscousTensor&  a_current);
  
  
  Real m_convergenceMetric;


public:
  virtual ~JFNKSolver() {
 #ifdef CH_USE_PETSC
    if (m_petscSolver) delete m_petscSolver;
#endif
  }

  JFNKSolver() 
#ifdef CH_USE_PETSC
    : m_petscSolver(NULL)
#endif
  {

  }


  virtual void define(const ProblemDomain& a_coarseDomain,
                      ConstitutiveRelation* a_constRel,
                      BasalFrictionRelation* a_basalFrictionRel,
                      const Vector<DisjointBoxLayout>& a_vectGrids,
                      const Vector<int>& a_vectRefRatio,
                      const RealVect& a_dxCrse,
                      IceThicknessIBC* a_bc,
                      int a_numLevels);

  
  //IceVelocitySolver non-linear solve
  virtual int solve(Vector<LevelData<FArrayBox>* >& a_horizontalVel,
                    Vector<LevelData<FArrayBox>* >& a_calvedIce,
                    Vector<LevelData<FArrayBox>* >& a_addedIce,
                    Vector<LevelData<FArrayBox>* >& a_removedIce,
                    Real& a_initialResidualNorm, Real& a_finalResidualNorm,
                    const Real a_convergenceMetric,
                    const Vector<LevelData<FArrayBox>* >& a_rhs,
                    const Vector<LevelData<FArrayBox>* >& a_C,
                    const Vector<LevelData<FArrayBox>* >& a_C0,
                    const Vector<LevelData<FArrayBox>* >& a_A,
                    const Vector<LevelData<FArrayBox>* >& a_muCoef,
                    Vector<RefCountedPtr<LevelSigmaCS > >& a_coordSys,
                    Real a_time,
                    int a_lbase, int a_maxLevel) ; 


  virtual int solve(Vector<LevelData<FArrayBox>* >& a_horizontalVel,
                    Vector<LevelData<FArrayBox>* >& a_calvedIce,
                    Vector<LevelData<FArrayBox>* >& a_addedIce,
                    Vector<LevelData<FArrayBox>* >& a_removedIce,
                    Real& a_initialResidualNorm, Real& a_finalResidualNorm,
                    const Real a_convergenceMetric,
                    const bool a_linear,
                    const Vector<LevelData<FArrayBox>* >& a_rhs,
                    const Vector<LevelData<FArrayBox>* >& a_C,
                    const Vector<LevelData<FArrayBox>* >& a_C0,
                    const Vector<LevelData<FArrayBox>* >& a_A,
                    const Vector<LevelData<FArrayBox>* >& a_muCoef,
                    Vector<RefCountedPtr<LevelSigmaCS > >& a_coordSys,
                    Real a_time,
                    int a_lbase, int a_maxLevel);

  
private:


  
  void imposeMaxAbs(Vector<LevelData<FArrayBox>*>& a_u,
                  Real a_limit);
  
};







#include "NamespaceFooter.H"
#endif