PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
Blatter.hh
Go to the documentation of this file.
1 /* Copyright (C) 2020, 2021, 2022, 2023 PISM Authors
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 PISM_BLATTER_H
20 #define PISM_BLATTER_H
21 
22 #include "pism/stressbalance/ShallowStressBalance.hh"
23 #include "pism/util/petscwrappers/SNES.hh"
24 #include "pism/util/petscwrappers/DM.hh"
25 #include "pism/util/petscwrappers/Vec.hh"
26 #include "pism/util/fem/FEM.hh"
27 
28 namespace pism {
29 
30 namespace fem {
31 class Element3;
32 class Q1Element3Face;
33 } // end of namespace fem
34 
35 namespace stressbalance {
36 
37 class Blatter : public ShallowStressBalance {
38 public:
39  Blatter(std::shared_ptr<const Grid> grid, int Mz, int coarsening_factor);
40  virtual ~Blatter() = default;
41 
42  void update(const Inputs &inputs, bool);
43 
44  std::shared_ptr<array::Array3D> velocity_u_sigma() const;
45  std::shared_ptr<array::Array3D> velocity_v_sigma() const;
46 
47  /*!
48  * 2D input parameters
49  */
50  struct Parameters {
51  // elevation (z coordinate) of the bottom domain boundary
52  double bed;
53  // thickness of the domain
54  double thickness;
55  // NodeType stored as double
56  double node_type;
57  // basal yield stress
58  double tauc;
59  // sea level elevation (used to determine if a location is grounded)
60  double sea_level;
61  // floatation function (positive where floating, zero or negative where grounded)
62  double floatation;
63  // FIXME: Add Dirichlet BC at a map plane location.
64  };
65 
66 protected:
67  // u and v components of ice velocity on the fine sigma grid
68  std::shared_ptr<array::Array3D> m_u_sigma, m_v_sigma;
69 
70  // 3D dof=2 DM used by m_snes
72  // storage for the solution
74  // storage for the old solution during parameter continuation
76  // solver
78 
80 
81  // Scaling of quadrature weights (note: this does not seem to matter).
82  double m_scaling;
83 
84  // Ice density times g
85  double m_rho_ice_g;
86 
87  // Water density times g
88  double m_rho_ocean_g;
89 
90  // The flow law Glen exponent
92 
93  // Use the eta-transformation to compute surface gradient
95 
96  // Viscosity regularization constant
98 
99  // Enhancement factor for ice viscosity.
101 
102  // True if the Eisenstat-Walker method of adjusting linear solver tolerances is enabled.
104 
105  static const int m_Nq = 100;
106  static const int m_n_work = 9;
107 
108  double m_work[m_n_work][m_Nq];
109 
111 
114 
115  void init_impl();
116 
117  void define_model_state_impl(const File &output) const;
118 
119  void write_model_state_impl(const File &output) const;
120 
121  static bool exterior_element(const int *node_type);
122 
123  static bool grounding_line(const double *F);
124 
125  static bool partially_submerged_face(int face, const double *z, const double *sea_level);
126 
127  void compute_node_type(double min_thickness);
128 
129  virtual void nodal_parameter_values(const fem::Q1Element3 &element,
130  Parameters **P,
131  int i,
132  int j,
133  int *node_type,
134  double *bottom,
135  double *thickness,
136  double *surface,
137  double *sea_level) const;
138 
139  virtual bool marine_boundary(int face,
140  const int *node_type,
141  const double *ice_bottom,
142  const double *sea_level);
143 
144  virtual bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex& I);
145 
146  virtual Vector2d u_bc(double x, double y, double z) const;
147 
148  void compute_jacobian(DMDALocalInfo *info, const Vector2d ***x, Mat A, Mat J);
149 
150  void jacobian_dirichlet(const DMDALocalInfo &info, Parameters **P, Mat J);
151 
152  virtual void jacobian_f(const fem::Q1Element3 &element,
153  const Vector2d *u_nodal,
154  const double *B_nodal,
155  double K[2 * fem::q13d::n_chi][2 * fem::q13d::n_chi]);
156 
157  virtual void jacobian_basal(const fem::Q1Element3Face &face,
158  const double *tauc_nodal,
159  const double *f_nodal,
160  const Vector2d *u_nodal,
161  double K[2 * fem::q13d::n_chi][2 * fem::q13d::n_chi]);
162 
163  void compute_residual(DMDALocalInfo *info, const Vector2d ***X, Vector2d ***R);
164 
165  void residual_dirichlet(const DMDALocalInfo &info,
166  Parameters **P,
167  const Vector2d ***x,
168  Vector2d ***R);
169 
170  virtual void residual_f(const fem::Q1Element3 &element,
171  const Vector2d *u_nodal,
172  const double *B_nodal,
173  Vector2d *residual);
174 
175  virtual void residual_source_term(const fem::Q1Element3 &element,
176  const double *surface,
177  const double *bed,
178  Vector2d *residual);
179 
180  virtual void residual_basal(const fem::Q1Element3 &element,
181  const fem::Q1Element3Face &face,
182  const double *tauc_nodal,
183  const double *f_nodal,
184  const Vector2d *u_nodal,
185  Vector2d *residual);
186 
187  virtual void residual_surface(const fem::Q1Element3 &element,
188  const fem::Q1Element3Face &face,
189  Vector2d *residual);
190 
191  virtual void residual_lateral(const fem::Q1Element3 &element,
192  const fem::Q1Element3Face &face,
193  const double *surface_nodal,
194  const double *z_nodal,
195  const double *sl_nodal,
196  Vector2d *residual);
197 
198  static PetscErrorCode jacobian_callback(DMDALocalInfo *info,
199  const Vector2d ***x,
200  Mat A, Mat J,
201  Blatter *solver);
202 
203  static PetscErrorCode function_callback(DMDALocalInfo *info,
204  const Vector2d ***x, Vector2d ***f,
205  Blatter *solver);
206 
207  virtual void init_2d_parameters(const Inputs &inputs);
208 
209  void init_ice_hardness(const Inputs &inputs, const petsc::DM &da);
210 
211  // Guts of the constructor. This method wraps PETSc calls to simplify error checking.
212  PetscErrorCode setup(DM pism_da, grid::Periodicity p, int Mz, int coarsening_factor,
213  const std::string &prefix);
214 
215  void set_initial_guess(const array::Array3D &u_sigma, const array::Array3D &v_sigma);
216 
217  void copy_solution();
218 
220 
221  void get_basal_velocity(array::Vector &result);
222 
223  void report_mesh_info();
224 
225  struct SolutionInfo {
227  SNESConvergedReason snes_reason;
228  KSPConvergedReason ksp_reason;
229  };
230 
233 };
234 
235 } // end of namespace stressbalance
236 } // end of namespace pism
237 
238 #endif /* PISM_BLATTER_H */
std::shared_ptr< const Grid > grid() const
Definition: Component.cc:105
High-level PISM I/O class.
Definition: File.hh:56
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2d.hh:29
A storage vector combining related fields in a struct.
Definition: Array2D.hh:32
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: Array3D.hh:33
3D Q1 element
Definition: Element.hh:351
void jacobian_dirichlet(const DMDALocalInfo &info, Parameters **P, Mat J)
Definition: jacobian.cc:168
void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Blatter.cc:683
void get_basal_velocity(array::Vector &result)
Definition: Blatter.cc:1118
virtual void residual_f(const fem::Q1Element3 &element, const Vector2d *u_nodal, const double *B_nodal, Vector2d *residual)
Definition: residual.cc:36
virtual Vector2d u_bc(double x, double y, double z) const
Definition: Blatter.cc:132
static const int m_n_work
Definition: Blatter.hh:106
virtual void residual_source_term(const fem::Q1Element3 &element, const double *surface, const double *bed, Vector2d *residual)
Definition: residual.cc:94
void compute_node_type(double min_thickness)
Definition: Blatter.cc:56
virtual void residual_basal(const fem::Q1Element3 &element, const fem::Q1Element3Face &face, const double *tauc_nodal, const double *f_nodal, const Vector2d *u_nodal, Vector2d *residual)
Definition: residual.cc:163
void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Blatter.cc:688
void compute_residual(DMDALocalInfo *info, const Vector2d ***X, Vector2d ***R)
Definition: residual.cc:311
virtual bool marine_boundary(int face, const int *node_type, const double *ice_bottom, const double *sea_level)
Definition: Blatter.cc:225
fem::Q1Element3Face m_face4
Definition: Blatter.hh:112
static bool partially_submerged_face(int face, const double *z, const double *sea_level)
Definition: Blatter.cc:192
virtual void residual_lateral(const fem::Q1Element3 &element, const fem::Q1Element3Face &face, const double *surface_nodal, const double *z_nodal, const double *sl_nodal, Vector2d *residual)
Definition: residual.cc:217
virtual void init_2d_parameters(const Inputs &inputs)
Definition: Blatter.cc:524
std::shared_ptr< array::Array3D > m_u_sigma
Definition: Blatter.hh:68
std::shared_ptr< array::Array3D > velocity_v_sigma() const
Definition: Blatter.cc:1197
SolutionInfo parameter_continuation()
Definition: Blatter.cc:807
void update(const Inputs &inputs, bool)
Definition: Blatter.cc:974
PetscErrorCode setup(DM pism_da, grid::Periodicity p, int Mz, int coarsening_factor, const std::string &prefix)
Definition: Blatter.cc:382
double m_work[m_n_work][m_Nq]
Definition: Blatter.hh:108
void init_ice_hardness(const Inputs &inputs, const petsc::DM &da)
Definition: Blatter.cc:571
static PetscErrorCode jacobian_callback(DMDALocalInfo *info, const Vector2d ***x, Mat A, Mat J, Blatter *solver)
Definition: jacobian.cc:354
void residual_dirichlet(const DMDALocalInfo &info, Parameters **P, const Vector2d ***x, Vector2d ***R)
Definition: residual.cc:262
static PetscErrorCode function_callback(DMDALocalInfo *info, const Vector2d ***x, Vector2d ***f, Blatter *solver)
Definition: residual.cc:457
static const int m_Nq
Definition: Blatter.hh:105
fem::Q1Element3Face m_face100
Definition: Blatter.hh:113
void compute_averaged_velocity(array::Vector &result)
Definition: Blatter.cc:1158
std::shared_ptr< array::Array3D > velocity_u_sigma() const
Definition: Blatter.cc:1193
Blatter(std::shared_ptr< const Grid > grid, int Mz, int coarsening_factor)
Definition: Blatter.cc:258
static bool exterior_element(const int *node_type)
Definition: Blatter.cc:146
virtual void jacobian_f(const fem::Q1Element3 &element, const Vector2d *u_nodal, const double *B_nodal, double K[2 *fem::q13d::n_chi][2 *fem::q13d::n_chi])
Definition: jacobian.cc:35
static bool grounding_line(const double *F)
Definition: Blatter.cc:166
virtual void jacobian_basal(const fem::Q1Element3Face &face, const double *tauc_nodal, const double *f_nodal, const Vector2d *u_nodal, double K[2 *fem::q13d::n_chi][2 *fem::q13d::n_chi])
Definition: jacobian.cc:122
void compute_jacobian(DMDALocalInfo *info, const Vector2d ***x, Mat A, Mat J)
Definition: jacobian.cc:199
array::Array2D< Parameters > m_parameters
Definition: Blatter.hh:79
virtual bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex &I)
Definition: Blatter.cc:124
virtual void residual_surface(const fem::Q1Element3 &element, const fem::Q1Element3Face &face, Vector2d *residual)
Definition: residual.cc:200
std::shared_ptr< array::Array3D > m_v_sigma
Definition: Blatter.hh:68
virtual ~Blatter()=default
Vector2d m_work2[m_n_work][m_Nq]
Definition: Blatter.hh:110
virtual void nodal_parameter_values(const fem::Q1Element3 &element, Parameters **P, int i, int j, int *node_type, double *bottom, double *thickness, double *surface, double *sea_level) const
Definition: Blatter.cc:627
void set_initial_guess(const array::Array3D &u_sigma, const array::Array3D &v_sigma)
Definition: Blatter.cc:1134
Shallow stress balance (such as the SSA).
const int n_chi
Number of shape functions on a Q1 element.
Definition: FEM.hh:213
Periodicity
Definition: Grid.hh:51
static double K(double psi_x, double psi_y, double speed, double epsilon)
static double F(double SL, double B, double H, double alpha)
const int J[]
Definition: ssafd_code.cc:34
const int I[]
Definition: ssafd_code.cc:24