SSA.hh (6181B)
1 // Copyright (C) 2004--2019 Jed Brown, Ed Bueler and Constantine Khroulev 2 // 3 // This file is part of PISM. 4 // 5 // PISM is free software; you can redistribute it and/or modify it under the 6 // terms of the GNU General Public License as published by the Free Software 7 // Foundation; either version 3 of the License, or (at your option) any later 8 // version. 9 // 10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY 11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 13 // details. 14 // 15 // You should have received a copy of the GNU General Public License 16 // along with PISM; if not, write to the Free Software 17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 18 19 #ifndef _SSA_H_ 20 #define _SSA_H_ 21 22 #include "pism/stressbalance/ShallowStressBalance.hh" 23 #include "pism/util/IceModelVec2CellType.hh" 24 25 namespace pism { 26 27 class Geometry; 28 29 namespace stressbalance { 30 31 //! Gives an extension coefficient to maintain ellipticity of SSA where ice is thin. 32 /*! 33 The SSA is basically a nonlinear elliptic, but vector-valued, equation which 34 determines the ice velocity field from the driving stress, the basal shear 35 stress, the ice hardness, and some boundary conditions. The problem loses 36 ellipticity (coercivity) if the thickness actually goes to zero. This class 37 provides an extension coefficient to maintain ellipticity. 38 39 More specifically, the SSA equations are 40 \f[ 41 \def\ddx#1{\ensuremath{\frac{\partial #1}{\partial x}}} 42 \def\ddy#1{\ensuremath{\frac{\partial #1}{\partial y}}} 43 - 2 \ddx{}\left[\nu H \left(2 \ddx{u} + \ddy{v}\right)\right] 44 - \ddy{}\left[\nu H \left(\ddy{u} + \ddx{v}\right)\right] 45 + \tau_{(b)x} = - \rho g H \ddx{h}, 46 \f] 47 and another similar equation for the \f$y\f$-component. Schoof 48 \ref SchoofStream shows that these PDEs are the variational equations for a 49 coercive functional, thus (morally) elliptic. 50 51 The quantity \f$\nu H\f$ is the nonlinear coefficient, and conceptually it is a 52 membrane strength. This class extends \f$\nu H\f$ to have a minimum value 53 at all points. It is a class, and not just a configuration constant, because 54 setting both the thickness \f$H\f$ and the value \f$\nu H\f$ are allowed, and 55 setting each of these does not affect the value of the other. 56 */ 57 class SSAStrengthExtension { 58 public: 59 SSAStrengthExtension(const Config &c); 60 61 void set_notional_strength(double my_nuH); 62 void set_min_thickness(double my_min_thickness); 63 double get_notional_strength() const; 64 double get_min_thickness() const; 65 private: 66 double m_min_thickness, m_constant_nu; 67 }; 68 69 //! Callback for constructing a new SSA subclass. The caller is 70 //! responsible for deleting the newly constructed SSA. 71 /*! The factory idiom gives a way to implement runtime polymorphism for the 72 choice of SSA algorithm. The factory is a function pointer that takes 73 all the arguments of an SSA constructor and returns a newly constructed instance. 74 Subclasses of SSA should provide an associated function pointer matching the 75 SSAFactory typedef */ 76 class SSA; 77 typedef SSA * (*SSAFactory)(IceGrid::ConstPtr); 78 79 80 //! PISM's SSA solver. 81 /*! 82 An object of this type solves equations for the vertically-constant horizontal 83 velocity of ice that is sliding over land or is floating. The equations are, in 84 their clearest divergence form 85 \f[ - \frac{\partial T_{ij}}{\partial x_j} - \tau_{(b)i} = f_i \f] 86 where \f$i,j\f$ range over \f$x,y\f$, \f$T_{ij}\f$ is a depth-integrated viscous 87 stress tensor (%i.e. equation (2.6) in [\ref SchoofStream]). 88 These equations determine velocity in a more-or-less elliptic manner. 89 Here \f$\tau_{(b)i}\f$ are the components of the basal shear stress applied to 90 the base of the ice. The right-hand side \f$f_i\f$ is the driving shear stress, 91 \f[ f_i = - \rho g H \frac{\partial h}{\partial x_i}. \f] 92 Here \f$H\f$ is the ice thickness and \f$h\f$ is the elevation of the surface of 93 the ice. More concretely, the SSA equations are 94 \f{align*} 95 - 2 \left[\nu H \left(2 u_x + v_y\right)\right]_x 96 - \left[\nu H \left(u_y + v_x\right)\right]_y 97 - \tau_{(b)1} &= - \rho g H h_x, \\ 98 - \left[\nu H \left(u_y + v_x\right)\right]_x 99 - 2 \left[\nu H \left(u_x + 2 v_y\right)\right]_y 100 - \tau_{(b)2} &= - \rho g H h_y, 101 \f} 102 where \f$u\f$ is the \f$x\f$-component of the velocity and \f$v\f$ is the 103 \f$y\f$-component of the velocity [\ref MacAyeal, \ref Morland, \ref WeisGreveHutter]. 104 105 Derived classes actually implement numerical methods to solve these equations. 106 This class is virtual, but it actually implements some helper functions believed 107 to be common to all implementations (%i.e. regular grid implementations) and it 108 provides the basic fields. 109 */ 110 class SSA : public ShallowStressBalance { 111 public: 112 SSA(IceGrid::ConstPtr g); 113 virtual ~SSA(); 114 115 SSAStrengthExtension *strength_extension; 116 117 virtual void update(const Inputs &inputs, bool full_update); 118 119 void set_initial_guess(const IceModelVec2V &guess); 120 121 virtual std::string stdout_report() const; 122 123 const IceModelVec2V& driving_stress() const; 124 protected: 125 virtual void define_model_state_impl(const File &output) const; 126 virtual void write_model_state_impl(const File &output) const; 127 128 virtual void init_impl(); 129 130 virtual DiagnosticList diagnostics_impl() const; 131 132 virtual void compute_driving_stress(const IceModelVec2S &ice_thickness, 133 const IceModelVec2S &surface_elevation, 134 const IceModelVec2CellType &cell_type, 135 const IceModelVec2Int *no_model_mask, 136 IceModelVec2V &result) const; 137 138 virtual void solve(const Inputs &inputs) = 0; 139 140 IceModelVec2CellType m_mask; 141 IceModelVec2V m_taud; 142 143 std::string m_stdout_ssa; 144 145 // objects used by the SSA solver (internally) 146 petsc::DM::Ptr m_da; // dof=2 DA 147 IceModelVec2V m_velocity_global; // global vector for solution 148 149 // profiling 150 int m_event_ssa; 151 }; 152 153 } // end of namespace stressbalance 154 } // end of namespace pism 155 156 #endif /* _SSA_H_ */