Changeset 791

Timestamp:

12/05/18 11:59:29 (5 years ago)

Author:

dubos

Message:

devel/Python : moved kernels (grad ...) and Davies relaxation as modules

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/LAM.py (added)
• dynamico/kernels.py (added)
• test/py/Baroclinic_3D_ullrich.py (modified) (2 diffs)

codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py

@@ -26 +26 @@
 from dynamico import precision as prec
 from dynamico.meshes import Cartesian_mesh as Mesh
+from dynamico.kernels import grad, curl, div, KE
+from dynamico.LAM import Davies
 
 import math as math
@@ -140 +142 @@
     gas = thermo.set_vs(v,s)
     return gas, w, z
-
-def grad(mesh, s_ik):
-    llm=s_ik.shape[1]
-    edge_num, left, right = mesh.edge_num, mesh.left, mesh.right
-    grad_ek = prec.zeros((edge_num, llm))
-    for e in range(edge_num):
-        left_e = left[e]-1 # Fortran => Python indexing
-        right_e = right[e]-1 # Fortran => Python indexing
-        grad_ek[e,:] = s_ik[right_e,:]-s_ik[left_e,:]
-    return grad_ek
-
-def div(mesh, u_ek):
-    llm=u_ek.shape[1]
-    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
-    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
-    div_ik = prec.zeros((primal_num, llm))
-    for i in range(primal_num):
-        div_i = 0.
-        for iedge in range(primal_deg[i]):
-            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
-            div_i = div_i + u_ek[e,:]*primal_ne[i,iedge]*le_de[e]
-        div_ik[i,:] = div_i / Ai[i]
-    return div_ik
-
-def KE(mesh, u_ek):
-    llm=u_ek.shape[1]
-    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
-    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
-    div_ik = prec.zeros((primal_num, llm))
-    for i in range(primal_num):
-        div_i = 0.
-        for iedge in range(primal_deg[i]):
-            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
-            div_i = div_i + le_de[e]*(u_ek[e,:]**2)
-        div_ik[i,:] = div_i * (.25/ Ai[i])
-    return div_ik
-
-def curl(mesh, u_ek):
-    llm=u_ek.shape[1]
-    dual_num, dual_deg, dual_edge, dual_ne, Av = mesh.dual_num, mesh.dual_deg, mesh.dual_edge, mesh.dual_ne, mesh.Av
-    curl_vk = prec.zeros((dual_num, llm))
-    for v in range(dual_num):
-        curl_v = 0.
-        for iedge in range(dual_deg[v]):
-            e = dual_edge[v,iedge]-1 # Fortran => Python indexing
-            curl_v = curl_v + u_ek[e,:]*dual_ne[v,iedge]
-        curl_vk[v,:] = curl_v / Av[v]
-    return curl_vk
-
-class Davies:
-    def __init__(self,N1,N2,x_i,y_i,x_e,y_e):
-        self.N1, self.N2 = N1, N2
-        self.beta_i = self.mask(x_i,y_i)
-        self.beta_e = self.mask(x_e,y_e)
-        self.x_e,self.y_e = x_e,y_e
-    def mask0(self,c,c0,delta): # 1D building block for Davies relaxation
-        N1, N2 = self.N1, self.N2
-        m = np.zeros(c.size)
-        c1,c2 = c0-N1*delta, c0-(N1+N2)*delta
-        for i in range(c.size):
-            ci=c[i]
-            m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0
-            if ci < c2 : m[i]=1.
-            if ci > c1 : m[i]=0.
-        return m
-    def relax(self, llm, step, flow):
-        beta_i, beta_e = self.beta_i, self.beta_e
-        m,S,u,Phi,W=flow
-        for l in range(llm):
-            step.mass[:,l]        = beta_i*step.mass[:,l]        + (1.-beta_i)*m[:,l]
-            step.theta_rhodz[:,l] = beta_i*step.theta_rhodz[:,l] + (1.-beta_i)*S[:,l]
-            step.u[:,l]           = beta_e*step.u[:,l]           + (1.-beta_e)*u[:,l]
-        for l in range(llm+1):
-            step.geopot[:,l]      = beta_i*step.geopot[:,l]      + (1.-beta_i)*Phi[:,l]
-            step.W[:,l]           = beta_i*step.W[:,l]           + (1.-beta_i)*W[:,l]
 
 class myDavies(Davies):