assemble_R checkReadyToSolve compute_enthalpy_CTS compute_lambda Enth Enth_s enthSystemCtx init k_from_T lambda m_B0 m_B_ks m_D0 m_D_ks m_E_e m_E_ij m_E_n m_E_s m_E_w m_EC m_Enth m_Enth3 m_Enth_s m_ice_c m_ice_density m_ice_k m_ice_thickness m_k_depends_on_T m_L_ks m_lambda m_margin_exclude_horizontal_advection m_margin_exclude_strain_heat m_margin_exclude_vertical_advection m_marginal m_nu m_p_air m_R m_R_cold m_R_factor m_R_temp m_strain_heating m_strain_heating3 m_U0 m_U_ks save_system set_basal_dirichlet_bc set_basal_heat_flux set_basal_neumann_bc set_surface_dirichlet_bc set_surface_heat_flux set_surface_neumann_bc solve ~enthSystemCtx

solve() void pism::energy::enthSystemCtx::solve ( std::vector< double > &  result ) Solve the tridiagonal system, in a single column, which determines the new values of the ice enthalpy. We are solving a convection-diffusion equation, treating the \( z \) direction implicitly and \( x, y \) directions explicitly. See bombproofenth for the documentation of the treatment of convection terms. The notes below document the way we treat diffusion using a simplified PDE. To discretize \[ \diff{}{z} \left( K(E) \diff{E}{z}\right) = \diff{E}{t} \] at a location \( k \) of the vertical grid, we use centered finite differences in space, backward differences in time, and evaluate \( K(E) \) at staggered-grid locations: \[ \frac{K_{k+\frac12}\frac{E_{k+1} - E_{k}}{\dz} - K_{k-\frac12}\frac{E_{k} - E_{k-1}}{\dz}}{\dz} = \frac{E_{k} - E^{n-1}_{k}}{\dt}, \] where \( E = E^{n} \) for clarity and \( K_{k\pm \frac12} = K(E^{n-1}_{k\pm \frac12}) \), i.e. we linearize this PDE by evaluating \( K(E) \) at the previous time-step. We define \[ R_i = \frac{\dt\, K_i}{\dz^2}, \] and the discretization takes form \[ -R_{k-\frac12} E_{k-1} + \left( 1 + R_{k-\frac12} + R_{k+\frac12} \right) E_{k} - R_{k+\frac12} E_{k+1} = E^{n-1}_{k}. \] In the assembly of the tridiagonal system this corresponds to \begin{align*} L_i &= - \frac12 (R_{i} + R_{i-1}),\\ D_i &= 1 + \frac12 (R_{i} + R_{i-1}) + \frac12 (R_{i} + R_{i+1}),\\ U_i &= - \frac12 (R_{i} + R_{i+1}),\\ b_i &= E^{n-1}_{i}, \end{align*} where \( L_i, D_i, U_i \) are lower-diagonal, diagonal, and upper-diagonal entries corresponding to an equation \( i \) and \( b_i \) is the corresponding right-hand side. (Staggered-grid values are approximated by interpolating from the regular grid). This method is unconditionally stable and has a maximum principle (see [MortonMayers, section 2.11]). Definition at line 480 of file enthSystem.cc. References pism::RuntimeError::add_context(), checkReadyToSolve(), pism::k, m_B0, m_B_ks, m_D0, m_D_ks, pism::columnSystemCtx::m_dt, pism::columnSystemCtx::m_dx, pism::columnSystemCtx::m_dy, m_E_e, m_E_n, m_E_s, m_E_w, m_Enth, pism::columnSystemCtx::m_i, m_ice_density, pism::columnSystemCtx::m_j, pism::columnSystemCtx::m_ks, m_L_ks, m_lambda, m_margin_exclude_horizontal_advection, m_margin_exclude_strain_heat, m_marginal, m_nu, m_R, pism::columnSystemCtx::m_solver, m_strain_heating, pism::columnSystemCtx::m_u, m_U0, m_U_ks, pism::columnSystemCtx::m_v, pism::columnSystemCtx::m_w, pism::columnSystemCtx::m_z, PISM_ERROR_LOCATION, pism::columnSystemCtx::reportColumnZeroPivotErrorMFile(), S(), and pism::energy::upwind(). Referenced by pism::energy::CHSystem::update_impl(), and pism::energy::EnthalpyModel::update_impl().