assemble_matrix assemble_rhs compute_hardav_staggered compute_nuH_norm compute_nuH_staggered compute_nuH_staggered_cfbc diagnostics_impl fracture_induced_softening init_impl integrated_viscosity is_marginal m_A m_b m_default_pc_failure_count m_default_pc_failure_max_count m_hardness m_KSP m_nuH m_nuH_old m_nuh_viewer m_nuh_viewer_size m_scaling m_velocity_old m_view_nuh m_work pc_setup_asm pc_setup_bjacobi picard_iteration picard_manager picard_strategy_regularization solve SSAFD update_nuH_viewers write_system_petsc ~SSAFD

solve() void pism::stressbalance::SSAFD::solve ( const Inputs &  inputs ) protectedvirtual Compute the vertically-averaged horizontal velocity from the shallow shelf approximation. This is the main procedure in the SSAFD. It manages the nonlinear solve process and the Picard iteration. The outer loop (over index k) is the nonlinear iteration. In this loop the effective viscosity is computed by compute_nuH_staggered() and then the linear system is set up and solved. Specifically, we call the PETSc procedure KSPSolve() to solve the linear system. Solving the linear system is also a loop, an iteration, but it occurs inside KSPSolve(). The user has full control of the KSP solve through the PETSc interface. The default choicess for KSP type -ksp_type and preconditioner type -pc_type are GMRES(30) for the former and block Jacobi with ILU on the blocks for the latter. The defaults usually work. These choices are important but poorly understood. The eigenvalues of the linearized SSA vary with ice sheet geometry and temperature in ways that are not well-studied. Nonetheless these eigenvalues determine the convergence of this (inner) linear iteration. A well-chosen preconditioner can put the eigenvalues in the right place so that the KSP can converge quickly. The preconditioner will behave differently on different numbers of processors. If the user wants the results of SSA calculations to be independent of the number of processors, then -pc_type none could be used, but performance will be poor. If you want to test different KSP methods, it may be helpful to see how many iterations were necessary. Use -ksp_monitor. Initial testing implies that CGS takes roughly half the iterations of GMRES(30), but is not significantly faster because the iterations are each roughly twice as slow. The outputs of PETSc options -ksp_monitor_singular_value, -ksp_compute_eigenvalues and -ksp_plot_eigenvalues -draw_pause N (the last holds plots for N seconds) may be useful to diagnose. The outer loop terminates when the effective viscosity times thickness is no longer changing much, according to the tolerance set by the option -ssafd_picard_rtol. The outer loop also terminates when a maximum number of iterations is exceeded. We save the velocity from the last time step in order to have a better estimate of the effective viscosity than the u=v=0 result. In truth there is an "outer outer" loop (over index l). This attempts to over-regularize the effective viscosity if the nonlinear iteration (the "outer" loop over k) is not converging with the default regularization. The same over-regularization is attempted if the KSP object reports that it has not converged. (An alternative for recovery in the KSP diverged case, suggested by Jed, is to revert to a direct linear solve, either for the whole domain (not scalable) or on the subdomains. This recovery alternative requires a more nontrivial choice but it may be worthwhile. Note the user can already do -pc_type asm -sub_pc_type lu at the command line, forcing subdomain direct solves.) FIXME: update this doxygen comment Implements pism::stressbalance::SSA. Definition at line 940 of file SSAFD.cc. References assemble_rhs(), pism::Geometry::cell_type, compute_hardav_staggered(), pism::array::Array2D< T >::copy_from(), pism::stressbalance::SSA::extrapolate_velocity(), pism::stressbalance::Inputs::geometry, pism::k, pism::Component::m_config, pism::Component::m_log, pism::stressbalance::ShallowStressBalance::m_velocity, m_velocity_old, picard_iteration(), picard_strategy_regularization(), PISM_ERROR_LOCATION, pism::array::Array::scale(), pism::array::Array::update_ghosts(), and write_system_petsc().