src.spyice.models.advection_diffusion#

This module implements a 1D implicit finite-difference solver for the advection–diffusion equation:

\[a \frac{\partial U}{\partial t} + b \frac{\partial U}{\partial z} + \frac{\partial}{\partial z} \left(c \frac{\partial U}{\partial z}\right) + d \frac{\partial W}{\partial t} = 0\]

The model supports both temperature and salinity evolution and includes optional boundary-condition treatments (Stefan, Buffo, Voller).

—

Class: AdvectionDiffusion#

class AdvectionDiffusion(argument, X, source, X_initial, W, W_initial, thickness_index, w, dt, dz, nz, t_passed, S_IC, Stefan=False, Buffo=False, Voller=False, bc_neumann=None)#

Class representing a 1D advection–diffusion model.

Parameters:

  • argument (str) – Either "temperature" or "salinity".

  • X (ndarray) – Current field values.

  • source (float) – Source term.

  • X_initial (ndarray) – Field values from previous timestep.

  • W (ndarray) – Current liquid fraction.

  • W_initial (ndarray) – Previous liquid fraction.

  • thickness_index (int) – Index of ice thickness boundary.

  • w (ndarray) – Convective brine velocity.

  • dt (float) – Time step.

  • dz (float) – Spatial step.

  • nz (int) – Number of spatial grid cells.

  • t_passed (float) – Simulation time.

  • S_IC (float) – Initial salinity condition.

  • Stefan (bool) – Enable Stefan boundary formulation.

  • Buffo (bool) – Enable Buffo boundary formulation.

  • Voller (bool) – Enable Voller enthalpy formulation.

  • bc_neumann (float) – Optional Neumann boundary gradient.

Raises:

  • AssertionError – If invalid argument is provided.

  • AssertionError – If both Stefan and Buffo are enabled.

—

Methods#

TDMAsolver#

TDMAsolver(a, b, c, d)#

Solve a tridiagonal system using the Thomas algorithm.

Parameters:

  • a (array) – Sub-diagonal.

  • b (array) – Main diagonal.

  • c (array) – Super-diagonal.

  • d (array) – RHS vector.

Returns:

Solution vector.

Buffosolver#

Buffosolver(a, b, c, f)#

Modified tridiagonal solver enforcing diagonal dominance.

Parameters:

  • a (ndarray) – Lower diagonal.

  • b (ndarray) – Upper diagonal.

  • c (ndarray) – Main diagonal.

  • f (ndarray) – RHS vector.

Returns:

Solution vector.

Raises:

AssertionError – If matrix is not diagonally dominant.

set_up_tridiagonal#

set_up_tridiagonal()#

Construct diagonal vectors for the implicit scheme matrix. Handles separate logic for temperature and salinity, including upwinding and boundary adjustments.

assemble_tridiagonal#

assemble_tridiagonal()#

Assemble full matrix from diagonal vectors.

Returns:

Assembled system matrix.

modify_tridiagonal_voller_scheme#

modify_tridiagonal_voller_scheme()#

Modify diagonal dominance for mushy-layer cells when using the Voller enthalpy scheme.

voller_X_array_set_zero_to_melt#

voller_X_array_set_zero_to_melt(X, X_boundary)#

Force solution values above melting threshold to boundary value.

set_source_term#

set_source_term()#

Compute discretized source term depending on field type.

Returns:

Source contribution.

unknowns_matrix#

unknowns_matrix(temperature_melt, non_constant_physical_properties=False)#

Solve the implicit linear system for the updated field.

Parameters:

temperature_melt (float) – Melting temperature.

Returns:

  • X_new (ndarray) — Updated solution

  • X_wind (ndarray) — Wind-corrected solution

  • A_before_correction (ndarray) — System matrix

factor_1#

factor_1(argument, a, c, dt, dz, nz)#

Compute diffusion discretization factor.

Returns:

  • ndarray for temperature

  • tuple of ndarrays for salinity (central, plus, minus)

factor_2#

factor_2(a, b, dt, dz, nz)#

Compute advection discretization factor.

Returns:

ndarray

factor_3#

factor_3(a, d, nz)#

Compute phase-coupling factor.

Returns:

ndarray

—

Numerical Characteristics#

  • Implicit time discretization

  • Finite difference spatial scheme

  • Tridiagonal linear system

  • Upwind advection (salinity)

  • Optional enthalpy formulation (Voller)

  • Stefan and Buffo boundary options

—

Dependencies#

  • spyice.parameters.user_input

  • spyice.coefficients.update_coefficients

  • spyice.rhs.apply_boundary_condition

  • spyice.rhs.correct_for_brine_movement

—

Notes#

Only one of Stefan or Buffo may be active at a time.

Diagonal dominance is required for the Buffo solver.

Zero-division protection is implemented in factor calculations.