Diffusion

Inheritance diagram of climlab.dynamics.Diffusion

General solver of the 1D diffusion equation:

\[\begin{split}\frac{\partial}{\partial t} \Psi(x,t) &= -\frac{1}{w(x)} \frac{\partial}{\partial x} \left[ w(x) ~ F(x,t) \right] \\ F &= -K ~ \frac{\partial \Psi}{\partial x}\end{split}\]

for a state variable \(\Psi(x,t)\) and arbitrary diffusivity \(K(x,t)\) in units of \(x^2 ~ t^{-1}\).

\(w(x)\) is an optional weighting function for the divergence operator on curvilinear grids.

The diffusivity \(K\) can be a single scalar, or optionally a vector specified at grid cell boundaries (so its length must be exactly 1 greater than the length of \(x\)).

\(K\) can be modified by the user at any time (e.g., after each timestep, if it depends on other state variables).

A fully implicit timestep is used for computational efficiency. Thus the computed tendency \(\frac{\partial \Psi}{\partial t}\) will depend on the timestep.

In addition to the tendency over the implicit timestep, the solver also calculates two diagnostics from the updated state:

  • diffusive_flux given by \(F(x)\) in units of \([\Psi]~[x]\)/s

  • diffusive_flux_convergence given by the right hand side of the first equation above, in units of \([\Psi]\)/s

This base class can be used without modification for diffusion in Cartesian coordinates (\(w=1\)) on a regularly spaced grid.

The state variable \(\Psi\) may be multi-dimensional, but the diffusion will operate along a single dimension only.

Other classes implement the weighting for spherical geometry.

class climlab.dynamics.diffusion.Diffusion(K=None, diffusion_axis=None, use_banded_solver=False, **kwargs)[source]

Bases: climlab.process.implicit.ImplicitProcess

A parent class for one dimensional implicit diffusion modules.

Initialization parameters

Parameters
  • K (float) – the diffusivity parameter in units of \(\frac{[\textrm{length}]^2}{\textrm{time}}\) where length is the unit of the spatial axis on which the diffusion is occuring.

  • diffusion_axis (str) – dictionary key for axis on which the diffusion is occuring in process’s domain axes dictionary

  • use_banded_solver (bool) – input flag, whether to use scipy.linalg.solve_banded() instead of numpy.linalg.solve() [default: False]

Note

The banded solver scipy.linalg.solve_banded() is faster than numpy.linalg.solve() but only works for one dimensional diffusion.

Object attributes

Additional to the parent class ImplicitProcess following object attributes are generated or modified during initialization:

Variables
  • param (dict) – parameter dictionary is extended by diffusivity parameter K (unit: \(\frac{[\textrm{length}]^2}{\textrm{time}}\))

  • use_banded_solver (bool) – input flag specifying numerical solving method (given during initialization)

  • diffusion_axis (str) – dictionary key for axis where diffusion is occuring: specified during initialization or output of method _guess_diffusion_axis()

  • _K_dimensionless (array) – diffusion parameter K multiplied by the timestep and divided by mean of diffusion axis delta in the power of two. Array has the size of diffusion axis bounds. \(K_{\textrm{dimensionless}}[i]= K \frac{\Delta t}{ \left(\overline{\Delta \textrm{bounds}} \right)^2}\)

  • _diffTriDiag (array) – tridiagonal diffusion matrix made by _make_diffusion_matrix() with input self._K_dimensionless

Example

Here is an example showing implementation of a vertical diffusion. It shows that a subprocess can work on just a subset of the parent process state variables.

import climlab
from climlab.dynamics.diffusion import Diffusion
import matplotlib.pyplot as plt

c = climlab.GreyRadiationModel()
K = 0.5
d = Diffusion(K=K, state = {'Tatm':c.state['Tatm']}, **c.param)

c.add_subprocess('diffusion',d)

### Integrate & Plot ###

fig = plt.figure( figsize=(6,4))
ax = fig.add_subplot(111)

ax.plot(c.lev, c.state['Tatm'], label='step 0')
c.step_forward()
ax.plot(c.lev, c.state['Tatm'], label='step 1')

ax.invert_xaxis()
ax.set_title('Diffusion subprocess')
ax.set_xlabel('level (mb)')
#ax.set_xticks([])
ax.set_ylabel('temperature (K)')
ax.legend(loc='best')
plt.show()

(Source code, png, hires.png, pdf)

../_images/example_diffusion.png
Attributes
K
depth

Depth at grid centers (m)

depth_bounds

Depth at grid interfaces (m)

diagnostics

Dictionary access to all diagnostic variables

input

