Source code for aotensor

"""
    Tensor computation module
    =========================

    The equation tensor for the coupled ocean-atmosphere model
    with temperature which allows for an extensible set of modes
    in the ocean and in the atmosphere.

    .. note :: These are calculated using the analytical expressions from De \
    Cruz, L., Demaeyer, J. and Vannitsem, S.: A modular arbitrary-order \
    ocean-atmosphere model: MAOOAM v1.0, Geosci. Model Dev. Discuss.\
            And the `Fortran Code <https://github.com/Climdyn/MAOOAM>`_

    .. note :: The python code is available here : \
    `aotensor.py <../_modules/aotensor.html>`_ .

    :Example:

    >>> aotensor, Li, Lj, Lk, Lv =aotensor.init_aotensor()

    Help Functions
    -------------------
    There are ndim coordinates that correspond to 4 physical quantities.
    These functions help to have the i-th coordinates of a quantity.

    * psi(i) -> i
    * theta(i) -> i + natm
    * A(i) -> i + 2*natm
    * T(i) -> i + 2*natm + noc
    * kdelta(i,j) -> (i==j)

    Global variables
    -------------------

    * real_eps = 2.2204460492503131e-16
    * t=np.zeros( ((ndim+1),(ndim+1),(ndim+1)),dtype=float)

    Dependencies
    ------------------------

    >>> from params_maooam import *
    >>> from inprod_analytic import *
    >>> from scipy.sparse import dok_matrix
    >>> from scipy.sparse import csr_matrix
    >>> import os

    Functions
    ------------------------

    * compute_aotensor
    * coeff(i, j, k, v)
    * simplify
    * init_aotensor


"""

import numpy as np
from scipy.sparse import dok_matrix
from scipy.sparse import csr_matrix

from params_maooam import *
from inprod_analytic import init_inprod
from inprod_analytic import atmos
from inprod_analytic import ocean

real_eps = 2.2204460492503131e-16

t = np.zeros(((ndim+1), (ndim+1), (ndim+1)), dtype=float)


def psi(i):
    return i


def theta(i):
    return i + natm


def A(i):
    return i + 2*natm


def T(i):
    return i + 2*natm + noc


def kdelta(i, j):

    if i == j:
        return 1

    else:
        return 0


