tSSAFEM.cc - 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
---
tSSAFEM.cc (42154B)
---
1 // Copyright (C) 2009--2019 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 #include "pism/util/IceGrid.hh"
20 #include "SSAFEM.hh"
21 #include "pism/util/FETools.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/basalstrength/basal_resistance.hh"
24 #include "pism/rheology/FlowLaw.hh"
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/Vars.hh"
28 #include "pism/stressbalance/StressBalance.hh"
29 #include "pism/geometry/Geometry.hh"
30
31 #include "pism/util/node_types.hh"
32
33 namespace pism {
34 namespace stressbalance {
35
36 /** The Q1 finite element SSA solver.
37 *
38 *
39 *
40 */
41 SSAFEM::SSAFEM(IceGrid::ConstPtr g)
42 : SSA(g),
43 m_bc_mask(NULL),
44 m_bc_values(NULL),
45 m_gc(*m_config),
46 m_element_index(*g),
47 m_element(*g),
48 m_quadrature(g->dx(), g->dy(), 1.0) {
49
50 const double ice_density = m_config->get_number("constants.ice.density");
51 m_alpha = 1 - ice_density / m_config->get_number("constants.sea_water.density");
52 m_rho_g = ice_density * m_config->get_number("constants.standard_gravity");
53
54 m_driving_stress_x = NULL;
55 m_driving_stress_y = NULL;
56
57 PetscErrorCode ierr;
58
59 m_dirichletScale = 1.0;
60 m_beta_ice_free_bedrock = m_config->get_number("basal_resistance.beta_ice_free_bedrock");
61
62 ierr = SNESCreate(m_grid->com, m_snes.rawptr());
63 PISM_CHK(ierr, "SNESCreate");
64
65 // Set the SNES callbacks to call into our compute_local_function and compute_local_jacobian.
66 m_callback_data.da = *m_da;
67 m_callback_data.ssa = this;
68
69 ierr = DMDASNESSetFunctionLocal(*m_da, INSERT_VALUES,
70 (DMDASNESFunction)function_callback,
71 &m_callback_data);
72 PISM_CHK(ierr, "DMDASNESSetFunctionLocal");
73
74 ierr = DMDASNESSetJacobianLocal(*m_da,
75 (DMDASNESJacobian)jacobian_callback,
76 &m_callback_data);
77 PISM_CHK(ierr, "DMDASNESSetJacobianLocal");
78
79 ierr = DMSetMatType(*m_da, "baij");
80 PISM_CHK(ierr, "DMSetMatType");
81
82 ierr = DMSetApplicationContext(*m_da, &m_callback_data);
83 PISM_CHK(ierr, "DMSetApplicationContext");
84
85 ierr = SNESSetDM(m_snes, *m_da);
86 PISM_CHK(ierr, "SNESSetDM");
87
88 // Default of maximum 200 iterations; possibly overridden by command line options
89 int snes_max_it = 200;
90 ierr = SNESSetTolerances(m_snes, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT,
91 snes_max_it, PETSC_DEFAULT);
92 PISM_CHK(ierr, "SNESSetTolerances");
93
94 ierr = SNESSetFromOptions(m_snes);
95 PISM_CHK(ierr, "SNESSetFromOptions");
96
97 // Allocate m_coefficients, which contains coefficient data at the nodes of all the elements.
98 m_coefficients.create(m_grid, "ssa_coefficients", WITH_GHOSTS, 1);
99
100 m_node_type.create(m_grid, "node_type", WITH_GHOSTS, 1);
101 m_node_type.set_attrs("internal", // intent
102 "node types: interior, boundary, exterior", // long name
103 "", "", "", 0); // no units or standard name
104
105 // ElementMap::nodal_values() expects a ghosted IceModelVec2S. Ghosts if this field are never
106 // assigned to and not communocated, though.
107 m_boundary_integral.create(m_grid, "boundary_integral", WITH_GHOSTS, 1);
108 m_boundary_integral.set_attrs("internal", // intent
109 "residual contribution from lateral boundaries", // long name
110 "", "", "", 0); // no units or standard name
111 }
112
113 SSA* SSAFEMFactory(IceGrid::ConstPtr g) {
114 return new SSAFEM(g);
115 }
116
117 SSAFEM::~SSAFEM() {
118 // empty
119 }
120
121 // Initialize the solver, called once by the client before use.
122 void SSAFEM::init_impl() {
123
124 SSA::init_impl();
125
126 // Use explicit driving stress if provided.
127 if (m_grid->variables().is_available("ssa_driving_stress_x") and
128 m_grid->variables().is_available("ssa_driving_stress_y")) {
129 m_driving_stress_x = m_grid->variables().get_2d_scalar("ssa_driving_stress_x");
130 m_driving_stress_y = m_grid->variables().get_2d_scalar("ssa_driving_stress_y");
131 }
132
133 m_log->message(2,
134 " [using the SNES-based finite element method implementation]\n");
135
136 // process command-line options
137 {
138 m_dirichletScale = 1.0e9;
139 m_dirichletScale = options::Real("-ssa_fe_dirichlet_scale",
140 "Enforce Dirichlet conditions with this additional scaling",
141 m_dirichletScale);
142
143 }
144
145 // On restart, SSA::init() reads the SSA velocity from a PISM output file
146 // into IceModelVec2V "velocity". We use that field as an initial guess.
147 // If we are not restarting from a PISM file, "velocity" is identically zero,
148 // and the call below clears m_velocity_global.
149
150 m_velocity_global.copy_from(m_velocity);
151 }
152
153 /** Solve the SSA system of equations.
154
155 The FEM solver computes values of the various coefficients (most notably: the vertically-averaged
156 ice hardness) at each of the element nodes *only once*.
157
158 When running in an ice-model context, at each time step, SSA::update() is called, which calls
159 SSAFEM::solve(). Since coefficients have generally changed between time steps, we need to recompute
160 coefficients. On the other hand, in the context of inversion, coefficients will not change between
161 iteration and there is no need to recompute their values.
162
163 So there are two different solve methods, SSAFEM::solve() and SSAFEM::solve_nocache(). The only
164 difference is that SSAFEM::solve() recomputes the cached values of the coefficients before calling
165 SSAFEM::solve_nocache().
166 */
167 void SSAFEM::solve(const Inputs &inputs) {
168
169 TerminationReason::Ptr reason = solve_with_reason(inputs);
170 if (reason->failed()) {
171 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "SSAFEM solve failed to converge (SNES reason %s)",
172 reason->description().c_str());
173 } else if (m_log->get_threshold() > 2) {
174 m_stdout_ssa += "SSAFEM converged (SNES reason " + reason->description() + ")";
175 }
176 }
177
178 TerminationReason::Ptr SSAFEM::solve_with_reason(const Inputs &inputs) {
179
180 // Set up the system to solve.
181 cache_inputs(inputs);
182
183 return solve_nocache();
184 }
185
186 //! Solve the SSA without first recomputing the values of coefficients at quad
187 //! points. See the disccusion of SSAFEM::solve for more discussion.
188 TerminationReason::Ptr SSAFEM::solve_nocache() {
189 PetscErrorCode ierr;
190
191 m_epsilon_ssa = m_config->get_number("stress_balance.ssa.epsilon");
192
193 options::String filename("-ssa_view", "");
194 if (filename.is_set()) {
195 petsc::Viewer viewer;
196 ierr = PetscViewerASCIIOpen(m_grid->com, filename->c_str(), viewer.rawptr());
197 PISM_CHK(ierr, "PetscViewerASCIIOpen");
198
199 ierr = PetscViewerASCIIPrintf(viewer, "SNES before SSASolve_FE\n");
200 PISM_CHK(ierr, "PetscViewerASCIIPrintf");
201
202 ierr = SNESView(m_snes, viewer);
203 PISM_CHK(ierr, "SNESView");
204
205 ierr = PetscViewerASCIIPrintf(viewer, "solution vector before SSASolve_FE\n");
206 PISM_CHK(ierr, "PetscViewerASCIIPrintf");
207
208 ierr = VecView(m_velocity_global.vec(), viewer);
209 PISM_CHK(ierr, "VecView");
210 }
211
212 m_stdout_ssa.clear();
213 if (m_log->get_threshold() >= 2) {
214 m_stdout_ssa = " SSA: ";
215 }
216
217 // Solve:
218 ierr = SNESSolve(m_snes, NULL, m_velocity_global.vec());
219 PISM_CHK(ierr, "SNESSolve");
220
221 // See if it worked.
222 SNESConvergedReason snes_reason;
223 ierr = SNESGetConvergedReason(m_snes, &snes_reason); PISM_CHK(ierr, "SNESGetConvergedReason");
224
225 TerminationReason::Ptr reason(new SNESTerminationReason(snes_reason));
226 if (not reason->failed()) {
227
228 // Extract the solution back from SSAX to velocity and communicate.
229 m_velocity.copy_from(m_velocity_global);
230 m_velocity.update_ghosts();
231
232 bool view_solution = options::Bool("-ssa_view_solution", "view solution of the SSA system");
233 if (view_solution) {
234 petsc::Viewer viewer;
235 ierr = PetscViewerASCIIOpen(m_grid->com, filename->c_str(), viewer.rawptr());
236 PISM_CHK(ierr, "PetscViewerASCIIOpen");
237
238 ierr = PetscViewerASCIIPrintf(viewer, "solution vector after SSASolve\n");
239 PISM_CHK(ierr, "PetscViewerASCIIPrintf");
240
241 ierr = VecView(m_velocity_global.vec(), viewer);
242 PISM_CHK(ierr, "VecView");
243 }
244
245 }
246
247 return reason;
248 }
249
250 //! Initialize stored data from the coefficients in the SSA. Called by SSAFEM::solve.
251 /**
252 This method is should be called after SSAFEM::init and whenever any geometry or temperature
253 related coefficients have changed. The method stores the values of the coefficients the nodes of
254 each element so that these are available to the residual and Jacobian evaluation methods.
255
256 We store the vertical average of the ice hardness to avoid re-doing this computation every
257 iteration.
258
259 In addition to coefficients at element nodes we store "node types" used to identify interior
260 elements, exterior elements, and boundary faces.
261 */
262 void SSAFEM::cache_inputs(const Inputs &inputs) {
263
264 // Hold on to pointers to the B.C. mask and values: they are needed in SNES callbacks and
265 // inputs.bc_{mask,values} are not available there.
266 m_bc_mask = inputs.bc_mask;
267 m_bc_values = inputs.bc_values;
268
269 const std::vector<double> &z = m_grid->z();
270
271 IceModelVec::AccessList list{&m_coefficients,
272 inputs.enthalpy,
273 &inputs.geometry->ice_thickness,
274 &inputs.geometry->bed_elevation,
275 &inputs.geometry->sea_level_elevation,
276 inputs.basal_yield_stress};
277
278 bool use_explicit_driving_stress = (m_driving_stress_x != NULL) && (m_driving_stress_y != NULL);
279 if (use_explicit_driving_stress) {
280 list.add({m_driving_stress_x, m_driving_stress_y});
281 }
282
283 ParallelSection loop(m_grid->com);
284 try {
285 for (Points p(*m_grid); p; p.next()) {
286 const int i = p.i(), j = p.j();
287
288 double thickness = inputs.geometry->ice_thickness(i, j);
289
290 Vector2 tau_d;
291 if (use_explicit_driving_stress) {
292 tau_d.u = (*m_driving_stress_x)(i, j);
293 tau_d.v = (*m_driving_stress_y)(i, j);
294 } else {
295 // tau_d above is set to zero by the Vector2
296 // constructor, but is not used
297 }
298
299 const double *enthalpy = inputs.enthalpy->get_column(i, j);
300 double hardness = rheology::averaged_hardness(*m_flow_law, thickness,
301 m_grid->kBelowHeight(thickness),
302 &z[0], enthalpy);
303
304 Coefficients c;
305 c.thickness = thickness;
306 c.bed = inputs.geometry->bed_elevation(i, j);
307 c.sea_level = inputs.geometry->sea_level_elevation(i, j);
308 c.tauc = (*inputs.basal_yield_stress)(i, j);
309 c.hardness = hardness;
310 c.driving_stress = tau_d;
311
312 m_coefficients(i, j) = c;
313 } // loop over owned grid points
314 } catch (...) {
315 loop.failed();
316 }
317 loop.check();
318
319 m_coefficients.update_ghosts();
320
321 const bool use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
322 if (use_cfbc) {
323 // Note: the call below uses ghosts of inputs.geometry->ice_thickness.
324 compute_node_types(inputs.geometry->ice_thickness,
325 m_config->get_number("stress_balance.ice_free_thickness_standard"),
326 m_node_type);
327 } else {
328 m_node_type.set(NODE_INTERIOR);
329 }
330
331 cache_residual_cfbc(inputs);
332
333 }
334
335 //! Compute quadrature point values of various coefficients given a quadrature `Q` and nodal values.
336 void SSAFEM::quad_point_values(const fem::Quadrature &Q,
337 const Coefficients *x,
338 int *mask,
339 double *thickness,
340 double *tauc,
341 double *hardness) const {
342 const fem::Germs *test = Q.test_function_values();
343 const unsigned int n = Q.n();
344
345 for (unsigned int q = 0; q < n; q++) {
346 double
347 bed = 0.0,
348 sea_level = 0.0;
349
350 thickness[q] = 0.0;
351 tauc[q] = 0.0;
352 hardness[q] = 0.0;
353
354 for (unsigned int k = 0; k < fem::q1::n_chi; k++) {
355 const fem::Germ &psi = test[q][k];
356
357 thickness[q] += psi.val * x[k].thickness;
358 bed += psi.val * x[k].bed;
359 sea_level += psi.val * x[k].sea_level;
360 tauc[q] += psi.val * x[k].tauc;
361 hardness[q] += psi.val * x[k].hardness;
362 }
363
364 mask[q] = m_gc.mask(sea_level, bed, thickness[q]);
365 }
366 }
367
368 //! Compute gravitational driving stress at quadrature points.
369 //! Uses explicitly-provided nodal values.
370 void SSAFEM::explicit_driving_stress(const fem::Quadrature &Q,
371 const Coefficients *x,
372 Vector2 *result) const {
373 const fem::Germs *test = Q.test_function_values();
374 const unsigned int n = Q.n();
375
376 for (unsigned int q = 0; q < n; q++) {
377 result[q] = 0.0;
378
379 for (unsigned int k = 0; k < fem::q1::n_chi; k++) {
380 const fem::Germ &psi = test[q][k];
381 result[q] += psi.val * x[k].driving_stress;
382 }
383 }
384 }
385
386 /** Compute the the driving stress at quadrature points.
387 The surface elevation @f$ h @f$ is defined by
388 @f{equation}{
389 h =
390 \begin{cases}
391 b + H, & \mathrm{grounded},\\
392 z_{sl} + \alpha H, & \mathrm{shelf},
393 \end{cases}
394 @f}
395 where @f$ b @f$ is the bed elevation, @f$ H @f$ is the ice thickness, and @f$ z_{sl} @f$ is the
396 sea level, and @f$ \alpha @f$ is the ratio of densities of ice and sea water.
397
398 Note that if @f$ b @f$, @f$ H @f$, and @f$ z_{sl} @f$ are bilinear (or @f$ C^{1} @f$ in
399 general), @f$ h @f$ is continuous but not continuously differentiable. In
400 other words, the gradient of @f$ h @f$ is not continuous, so near the
401 grounding line we cannot compute the gravitational driving stress
402 @f$ \tau_{d} = - \rho g H \nabla h @f$ using the @f$ Q_1 @f$ or @f$ P_1 @f$ FE basis
403 expansion of @f$ h @f$.
404
405 We differentiate the equation above instead, getting
406 @f{equation}{
407 \nabla h =
408 \begin{cases}
409 \nabla b + \nabla H, & \mathrm{grounded},\\
410 \nabla z_{sl} + \alpha \nabla H, & \mathrm{shelf}.
411 \end{cases}
412 @f}
413
414 and use basis expansions of @f$ \nabla b @f$ and @f$ \nabla H @f$.
415
416 Overall, we have
417
418 @f{align*}{
419 \tau_{d} &= \rho g H \nabla h\\
420 &=
421 - \rho g H
422 \begin{cases}
423 \nabla b + \nabla H, & \mathrm{grounded},\\
424 \alpha \nabla H, & \mathrm{shelf},
425 \end{cases}
426 @f}
427
428 because @f$ z = z_{sl} @f$ defines the geoid surface and so *its gradient
429 does not contribute to the driving stress*.
430 */
431 void SSAFEM::driving_stress(const fem::Quadrature &Q,
432 const Coefficients *x,
433 Vector2 *result) const {
434 const fem::Germs *test = Q.test_function_values();
435 const unsigned int n = Q.n();
436
437 for (unsigned int q = 0; q < n; q++) {
438 double
439 sea_level = 0.0,
440 b = 0.0,
441 b_x = 0.0,
442 b_y = 0.0,
443 H = 0.0,
444 H_x = 0.0,
445 H_y = 0.0;
446
447 result[q] = 0.0;
448
449 for (unsigned int k = 0; k < fem::q1::n_chi; k++) {
450 const fem::Germ &psi = test[q][k];
451
452 b += psi.val * x[k].bed;
453 b_x += psi.dx * x[k].bed;
454 b_y += psi.dy * x[k].bed;
455
456 H += psi.val * x[k].thickness;
457 H_x += psi.dx * x[k].thickness;
458 H_y += psi.dy * x[k].thickness;
459
460 sea_level += psi.val * x[k].sea_level;
461 }
462
463 const int M = m_gc.mask(sea_level, b, H);
464 const bool grounded = mask::grounded(M);
465
466 const double
467 pressure = m_rho_g * H,
468 h_x = grounded ? b_x + H_x : m_alpha * H_x,
469 h_y = grounded ? b_y + H_y : m_alpha * H_y;
470
471 result[q].u = - pressure * h_x;
472 result[q].v = - pressure * h_y;
473 }
474 }
475
476
477 /** @brief Compute the "(regularized effective viscosity) x (ice thickness)" and effective viscous
478 * bed strength from the current solution, at a single quadrature point.
479 *
480 * @param[in] thickness ice thickness
481 * @param[in] hardness ice hardness
482 * @param[in] mask cell type mask
483 * @param[in] tauc basal yield stress
484 * @param[in] U the value of the solution
485 * @param[in] U_x x-derivatives of velocity components
486 * @param[in] U_y y-derivatives of velocity components
487 * @param[out] nuH product of the ice viscosity and thickness @f$ \nu H @f$
488 * @param[out] dnuH derivative of @f$ \nu H @f$ with respect to the
489 * second invariant @f$ \gamma @f$. Set to NULL if
490 * not desired.
491 * @param[out] beta basal drag coefficient @f$ \beta @f$
492 * @param[out] dbeta derivative of @f$ \beta @f$ with respect to the
493 * second invariant @f$ \gamma @f$. Set to NULL if
494 * not desired.
495 */
496 void SSAFEM::PointwiseNuHAndBeta(double thickness,
497 double hardness,
498 int mask,
499 double tauc,
500 const Vector2 &U,
501 const Vector2 &U_x,
502 const Vector2 &U_y,
503 double *nuH, double *dnuH,
504 double *beta, double *dbeta) {
505
506 if (thickness < strength_extension->get_min_thickness()) {
507 *nuH = strength_extension->get_notional_strength();
508 if (dnuH) {
509 *dnuH = 0;
510 }
511 } else {
512 m_flow_law->effective_viscosity(hardness, secondInvariant_2D(U_x, U_y),
513 nuH, dnuH);
514
515 *nuH = m_epsilon_ssa + *nuH * thickness;
516
517 if (dnuH) {
518 *dnuH *= thickness;
519 }
520 }
521
522 if (mask::grounded_ice(mask)) {
523 m_basal_sliding_law->drag_with_derivative(tauc, U.u, U.v, beta, dbeta);
524 } else {
525 *beta = 0;
526
527 if (mask::ice_free_land(mask)) {
528 *beta = m_beta_ice_free_bedrock;
529 }
530
531 if (dbeta) {
532 *dbeta = 0;
533 }
534 }
535 }
536
537
538 //! Compute and cache residual contributions from the integral over the lateral boundary.
539 /**
540
541 This method computes FIXME.
542
543 */
544 void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {
545
546 fem::BoundaryQuadrature2 bq(m_grid->dx(), m_grid->dy(), 1.0);
547
548 const unsigned int Nk = fem::q1::n_chi;
549 const unsigned int Nq = bq.n();
550 const unsigned int n_sides = fem::q1::n_sides;
551 using mask::ocean;
552
553 const Vector2 *outward_normal = fem::q1::outward_normals();
554
555 const bool
556 use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc"),
557 is_dry_simulation = m_config->get_flag("ocean.always_grounded");
558
559 const double
560 ice_density = m_config->get_number("constants.ice.density"),
561 ocean_density = m_config->get_number("constants.sea_water.density"),
562 standard_gravity = m_config->get_number("constants.standard_gravity");
563
564 // Reset the boundary integral so that all values are overwritten.
565 m_boundary_integral.set(0.0);
566
567 if (not use_cfbc) {
568 // If CFBC is not used then we're done.
569 return;
570 }
571
572 IceModelVec::AccessList list{&m_node_type,
573 &inputs.geometry->ice_thickness,
574 &inputs.geometry->bed_elevation,
575 &inputs.geometry->sea_level_elevation,
576 &m_boundary_integral};
577
578 // Iterate over the elements.
579 const int
580 xs = m_element_index.xs,
581 xm = m_element_index.xm,
582 ys = m_element_index.ys,
583 ym = m_element_index.ym;
584
585 ParallelSection loop(m_grid->com);
586 try {
587 for (int j = ys; j < ys + ym; j++) {
588 for (int i = xs; i < xs + xm; i++) {
589 // Initialize the map from global to element degrees of freedom.
590 m_element.reset(i, j);
591
592 int node_type[Nk];
593 m_element.nodal_values(m_node_type, node_type);
594
595 // an element is "interior" if all its nodes are interior or boundary
596 const bool interior_element = (node_type[0] < NODE_EXTERIOR and
597 node_type[1] < NODE_EXTERIOR and
598 node_type[2] < NODE_EXTERIOR and
599 node_type[3] < NODE_EXTERIOR);
600
601 // residual contributions at element nodes
602 std::vector<Vector2> I(Nk);
603
604 if (not interior_element) {
605 // not an interior element: the contribution is zero
606 continue;
607 }
608
609 double H_nodal[Nk];
610 m_element.nodal_values(inputs.geometry->ice_thickness, H_nodal);
611
612 double b_nodal[Nk];
613 m_element.nodal_values(inputs.geometry->bed_elevation, b_nodal);
614
615 double sl_nodal[Nk];
616 m_element.nodal_values(inputs.geometry->sea_level_elevation, sl_nodal);
617
618 // storage for test function values psi[0] for the first "end" of a side, psi[1] for the
619 // second
620 double psi[2] = {0.0, 0.0};
621
622 // loop over element sides
623 for (unsigned int side = 0; side < n_sides; ++side) {
624
625 // nodes incident to the current side
626 const int
627 n0 = fem::q1::incident_nodes[side][0],
628 n1 = fem::q1::incident_nodes[side][1];
629
630 if (not (node_type[n0] == NODE_BOUNDARY and node_type[n1] == NODE_BOUNDARY)) {
631 // not a boundary side; skip it
632 continue;
633 }
634
635 // in our case (i.e. uniform spacing in x and y directions) W is the same at all
636 // quadrature points along a side.
637 const double W = bq.weights(side);
638
639 for (unsigned int q = 0; q < Nq; ++q) {
640
641 // test functions at nodes incident to the current side, evaluated at the quadrature point
642 // q
643 psi[0] = bq.germ(side, q, n0).val;
644 psi[1] = bq.germ(side, q, n1).val;
645
646 // Compute ice thickness and bed elevation at a quadrature point. This uses a 1D basis
647 // expansion on the side.
648 const double
649 H = H_nodal[n0] * psi[0] + H_nodal[n1] * psi[1],
650 bed = b_nodal[n0] * psi[0] + b_nodal[n1] * psi[1],
651 sea_level = sl_nodal[n0] * psi[0] + sl_nodal[n1] * psi[1];
652
653 const bool floating = ocean(m_gc.mask(sea_level, bed, H));
654
655 // ocean pressure difference at a quadrature point
656 const double dP = margin_pressure_difference(floating, is_dry_simulation,
657 H, bed, sea_level,
658 ice_density, ocean_density,
659 standard_gravity);
660
661 // This integral contributes to the residual at 2 nodes (the ones incident to the
662 // current side). This is is written in a way that allows *adding* (... += ...) the
663 // boundary contribution in the residual computation.
664 I[n0] += W * (- psi[0] * dP) * outward_normal[side];
665 I[n1] += W * (- psi[1] * dP) * outward_normal[side];
666 // FIXME: I need to include the special case corresponding to ice margins next
667 // to fjord walls, nunataks, etc. In this case dP == 0.
668 //
669 // FIXME: set pressure difference to zero at grounded locations at domain
670 // boundaries.
671 } // q-loop
672
673 } // loop over element sides
674
675 m_element.add_contribution(&I[0], m_boundary_integral);
676
677 } // i-loop
678 } // j-loop
679 } catch (...) {
680 loop.failed();
681 }
682 loop.check();
683 }
684
685
686 //! Implements the callback for computing the residual.
687 /*!
688 * Compute the residual \f[r_{ij}= G(x, \psi_{ij}) \f] where \f$G\f$
689 * is the weak form of the SSA, \f$x\f$ is the current approximate
690 * solution, and the \f$\psi_{ij}\f$ are test functions.
691 *
692 * The weak form of the SSA system is
693 */
694 void SSAFEM::compute_local_function(Vector2 const *const *const velocity_global,
695 Vector2 **residual_global) {
696
697 const bool use_explicit_driving_stress = (m_driving_stress_x != NULL) && (m_driving_stress_y != NULL);
698
699 const bool use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
700
701 const unsigned int Nk = fem::q1::n_chi;
702 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
703
704 IceModelVec::AccessList list{&m_node_type, &m_coefficients, &m_boundary_integral};
705
706 // Set the boundary contribution of the residual. This is computed at the nodes, so we don't want
707 // to set it using ElementMap::add_contribution() because that would lead to
708 // double-counting. Also note that without CFBC m_boundary_integral is exactly zero.
709 for (Points p(*m_grid); p; p.next()) {
710 const int i = p.i(), j = p.j();
711
712 residual_global[j][i] = m_boundary_integral(i, j);
713 }
714
715 // Start access to Dirichlet data if present.
716 fem::DirichletData_Vector dirichlet_data(m_bc_mask, m_bc_values, m_dirichletScale);
717
718 // Storage for the current solution and its derivatives at quadrature points.
719 Vector2 U[Nq_max], U_x[Nq_max], U_y[Nq_max];
720
721 // Iterate over the elements.
722 const int
723 xs = m_element_index.xs,
724 xm = m_element_index.xm,
725 ys = m_element_index.ys,
726 ym = m_element_index.ym;
727
728 ParallelSection loop(m_grid->com);
729 try {
730 for (int j = ys; j < ys + ym; j++) {
731 for (int i = xs; i < xs + xm; i++) {
732 // Initialize the map from global to element degrees of freedom.
733 m_element.reset(i, j);
734 int node_type[Nk];
735 m_element.nodal_values(m_node_type, node_type);
736
737 // an element is "interior" if all its nodes are interior or boundary
738 const bool interior_element = (node_type[0] < NODE_EXTERIOR and
739 node_type[1] < NODE_EXTERIOR and
740 node_type[2] < NODE_EXTERIOR and
741 node_type[3] < NODE_EXTERIOR);
742
743 if (use_cfbc and (not interior_element)) {
744 // an exterior element in the CFBC case
745 continue;
746 }
747
748 // Note: without CFBC all elements are "interior".
749
750 fem::Quadrature &Q = m_quadrature;
751
752 // Number of quadrature points.
753 const unsigned int Nq = Q.n();
754
755 // An Nq by Nk array of test function values.
756 const fem::Germs *test = Q.test_function_values();
757
758 // Jacobian times weights for quadrature.
759 const double* W = Q.weights();
760
761 // Storage for the solution and residuals at element nodes.
762 Vector2 residual[Nk];
763
764 int mask[Nq_max];
765 double thickness[Nq_max];
766 double tauc[Nq_max];
767 double hardness[Nq_max];
768 Vector2 tau_d[Nq_max];
769
770 {
771 Coefficients coeffs[Nk];
772 m_element.nodal_values(m_coefficients, coeffs);
773
774 quad_point_values(Q, coeffs, mask, thickness, tauc, hardness);
775
776 if (use_explicit_driving_stress) {
777 explicit_driving_stress(Q, coeffs, tau_d);
778 } else {
779 driving_stress(Q, coeffs, tau_d);
780 }
781 }
782
783 {
784 // Obtain the value of the solution at the nodes adjacent to the element.
785 Vector2 velocity_nodal[Nk];
786 m_element.nodal_values(velocity_global, velocity_nodal);
787
788 // These values now need to be adjusted if some nodes in the element have Dirichlet data.
789 if (dirichlet_data) {
790 // Set elements of velocity_nodal that correspond to Dirichlet nodes to prescribed
791 // values.
792 dirichlet_data.enforce(m_element, velocity_nodal);
793 // mark Dirichlet nodes in m_element so that they are not touched by
794 // add_contribution() below
795 dirichlet_data.constrain(m_element);
796 }
797
798 // Compute the solution values and its gradient at the quadrature points.
799 quadrature_point_values(Q, velocity_nodal, // input
800 U, U_x, U_y); // outputs
801 }
802
803 // Zero out the element-local residual in preparation for updating it.
804 for (unsigned int k = 0; k < Nk; k++) {
805 residual[k].u = 0;
806 residual[k].v = 0;
807 }
808
809 // loop over quadrature points:
810 for (unsigned int q = 0; q < Nq; q++) {
811
812 double eta = 0.0, beta = 0.0;
813 PointwiseNuHAndBeta(thickness[q], hardness[q], mask[q], tauc[q],
814 U[q], U_x[q], U_y[q], // inputs
815 &eta, NULL, &beta, NULL); // outputs
816
817 // The next few lines compute the actual residual for the element.
818 const Vector2 tau_b = U[q] * (- beta); // basal shear stress
819
820 const double
821 jw = W[q],
822 u_x = U_x[q].u,
823 v_y = U_y[q].v,
824 u_y_plus_v_x = U_y[q].u + U_x[q].v;
825
826 // Loop over test functions.
827 for (unsigned int k = 0; k < Nk; k++) {
828 const fem::Germ &psi = test[q][k];
829
830 residual[k].u += jw * (eta * (psi.dx * (4.0 * u_x + 2.0 * v_y) + psi.dy * u_y_plus_v_x)
831 - psi.val * (tau_b.u + tau_d[q].u));
832 residual[k].v += jw * (eta * (psi.dx * u_y_plus_v_x + psi.dy * (2.0 * u_x + 4.0 * v_y))
833 - psi.val * (tau_b.v + tau_d[q].v));
834 } // k (test functions)
835 } // q (quadrature points)
836
837 m_element.add_contribution(residual, residual_global);
838 } // i-loop
839 } // j-loop
840 } catch (...) {
841 loop.failed();
842 }
843 loop.check();
844
845 // Until now we have not touched rows in the residual corresponding to Dirichlet data.
846 // We fix this now.
847 if (dirichlet_data) {
848 dirichlet_data.fix_residual(velocity_global, residual_global);
849 }
850
851 if (use_cfbc) {
852 // Prescribe homogeneous Dirichlet B.C. at ice-free nodes. This uses the fact that m_node_type
853 // can be used as a "Dirichlet B.C. mask", i.e. ice-free nodes (and only ice-free nodes) are
854 // marked with ones.
855 fem::DirichletData_Vector dirichlet_ice_free(&m_node_type, NULL, m_dirichletScale);
856 dirichlet_ice_free.fix_residual_homogeneous(residual_global);
857 }
858
859 monitor_function(velocity_global, residual_global);
860 }
861
862 void SSAFEM::monitor_function(Vector2 const *const *const velocity_global,
863 Vector2 const *const *const residual_global) {
864 PetscErrorCode ierr;
865 bool monitorFunction = options::Bool("-ssa_monitor_function", "monitor the SSA residual");
866 if (not monitorFunction) {
867 return;
868 }
869
870 ierr = PetscPrintf(m_grid->com,
871 "SSA Solution and Function values (pointwise residuals)\n");
872 PISM_CHK(ierr, "PetscPrintf");
873
874 ParallelSection loop(m_grid->com);
875 try {
876 for (Points p(*m_grid); p; p.next()) {
877 const int i = p.i(), j = p.j();
878
879 ierr = PetscSynchronizedPrintf(m_grid->com,
880 "[%2d, %2d] u=(%12.10e, %12.10e) f=(%12.4e, %12.4e)\n",
881 i, j,
882 velocity_global[j][i].u, velocity_global[j][i].v,
883 residual_global[j][i].u, residual_global[j][i].v);
884 PISM_CHK(ierr, "PetscSynchronizedPrintf");
885 }
886 } catch (...) {
887 loop.failed();
888 }
889 loop.check();
890
891 ierr = PetscSynchronizedFlush(m_grid->com, NULL);
892 PISM_CHK(ierr, "PetscSynchronizedFlush");
893 }
894
895
896 //! Implements the callback for computing the Jacobian.
897 /*!
898 Compute the Jacobian
899
900 @f[ J_{ij}{kl} \frac{d r_{ij}}{d x_{kl}}= G(x, \psi_{ij}) @f]
901
902 where \f$G\f$ is the weak form of the SSA, \f$x\f$ is the current
903 approximate solution, and the \f$\psi_{ij}\f$ are test functions.
904
905 */
906 void SSAFEM::compute_local_jacobian(Vector2 const *const *const velocity_global, Mat Jac) {
907
908 const unsigned int Nk = fem::q1::n_chi;
909 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
910
911 const bool use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
912
913 // Zero out the Jacobian in preparation for updating it.
914 PetscErrorCode ierr = MatZeroEntries(Jac);
915 PISM_CHK(ierr, "MatZeroEntries");
916
917 IceModelVec::AccessList list{&m_node_type, &m_coefficients};
918
919 // Start access to Dirichlet data if present.
920 fem::DirichletData_Vector dirichlet_data(m_bc_mask, m_bc_values, m_dirichletScale);
921
922 // Storage for the current solution at quadrature points.
923 Vector2 U[Nq_max], U_x[Nq_max], U_y[Nq_max];
924
925 // Loop through all the elements.
926 int
927 xs = m_element_index.xs,
928 xm = m_element_index.xm,
929 ys = m_element_index.ys,
930 ym = m_element_index.ym;
931
932 ParallelSection loop(m_grid->com);
933 try {
934 for (int j = ys; j < ys + ym; j++) {
935 for (int i = xs; i < xs + xm; i++) {
936 // Initialize the map from global to element degrees of freedom.
937 m_element.reset(i, j);
938
939 int node_type[Nk];
940 m_element.nodal_values(m_node_type, node_type);
941 // an element is "interior" if all its nodes are interior or boundary
942 const bool interior_element = (node_type[0] < NODE_EXTERIOR and
943 node_type[1] < NODE_EXTERIOR and
944 node_type[2] < NODE_EXTERIOR and
945 node_type[3] < NODE_EXTERIOR);
946
947 if (use_cfbc and (not interior_element)) {
948 // an exterior element in the CFBC case
949 continue;
950 }
951
952 fem::Quadrature &Q = m_quadrature;
953
954 // Number of quadrature points.
955 const unsigned int Nq = Q.n();
956
957 // Jacobian times weights for quadrature.
958 const double* W = Q.weights();
959
960 // Values of the finite element test functions at the quadrature points.
961 // This is an Nq by Nk array of function germs
962 const fem::Germs *test = Q.test_function_values();
963
964 int mask[Nq_max];
965 double thickness[Nq_max];
966 double tauc[Nq_max];
967 double hardness[Nq_max];
968
969 {
970 Coefficients coeffs[Nk];
971 m_element.nodal_values(m_coefficients, coeffs);
972
973 quad_point_values(Q, coeffs,
974 mask, thickness, tauc, hardness);
975 }
976
977 {
978 // Values of the solution at the nodes of the current element.
979 Vector2 velocity_nodal[Nk];
980 // Obtain the value of the solution at the adjacent nodes to the element.
981 m_element.nodal_values(velocity_global, velocity_nodal);
982
983 // These values now need to be adjusted if some nodes in the element have
984 // Dirichlet data.
985 if (dirichlet_data) {
986 dirichlet_data.enforce(m_element, velocity_nodal);
987 dirichlet_data.constrain(m_element);
988 }
989 // Compute the values of the solution at the quadrature points.
990 quadrature_point_values(Q, velocity_nodal, U, U_x, U_y);
991 }
992
993 // Element-local Jacobian matrix (there are Nk vector valued degrees
994 // of freedom per element, for a total of (2*Nk)*(2*Nk) = 16
995 // entries in the local Jacobian.
996 double K[2*Nk][2*Nk];
997 // Build the element-local Jacobian.
998 ierr = PetscMemzero(K, sizeof(K));
999 PISM_CHK(ierr, "PetscMemzero");
1000
1001 for (unsigned int q = 0; q < Nq; q++) {
1002 const double
1003 jw = W[q],
1004 u = U[q].u,
1005 v = U[q].v,
1006 u_x = U_x[q].u,
1007 v_y = U_y[q].v,
1008 u_y_plus_v_x = U_y[q].u + U_x[q].v;
1009
1010 double eta = 0.0, deta = 0.0, beta = 0.0, dbeta = 0.0;
1011 PointwiseNuHAndBeta(thickness[q], hardness[q], mask[q], tauc[q],
1012 U[q], U_x[q], U_y[q],
1013 &eta, &deta, &beta, &dbeta);
1014
1015 for (unsigned int l = 0; l < Nk; l++) { // Trial functions
1016
1017 // Current trial function and its derivatives:
1018 const fem::Germ &phi = test[q][l];
1019
1020 // Derivatives of \gamma with respect to u_l and v_l:
1021 const double
1022 gamma_u = (2.0 * u_x + v_y) * phi.dx + 0.5 * u_y_plus_v_x * phi.dy,
1023 gamma_v = 0.5 * u_y_plus_v_x * phi.dx + (u_x + 2.0 * v_y) * phi.dy;
1024
1025 // Derivatives of \eta = \nu*H with respect to u_l and v_l:
1026 const double
1027 eta_u = deta * gamma_u,
1028 eta_v = deta * gamma_v;
1029
1030 // Derivatives of the basal shear stress term (\tau_b):
1031 const double
1032 taub_xu = -dbeta * u * u * phi.val - beta * phi.val, // x-component, derivative with respect to u_l
1033 taub_xv = -dbeta * u * v * phi.val, // x-component, derivative with respect to u_l
1034 taub_yu = -dbeta * v * u * phi.val, // y-component, derivative with respect to v_l
1035 taub_yv = -dbeta * v * v * phi.val - beta * phi.val; // y-component, derivative with respect to v_l
1036
1037 for (unsigned int k = 0; k < Nk; k++) { // Test functions
1038
1039 // Current test function and its derivatives:
1040
1041 const fem::Germ &psi = test[q][k];
1042
1043 if (eta == 0) {
1044 ierr = PetscPrintf(PETSC_COMM_SELF, "eta=0 i %d j %d q %d k %d\n", i, j, q, k);
1045 PISM_CHK(ierr, "PetscPrintf");
1046 }
1047
1048 // u-u coupling
1049 K[k*2 + 0][l*2 + 0] += jw * (eta_u * (psi.dx * (4 * u_x + 2 * v_y) + psi.dy * u_y_plus_v_x)
1050 + eta * (4 * psi.dx * phi.dx + psi.dy * phi.dy) - psi.val * taub_xu);
1051 // u-v coupling
1052 K[k*2 + 0][l*2 + 1] += jw * (eta_v * (psi.dx * (4 * u_x + 2 * v_y) + psi.dy * u_y_plus_v_x)
1053 + eta * (2 * psi.dx * phi.dy + psi.dy * phi.dx) - psi.val * taub_xv);
1054 // v-u coupling
1055 K[k*2 + 1][l*2 + 0] += jw * (eta_u * (psi.dx * u_y_plus_v_x + psi.dy * (2 * u_x + 4 * v_y))
1056 + eta * (psi.dx * phi.dy + 2 * psi.dy * phi.dx) - psi.val * taub_yu);
1057 // v-v coupling
1058 K[k*2 + 1][l*2 + 1] += jw * (eta_v * (psi.dx * u_y_plus_v_x + psi.dy * (2 * u_x + 4 * v_y))
1059 + eta * (psi.dx * phi.dx + 4 * psi.dy * phi.dy) - psi.val * taub_yv);
1060
1061 } // l
1062 } // k
1063 } // q
1064 m_element.add_contribution(&K[0][0], Jac);
1065 } // j
1066 } // i
1067 } catch (...) {
1068 loop.failed();
1069 }
1070 loop.check();
1071
1072
1073 // Until now, the rows and columns correspoinding to Dirichlet data
1074 // have not been set. We now put an identity block in for these
1075 // unknowns. Note that because we have takes steps to not touching
1076 // these columns previously, the symmetry of the Jacobian matrix is
1077 // preserved.
1078 if (dirichlet_data) {
1079 dirichlet_data.fix_jacobian(Jac);
1080 }
1081
1082 if (use_cfbc) {
1083 // Prescribe homogeneous Dirichlet B.C. at ice-free nodes. This uses the fact that m_node_type
1084 // can be used as a "Dirichlet B.C. mask", i.e. ice-free nodes (and only ice-free nodes) are
1085 // marked with ones.
1086 fem::DirichletData_Vector dirichlet_ice_free(&m_node_type, NULL, m_dirichletScale);
1087 dirichlet_ice_free.fix_jacobian(Jac);
1088 }
1089
1090 ierr = MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
1091 PISM_CHK(ierr, "MatAssemblyBegin");
1092
1093 ierr = MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
1094 PISM_CHK(ierr, "MatAssemblyEnd");
1095
1096 ierr = MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
1097 PISM_CHK(ierr, "MatSetOption");
1098
1099 ierr = MatSetOption(Jac, MAT_SYMMETRIC, PETSC_TRUE);
1100 PISM_CHK(ierr, "MatSetOption");
1101
1102 monitor_jacobian(Jac);
1103 }
1104
1105 void SSAFEM::monitor_jacobian(Mat Jac) {
1106 PetscErrorCode ierr;
1107 bool mon_jac = options::Bool("-ssa_monitor_jacobian", "monitor the SSA Jacobian");
1108
1109 if (not mon_jac) {
1110 return;
1111 }
1112
1113 // iter has to be a PetscInt because it is used in the
1114 // SNESGetIterationNumber() call below.
1115 PetscInt iter = 0;
1116 ierr = SNESGetIterationNumber(m_snes, &iter);
1117 PISM_CHK(ierr, "SNESGetIterationNumber");
1118
1119 char file_name[PETSC_MAX_PATH_LEN];
1120 snprintf(file_name, PETSC_MAX_PATH_LEN, "PISM_SSAFEM_J%d.m", (int)iter);
1121
1122 m_log->message(2,
1123 "writing Matlab-readable file for SSAFEM system A xsoln = rhs to file `%s' ...\n",
1124 file_name);
1125
1126 petsc::Viewer viewer(m_grid->com);
1127
1128 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);
1129 PISM_CHK(ierr, "PetscViewerSetType");
1130
1131 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1132 PISM_CHK(ierr, "PetscViewerPushFormat");
1133
1134 ierr = PetscViewerFileSetName(viewer, file_name);
1135 PISM_CHK(ierr, "PetscViewerFileSetName");
1136
1137 ierr = PetscObjectSetName((PetscObject) Jac, "A");
1138 PISM_CHK(ierr, "PetscObjectSetName");
1139
1140 ierr = MatView(Jac, viewer);
1141 PISM_CHK(ierr, "MatView");
1142
1143 ierr = PetscViewerPopFormat(viewer);
1144 PISM_CHK(ierr, "PetscViewerPopFormat");
1145 }
1146
1147 //!
1148 PetscErrorCode SSAFEM::function_callback(DMDALocalInfo *info,
1149 Vector2 const *const *const velocity,
1150 Vector2 **residual,
1151 CallbackData *fe) {
1152 try {
1153 (void) info;
1154 fe->ssa->compute_local_function(velocity, residual);
1155 } catch (...) {
1156 MPI_Comm com = MPI_COMM_SELF;
1157 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)fe->da, &com); CHKERRQ(ierr);
1158 handle_fatal_errors(com);
1159 SETERRQ(com, 1, "A PISM callback failed");
1160 }
1161 return 0;
1162 }
1163
1164 PetscErrorCode SSAFEM::jacobian_callback(DMDALocalInfo *info,
1165 Vector2 const *const *const velocity,
1166 Mat A, Mat J, CallbackData *fe) {
1167 try {
1168 (void) A;
1169 (void) info;
1170 fe->ssa->compute_local_jacobian(velocity, J);
1171 } catch (...) {
1172 MPI_Comm com = MPI_COMM_SELF;
1173 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)fe->da, &com); CHKERRQ(ierr);
1174 handle_fatal_errors(com);
1175 SETERRQ(com, 1, "A PISM callback failed");
1176 }
1177 return 0;
1178 }
1179
1180 } // end of namespace stressbalance
1181 } // end of namespace pism