Dictionary access to all input variables

lat

Latitude of grid centers (degrees North)

lat_bounds

Latitude of grid interfaces (degrees North)

lev

Pressure levels at grid centers (hPa or mb)

lev_bounds

Pressure levels at grid interfaces (hPa or mb)

lon

Longitude of grid centers (degrees)

lon_bounds

Longitude of grid interfaces (degrees)

timestep

The amount of time over which step_forward() is integrating in unit seconds.

Methods

add_diagnostic(name[, value])

Create a new diagnostic variable called name for this process and initialize it with the given value.

add_input(name[, value])

Create a new input variable called name for this process and initialize it with the given value.

add_subprocess(name, proc)

Adds a single subprocess to this process.

add_subprocesses(procdict)

Adds a dictionary of subproceses to this process.

compute()

Computes the tendencies for all state variables given current state and specified input.

compute_diagnostics([num_iter])

Compute all tendencies and diagnostics, but don’t update model state.

declare_diagnostics(diaglist)

Add the variable names in inputlist to the list of diagnostics.

declare_input(inputlist)

Add the variable names in inputlist to the list of necessary inputs.

integrate_converge([crit, verbose])

Integrates the model until model states are converging.

integrate_days([days, verbose])

Integrates the model forward for a specified number of days.

integrate_years([years, verbose])

Integrates the model by a given number of years.

remove_diagnostic(name)

Removes a diagnostic from the process.diagnostic dictionary and also delete the associated process attribute.

remove_subprocess(name[, verbose])

Removes a single subprocess from this process.

set_state(name, value)

Sets the variable name to a new state value.

set_timestep([timestep, num_steps_per_year])

Calculates the timestep in unit seconds and calls the setter function of timestep()

step_forward()

Updates state variables with computed tendencies.

to_xarray([diagnostics])

Convert process variables to xarray.Dataset format.

K
_implicit_solver()[source]

Invertes and solves the matrix problem for diffusion matrix and temperature T.

The method is called by the _compute() function of the ImplicitProcess class and solves the matrix problem

\[A \cdot T_{\textrm{new}} = T_{\textrm{old}}\]

for diffusion matrix A and corresponding temperatures. \(T_{\textrm{old}}\) is in this case the current state variable which already has been adjusted by the explicit processes. \(T_{\textrm{new}}\) is the new state of the variable. To derive the temperature tendency of the diffusion process the adjustment has to be calculated and muliplied with the timestep which is done by the _compute() function of the ImplicitProcess class.

This method calculates the matrix inversion for every state variable and calling either solve_implicit_banded() or numpy.linalg.solve() dependent on the flag self.use_banded_solver.

Variables
  • state (dict) – method uses current state variables but does not modify them

  • use_banded_solver (bool) – input flag whether to use _solve_implicit_banded() or numpy.linalg.solve() to do the matrix inversion

  • _diffTriDiag (array) – the diffusion matrix which is given with the current state variable to the method solving the matrix problem

_update_diagnostics(newstate)[source]

This method is called each timestep after the new state is computed with the implicit solver. Daughter classes can implement this method to compute any diagnostic quantities using the new state.

climlab.dynamics.diffusion._guess_diffusion_axis(process_or_domain)[source]

Scans given process, domain or dictionary of domains for a diffusion axis and returns appropriate name.

In case only one axis with length > 1 in the process or set of domains exists, the name of that axis is returned. Otherwise an error is raised.

Parameters

process_or_domain (Process, _Domain or dict of domains) – input from where diffusion axis should be guessed

Raises

ValueError if more than one diffusion axis is possible.

Returns

name of the diffusion axis

Return type

str

climlab.dynamics.diffusion._make_diffusion_matrix(K, weight1=None, weight2=None)[source]

Builds the general diffusion matrix with dimension nxn.

Note

\(n\) = number of points of diffusion axis \(n+1\) = number of bounts of diffusion axis

Function-all argument

Parameters
  • K (array) – dimensionless diffusivities at cell boundaries (size: 1xn+1)

  • weight1 (array) – weight_1 (size: 1xn+1)

  • weight2 (array) – weight_2 (size: 1xn)

Returns

completely listed tridiagonal diffusion matrix (size: nxn)

Return type

array

Note

The elements of array K are acutally dimensionless:

\[K[i] = K_{\textrm{physical}} \frac{\Delta t}{(\Delta y)^2}\]

where \(K_{\textrm{physical}}\) is in unit \(\frac{\textrm{length}^2}{\textrm{time}}\)

The diffusion matrix is build like the following

