advective_velocity basal_frictional_heating compute_vertical_velocity compute_volumetric_strain_heating define_model_state_impl diagnostics_impl diffusive_flux init m_cfl_2d m_cfl_3d m_modifier m_shallow_stress_balance m_strain_heating m_w max_diffusivity max_timestep_cfl_2d max_timestep_cfl_3d modifier shallow stdout_report StressBalance ts_diagnostics_impl update velocity_u velocity_v velocity_w volumetric_strain_heating write_model_state_impl ~StressBalance

compute_volumetric_strain_heating() void pism::stressbalance::StressBalance::compute_volumetric_strain_heating ( const Inputs &  inputs ) protectedvirtual Computes the volumetric strain heating using horizontal velocity. Following the notation used in [BBssasliding], let \(u\) be a three-dimensional vector velocity field. Then the strain rate tensor \(D_{ij}\) is defined by \[ D_{ij} = \frac 12 \left(\diff{u_{i}}{x_{j}} + \diff{u_{j}}{x_{i}} \right), \] Where \(i\) and \(j\) range from \(1\) to \(3\). The flow law in the viscosity form states \[ \tau_{ij} = 2 \eta D_{ij}, \] and the nonlinear ice viscosity satisfies \[ 2 \eta = B(T) D^{(1/n) - 1}. \] Here \(D^{2}\) is defined by \(2D^{2} = D_{ij}D_{ij}\) (using the summation convention) and \(B(T) = A(T)^{-1/n}\) is the ice hardness. Now the volumetric strain heating is \[ \Sigma = \sum_{i,j=1}^{3}D_{ij}\tau_{ij} = 2 B(T) D^{(1/n) + 1}. \] We use an approximation of \(D_{ij}\) common in shallow ice models: we assume that horizontal derivatives of the vertical velocity are much smaller than \(z\) derivatives horizontal velocity components \(u\) and \(v\). (We drop \(w_x\) and \(w_y\) terms in \(D_{ij}\).) we use the incompressibility of ice to approximate \(w_z\): \[ w_z = - (u_x + v_y). \] Requires ghosts of u and v velocity components and uses the fact that u and v above the ice are filled using constant extrapolation. Resulting field does not have ghosts. Below is the Maxima code that produces the expression evaluated by D2(). derivabbrev : true; U : [u, v, w]; X : [x, y, z]; depends(U, X); gradef(w, x, 0); gradef(w, y, 0); gradef(w, z, -(diff(u, x) + diff(v, y))); d[i,j] := 1/2 * (diff(U[i], X[j]) + diff(U[j], X[i])); D : genmatrix(d, 3, 3), ratsimp, factor; tex('D = D); tex('D^2 = 1/2 * mat_trace(D . D)); Returns 0 on success Definition at line 501 of file StressBalance.cc. References pism::Geometry::cell_type, pism::ParallelSection::check(), pism::stressbalance::D2(), pism::stressbalance::Inputs::enthalpy, pism::rheology::FlowLaw::exponent(), pism::ParallelSection::failed(), pism::stressbalance::Inputs::geometry, pism::array::Array3D::get_column(), pism::rheology::FlowLaw::hardness_n(), pism::Geometry::ice_thickness, pism::k, pism::Component::m_grid, m_modifier, m_shallow_stress_balance, m_strain_heating, n, and PISM_CHK. Referenced by update().