tSSAFEM.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tSSAFEM.hh (5754B)
---
1 // Copyright (C) 2009--2017 Jed Brown and Ed Bueler and Constantine Khroulev and David Maxwell
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 _SSAFEM_H_
20 #define _SSAFEM_H_
21
22 #include "SSA.hh"
23 #include "pism/util/FETools.hh"
24 #include "pism/util/petscwrappers/SNES.hh"
25 #include "pism/util/TerminationReason.hh"
26 #include "pism/util/Mask.hh"
27
28 namespace pism {
29
30 namespace stressbalance {
31
32 //! Factory function for constructing a new SSAFEM.
33 SSA * SSAFEMFactory(IceGrid::ConstPtr grid);
34
35 //! PISM's SSA solver: the finite element method implementation written by Jed and David
36 /*!
37 Jed's original code is in rev 831: src/base/ssaJed/...
38 The SSAFEM duplicates most of the functionality of SSAFD, using the finite element method.
39 */
40 class SSAFEM : public SSA {
41 public:
42 SSAFEM(IceGrid::ConstPtr g);
43
44 virtual ~SSAFEM();
45
46 protected:
47 virtual void init_impl();
48 void cache_inputs(const Inputs &inputs);
49
50 //! Storage for SSA coefficients at element nodes.
51 //!
52 //! All fields must be "double" or structures containing "double"
53 //! for IceModelVec2Fat to work correctly.
54 struct Coefficients {
55 //! ice thickness
56 double thickness;
57 //! bed elevation
58 double bed;
59 //! sea level
60 double sea_level;
61 //! basal yield stress
62 double tauc;
63 //! ice hardness
64 double hardness;
65 //! prescribed gravitational driving stress
66 Vector2 driving_stress;
67 };
68
69 const IceModelVec2Int *m_bc_mask;
70 const IceModelVec2V *m_bc_values;
71
72 GeometryCalculator m_gc;
73 double m_alpha;
74 double m_rho_g;
75
76 IceModelVec2Fat<Coefficients> m_coefficients;
77
78 void quad_point_values(const fem::Quadrature &Q,
79 const Coefficients *x,
80 int *mask,
81 double *thickness,
82 double *tauc,
83 double *hardness) const;
84
85 void explicit_driving_stress(const fem::Quadrature &Q,
86 const Coefficients *x,
87 Vector2 *driving_stress) const;
88
89 void driving_stress(const fem::Quadrature &Q,
90 const Coefficients *x,
91 Vector2 *driving_stress) const;
92
93 void PointwiseNuHAndBeta(double thickness,
94 double hardness,
95 int mask,
96 double tauc,
97 const Vector2 &U,
98 const Vector2 &U_x,
99 const Vector2 &U_y,
100 double *nuH, double *dnuH,
101 double *beta, double *dbeta);
102
103 void compute_local_function(Vector2 const *const *const velocity,
104 Vector2 **residual);
105
106 void compute_local_jacobian(Vector2 const *const *const velocity, Mat J);
107
108 virtual void solve(const Inputs &inputs);
109
110 TerminationReason::Ptr solve_with_reason(const Inputs &inputs);
111
112 TerminationReason::Ptr solve_nocache();
113
114 //! Adaptor for gluing SNESDAFormFunction callbacks to an SSAFEM.
115 /* The callbacks from SNES are mediated via SNESDAFormFunction, which has the
116 convention that its context argument is a pointer to a struct
117 having a DA as its first entry. The CallbackData fulfills
118 this requirement, and allows for passing the callback on to an honest
119 SSAFEM object. */
120 struct CallbackData {
121 DM da;
122 SSAFEM *ssa;
123 };
124
125 // objects used internally
126 CallbackData m_callback_data;
127
128 petsc::SNES m_snes;
129
130 //! Storage for node types (interior, boundary, exterior).
131 IceModelVec2Int m_node_type;
132 //! Boundary integral (CFBC contribution to the residual).
133 IceModelVec2V m_boundary_integral;
134
135 double m_dirichletScale;
136 double m_beta_ice_free_bedrock;
137 double m_epsilon_ssa;
138
139 fem::ElementIterator m_element_index;
140 fem::ElementMap m_element;
141 fem::Q1Quadrature4 m_quadrature;
142
143 // Support for direct specification of driving stress to the FEM SSA solver. This helps
144 // with certain test cases where the grid is periodic but the driving stress cannot be the
145 // gradient of a periodic function. (See commit ffb4be16.)
146 const IceModelVec2S *m_driving_stress_x;
147 const IceModelVec2S *m_driving_stress_y;
148 private:
149 void cache_residual_cfbc(const Inputs &inputs);
150 void monitor_jacobian(Mat Jac);
151 void monitor_function(Vector2 const *const *const velocity_global,
152 Vector2 const *const *const residual_global);
153
154 //! SNES callbacks.
155 /*! These simply forward the call on to the SSAFEM member of the CallbackData */
156 static PetscErrorCode function_callback(DMDALocalInfo *info,
157 Vector2 const *const *const velocity,
158 Vector2 **residual,
159 CallbackData *fe);
160 static PetscErrorCode jacobian_callback(DMDALocalInfo *info,
161 Vector2 const *const *const xg,
162 Mat A, Mat J, CallbackData *fe);
163 };
164
165
166 } // end of namespace stressbalance
167 } // end of namespace pism
168
169 #endif /* _SSAFEM_H_ */