Source code for inprod_analytic

"""
    Inner products module
    ======================

    Inner products between the truncated set of basis functions for the \

    ocean and atmosphere streamfunction fields.


    .. 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 from \
	Cehelsky, P., & Tung, K. K. : Theories of multiple equilibria and \
	weather regimes-A critical reexamination. Part II: Baroclinic two-layer \
	models. Journal of the atmospheric sciences, 44(21), 3282-3303, 1987.

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

    :Example:

    >>> from inprod_analytic import init_inprod
    >>> init_inprod()

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

    >>> awavenum = np.empty(natm, dtype=object)
    >>> owavenum = np.empty(noc, dtype=object)
    >>> atmos = atm_tensors(natm)
    >>> ocean = ocean_tensors(noc)

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

    it uses the modules :

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

    >>> from params_maooam import nbatm
    >>> from params_maooam import nboc
    >>> from params_maooam import natm
    >>> from params_maooam import noc
    >>> from params_maooam import n
    >>> from params_maooam import oms
    >>> from params_maooam import ams
    >>> from params_maooam import pi

    Classes
    -------

    * atm_wavenum(typ,P,N,H,Nx,Ny)
    * ocean_wavenum(P,H,Nx,Ny)
    * atm_tensors(natm)
    * ocean_tensors(noc)
"""

import numpy as np
from scipy.sparse import csr_matrix

from params_maooam import nbatm
from params_maooam import nboc
from params_maooam import natm
from params_maooam import noc
from params_maooam import n
from params_maooam import oms
from params_maooam import ams
from params_maooam import pi

sqrtof2 = 1.4142135381698608398437500000