\[\begin{split}\textrm{diffTriDiag}= \left[ \begin{array}{cccccc} 1+\frac{s_1 }{w_{2,0}} & -\frac{s_1}{w_{2,0}} & 0 & & ... & 0 \\ -\frac{s_1}{w_{2,1}} & 1+\frac{s_1 + s_2}{w_{2,1}} & -\frac{s_2}{w_{2,1}} & 0 & ... & 0 \\ 0 & -\frac{s_2}{w_{2,2}} & 1+\frac{s_2 + s_3}{w_{2,2}} & -\frac{s_3}{w_{2,2}} &... & 0 \\ & & \ddots & \ddots & \ddots & \\ 0 & 0 & ... & -\frac{s_{n-2}}{w_{2,n-2}} & 1+\frac{s_{n-2} + s_{n-1}}{w_{2,{n-2}}} & -\frac{s_{n-1}}{w_{2,{n-2}}} \\ 0 & 0 & ... & 0 & -\frac{s_{n-1}}{w_{2,n-1}} & 1+\frac{s_{n-1}}{w_{2,n-1}} \\ \end{array} \right]\end{split}\]

where

\[\begin{split}\begin{array}{lllllll} K &= [K_0, &K_1, &K_2, &...,&K_{n-1}, &K_{n}] \\ w_1 &= [w_{1,0}, &w_{1,1},&w_{1,2},&...,&w_{1,n-1},&w_{1,n}] \\ w_2 &= [w_{2,0}, &w_{2,1},&w_{2,2},&...,&w_{2,n-1}] \end{array}\end{split}\]

and following subsitute:

\[s_i = w_{1,i} K_i\]
climlab.dynamics.diffusion._make_meridional_diffusion_matrix(K, lataxis)[source]

Calls _make_diffusion_matrix() with appropriate weights for the meridional diffusion case.

Parameters
  • K (array) – dimensionless diffusivities at cell boundaries of diffusion axis lataxis

  • lataxis (axis) – latitude axis where diffusion is occuring

Weights are computed as the following:

\[\begin{split}\begin{array}{ll} w_1 &= \cos(\textrm{bounds}) \\ &= \left[ \cos(b_0), \cos(b_1), \cos(b_2), \ ... \ , \cos(b_{n-1}), \cos(b_n) \right] \\ w_2 &= \cos(\textrm{points}) \\ &= \left[ \cos(p_0), \cos(p_1), \cos(p_2), \ ... \ , \cos(p_{n-1}) \right] \end{array}\end{split}\]

when bounds and points from lataxis are written as

\[\begin{split}\begin{array}{ll} \textrm{bounds} &= [b_0, b_1, b_2, \ ... \ , b_{n-1}, b_{n}] \\ \textrm{points} &= [p_0, p_1, p_2, \ ... \ , p_{n-1}] \end{array}\end{split}\]

Giving this input to _make_diffusion_matrix() results in a matrix like:

\[\begin{split}\textrm{diffTriDiag}= \left[ \begin{array}{cccccc} 1+\frac{u_1 }{\cos(p_0)} & -\frac{u_1}{\cos(p_0)} & 0 & & ... & 0 \\ -\frac{u_1}{\cos(p_1)} & 1+\frac{u_1 + u_2}{\cos(p_1)} & -\frac{u_2}{\cos(b_1)} & 0 & ... & 0 \\ 0 & -\frac{u_2}{\cos(p_2)} & 1+\frac{u_2 + u_3}{\cos(p_2)} & -\frac{u_3}{\cos(p_2)} &... & 0 \\ & & \ddots & \ddots & \ddots & \\ 0 & 0 & ... & -\frac{u_{n-2}}{\cos(p_{n-2})} & 1+\frac{u_{n-2} + u_{n-1}}{\cos(p_{n-2})} & -\frac{u_{n-1}}{\cos(p_{n-2})} \\ 0 & 0 & ... & 0 & -\frac{u_{n-1}}{\cos(p_{n-1})} & 1+\frac{u_{n-1}}{\cos(p_{n-1})} \\ \end{array} \right]\end{split}\]

with the substitue of:

\[u_i = \cos(b_i) K_i\]
climlab.dynamics.diffusion._solve_implicit_banded(current, banded_matrix)[source]

Uses a banded solver for matrix inversion of a tridiagonal matrix.

Converts the complete listed tridiagonal matrix (nxn) into a three row matrix (3xn) and calls scipy.linalg.solve_banded().

Parameters
  • current (array) – the current state of the variable for which matrix inversion should be computed

  • banded_matrix (array) – complete diffusion matrix (dimension: nxn)

Returns

output of scipy.linalg.solve_banded()

Return type

array