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