Merge pull request #111 from brian-rose/advdiff

New MeridionalAdvectionDiffusion class

climlab/__init__.py

@@ -7,7 +7,7 @@
 has been documented for a comprehensive understanding and traceability.
 '''
 
-__version__ = '0.7.4.dev1'
+__version__ = '0.7.4.dev2'
 
 # this should ensure that we can still import constants.py as climlab.constants
 from .utils import constants, thermo, legendre

climlab/dynamics/__init__.py

@@ -3,10 +3,15 @@
 
 :class:`~climlab.dynamics.BudykoTransport` is a relaxation to global mean.
 
-Other modules are 1D advection- diffusion solvers (implemented using implicit timestepping).
+Other modules are 1D advection-diffusion solvers (implemented using implicit timestepping).
 
 :class:`~climlab.dynamics.AdvectionDiffusion` is a general-purpose 1D
-advection-diffusion process.
+advection-diffusion process. It can be used out-of-the-box for models with
+Cartesian grid geometry, but also accepts weighting functions for the
+divergence operator on curvilinear grids.
+
+:class:`~climlab.dynamics.MeridionalAdvectionDiffusion` implements the
+1D advection-diffusion process on the sphere (flux in the north-south direction).
 
 Subclass :class:`~climlab.dynamics.MeridionalHeatDiffusion` is the appropriate class
 for the traditional diffusive EBM, in which transport is parameterized as a
@@ -19,6 +24,6 @@
 from __future__ import absolute_import
 from .budyko_transport import BudykoTransport
 from .advection_diffusion import AdvectionDiffusion, Diffusion
-from .meridional_diffusion import MeridionalDiffusion
+from .meridional_advection_diffusion import MeridionalAdvectionDiffusion, MeridionalDiffusion
 from .meridional_heat_diffusion import MeridionalHeatDiffusion
 from .meridional_moist_diffusion import MeridionalMoistDiffusion

climlab/dynamics/adv_diff_numerics.py

@@ -206,6 +206,31 @@
 Solving for the future value :math:`\boldsymbol{\psi}^{n+1}` is then accomplished
 by solving the :math:`J \times J` tridiagonal linear system using standard routines.
 
+Analytical benchmark
+--------------------
+
+Here is an analytical case to be used for testing purposes to validate the numerical code.
+This is implemented in the CLIMLAB test suite.
+
+- :math:`K=K_0` is constant
+- :math:`w(x) = 1` everywhere (Cartesian coordinates)
+- :math:`F = 0` everywhere
+- :math:`\psi(x,0) = \psi_0 \sin^2\left(\frac{\pi x}{L}\right)`
+- :math:`u(x) = U_0 \sin\left(\frac{\pi x}{L}\right)`
+for a domain with endpoints at :math:`x=0` and :math:`x=L`.
+
+The analytical solution is
+
+.. math::
+
+    \begin{align}
+    \mathcal{F} &= \psi_0 \sin\left(\frac{\pi x}{L}\right) \left[U_0 \sin^2\left(\frac{\pi x}{L}\right) - 2K \frac{\pi}{L} \cos\left(\frac{\pi x}{L}\right) \right]  \\
+    \frac{\partial \psi}{\partial t} &= -\psi_0 \frac{\pi}{L} \left\{ 3 U_0 \sin^2\left(\frac{\pi x}{L}\right) \cos\left(\frac{\pi x}{L}\right) -2K\frac{\pi}{L} \left[\cos^2\left(\frac{\pi x}{L}\right) -\sin^2\left(\frac{\pi x}{L}\right) \right] \right\}
+    \end{align}
+
+which satisfies the boundary condition :math:`\mathcal{F} = 0` at :math:`x=0` and :math:`x=L`.
+
+
 Module function reference
 -------------------------
 

climlab/dynamics/meridional_advection_diffusion.py

@@ -0,0 +1,81 @@
+r"""General solver of the 1D meridional advection-diffusion equation on the sphere:
+
+.. math::
+
+    \frac{\partial}{\partial t} \psi(\phi,t) &= -\frac{1}{a \cos\phi} \frac{\partial}{\partial \phi} \left[ \cos\phi ~ F(\phi,t) \right] \\
+    F &= U(\phi) \psi(\phi) -\frac{K(\phi)}{a} ~ \frac{\partial \psi}{\partial \phi}
+
+for a state variable :math:`\psi(\phi,t)`, arbitrary diffusivity :math:`K(\phi)`
+in units of :math:`x^2 ~ t^{-1}`, and advecting velocity :math:`U(\phi)`.
+:math:`\phi` is latitude and :math:`a` is the Earth's radius (in meters).
+
+:math:`K` and :math:`U` can be scalars,
+or optionally vector *specified at grid cell boundaries*
+(so their lengths must be exactly 1 greater than the length of :math:`\phi`).
+
+:math:`K` and :math:`U` can be modified by the user at any time
+(e.g., after each timestep, if they depend on other state variables).
+
+A fully implicit timestep is used for computational efficiency. Thus the computed
+tendency :math:`\frac{\partial \psi}{\partial t}` will depend on the timestep.
+
+In addition to the tendency over the implicit timestep,
+the solver also calculates several diagnostics from the updated state:
+
+- ``diffusive_flux`` given by :math:`-\frac{K(\phi)}{a} ~ \frac{\partial \psi}{\partial \phi}` in units of :math:`[\psi]~[x]`/s
+- ``advective_flux`` given by :math:`U(\phi) \psi(\phi)` (same units)
+- ``total_flux``, the sum of advective, diffusive and prescribed fluxes
+- ``flux_convergence`` (or instantanous scalar tendency) given by the right hand side of the first equation above, in units of :math:`[\psi]`/s
+
+Non-uniform grid spacing is supported.
+
+The state variable :math:`\psi` may be multi-dimensional, but the diffusion
+will operate along the latitude dimension only.
+"""
+from __future__ import division
+import numpy as np
+from .advection_diffusion import AdvectionDiffusion, Diffusion
+from climlab import constants as const
+
+
+class MeridionalAdvectionDiffusion(AdvectionDiffusion):
+    """A parent class for meridional advection-diffusion processes.
+    """
+    def __init__(self,
+                 K=0.,
+                 U=0.,
+                 use_banded_solver=False,
+                 prescribed_flux=0.,
+                 **kwargs):
+        super(MeridionalAdvectionDiffusion, self).__init__(K=K, U=U,
+                diffusion_axis='lat', use_banded_solver=use_banded_solver, **kwargs)
+        # Conversion of delta from degrees (grid units) to physical length units
+        phi_stag = np.deg2rad(self.lat_bounds)
+        phi = np.deg2rad(self.lat)
+        self._Xcenter[...,:] = phi*const.a
+        self._Xbounds[...,:] = phi_stag*const.a
+        self._weight_bounds[...,:] = np.cos(phi_stag)
+        self._weight_center[...,:] = np.cos(phi)
+        #  Now properly compute the weighted diffusion matrix
+        self.prescribed_flux = prescribed_flux
+        self.K = K
+        self.U = U
+
+
+class MeridionalDiffusion(MeridionalAdvectionDiffusion):
+    """A parent class for meridional diffusion-only processes,
+    with advection set to zero.
+
+    Otherwise identical to the parent class.
+    """
+    def __init__(self,
+                 K=0.,
+                 use_banded_solver=False,
+                 prescribed_flux=0.,
+                 **kwargs):
+        # Just initialize the AdvectionDiffusion class with U=0
+        super(MeridionalDiffusion, self).__init__(
+            U=0.,
+            K=K,
+            prescribed_flux=prescribed_flux,
+            use_banded_solver=use_banded_solver, **kwargs)

climlab/dynamics/meridional_heat_diffusion.py

@@ -37,7 +37,7 @@
 """
 from __future__ import division
 import numpy as np
-from .meridional_diffusion import MeridionalDiffusion
+from .meridional_advection_diffusion import MeridionalDiffusion
 from climlab import constants as const
 
 
@@ -46,22 +46,22 @@ class MeridionalHeatDiffusion(MeridionalDiffusion):
 
     Solves the meridional heat diffusion equation
 
-    $$ C \frac{\partial T}{\partial t} = -\frac{1}{\cos\phi} \frac{\partial}{\partial \phi} \left[ -D \cos\phi \frac{\partial T}{\partial \phi} \right]$$
+    .. math::
 
-    on an evenly-spaced latitude grid, with a state variable $T$, a heat capacity $C$ and diffusivity $D$.
+        C \frac{\partial T}{\partial t} = -\frac{1}{\cos\phi} \frac{\partial}{\partial \phi} \left[ -D \cos\phi \frac{\partial T}{\partial \phi} \right]
 
-    Assuming $T$ is a temperature in $K$ or $^\circ$C, then the units are:
+    on an evenly-spaced latitude grid, with a state variable :math:`T`,
+    a heat capacity :math:`C` and diffusivity :math:`D`.
 
-    - $D$ in W m$^{-2}$ K$^{-1}$
-    - $C$ in J m$^{-2}$ K$^{-1}$
+    Assuming :math:`T` is a temperature in K or degC, then the units are:
 
-    If the state variable has other units, then $D$ and $C$ should be expressed
-    per state variabe unit.
+    - :math:`D` in W m-2 K-1
+    - :math:`C` in J m-2 K-1
 
-    $D$ is provided as input, and can be either scalar
-    or vector defined at latitude boundaries (length).
+    :math:`D` is provided as input, and can be either scalar
+    or vector defined at latitude boundaries.
 
-    $C$ is normally handled automatically for temperature state variables in CLIMLAB.
+    :math:`C` is normally handled automatically for temperature state variables in CLIMLAB.
     '''
     def __init__(self,
                  D=0.555,  # in W / m^2 / degC

conda-recipe/meta.yaml

@@ -1,4 +1,4 @@
-{% set version = "0.7.4.dev1" %}
+{% set version = "0.7.4.dev2" %}
 
 package:
   name: climlab

docs/source/api/climlab.dynamics.MeridionalAdvectionDiffusion.rst

@@ -0,0 +1,12 @@
+MeridionalAdvectionDiffusion
+---------
+
+.. inheritance-diagram:: climlab.dynamics.MeridionalAdvectionDiffusion
+   :parts: 1
+   :private-bases:
+
+.. automodule:: climlab.dynamics.meridional_advection_diffusion
+    :members:
+    :private-members:
+    :undoc-members:
+    :show-inheritance:

docs/source/api/climlab.dynamics.rst

@@ -11,8 +11,7 @@ climlab.dynamics
 
     climlab.dynamics.BudykoTransport
     climlab.dynamics.AdvectionDiffusion
-    climlab.dynamics.Diffusion
-    climlab.dynamics.MeridionalDiffusion
+    climlab.dynamics.MeridionalAdvectionDiffusion
     climlab.dynamics.MeridionalHeatDiffusion
     climlab.dynamics.MeridionalMoistDiffusion
     climlab.dynamics.adv_diff_numerics

setup.py

@@ -1,7 +1,7 @@
 import os, sys
 import textwrap
 
-VERSION = '0.7.4.dev1'
+VERSION = '0.7.4.dev2'
 
 # BEFORE importing setuptools, remove MANIFEST. Otherwise it may not be
 # properly updated when the contents of directories change (true for distutils,