m_equations.py

class m_equations.Equations(MeshClass, ElemClass, TracersClass, MeltingClass, FilesClass, t)

Bases: object

compute_u()

Var:

Updates the heat fluxes \(\boldsymbol{q}_{\rm ice}\) and \(\boldsymbol{q}_{\rm ocean}\) and computes new phase transition velocity \(\boldsymbol{u}\).

equation_Stokes()

Var:

Defines the weak form of the Stokes problem

Returns:

Weak form of the continuity equation

\[0 = \int_\Omega \nabla\cdot \boldsymbol{v}\ \hat{p}\ \textrm{d}\Omega \]

Weak form of the momentum equation if elasticity = False

\[\begin{split}0 = \int_\Omega \hat{p}\nabla\cdot \hat{\boldsymbol{v}} - \rho g (\boldsymbol{e}_z \cdot \hat{\boldsymbol{v}}) -\eta\ (\nabla \boldsymbol{v} + \nabla^T \boldsymbol{v}): \nabla^T \hat{\boldsymbol{v}}\ \textrm{d}\Omega\\ \end{split}\]

Weak form of the momentum equation if elasticity = True

\[\begin{split}0 = \int_\Omega \hat{p}\nabla\cdot \hat{\boldsymbol{v}} - \rho g (\boldsymbol{e}_z \cdot \hat{\boldsymbol{v}}) - Z\eta\ (\nabla \boldsymbol{v} + \nabla^T \boldsymbol{v}): \nabla^T \hat{\boldsymbol{v}} - (1-Z) (\boldsymbol{\sigma} : \nabla^T \hat{\boldsymbol{v}})\ \textrm{d}\Omega\\ \end{split}\]

In case of a free surface condition at the top boundary BC_Stokes_problem[0][0] = "free_surface", the ‘drunken sailor’ stabilization term is added

\[- \int_{\Gamma_{\rm top}} \rho_s g\lambda\Delta t (\boldsymbol{v} \cdot \boldsymbol{n})(\boldsymbol{e}_z \cdot \hat{\boldsymbol{v}}) \ \textrm{d}\Gamma_{\rm top}\]

In case of a free surface condition at the bottom boundary BC_Stokes_problem[1][0] = "free_surface", the hydrostatic pressure and the ‘drunken sailor’ stabilization term are added

\[+ \int_{\Gamma_{\rm bot}} (h_{\rm bot}\rho_l - h\rho_s)(\boldsymbol{n} \cdot \hat{\boldsymbol{v}}) + (\rho_l - \rho_s) g\lambda\Delta t [(\boldsymbol{v} + \boldsymbol{u}) \cdot \boldsymbol{n}](\boldsymbol{e}_z \cdot \hat{\boldsymbol{v}}) \ \textrm{d}\Gamma_{\rm bot}\]

equation_free_surface_bottom()

Var:

Defines the weak form, LinearVariationalProblem and LinearVariationalSolver for the bottom surface evolution equation.

Returns:

\[\begin{split}\begin{align} 0 = &\int_\Omega \nabla h_{\rm bot} \cdot \nabla \hat{h}_{\rm bot}\ \textrm{d}\Omega\\ &- \int_{\Gamma_{\rm bot}} \frac{\partial h_{\rm bot}}{\partial \boldsymbol{n}} \hat{h}_{\rm bot}\ \textrm{d}\Gamma_{\rm bot}\\ &- \int_{\Gamma_{\rm bot}} \left[h_{\rm bot} - h_{\rm bot}^k - \Delta t \left(\frac{\partial h_{\rm bot}}{\partial x} (v_x + u_x) - (v_z + u_z) \right)\right] \left(\frac{\partial \hat{h}_{\rm bot}}{\partial \boldsymbol{n}} - \frac{\hat{h}_{\rm bot}}{\gamma h_e} \right)\ \textrm{d}\Gamma_{\rm bot} \end{align}\end{split}\]

where

• \(h_{\rm bot}\) is the bottom topography trial function

• \(h_{\rm bot}^k\) is the bottom topography function from the previous time step

• \(\hat{h}_{\rm bot}\) is the bottom topography test function

• \(v_x\) and \(v_z\) are the horizontal and vertical components of the velocity

• \(u_x\) and \(u_z\) are the horizontal and vertical components of the phase change velocity

• \(\Delta t\) is the time step length

• \(\gamma\) is a stabilization parameter

• \(h_e\) is the element size at the bottom boundary

equation_free_surface_top()

Var:

Defines the weak form, LinearVariationalProblem and LinearVariationalSolver for the top surface evolution equation.

Returns:

\[\begin{split}\begin{align} 0 = &\int_\Omega \nabla h_{\rm top} \cdot \nabla \hat{h}_{\rm top}\ \textrm{d}\Omega\\ &- \int_{\Gamma_{\rm top}} \frac{\partial h_{\rm top}}{\partial \boldsymbol{n}} \hat{h}_{\rm top}\ \textrm{d}\Gamma_{\rm top}\\ &- \int_{\Gamma_{\rm top}} \left[h_{\rm top} - h_{\rm top}^k - \Delta t \left(\frac{\partial h_{\rm top}}{\partial x} v_x - v_z \right)\right] \left(\frac{\partial \hat{h}_{\rm top}}{\partial \boldsymbol{n}} - \frac{\hat{h}_{\rm top}}{\gamma h_e} \right)\ \textrm{d}\Gamma_{\rm top} \end{align}\end{split}\]

