PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
|
◆ compute_diffusivity()
Compute the SIA diffusivity. If full_update, also store delta on the staggered grid. Recall that Q = -D \nabla h is the diffusive flux in the mass-continuity equation \frac{\partial H}{\partial t} = M - S - \nabla \cdot (Q + \mathbf{U}_b H), where h is the ice surface elevation, M is the top surface accumulation/ablation rate, S is the basal mass balance and \mathbf{U}_b is the thickness-advective (in PISM: usually SSA) ice velocity. Recall also that at any particular point in the map-plane (i.e. if x and y are fixed) D = 2\int_b^h F(z)P(z)(h-z)dz, where F(z) is a constitutive function and P(z) is the pressure at a level z. By defining \delta(z) = 2F(z)P(z) one can write D = \int_b^h\delta(z)(h-z)dz. The advantage is that it is then possible to avoid re-evaluating F(z) (which is computationally expensive) in the horizontal ice velocity (see compute_3d_horizontal_velocity()) computation. This method computes D and stores \delta in delta[0,1] if full_update is true. The trapezoidal rule is used to approximate the integral.
Definition at line 529 of file SIAFD.cc. References pism::Geometry::cell_type, pism::ParallelSection::check(), pism::Time::current(), pism::ParallelSection::failed(), pism::RuntimeError::formatted(), pism::array::Array3D::get_column(), pism::GlobalMax(), pism::GlobalSum(), pism::Geometry::ice_surface_elevation, pism::Geometry::ice_thickness, interglacial(), pism::k, m_bed_smoother, pism::Component::m_config, pism::stressbalance::SSB_Modifier::m_D_max, m_delta_0, m_delta_1, m_e_factor, m_e_factor_interglacial, pism::stressbalance::SSB_Modifier::m_EC, pism::stressbalance::SSB_Modifier::m_flow_law, pism::Component::m_grid, pism::Component::m_log, m_seconds_per_year, m_work_2d_0, m_work_2d_1, PISM_ERROR_LOCATION, pism::array::Array::set(), pism::array::Array3D::set_column(), pism::stressbalance::BedSmoother::smoothed_thk(), pism::array::Array::stencil_width(), pism::stressbalance::BedSmoother::theta(), pism::Component::time(), pism::grid::X_PERIODIC, and pism::grid::Y_PERIODIC. Referenced by update(). |