PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800

◆ 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::array::max(), 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().