tSSA.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tSSA.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_ */