InverseVerticallyIntegratedVelocitySolver.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 _INVERSE_VI_VEL_SOLVER_H_
#define _INVERSE_VI_VEL_SOLVER_H_

#include "JFNKSolver.H"
#include "LevelDataOps.H"
#include "SurfaceFlux.H"
#include "BasalFriction.H"
#include "MuCoefficient.H"
#include "IceConstants.H"
#include "LevelDataBasalFriction.H"
#include "CGOptimize.H"
#include "NamespaceHeader.H"



class InverseIceVelocitySolver :  public IceVelocitySolver
{

  static void copyAMR ( Vector<RefCountedPtr<LevelData<FArrayBox> > >&a_copy, 
                        const Vector<LevelData<FArrayBox>*>& a_orig )
  {
    a_copy.resize(0);

    for (int lev = 0; lev < a_orig.size(); lev++)
      {
        if (a_orig[lev])
          {
            const DisjointBoxLayout& dbl = a_orig[lev]->disjointBoxLayout();
            LevelData<FArrayBox>* lptr = new LevelData<FArrayBox>
              (dbl, a_orig[lev]->nComp(),a_orig[lev]->ghostVect());
            for (DataIterator dit(dbl); dit.ok(); ++dit)
              {
                (*lptr)[dit].copy( (*a_orig[lev])[dit]);
              }
            a_copy.push_back(RefCountedPtr<LevelData<FArrayBox> >( lptr));
          }     
      }
  }


public:
  
  virtual void setPreviousTime(Real a_time) =0;
  
  static BasalFriction* basalFriction
  (const Vector<LevelData<FArrayBox>*>& a_C, const Vector<int>& a_ratio, const RealVect& a_dxCrse )
  {
    Vector<RefCountedPtr<LevelData<FArrayBox> > > C;
    copyAMR(C, a_C);
    MultiLevelDataBasalFriction* ptr = new  MultiLevelDataBasalFriction(C,a_dxCrse,a_ratio);
    return ptr;
  }

  static MuCoefficient* muCoefficient
  (const Vector<LevelData<FArrayBox>*>& a_mu, const Vector<int>& a_ratio, const RealVect& a_dxCrse )
  {
    Vector<RefCountedPtr<LevelData<FArrayBox> > > mu;
    copyAMR(mu, a_mu);
    MultiLevelDataMuCoefficient* ptr = new  MultiLevelDataMuCoefficient(mu,a_dxCrse,a_ratio);
    return ptr;
  }

  virtual BasalFriction* basalFriction() = 0;
  
  virtual MuCoefficient* muCoefficient() = 0;


  
  
};



