MIT License
Copyright (c) 2020 computational-sediment-hyd
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
import numpy as np
def calupre(p,u,v,w,Vdir,Vbound,alphaU,dx,dy,dz,nu,axis):
# x: axis = 0, y: axis = 1, z: axis = 2
zero = float(0.0)
V0 = np.copy(Vdir)
nx, ny, nz = p.shape
ac = np.zeros_like(Vdir)
axp = np.zeros_like(Vdir)
axm = np.zeros_like(Vdir)
ayp = np.zeros_like(Vdir)
aym = np.zeros_like(Vdir)
azp = np.zeros_like(Vdir)
azm = np.zeros_like(Vdir)
b = np.zeros_like(Vdir)
act = np.zeros_like(Vdir)
if axis == 0:
xs, ys, zs = 1, 0, 0
xmb, xpb, ymb, ypb, zmb, zpb = -1, -1, 0, ny-1, 0, nz-1
elif axis == 1:
xs, ys, zs = 0, 1, 0
xmb, xpb, ymb, ypb, zmb, zpb = 0, nx-1, -1, -1, 0, nz-1
elif axis == 2:
xs, ys, zs = 0, 0, 1
xmb, xpb, ymb, ypb, zmb, zpb = 0, nx-1, 0, ny-1, -1, -1
for i in range(xs, nx):
for j in range(ys, ny):
for k in range(zs, nz):
c = (i,j,k)
if axis == 0:
uhxp = 0.5*u[i,j,k] + 0.5*u[i+1,j,k]
uhxm = 0.5*u[i,j,k] + 0.5*u[i-1,j,k]
vhyp = 0.5*v[i,j+1,k] + 0.5*v[i-1,j+1,k]
vhym = 0.5*v[i,j ,k] + 0.5*v[i-1,j ,k]
whzp = 0.5*w[i,j,k+1] + 0.5*w[i-1,j,k+1]
whzm = 0.5*w[i,j ,k] + 0.5*w[i-1,j,k]
elif axis == 1:
uhxp = 0.5*u[i+1,j,k] + 0.5*u[i+1,j-1,k]
uhxm = 0.5*u[i,j,k] + 0.5*u[i,j-1,k]
vhyp = 0.5*v[i,j,k] + 0.5*v[i,j+1,k]
vhym = 0.5*v[i,j,k] + 0.5*v[i,j-1,k]
whzp = 0.5*w[i,j,k+1] + 0.5*w[i,j-1,k+1]
whzm = 0.5*w[i,j ,k] + 0.5*w[i,j-1,k]
elif axis == 2:
uhxp = 0.5*u[i+1,j,k] + 0.5*u[i+1,j,k-1]
uhxm = 0.5*u[i,j,k] + 0.5*u[i,j,k-1]
vhyp = 0.5*v[i,j+1,k] + 0.5*v[i,j+1,k-1]
vhym = 0.5*v[i,j,k] + 0.5*v[i,j,k-1]
whzp = 0.5*w[i,j,k] + 0.5*w[i,j,k+1]
whzm = 0.5*w[i,j,k] + 0.5*w[i,j,k-1]
axp[c] = (max([-uhxp,0.0]) + nu/dx)*dy*dz
axm[c] = (max([ uhxm,0.0]) + nu/dx)*dy*dz
ayp[c] = (max([-vhyp,0.0]) + nu/dy)*dx*dz
aym[c] = (max([ vhym,0.0]) + nu/dy)*dx*dz
azp[c] = (max([-whzp,0.0]) + nu/dz)*dx*dy
azm[c] = (max([ whzm,0.0]) + nu/dz)*dx*dy
ac[c] = axp[c]+axm[c]+ayp[c]+aym[c]+azp[c]+azm[c]
if axis == 0:
b[c] = -dy*dz*( p[i,j,k] - p[i-1,j,k] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c]
elif axis == 1:
b[c] = -dx*dz*( p[i,j,k] - p[i,j-1,k] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c]
elif axis == 2:
b[c] = -dx*dy*( p[i,j,k] - p[i,j,k-1] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c]
ac[c] /= alphaU
act[c] = ac[c]
if i == xmb :
ac[c] += axm[c]
elif i == xpb :
ac[c] += axp[c]
if j == ymb :
ac[c] += aym[c]
elif j == ypb :
ac[c] += ayp[c]
if k == zmb :
ac[c] += azm[c]
elif k == zpb :
ac[c] += azp[c]
b[c] += 2.0*azp[c]*Vbound
def sor(u):
u0 = np.copy(u) # deep copy
for i in range(xs, nx):
for j in range(ys, ny):
for k in range(zs, nz):
c = (i,j,k)
xp,xm,yp,ym,zp,zm = (i+1,j,k),(i-1,j,k),(i,j+1,k),(i,j-1,k),(i,j,k+1),(i,j,k-1)
if i == xmb :
uxp ,uxm = u0[xp], zero
elif i == xpb :
uxp ,uxm = zero, u0[xm]
else:
uxp ,uxm = u0[xp], u0[xm]
if j == ymb :
uyp ,uym = u0[yp], zero
elif j == ypb :
uyp ,uym = zero , u0[ym]
else:
uyp ,uym = u0[yp], u0[ym]
if k == zmb :
uzp, uzm = u0[zp], zero
elif k == zpb :
uzp, uzm = zero, u0[zm]
else:
uzp ,uzm = u0[zp], u0[zm]
ut = (axp[c]*uxp + axm[c]*uxm \
+ ayp[c]*uyp + aym[c]*uym \
+ azp[c]*uzp + azm[c]*uzm + b[c])/ac[c]
u0[c] = ut
return u0
unew = np.copy(Vdir) # deep copy
for _ in range(100):
utmp = sor(unew)
err = np.abs(unew - utmp).max()
unew = np.copy(utmp)
if err < 10**(-3) : break
# print(err)
return unew, act
def caldp(p,u,v,w,acu,acv,acw,dx,dy,dz):
zero = float(0.0)
dp = np.zeros_like(p)
ac = np.zeros_like(p)
axp = np.zeros_like(p)
axm = np.zeros_like(p)
ayp = np.zeros_like(p)
aym = np.zeros_like(p)
azp = np.zeros_like(p)
azm = np.zeros_like(p)
b = np.zeros_like(p)
nx, ny, nz = p.shape
for i in range(nx):
for j in range(ny):
for k in range(nz):
c = (i, j, k)
axp[c] = zero if np.abs( acu[i+1,j, k ] ) < 10**(-8) else dy**2*dz**2/acu[i+1,j, k ]
axm[c] = zero if np.abs( acu[i ,j, k ] ) < 10**(-8) else dy**2*dz**2/acu[i ,j, k ]
ayp[c] = zero if np.abs( acv[i ,j+1,k ] ) < 10**(-8) else dx**2*dz**2/acv[i ,j+1,k ]
aym[c] = zero if np.abs( acv[i ,j, k ] ) < 10**(-8) else dx**2*dz**2/acv[i ,j, k ]
azp[c] = zero if np.abs( acw[i ,j, k+1] ) < 10**(-8) else dx**2*dy**2/acw[i ,j, k+1]
azm[c] = zero if np.abs( acw[i ,j, k ] ) < 10**(-8) else dx**2*dy**2/acw[i ,j, k ]
ac[c] = axp[c]+axm[c]+ayp[c]+aym[c]+azp[c]+azm[c]
b[c] = - (u[i+1,j,k] - u[i,j,k])*dy*dz \
- (v[i,j+1,k] - v[i,j,k])*dx*dz \
- (w[i,j,k+1] - w[i,j,k])*dx*dy
def sor(dp0, omega):
p0 = np.copy(dp0) # deep copy
for i in range(nx):
for j in range(ny):
for k in range(nz):
c = (i,j,k)
xp,xm,yp,ym,zp,zm = (i+1,j,k),(i-1,j,k),(i,j+1,k),(i,j-1,k),(i,j,k+1),(i,j,k-1)
if i == 0:
pxp ,pxm = p0[xp], zero
elif i == nx-1:
pxp ,pxm = zero, p0[xm]
else:
pxp ,pxm = p0[xp], p0[xm]
if j == 0:
pyp ,pym = p0[yp], zero
elif j == ny-1:
pyp ,pym = zero , p0[ym]
else:
pyp ,pym = p0[yp], p0[ym]
if k == 0:
pzp, pzm = p0[zp], zero
elif k == nz-1:
pzp, pzm = zero, p0[zm]
else:
pzp ,pzm = p0[zp], p0[zm]
pt = (axp[c]*pxp + axm[c]*pxm \
+ ayp[c]*pyp + aym[c]*pym \
+ azp[c]*pzp + azm[c]*pzm + b[c])/ac[c]
p0[c] += omega*(pt - p0[c])
# p0[c] = ut
return p0
dpnew = np.copy(dp) # deep copy
for nn in range(1000):
dptmp = sor(dpnew, float(1.8))
err = np.abs(dpnew - dptmp).max()
dpnew = np.copy(dptmp)
if err < 10**(-6) :
print('SOR N =', nn)
break
return dpnew
def SIMPLE(p,u,v,w,Uslip,Vslip,alphaU,alphaP,dx,dy,dz,nu):
zero = float(0.0)
nx, ny, nz = p.shape
for nn in range(1000):
u, acu = calupre(p,u,v,w,u,Uslip,alphaU,dx,dy,dz,nu,axis=0)
v, acv = calupre(p,u,v,w,v,Vslip,alphaU,dx,dy,dz,nu,axis=1)
w, acw = calupre(p,u,v,w,w,zero,alphaU,dx,dy,dz ,nu,axis=2)
dp = caldp(p,u,v,w,acu,acv,acw,dx,dy,dz)
p += alphaP*dp
unew = np.copy(u) # deep copy
vnew = np.copy(v) # deep copy
wnew = np.copy(w) # deep copy
for i in range(nx):
for j in range(ny):
for k in range(nz):
if i !=0 : unew[i,j,k] = u[i,j,k] + (dp[i-1,j,k] - dp[i,j,k])*dy*dz/acu[i,j,k]
if j !=0 : vnew[i,j,k] = v[i,j,k] + (dp[i,j-1,k] - dp[i,j,k])*dx*dz/acv[i,j,k]
if k !=0 : wnew[i,j,k] = w[i,j,k] + (dp[i,j,k-1] - dp[i,j,k])*dx*dy/acw[i,j,k]
err = max([np.abs(unew - u).max(), np.abs(vnew - v).max(), np.abs(wnew - w).max()])
print('N', nn, 'error =', err)
u = np.copy(unew) # deep copy
v = np.copy(vnew) # deep copy
w = np.copy(wnew) # deep copy
if err < 10**(-6) : break
return p,u,v,w
%%time
dx,dy,dz = 0.05,0.05,0.05
nx,ny,nz = 20,20,20
nu = 0.01
alphaU = 0.5
alphaP = 0.8
Uslip = 1.0
Vslip = 0.0
p = np.zeros((nx,ny,nz))
u = np.zeros((nx+1,ny,nz))
v = np.zeros((nx,ny+1,nz))
w = np.zeros((nx,ny,nz+1))
p,u,v,w = SIMPLE(p,u,v,w,Uslip,Vslip,alphaU,alphaP,dx,dy,dz,nu)
SOR N = 125 N 0 error = 0.1054172791001945 SOR N = 120 N 1 error = 0.026840052366420453 SOR N = 109 N 2 error = 0.009562053929592945 SOR N = 106 N 3 error = 0.005360998254980088 SOR N = 93 N 4 error = 0.004474142497522898 SOR N = 94 N 5 error = 0.0032290649302523233 SOR N = 69 N 6 error = 0.002666631525062191 SOR N = 84 N 7 error = 0.002131719922800368 SOR N = 67 N 8 error = 0.001777679411929134 SOR N = 76 N 9 error = 0.0014729779128486165 SOR N = 60 N 10 error = 0.0012394939779675207 SOR N = 73 N 11 error = 0.0010233730467733015 SOR N = 63 N 12 error = 0.0009551915739345884 SOR N = 66 N 13 error = 0.0008101360144389114 SOR N = 64 N 14 error = 0.0007446578525629466 SOR N = 62 N 15 error = 0.0006630197513536884 SOR N = 60 N 16 error = 0.0006007713345940191 SOR N = 60 N 17 error = 0.0005408998700914824 SOR N = 60 N 18 error = 0.0004853426230619351 SOR N = 54 N 19 error = 0.00045843700257259916 SOR N = 60 N 20 error = 0.0003988421881952475 SOR N = 53 N 21 error = 0.00043021171823730275 SOR N = 54 N 22 error = 0.000336056822965225 SOR N = 54 N 23 error = 0.00034334116655392044 SOR N = 53 N 24 error = 0.00029647191465407077 SOR N = 53 N 25 error = 0.00027364370223947887 SOR N = 52 N 26 error = 0.0002676683707877048 SOR N = 51 N 27 error = 0.0002592662461848305 SOR N = 50 N 28 error = 0.00024920470816780504 SOR N = 50 N 29 error = 0.00023910012218467114 SOR N = 50 N 30 error = 0.00022888050150252082 SOR N = 50 N 31 error = 0.00021876325415018383 SOR N = 56 N 32 error = 0.00018654047333169221 SOR N = 51 N 33 error = 0.00018874049753789257 SOR N = 45 N 34 error = 0.00017905529517292518 SOR N = 48 N 35 error = 0.00015435267801000574 SOR N = 43 N 36 error = 0.00015495249139380052 SOR N = 43 N 37 error = 0.0001500953516865855 SOR N = 42 N 38 error = 0.0001435946207644645 SOR N = 42 N 39 error = 0.00013688311400816833 SOR N = 42 N 40 error = 0.00013036830247009634 SOR N = 41 N 41 error = 0.0001241165939866451 SOR N = 41 N 42 error = 0.00011816156683286394 SOR N = 41 N 43 error = 0.00011247299881905759 SOR N = 41 N 44 error = 0.00010622567236323599 SOR N = 41 N 45 error = 0.00010114366322544477 SOR N = 41 N 46 error = 9.635339645153174e-05 SOR N = 40 N 47 error = 9.179793736047159e-05 SOR N = 40 N 48 error = 8.747896623534368e-05 SOR N = 40 N 49 error = 8.342524406629304e-05 SOR N = 40 N 50 error = 7.964053028247653e-05 SOR N = 40 N 51 error = 7.597562643929523e-05 SOR N = 40 N 52 error = 7.247874027371815e-05 SOR N = 40 N 53 error = 6.918152121720977e-05 SOR N = 40 N 54 error = 6.60482520006278e-05 SOR N = 39 N 55 error = 6.319257672732226e-05 SOR N = 39 N 56 error = 6.473462051934109e-05 SOR N = 39 N 57 error = 6.366427739659675e-05 SOR N = 39 N 58 error = 6.183994440921159e-05 SOR N = 39 N 59 error = 5.989779351682489e-05 SOR N = 39 N 60 error = 5.794540768976064e-05 SOR N = 38 N 61 error = 5.6018841795471563e-05 SOR N = 62 N 62 error = 0.00014012859791724674 SOR N = 70 N 63 error = 0.00030523211189419086 SOR N = 70 N 64 error = 0.00031769023621645853 SOR N = 70 N 65 error = 0.00032058495419517996 SOR N = 70 N 66 error = 0.0003027958441054607 SOR N = 65 N 67 error = 0.00016303298345070327 SOR N = 31 N 68 error = 2.7589976646716363e-05 SOR N = 34 N 69 error = 3.398227355037864e-05 SOR N = 33 N 70 error = 3.3124064515527296e-05 SOR N = 33 N 71 error = 3.224077203133058e-05 SOR N = 32 N 72 error = 3.141954219593179e-05 SOR N = 32 N 73 error = 3.0402919470054468e-05 SOR N = 32 N 74 error = 2.8088649561580636e-05 SOR N = 32 N 75 error = 2.897129306045354e-05 SOR N = 31 N 76 error = 2.858906798969274e-05 SOR N = 31 N 77 error = 2.7884534702316e-05 SOR N = 31 N 78 error = 2.715388787510875e-05 SOR N = 31 N 79 error = 2.6450167137037628e-05 SOR N = 30 N 80 error = 2.5925174391039363e-05 SOR N = 30 N 81 error = 2.599445992423899e-05 SOR N = 30 N 82 error = 2.5347989921209457e-05 SOR N = 30 N 83 error = 2.4645267609368915e-05 SOR N = 29 N 84 error = 2.3961220961987717e-05 SOR N = 29 N 85 error = 2.3306482079238355e-05 SOR N = 29 N 86 error = 2.267666541694302e-05 SOR N = 29 N 87 error = 2.206940890855935e-05 SOR N = 29 N 88 error = 2.1482939978306748e-05 SOR N = 28 N 89 error = 2.0918162684901986e-05 SOR N = 28 N 90 error = 2.036955693593412e-05 SOR N = 28 N 91 error = 1.9839187604175912e-05 SOR N = 28 N 92 error = 1.9324913766799456e-05 SOR N = 28 N 93 error = 1.8826719417108784e-05 SOR N = 28 N 94 error = 1.8343932562586707e-05 SOR N = 28 N 95 error = 1.7875704076159016e-05 SOR N = 27 N 96 error = 1.7421960362223876e-05 SOR N = 27 N 97 error = 1.698314896375619e-05 SOR N = 27 N 98 error = 1.6556941404843872e-05 SOR N = 27 N 99 error = 1.6144575745419276e-05 SOR N = 27 N 100 error = 1.5744326605116044e-05 SOR N = 26 N 101 error = 1.5356390333981507e-05 SOR N = 26 N 102 error = 1.4980510895518107e-05 SOR N = 26 N 103 error = 1.4614357406378398e-05 SOR N = 26 N 104 error = 1.4260768505386379e-05 SOR N = 26 N 105 error = 1.3917367757965149e-05 SOR N = 26 N 106 error = 1.3583933360872269e-05 SOR N = 26 N 107 error = 1.3260126400765904e-05 SOR N = 25 N 108 error = 1.2944433334938221e-05 SOR N = 25 N 109 error = 1.264212657456354e-05 SOR N = 25 N 110 error = 1.2343689665961222e-05 SOR N = 25 N 111 error = 1.2056259599324548e-05 SOR N = 25 N 112 error = 1.1776533205526407e-05 SOR N = 25 N 113 error = 1.1504563897341002e-05 SOR N = 25 N 114 error = 1.1240385126981556e-05 SOR N = 25 N 115 error = 1.0983425486971177e-05 SOR N = 25 N 116 error = 1.0733679783103689e-05 SOR N = 24 N 117 error = 1.0488599582136882e-05 SOR N = 24 N 118 error = 1.025664598860998e-05 SOR N = 24 N 119 error = 1.0024301006189562e-05 SOR N = 24 N 120 error = 9.800993931941004e-06 SOR N = 24 N 121 error = 9.58337486117733e-06 SOR N = 24 N 122 error = 9.371530223512003e-06 SOR N = 24 N 123 error = 9.165800510690936e-06 SOR N = 24 N 124 error = 8.964859223775656e-06 SOR N = 24 N 125 error = 8.76957528181399e-06 SOR N = 24 N 126 error = 8.578994010716157e-06 SOR N = 23 N 127 error = 8.391625714365691e-06 SOR N = 23 N 128 error = 8.214163884662229e-06 SOR N = 23 N 129 error = 8.03603539067943e-06 SOR N = 23 N 130 error = 7.864470786350664e-06 SOR N = 23 N 131 error = 7.696970618281673e-06 SOR N = 23 N 132 error = 7.533437274775956e-06 SOR N = 23 N 133 error = 7.37399148104112e-06 SOR N = 23 N 134 error = 7.218437776845832e-06 SOR N = 22 N 135 error = 7.06741559322599e-06 SOR N = 22 N 136 error = 6.919704353486322e-06 SOR N = 22 N 137 error = 6.773784997748944e-06 SOR N = 22 N 138 error = 6.633038374032063e-06 SOR N = 22 N 139 error = 6.495210704782206e-06 SOR N = 22 N 140 error = 6.3607407126131665e-06 SOR N = 22 N 141 error = 6.229707239652216e-06 SOR N = 22 N 142 error = 6.101586070533793e-06 SOR N = 21 N 143 error = 5.973629219918619e-06 SOR N = 21 N 144 error = 5.853906029584799e-06 SOR N = 21 N 145 error = 5.73379941773422e-06 SOR N = 21 N 146 error = 5.6174751538740075e-06 SOR N = 21 N 147 error = 5.503453911537282e-06 SOR N = 21 N 148 error = 5.390388731818518e-06 SOR N = 20 N 149 error = 5.2864563543086884e-06 SOR N = 20 N 150 error = 5.173697379928788e-06 SOR N = 20 N 151 error = 5.070093015757671e-06 SOR N = 20 N 152 error = 4.969354221368016e-06 SOR N = 20 N 153 error = 4.869661774098422e-06 SOR N = 20 N 154 error = 4.772086256932262e-06 SOR N = 20 N 155 error = 4.676430728212111e-06 SOR N = 19 N 156 error = 4.5836499956708465e-06 SOR N = 20 N 157 error = 4.4904264110967596e-06 SOR N = 19 N 158 error = 4.403071816150295e-06 SOR N = 19 N 159 error = 4.31709385506629e-06 SOR N = 19 N 160 error = 4.22919505516095e-06 SOR N = 19 N 161 error = 4.14623453118268e-06 SOR N = 19 N 162 error = 4.0643825724162586e-06 SOR N = 18 N 163 error = 3.9845962542806035e-06 SOR N = 18 N 164 error = 3.909632786375239e-06 SOR N = 18 N 165 error = 3.827390240906947e-06 SOR N = 18 N 166 error = 3.7520635074728137e-06 SOR N = 17 N 167 error = 3.668413285806693e-06 SOR N = 17 N 168 error = 3.6123403394650033e-06 SOR N = 17 N 169 error = 3.535180677832761e-06 SOR N = 17 N 170 error = 3.465793159046493e-06 SOR N = 16 N 171 error = 3.3992023011630845e-06 SOR N = 16 N 172 error = 3.3376731522483105e-06 SOR N = 16 N 173 error = 3.266856210965008e-06 SOR N = 16 N 174 error = 3.20231841649532e-06 SOR N = 16 N 175 error = 3.1397203980598754e-06 SOR N = 15 N 176 error = 3.0712207199357078e-06 SOR N = 16 N 177 error = 3.027014526146843e-06 SOR N = 15 N 178 error = 2.9520133110194635e-06 SOR N = 15 N 179 error = 2.9027338558917926e-06 SOR N = 15 N 180 error = 2.8461905576537827e-06 SOR N = 15 N 181 error = 2.790611308700619e-06 SOR N = 14 N 182 error = 2.7410191304755305e-06 SOR N = 15 N 183 error = 2.688334522854552e-06 SOR N = 14 N 184 error = 2.6389620582933926e-06 SOR N = 14 N 185 error = 2.577222071448171e-06 SOR N = 14 N 186 error = 2.530575231890486e-06 SOR N = 14 N 187 error = 2.4816696066465305e-06 SOR N = 13 N 188 error = 2.4391153694014456e-06 SOR N = 14 N 189 error = 2.3753895243561196e-06 SOR N = 13 N 190 error = 2.350730326894368e-06 SOR N = 13 N 191 error = 2.2904570768089716e-06 SOR N = 13 N 192 error = 2.250331456510324e-06 SOR N = 12 N 193 error = 2.2228703493720747e-06 SOR N = 13 N 194 error = 2.145835665184892e-06 SOR N = 12 N 195 error = 2.1459965736936315e-06 SOR N = 12 N 196 error = 2.074253453526742e-06 SOR N = 12 N 197 error = 2.0414637246501943e-06 SOR N = 11 N 198 error = 2.014235961300681e-06 SOR N = 12 N 199 error = 1.9514262204478605e-06 SOR N = 10 N 200 error = 1.950941234313275e-06 SOR N = 12 N 201 error = 1.8628834946021744e-06 SOR N = 9 N 202 error = 1.9052719819612207e-06 SOR N = 11 N 203 error = 1.76883247940407e-06 SOR N = 9 N 204 error = 1.8228045058632514e-06 SOR N = 11 N 205 error = 1.707919988103157e-06 SOR N = 8 N 206 error = 1.7580508588632693e-06 SOR N = 12 N 207 error = 1.6294349520218354e-06 SOR N = 7 N 208 error = 1.8017573641115892e-06 SOR N = 13 N 209 error = 1.701994712034749e-06 SOR N = 7 N 210 error = 1.8458912151714246e-06 SOR N = 13 N 211 error = 1.665398729283618e-06 SOR N = 7 N 212 error = 1.8073989873079732e-06 SOR N = 13 N 213 error = 1.643898838599861e-06 SOR N = 7 N 214 error = 1.7307470111666001e-06 SOR N = 12 N 215 error = 1.5541772441993174e-06 SOR N = 6 N 216 error = 1.6524432158090963e-06 SOR N = 13 N 217 error = 1.7998636709370963e-06 SOR N = 6 N 218 error = 1.5734240664116994e-06 SOR N = 12 N 219 error = 1.6644250930217264e-06 SOR N = 6 N 220 error = 1.5392734142929965e-06 SOR N = 12 N 221 error = 1.711829663561537e-06 SOR N = 6 N 222 error = 1.4371711831717704e-06 SOR N = 12 N 223 error = 1.6276829473069188e-06 SOR N = 6 N 224 error = 1.3932405514907598e-06 SOR N = 11 N 225 error = 1.5611043250804424e-06 SOR N = 6 N 226 error = 1.3042823898551381e-06 SOR N = 11 N 227 error = 1.556537582693418e-06 SOR N = 6 N 228 error = 1.2246238252841546e-06 SOR N = 10 N 229 error = 1.4791853819307033e-06 SOR N = 5 N 230 error = 1.0862523535382085e-06 SOR N = 11 N 231 error = 1.6347966111593393e-06 SOR N = 5 N 232 error = 1.0304977162561846e-06 SOR N = 11 N 233 error = 1.4893600061155476e-06 SOR N = 5 N 234 error = 1.0153743705970664e-06 SOR N = 10 N 235 error = 1.49247726410201e-06 SOR N = 5 N 236 error = 9.118841478861217e-07 Wall time: 10min 3s
import xarray as xr
uc = np.zeros((nx,ny,nz))
vc = np.zeros((nx,ny,nz))
wc = np.zeros((nx,ny,nz))
for i in range(nx) : uc[i,:,:] = u[i,:,:] + u[i+1,:,:]
for j in range(ny) : vc[:,j,:] = v[:,j,:] + v[:,j+1,:]
for k in range(nz) : wc[:,:,k] = w[:,:,k] + w[:,:,k+1]
dss = xr.Dataset( { 'p': (['x','y','z'], p), 'u': (['x','y','z'], uc), 'v': (['x','y','z'], vc), 'w': (['x','y','z'], wc) }
, coords={'x' : np.arange(nx)*dx + 0.5*dx , 'y' : np.arange(ny)*dy + 0.5*dy, 'z' : np.arange(nz)*dz + 0.5*dz }
)
dss.to_netcdf('output.nc')
dss.close()