Loading web-font TeX/Math/Italic
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages

◆ compute_diffusivity()

void pism::stressbalance::SIAFD::compute_diffusivity ( bool  full_update,
const Geometry geometry,
const array::Array3D enthalpy,
const array::Array3D age,
const array::Staggered1 h_x,
const array::Staggered1 h_y,
array::Staggered1 result 
)
protectedvirtual

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.

Parameters
[in]full_updatethe flag specitying if we're doing a "full" update.
[in]h_xx-component of the surface gradient, on the staggered grid
[in]h_yy-component of the surface gradient, on the staggered grid
[out]resultdiffusivity of the SIA flow

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().