class InverseVerticallyIntegratedVelocitySolver : public InverseIceVelocitySolver, 
                                                  public ObjectiveWithGradient<Vector<LevelData<FArrayBox>*> >
{

public:

  class Configuration
  {
  public: 

    Configuration();

    ~Configuration();


    void parse(const char* a_prefix = "control");


    int m_minLevelForOptimization;

    Real m_minTimeBetweenOptimizations;

    Real m_dtTypical;
    
    enum BoundMethod {none, projection, MAX_BOUND_METHOD};
    BoundMethod m_boundMethod;
    
    enum VelocityMisfitType {speed, velocity, log_speed, MAX_VELOCITY_MISFIT_TYPE};
    VelocityMisfitType m_velMisfitType;

    bool m_optimizeX0;
    bool m_optimizeX1;
    bool m_gradMuCoefShelfOnly;
    
    Real m_lowerX0;

    Real m_upperX0; 

    Real m_lowerX1;

    Real m_upperX1; 

    Real m_initialLowerMuCoef;

    Real m_initialUpperMuCoef;  

    Real m_initialLowerC;

    Real m_initialUpperC;

    Real m_velMisfitCoefficient;

    Real m_divuhMisfitCoefficient;

    Real m_divuhMisfitSmooth;

    Real m_gradCsqRegularization; 

    Real m_gradMuCoefsqRegularization; 

    

    
    Real m_gradX0sqRegularization; 

    Real m_gradX1sqRegularization;

    Real m_X0Regularization;

    Real m_X1Regularization;

    Real m_X0TimeRegularization;

    Real m_X1TimeRegularization;

    Real m_thicknessThreshold;

    Real m_CGmaxIter;

    Real m_CGtol;

    Real m_CGsecantParameter;

    Real m_CGsecantStepMaxGrow;

    int m_CGsecantMaxIter;

    Real m_CGsecantTol;

    Real m_CGhang;

    Real m_CGrestartInterval;

    SurfaceData* m_velObs_c;
    
    SurfaceData* m_velObs_x;
    
#if CH_SPACEDIM > 1
    SurfaceData* m_velObs_y;
#endif

    bool m_writeInnerSteps;

    std::string m_innerStepFileNameBase;

    std::string m_outerStepFileNameBase;

    SurfaceData* m_divuhObs_c;

    SurfaceData* m_divuhObs_a;    

    SurfaceData* m_gradientFactor;

    int m_outerCounter;

    

  };

private:

  Configuration m_config;

  Vector<LevelDataOps<FArrayBox> > m_vectOps;

  template <class T>
  void create(Vector<LevelData<T>*>& a_data, 
              const int a_ncomp, const IntVect& a_ghost)
  {
    if (a_data.size() <= m_finest_level)
      a_data.resize(m_finest_level+1);

    for (int lev =0; lev <= m_finest_level; lev++)
      {
        if (a_data[lev] != NULL)
          delete a_data[lev];
        a_data[lev] = new LevelData<T>(m_grids[lev], a_ncomp, a_ghost);
      }
  }

  template <class T>
  void free(Vector<LevelData<T>*>& a_data)
  {
    for (int lev =0; lev < a_data.size(); lev++)
      {
        if (a_data[lev] != NULL)
          {
            delete a_data[lev];
            a_data[lev] = NULL;
          }
      }
  }

  // set a_x = 0.0
  void setToZero(Vector<LevelData<FArrayBox>* >& a_a);


  void  mapX(const Vector<LevelData<FArrayBox>* >& a_x);

  void writeState
  (const std::string& a_file, int a_counter,
   const Vector<LevelData<FArrayBox>* >& a_x,
   const Vector<LevelData<FArrayBox>* >& a_g) const;


  void computeDivUH();

  //calculate the rhs to the adjoint equation, and the misfits
  void computeAdjointRhs();

  // add directional derivatives of unregularized problem to a_g
  void computeGradient(Vector<LevelData<FArrayBox>* >& a_g, 
                       const  Vector<LevelData<FArrayBox>* >& a_x);

  void regularizeGradient(Vector<LevelData<FArrayBox>* >& a_g, 
                          const  Vector<LevelData<FArrayBox>* >& a_x);

  void applyProjection(Vector<LevelData<FArrayBox>* >& a_g, 
                       const  Vector<LevelData<FArrayBox>* >& a_x);
    
  void solveStressEqn (Vector<LevelData<FArrayBox>* >& a_u,
                       const bool a_adjoint,
                       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);
  Real m_time; 

  Real m_prev_time;
  
  Real XnTimeCoef()
  {
    return  (m_time > m_prev_time + TINY_NORM)?( 1.0 / (m_time-m_prev_time)):0.0;
  }

  Real X0Regularization()
  {
    return m_config.m_X0Regularization +  m_config.m_X0TimeRegularization * XnTimeCoef();
  }

  Real X1Regularization()
  {
    return m_config.m_X1Regularization +  m_config.m_X1TimeRegularization * XnTimeCoef();
  }

  int m_innerCounter, m_outerCounter;
  
  bool m_velocityInitialised;

  int m_finest_level;

  Vector<DisjointBoxLayout> m_grids;
  
  Vector<int> m_refRatio;
  
  Vector<RealVect> m_dx;

  Vector<ProblemDomain> m_domain;

  std::string outerStateFile() const;

  std::string innerStateFile() const;

  Vector<RefCountedPtr<LevelSigmaCS > > m_coordSys;
  
  Vector<LevelData<FArrayBox>*> m_A;
  
  Vector<LevelData<FArrayBox>*> m_C0;

  // coefficients of basal friction.
  // m_C = exp(X[0]) * m_COrigin;
  Vector<LevelData<FArrayBox>*> m_COrigin; // original value of C
  Vector<LevelData<FArrayBox>*> m_C      ; // current value of C
  Vector<LevelData<FArrayBox>*> m_Cmasked ; // C with shelf values set 
  Vector<LevelData<FArrayBox>*> m_lapC;
  Vector<LevelData<FArrayBox>*> m_gradCSq;

  Vector<LevelData<FArrayBox>*> m_muCoefOrigin; // original value of muCoef
  Vector<LevelData<FArrayBox>*> m_muCoef      ; // current value of muCoef
  Vector<LevelData<FArrayBox>*> m_lapMuCoef;
  Vector<LevelData<FArrayBox>*> m_gradMuCoefSq;
  //Vector<LevelData<FluxBox>*> m_faceMuCoef;
  
  Vector<LevelData<FArrayBox>*> m_lapX;
  Vector<LevelData<FArrayBox>*> m_gradXSq;
  
  Vector<LevelData<FArrayBox>*> m_rhs;

  Vector<LevelData<FArrayBox>*> m_adjRhs;

  Vector<LevelData<FArrayBox>*> m_adjVel;
  
  Vector<LevelData<FArrayBox>*> m_velObs;

  Vector<LevelData<FArrayBox>*> m_velCoef;

  Vector<LevelData<FArrayBox>*> m_divuhObs;

  Vector<LevelData<FArrayBox>*> m_divuhCoef;

  Vector<LevelData<FArrayBox>*> m_velb;

  Vector<LevelData<FArrayBox>*> m_vels;

  Vector<LevelData<FArrayBox>*> m_divuh;
  
  Vector<LevelData<FArrayBox>*> m_velocityMisfit;
  
  Vector<LevelData<FArrayBox>*> m_divuhMisfit;

  //physics pointers
  BasalFrictionRelation* m_basalFrictionRelation;
  ConstitutiveRelation* m_constitutiveRelation;
  IceThicknessIBC* m_thicknessIBC;

  Real m_bestMisfit;
  Vector<LevelData<FArrayBox>*> m_bestVel;
  Vector<LevelData<FArrayBox>*> m_bestC;
  Vector<LevelData<FArrayBox>*> m_bestMuCoef;
  bool m_optimization_done; // don't want to provide the dat aabove unless it is meaningful
  
  const AmrIceBase* m_amrIce;

  // need to pass these to functions like JFNKSolver::solve
  Vector<LevelData<FArrayBox>* > m_calvedIce, m_addedIce, m_removedIce;
  
public:

  InverseVerticallyIntegratedVelocitySolver();

  ~InverseVerticallyIntegratedVelocitySolver();

  void setPreviousTime(Real a_time)
  {
    m_prev_time = a_time;
  }
  
  void define(const ProblemDomain& a_coarseDomain,
              ConstitutiveRelation* a_constRel,
              BasalFrictionRelation* a_FrictionRel,
              const Vector<DisjointBoxLayout>& a_grids,
              const Vector<int>& a_refRatio,
              const RealVect& a_dxCrse,
              IceThicknessIBC* a_bc,
              int a_numLevels)
  {
    MayDay::Error("not sure why we got here...");
  }

  void define(const AmrIceBase& a_amrIce,
              const ProblemDomain& a_coarseDomain,
              ConstitutiveRelation* a_constRel,
              BasalFrictionRelation* a_FrictionRel,
              const Vector<DisjointBoxLayout>& a_grids,
              const Vector<int>& a_refRatio,
              const RealVect& a_dxCrse,
              IceThicknessIBC* a_bc,
              int a_numLevels);
  
 
  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);




  //implementation of  ObjectiveWithGradient
  void computeObjectiveAndGradient
  (Real& a_fm, Real& a_fp,
   Vector<LevelData<FArrayBox>* >& a_g, 
   const  Vector<LevelData<FArrayBox>* >& a_x, 
   bool a_inner);
  
  void restart(){;}

  //allocate storage for a_ca, copy from a_b
  void create(Vector<LevelData<FArrayBox>* >& a_a,  const  Vector<LevelData<FArrayBox>* >& a_b);
  
  //apply preconditioner s = M^{-1}r
  void preCond(Vector<LevelData<FArrayBox>* >& a_s,  const Vector<LevelData<FArrayBox>* >& a_r);


  void free(Vector<LevelData<FArrayBox>*>& a_data)
  {
    free<FArrayBox>(a_data);
  }


  // add a scalar to a_a
  void plus(Vector<LevelData<FArrayBox>* >& a_a, Real a_r );

  // set a_x = s * a_x
  void scale(Vector<LevelData<FArrayBox>* >& a_x, 
             const  Real a_s);

  // set a_y = a_x
  void assign(Vector<LevelData<FArrayBox>* >& a_y, 
              const Vector<LevelData<FArrayBox>* >& a_x);
  
  // set a_y = a_y + a_s * a_x
  void incr(Vector<LevelData<FArrayBox>* >& y, 
            const Vector<LevelData<FArrayBox>* >& x, Real s);
  
  // return a_y.a_x
  Real dotProduct(Vector<LevelData<FArrayBox>* >& a_y, 
                  const Vector<LevelData<FArrayBox>* >& a_x);

  // no of degrees of freedom
  int nDoF(const Vector<LevelData<FArrayBox>* >& x);

  BasalFriction* basalFriction()
  {
    BasalFriction* ptr = NULL;
    if (m_optimization_done)
      ptr = InverseIceVelocitySolver::basalFriction(m_bestC, m_refRatio, m_dx[0]);
    return ptr;
  }
  
  MuCoefficient* muCoefficient()
  {
    MuCoefficient* ptr = NULL;
    if (m_optimization_done)
      ptr =  InverseIceVelocitySolver::muCoefficient(m_bestMuCoef, m_refRatio, m_dx[0]);
    return ptr;
  }
};



#include "NamespaceFooter.H"
#endif