License

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 module

In [1]:
import numpy as np
import numba

calculate a temporary velocity

In [2]:
@numba.jit(nopython=True, parallel=False)
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

calculate a pressure correction value

In [3]:
@numba.jit(nopython=True, parallel=False)
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

a convergent calculation of SIMPLE method

In [4]:
@numba.jit(nopython=True, parallel=False)
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

main

In [6]:
%%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: 6.34 s

export

In [7]:
import xarray as xr
In [8]:
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()