where

• \(h_{\rm top}\) is the top topography trial function

• \(h_{\rm top}^k\) is the top topography function from the previous time step

• \(\hat{h}_{\rm top}\) is the top topography test function

• \(v_x\) and \(v_z\) are the horizontal and vertical components of the velocity

• \(\Delta t\) is the time step length

• \(\gamma\) is a stabilization parameter

• \(h_e\) is the element size at the top boundary

equation_heat()

Var:

Defines the weak form of the time-dependent heat transfer equation.

Returns:

\[\begin{split}\begin{align} 0 = &\int_\Omega \rho_s c_p \frac{T - T^k}{\Delta t}\hat{T}\\ &+ 0.5[\rho_s c_p((\boldsymbol{v} - \boldsymbol{v}_{\rm mesh}) \cdot \nabla T)\ \hat{T} + k\ (\nabla T \cdot \nabla \hat{T})]\\ &+ 0.5[\rho_s c_p((\boldsymbol{v} - \boldsymbol{v}_{\rm mesh}) \cdot \nabla T^k)\ \hat{T} + k\ (\nabla T^k \cdot \nabla \hat{T})]\ \textrm{d}\Omega\\ \end{align}\end{split}\]

where

• \(T\) is the temperature trial function

• \(T^k\) is the temperature function from the previous time step

• \(\hat{T}\) is the temperature test function

• \(\boldsymbol{v}\) is the velocity

• \(\boldsymbol{v}_{\rm mesh}\) is the mesh velocity

• \(\Delta t\) is the time step length

If tidal_dissipation = True, the corresponding volumetric heating term defined in m_rheology.tidal_heating() is added

\[- \int_\Omega H_{\rm tidal}(\eta)\ \hat{T}\ \textrm{d}\Omega\]

In case of radiation boundary condition at the top boundary BC_T_top = "radiation", the corresponding boundary condition is added

\[- \int_{\Gamma_{\rm top}} (1 - A)I\hat{T}\ + T^4\varepsilon\sigma_{SB}\hat{T}\ \textrm{d}\Gamma_{\rm top}\]

where

• \(I\) is the insolation at the body’s surface (in m_parameters)

• \(A\) is the albedo of the body’s surface (in m_parameters)

• \(\varepsilon\) is the body’s surface emmisivity (in m_parameters)

• \(\sigma_{SB}\) is the Stefan-Boltzmann constant (in m_constants)

Note

If BC_heat_transfer[0][0] = "radiation", nonlinear solver will be automatically used.

equation_heat_ini()

Var:

Defines the weak form of the stationary heat transfer equation.

\[0 = \int_\Omega k\ (\nabla T \cdot \nabla \hat{T})\ \textrm{d}\Omega\]

If BC_heat_transfer[0][0] = "radiation", the corresponding boundary condition is added

\[- \int_{\Gamma_{\rm top}} I(1 - A)\hat{T}\ + T^4\varepsilon\sigma_{SB}\hat{T}\ \textrm{d}\Gamma_{\rm top}\]

where

• \(I\) is the insolation at the body’s surface

• \(A\) is the albedo of the body’s surface

• \(\varepsilon\) is the body’s surface emmisivity

• \(\sigma_{SB}\) is the Stefan-Boltzmann constant

Note

If BC_heat_transfer[0][0] = "radiation", nonlinear solver will be automatically used.

equation_mesh_displacement()

Var:

Defines the weak form, LinearVariationalProblem and LinearVariationalSolver for mesh displacement distribution equation.

Returns:

If mesh_displacement_laplace = "full" then

\[0 = \int_\Omega \nabla h_2 \cdot \nabla \hat{h}_2\ \textrm{d}\Omega\]

If mesh_displacement_laplace = "z-only" then

\[0 = \int_\Omega \frac{\partial h_2}{\partial z}\frac{\partial \hat{h}_2}{\partial z} \ \textrm{d}\Omega\]

where

• \(h_2\) is the mesh displacement trial function

• \(\hat{h}_2\) is the mesh displacement test function

rotate_and_interpolate_stress()

Computes the vorticity tensor \(\boldsymbol{W} = (\nabla \boldsymbol{v} - (\nabla \boldsymbol{v})^T)/2\), rotates the elastic stress and interpolates it at the mesh. See m_interpolation.stress_rotation().

solve_Stokes_problem(step, FilesClass)

Solves the Stokes problem.

solve_heat_equation()

solve_initial_heat_equation()

solve_topography_evolution()

thermal_perturbation()

Var:

Perturbs the initial temperature profile.

Parameters:

• perturb_ampl – amplitude of the thermal perturbation (\(A\))

• perturb_freq – spatial frequency of the thermal perturbation (\(f\))

Returns:

\[T^\prime = T - A\sin\left(\frac{\pi y}{h}\right) \cdot\cos\left(\frac{2\pi fx}{l}\right)\]

update_stress()

update_viscosity(step, Picard_iter)

Var:

Updates the viscosity fuction