Changeset 791 for codes/icosagcm
- Timestamp:
- 12/05/18 11:59:29 (5 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 2 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
r790 r791 26 26 from dynamico import precision as prec 27 27 from dynamico.meshes import Cartesian_mesh as Mesh 28 from dynamico.kernels import grad, curl, div, KE 29 from dynamico.LAM import Davies 28 30 29 31 import math as math … … 140 142 gas = thermo.set_vs(v,s) 141 143 return gas, w, z 142 143 def grad(mesh, s_ik):144 llm=s_ik.shape[1]145 edge_num, left, right = mesh.edge_num, mesh.left, mesh.right146 grad_ek = prec.zeros((edge_num, llm))147 for e in range(edge_num):148 left_e = left[e]-1 # Fortran => Python indexing149 right_e = right[e]-1 # Fortran => Python indexing150 grad_ek[e,:] = s_ik[right_e,:]-s_ik[left_e,:]151 return grad_ek152 153 def div(mesh, u_ek):154 llm=u_ek.shape[1]155 primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge156 primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai157 div_ik = prec.zeros((primal_num, llm))158 for i in range(primal_num):159 div_i = 0.160 for iedge in range(primal_deg[i]):161 e = primal_edge[i,iedge]-1 # Fortran => Python indexing162 div_i = div_i + u_ek[e,:]*primal_ne[i,iedge]*le_de[e]163 div_ik[i,:] = div_i / Ai[i]164 return div_ik165 166 def KE(mesh, u_ek):167 llm=u_ek.shape[1]168 primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge169 primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai170 div_ik = prec.zeros((primal_num, llm))171 for i in range(primal_num):172 div_i = 0.173 for iedge in range(primal_deg[i]):174 e = primal_edge[i,iedge]-1 # Fortran => Python indexing175 div_i = div_i + le_de[e]*(u_ek[e,:]**2)176 div_ik[i,:] = div_i * (.25/ Ai[i])177 return div_ik178 179 def curl(mesh, u_ek):180 llm=u_ek.shape[1]181 dual_num, dual_deg, dual_edge, dual_ne, Av = mesh.dual_num, mesh.dual_deg, mesh.dual_edge, mesh.dual_ne, mesh.Av182 curl_vk = prec.zeros((dual_num, llm))183 for v in range(dual_num):184 curl_v = 0.185 for iedge in range(dual_deg[v]):186 e = dual_edge[v,iedge]-1 # Fortran => Python indexing187 curl_v = curl_v + u_ek[e,:]*dual_ne[v,iedge]188 curl_vk[v,:] = curl_v / Av[v]189 return curl_vk190 191 class Davies:192 def __init__(self,N1,N2,x_i,y_i,x_e,y_e):193 self.N1, self.N2 = N1, N2194 self.beta_i = self.mask(x_i,y_i)195 self.beta_e = self.mask(x_e,y_e)196 self.x_e,self.y_e = x_e,y_e197 def mask0(self,c,c0,delta): # 1D building block for Davies relaxation198 N1, N2 = self.N1, self.N2199 m = np.zeros(c.size)200 c1,c2 = c0-N1*delta, c0-(N1+N2)*delta201 for i in range(c.size):202 ci=c[i]203 m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0204 if ci < c2 : m[i]=1.205 if ci > c1 : m[i]=0.206 return m207 def relax(self, llm, step, flow):208 beta_i, beta_e = self.beta_i, self.beta_e209 m,S,u,Phi,W=flow210 for l in range(llm):211 step.mass[:,l] = beta_i*step.mass[:,l] + (1.-beta_i)*m[:,l]212 step.theta_rhodz[:,l] = beta_i*step.theta_rhodz[:,l] + (1.-beta_i)*S[:,l]213 step.u[:,l] = beta_e*step.u[:,l] + (1.-beta_e)*u[:,l]214 for l in range(llm+1):215 step.geopot[:,l] = beta_i*step.geopot[:,l] + (1.-beta_i)*Phi[:,l]216 step.W[:,l] = beta_i*step.W[:,l] + (1.-beta_i)*W[:,l]217 144 218 145 class myDavies(Davies):
Note: See TracChangeset
for help on using the changeset viewer.