[docs]def compute_aotensor(): """ Computes the three-dimensional tensor t Takes the inner products of inprod_analytic \ and computes the tensor :param t: tensor t is a global variable of aotensor :type t: array((37,37,37),float) :return: change the global tensor :rtype: void :Example: >>> compute_aotensor() .. warning:: Needs the global variable aotensor and the global inner \ products to be initialized. """ t[theta(1), 0, 0] = (Cpa / (1 - atmos.a[0, 0] * sig0)) for i in range(1, natm+1): for j in range(1, natm+1): t[psi(i), psi(j), 0] = -((atmos.c[(i-1), (j-1)] * betp) / atmos.a[(i-1), (i-1)]) - (kd * kdelta((i-1), (j-1))) / 2 t[theta(i), psi(j), 0] = ((atmos.a[(i-1), (j-1)] * kd * sig0) / (-2 + 2*atmos.a[(i-1), (i-1)]*sig0)) t[psi(i), theta(j), 0] = (kd * kdelta((i-1), (j-1))) / 2 t[theta(i), theta(j), 0] = (-((sig0 * (2. * atmos.c[(i-1), (j-1)] \ * betp + atmos.a[(i-1), (j-1)] * (kd + 4. * kdp))))\ + 2.*(LSBpa + sc*Lpa) * kdelta((i-1), (j-1))) \ / (-2. + 2.*atmos.a[(i-1), (i-1)]*sig0) for k in range(1, natm+1): t[psi(i), psi(j), psi(k)] = -(atmos.b[(i-1), (j-1), (k-1)] \ / atmos.a[(i-1), (i-1)]) t[psi(i), theta(j), theta(k)] = -((atmos.b[(i-1), (j-1), (k-1)] \ / atmos.a[(i-1), (i-1)])) t[theta(i), psi(j), theta(k)] = (atmos.g[(i-1), (j-1), (k-1)] \ - atmos.b[(i-1), (j-1), (k-1)]*sig0) / \ (-1 + atmos.a[(i-1), (i-1)]*sig0) t[theta(i), theta(j), psi(k)] = (atmos.b[(i-1), (j-1), (k-1)] \ * sig0) / (1 - atmos.a[(i-1), (i-1)] * sig0) for j in range(1, noc+1): t[psi(i), A(j), 0] = kd * atmos.d[(i-1), (j-1)] / \ (2 * atmos.a[(i-1), (i-1)]) t[theta(i), A(j), 0] = kd * (atmos.d[(i-1), (j-1)] * sig0) \ / (2 - 2 * atmos.a[(i-1), (i-1)] * sig0) t[theta(i), T(j), 0] = atmos.s[(i-1), (j-1)] * (2 * LSBpo + Lpa) \ / (2 - 2 * atmos.a[(i-1), (i-1)] * sig0) for i in range(1, noc+1): for j in range(1, natm+1): t[A(i), psi(j), 0] = ocean.K[(i-1), (j-1)] * dp \ / (ocean.M[(i-1), (i-1)] + G) t[A(i), theta(j), 0] = -(ocean.K[(i-1), (j-1)]) * dp \ / (ocean.M[(i-1), (i-1)] + G) for j in range(1, noc+1): t[A(i), A(j), 0] = -((ocean.N[(i-1), (j-1)] * betp + \ ocean.M[(i-1), (i-1)] * (rp + dp) * kdelta((i-1), (j-1)))) \ / (ocean.M[(i-1), (i-1)] + G) for k in range(1, noc+1): t[A(i), A(j), A(k)] = -(ocean.C[(i-1), (j-1), (k-1)]) \ / (ocean.M[(i-1), (i-1)] + G) for i in range(1, noc+1): t[T(i), 0, 0] = Cpo * ocean.W[(i-1), 0] for j in range(1, natm+1): t[T(i), theta(j), 0] = ocean.W[(i-1), (j-1)] * (2 * sc * Lpo + sbpa) for j in range(1, noc+1): t[T(i), T(j), 0] = -((Lpo + sbpo)) * kdelta((i-1), (j-1)) for k in range(1, noc+1): t[T(i), A(j), T(k)] = -(ocean.O[(i-1), (j-1), (k-1)])
[docs]def coeff(i, j, k, v): """ Affects v for :math:`t_{i,j,k}` making that tensor[i] upper triangular. Used in compute_aotensor. :param i: first coordinates :type i: int in [1,37] :param j: second coodinates :type j: int in [1,37] :param k: third coordinates :type k: int in [1,37] :param v: value :type v: float :return: change the global tensor :rtype: void :Example: >>> coeff(i, j, k, v) """ if(abs(v) >= real_eps): if j <= k: t[i, j, k] = v else: t[i, k, j] = v
[docs]def simplify(): """ Make sure that tensor[i] is upper triangular. To do after compute_aotensor(). :param t: tensor t is a global variable of aotensor :type t: array((ndim+1,ndim+1,ndim+1),float) :return: change the global tensor :rtype: void :Example: >>> simplify() """ for i in range(1, ndim+1): for j in range(0, ndim+1): for k in range(0, j): if t[i, j, k] != 0: t[i, k, j] = t[i, j, k]+t[i, k, j] t[i, j, k] = 0.
[docs]def init_aotensor(): """ Initialize the tensor. :return: aotensor, Li, Lj, Lk, Lv :type aotensor: (int i, int j, int k, float v) in list :type Li: int in list :type Lj: int in list :type Lk: int in list :type Lv: float in list :Example: >>> aotensor, Li, Lj, Lk, Lv = init_aotensor() """ init_inprod() compute_aotensor() simplify() global tensor_liste global tensor tensor_liste = [] for i in range(0, ndim+1): Xbis = csr_matrix(t[i]) X = Xbis tensor_liste.append(X) tensor = np.array(tensor_liste) global aotensor aotensor = [] for i in range(1, ndim+1): X = tensor[i].nonzero() for m in range(0, tensor[i].nnz): j = X[0][m] k = X[1][m] aotensor.append((i, j, k, tensor[i][j, k])) Li = [] Lj = [] Lk = [] Lv = [] for i in range(0, ndim+1): M = t[i] Lj1 = M.nonzero()[0] Lk1 = M.nonzero()[1] for m in range(0, np.count_nonzero(M)): j = Lj1[m] k = Lk1[m] Li.append(i) Lj.append(j) Lk.append(k) Lv.append(M[j, k]) return aotensor, Li, Lj, Lk, Lv