[docs]class atm_wavenum: """ Class to define atmosphere wavenumbers. Attributes : * typ (char) = 'A','K' or 'L'. * M (int) * P (int) * H (int) * Nx (int) * Ny (int) """ def __init__(self, typ, P, M, H, Nx, Ny): """creates the wavenumbers""" self.typ = typ self.P = P self.M = M self.H = H self.Nx = Nx self.Ny = Ny def __repr__(self): return "P = {}, M= {},H={}, Nx= {}, Ny={}".format( self.P, self.M, self.H, self.Nx, self.Ny)
[docs]class ocean_wavenum: """ Class to define ocean wavenumbers Attributes : * P (int) * H (int) * Nx (int) * Ny (int) """ def __init__(self, P, H, Nx, Ny): """creates the wavenumbers""" self.P = P self.H = H self.Nx = Nx self.Ny = Ny def __repr__(self): return "P = {}, H= {}, Nx= {}, Ny={}".format( self.P, self.H, self.Nx, self.Ny)
[docs]class atm_tensors: """ Class which contains all the coefficients a,c,d,s,b,g needed for the tensor computation : Attributes : * :math:`a_{i, j}` * :math:`c_{i, j}` * :math:`d_{i, j}` * :math:`s_{i, j}` * :math:`b_{i, j, k}` * :math:`g_{i, j, k}` Return : * The object will be name `atmos`. """ def __init__(self, natm): self.a = np.zeros((natm, natm), dtype=float) self.c = np.zeros((natm, natm), dtype=float) self.d = np.zeros((natm, noc), dtype=float) self.s = np.zeros((natm, noc), dtype=float) self.b = np.zeros((natm, natm, natm)) self.g = np.zeros((natm, natm, natm)) def __repr__(self): return "matrice d : \n{} \n".format(self.d) # !--------------------------------------------------------! # ! 1. Inner products in the equations for the atmosphere ! # !--------------------------------------------------------!
[docs] def calculate_a(self): r""" .. math:: a_{i, j} = (F_i, {\nabla}^2 F_j) .. note:: Eigenvalues of the Laplacian (atmospheric) """ if natm == 0: exit("*** Problem with calculate_a : natm==0!***") for i in range(0, natm): ti = awavenum[i] self.a[i, i] = - (n**2) * ti.Nx**2 - ti.Ny**2
[docs] def calculate_b(self): r""" .. math:: b_{i, j, k} = (F_i, J(F_j, \nabla^2 F_k)) .. note:: Atmospheric g and a tensors must be computed before \ calling this routine. """ if natm == 0: exit("*** Problem with calculate_b : natm==0!***") for i in range(0, natm): for j in range(0, natm): for k in range(0, natm): val = self.a[k, k]*self.g[i, j, k] self.b[i, j, k] = val
[docs] def calculate_c_atm(self): """ .. math:: c_{i,j} = (F_i, \partial_x F_j) .. note:: Beta term for the atmosphere Strict function !! Only accepts KL type. For any other combination, it will not calculate anything. """ if natm == 0: exit("*** Problem with calculate_c_atm : natm==0!***") for i in range(0, natm): for j in range(0, natm): val = 0. Ti = awavenum[i] Tj = awavenum[j] if ((Ti.typ, Tj.typ) == ('K', 'L')): val = delta(Ti.M - Tj.H) * delta(Ti.P - Tj.P) val = n * Ti.M * val if (val != 0.): self.c[i, j] = val self.c[j, i] = -val
[docs] def calculate_d(self, ocean): r""" .. math:: d_{i,j} = (F_i, \nabla^2 \eta_j) .. note:: Forcing of the ocean on the atmosphere. Atmospheric s tensor and oceanic M tensor must be computed before calling this routine ! """ if natm == 0: exit("*** Problem with calculate_d : natm==0!***") for i in range(0, natm): for j in range(0, noc): self.d[i, j] = self.s[i, j] * ocean.M[j, j]
[docs] def calculate_g(self): """ .. math:: g_{i,j,k} = (F_i, J(F_j, F_k)) .. note:: This is a strict function: it only accepts AKL, KKL and LLL types. For any other combination, it will not calculate anything. """ if natm == 0: exit("*** Problem with calculate_h : natm==0!***") for i in range(0, natm): for j in range(0, natm): for k in range(0, natm): Ti = awavenum[i] Tj = awavenum[j] Tk = awavenum[k] val = 0. if (Ti.typ, Tj.typ, Tk.typ) == ('A', 'K', 'L'): vb1 = B1(Ti.P, Tj.P, Tk.P) vb2 = B2(Ti.P, Tj.P, Tk.P) val = -2 * (sqrtof2 / pi) * Tj.M * delta(Tj.M - Tk.H) \ * flambda(Ti.P + Tj.P + Tk.P) if (val != 0): val = val * (((vb1**2) / (vb1**2 - 1)) - ((vb2**2) / (vb2**2 - 1))) if ((Ti.typ, Tj.typ, Tk.typ) == ('K', 'K', 'L')): vs1 = S1(Tj.P, Tk.P, Tj.M, Tk.H) vs2 = S2(Tj.P, Tk.P, Tj.M, Tk.H) val = vs1 * (delta(Ti.M - Tk.H - Tj.M) \ * delta(Ti.P - Tk.P + Tj.P) \ - delta(Ti.M - Tk.H - Tj.M) \ * delta(Ti.P + Tk.P - Tj.P) \ + (delta(Tk.H - Tj.M + Ti.M) \ + delta(Tk.H - Tj.M - Ti.M)) \ * delta(Tk.P + Tj.P - Ti.P)) \ + vs2 * (delta(Ti.M - Tk.H - Tj.M) \ * delta(Ti.P - Tk.P - Tj.P) \ + (delta(Tk.H - Tj.M - Ti.M) \ + delta(Ti.M + Tk.H - Tj.M)) \ * (delta(Ti.P - Tk.P + Tj.P) \ - delta(Tk.P - Tj.P + Ti.P))) val = val*n if (val != 0.): self.g[i, j, k] = val self.g[j, k, i] = val self.g[k, i, j] = val self.g[i, k, j] = -val self.g[j, i, k] = -val self.g[k, j, i] = -val for i in range(0, natm): for j in range(i+1, natm): for k in range(j+1, natm): Ti = awavenum[i] Ti = awavenum[i] Tj = awavenum[j] Tk = awavenum[k] val = 0. if ((Ti.typ, Tj.typ, Tk.typ) == ('L', 'L', 'L')): vs3 = S3(Tj.P, Tk.P, Tj.H, Tk.H) vs4 = S4(Tj.P, Tk.P, Tj.H, Tk.H) val = vs3 * ((delta(Tk.H - Tj.H - Ti.H) \ - delta(Tk.H - Tj.H + Ti.H)) \ * delta(Tk.P + Tj.P - Ti.P) \ + delta(Tk.H + Tj.H - Ti.H) \ * (delta(Tk.P - Tj.P + Ti.P) \ - delta(Tk.P - Tj.P - Ti.P))) + vs4 \ * ((delta(Tk.H + Tj.H - Ti.H) \ * delta(Tk.P - Tj.P - Ti.P)) \ + (delta(Tk.H - Tj.H + Ti.H)\ - delta(Tk.H - Tj.H - Ti.H)) \ * (delta(Tk.P - Tj.P - Ti.P) \ - delta(Tk.P - Tj.P + Ti.P))) val = val*n if (val != 0.): atmos.g[i, j, k] = val atmos.g[j, k, i] = val atmos.g[k, i, j] = val atmos.g[i, k, j] = -val atmos.g[j, i, k] = -val atmos.g[k, j, i] = -val
[docs] def calculate_s(self): """ .. math:: s_{i,j} = (F_i, \eta_j) .. note:: Forcing (thermal) of the ocean on the atmosphere. """ if natm == 0: exit("*** Problem with calculate_s : natm==0!***") if noc == 0: exit("*** Problem with calculate_s : noc==0!***") for i in range(0, natm): for j in range(0, noc): Ti = awavenum[i] Dj = owavenum[j] val = 0. if (Ti.typ == 'A'): val = flambda(Dj.H) * flambda(Dj.P + Ti.P) if (val != 0.): val = val * 8 * sqrtof2 * Dj.P / \ (pi**2 * (Dj.P**2 - Ti.P**2) * Dj.H) if (Ti.typ == 'K'): val = flambda(2 * Ti.M + Dj.H) * delta(Dj.P - Ti.P) if (val != 0): val = val*4*Dj.H/(pi * (-4 * Ti.M**2 + Dj.H**2)) if (Ti.typ == 'L'): val = delta(Dj.P - Ti.P) * delta(2 * Ti.H - Dj.H) if (val != 0.): self.s[i, j] = val
[docs]class ocean_tensors: """ Class which contains all the coefficients k,m,n,w,o,c \ needed for the tensor computation : Attributes : * :math:`K_{i,j}` * :math:`M_{i,j}` * :math:`N_{i,j}` * :math:`W_{i,j}` * :math:`O_{i,j,k}` * :math:`C_{i,j,k}` Return : * The object will be name ocean """ def __init__(self, noc): self.K = np.zeros((noc, natm), dtype=float) self.M = np.zeros((noc, noc), dtype=float) self.N = np.zeros((noc, noc), dtype=float) self.W = np.zeros((noc, natm), dtype=float) self.O = np.zeros((noc, noc, noc), dtype=float) self.C = np.zeros((noc, noc, noc), dtype=float) def __repr__(self): return "matrice C : \n{}".format(self.C) # !--------------------------------------------------------! # ! 2. Inner products in the equations for the ocean ! # !--------------------------------------------------------!
[docs] def calculate_K(self, atmos): r""" Forcing of the atmosphere on the ocean. .. math:: K_{i,j} = (\eta_i, \nabla^2 F_j) .. note:: atmospheric a and s tensors must be computed before calling this function ! """ if noc == 0: exit("*** Problem with calculate_K : noc==0!***") for i in range(0, noc): for j in range(0, natm): self.K[i, j] = atmos.s[j, i] * atmos.a[j, j]
[docs] def calculate_M(self): r""" Forcing of the ocean fields on the ocean. .. math:: M_{i,j} = (\eta_i, \nabla^2 \eta_j) """ if noc == 0: exit("*** Problem with calculate_M : noc==0!***") for i in range(0, noc): Di = owavenum[i] self.M[i, i] = - (n**2) * Di.Nx**2 - Di.Ny**2
[docs] def calculate_N(self): """ Beta term for the ocean .. math:: N_{i,j} = (\eta_i, \partial_x \eta_j) """ if noc == 0: exit("*** Problem with calculate_N : noc==0!***") val = 0. for i in range(0, noc): for j in range(0, noc): Di = owavenum[i] Dj = owavenum[j] val = delta(Di.P - Dj.P) * flambda(Di.H + Dj.H) if (val != 0.): self.N[i, j] = val * (-2) * Dj.H * Di.H * n / \ ((Dj.H**2 - Di.H**2) * pi)
[docs] def calculate_O(self): """ Temperature advection term (passive scalar) .. math:: O_{i,j,k} = (\eta_i, J(\eta_j, \eta_k)) """ if noc == 0: exit("*** Problem with calculate_O : noc==0!***") val = 0. for i in range(0, noc): for j in range(i, noc): for k in range(i, noc): Di = owavenum[i] Dj = owavenum[j] Dk = owavenum[k] vs3 = S3(Dj.P, Dk.P, Dj.H, Dk.H) vs4 = S4(Dj.P, Dk.P, Dj.H, Dk.H) val = vs3*((delta(Dk.H - Dj.H - Di.H) \ - delta(Dk.H - Dj.H + Di.H)) \ * delta(Dk.P + Dj.P - Di.P) \ + delta(Dk.H + Dj.H - Di.H) \ * (delta(Dk.P - Dj.P + Di.P) \ - delta(Dk.P - Dj.P - Di.P))) \ + vs4 * ((delta(Dk.H + Dj.H - Di.H) \ * delta(Dk.P - Dj.P - Di.P)) \ + (delta(Dk.H - Dj.H + Di.H) \ - delta(Dk.H - Dj.H - Di.H)) \ * (delta(Dk.P - Dj.P - Di.P) \ - delta(Dk.P - Dj.P + Di.P))) \ val = val * n / 2 if (val != 0.): self.O[i, j, k] = val self.O[j, k, i] = val self.O[k, i, j] = val self.O[i, k, j] = -val self.O[j, i, k] = -val self.O[k, j, i] = -val
[docs] def calculate_C_oc(self): r""" .. math:: C_{i,j,k} = (\eta_i, J(\eta_j,\nabla^2 \eta_k)) .. note :: Requires :math:`O_{i,j,k}` \ and :math:`M_{i,j}` to be calculated beforehand. """ if noc == 0: exit("*** Problem with calculate_C : noc==0!***") val = 0. for i in range(0, noc): for j in range(0, noc): for k in range(0, noc): val = self.M[k, k] * self.O[i, j, k] if (val != 0.): self.C[i, j, k] = val
[docs] def calculate_W(self, atmos): """ Short-wave radiative forcing of the ocean. .. math:: W_{i,j} = (\eta_i, F_j) .. note :: atmospheric s tensor must be computed before calling this function ! """ if noc == 0: exit("*** Problem with calculate_W : noc==0!***") if natm == 0: exit("*** Problem with calculate_W : natm==0!***") for i in range(0, noc): for j in range(0, natm): self.W[i, j] = atmos.s[j, i]
# !-----------------------------------------------------! # ! ! # ! Definition of the Helper functions from Cehelsky ! # ! \ Tung ! # ! ! # !-----------------------------------------------------! def B1(Pi, Pj, Pk): return (Pk+Pj)/float(Pi) def B2(Pi, Pj, Pk): return (Pk-Pj)/float(Pi) def delta(r): if (r == 0): return 1. else: return 0. def flambda(r): if (r % 2 == 0): return 0. else: return 1. def S1(Pj, Pk, Mj, Hk): return -((Pk * Mj + Pj * Hk)) / 2. def S2(Pj, Pk, Mj, Hk): return (Pk * Mj - Pj * Hk) / 2. def S3(Pj, Pk, Hj, Hk): return (Pk * Hj + Pj * Hk) / 2. def S4(Pj, Pk, Hj, Hk): return (Pk * Hj - Pj * Hk) / 2. # initialization of the variables awavenum = np.empty(natm, dtype=object) owavenum = np.empty(noc, dtype=object) atmos = atm_tensors(natm) ocean = ocean_tensors(noc)
[docs]def init_inprod(): """creates and computes the inner products.""" j = -1 # awavenum definition # here is the constructor : def __init__(self,typ,P,M,H,Nx,Ny): for i in range(0, nbatm): if ams[i, 0] == 1: awavenum[j+1] = atm_wavenum( 'A', ams[i, 1], 0, 0, 0, ams[i, 1]) awavenum[j+2] = atm_wavenum( 'K', ams[i, 1], ams[i, 0], 0, ams[i, 0], ams[i, 1]) awavenum[j+3] = atm_wavenum( 'L', ams[i, 1], 0, ams[i, 0], ams[i, 0], ams[i, 1]) j = j+3 else: awavenum[j+1] = atm_wavenum( 'K', ams[i, 1], ams[i, 0], 0, ams[i, 0], ams[i, 1]) awavenum[j+2] = atm_wavenum( 'L', ams[i, 1], 0, ams[i, 0], ams[i, 0], ams[i, 1]) j = j+2 # owavenum definition # here is the constructor : def __init__(self,P,H,Nx,Ny): for i in range(0, nboc): owavenum[i] = ocean_wavenum( oms[i, 1], oms[i, 0], oms[i, 0]/2., oms[i, 1]) # Computation of the atmospheric inner products tensors atmos.calculate_a() atmos.calculate_g() atmos.calculate_s() atmos.calculate_b() atmos.calculate_c_atm() # Computation of the oceanic inner products tensors ocean.calculate_M() ocean.calculate_N() ocean.calculate_O() ocean.calculate_C_oc() ocean.calculate_W(atmos) ocean.calculate_K(atmos) # A last atmospheric one that needs ocean.M atmos.calculate_d(ocean)