Pressure Driven Uneven Solver
%precision 3
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as sl
import scipy.sparse as sp
import scipy.sparse.linalg as spla
# the following allows us to plot triangles indicating convergence order
from matplotlib import cm
# and we will create some animations!
import matplotlib.animation as animation
from IPython.display import HTML
from pprint import pprint
def pressure_poisson_jacobi(p, dx, dy, RHS, rtol = 1.e-5, logs = False):
""" Solve the pressure Poisson equation (PPE)
using Jacobi iteration assuming mesh spacing of
dx and dy (we assume at the moment that dx=dy)
and the RHS function given by RHS.
Assumes imposition of a Neumann BC on all boundaries.
Return the pressure field.
"""
# iterate
tol = 10.*rtol
it = 0
p_old = np.copy(p)
imax = len(p)
jmax = np.size(p)//imax
while tol > rtol:
it += 1
#this version is valid for dx!=dy
p[1:-1, 1:-1] = 1.0/(2.0+2.0*(dx**2)/(dy**2))*(p_old[2:, 1:-1] + p_old[:-2, 1:-1] +
(p_old[1:-1, 2:] + p_old[1:-1,:-2])*(dx**2)/(dy**2)
- (dx**2)*RHS[1:-1, 1:-1])
# apply zero gradient Neumann boundary conditions at the no slip walls
p[:, 0] = p[:, 1]
p[:, -1] = p[:, -2]
p[0, :-(jmax-jmax//2)] = p[1, :-(jmax-jmax//2)]
# relative change in pressure
tol = sl.norm(p - p_old)/np.maximum(1.0e-10,sl.norm(p))
#swap arrays without copying the data
temp = p_old
p_old = p
p = temp
if logs: print('pressure solve iterations = {:4d}'.format(it))
return p
def calculate_ppm_RHS_upwind(rho, u, v, RHS, dt, dx, dy):
""" Calculate the RHS of the
Poisson equation resulting from the projection method.
Use upwind differences for the first derivatives of u and v.
"""
RHS[1:-1, 1:-1] =rho*(np.select([u[1:-1, 1:-1] > 0, u[1:-1, 1:-1] <= 0],
[np.diff(u[:-1, 1:-1], n=1, axis=0)/dx,
np.diff(u[1:, 1:-1], n=1, axis=0)/dx]) +
np.select([v[1:-1, 1:-1] > 0, v[1:-1, 1:-1] <= 0],
[np.diff(v[1:-1, :-1], n=1, axis=1)/dy,
np.diff(v[1:-1, 1:], n=1, axis=1)/dy]))
return RHS
def calculate_ppm_RHS_central(rho, u, v, RHS, dt, dx, dy):
""" Calculate the RHS of the
Poisson equation resulting from the projection method.
Use central differences for the first derivatives of u and v.
"""
RHS[1:-1, 1:-1] = rho * (
(1.0/dt) * ( (u[2:, 1:-1] - u[:-2, 1:-1]) / (2.0*dx)
+ (v[1:-1, 2:] - v[1:-1, :-2]) / (2.0*dy) ))
return RHS
def project_velocity(rho, u, v, dt, dx, dy, p):
""" Update the velocity to be divergence free using the pressure.
"""
u[1:-1, 1:-1] = u[1:-1, 1:-1] - dt * (1./rho) * (
(p[2:, 1:-1] - p[:-2, 1:-1])/(2*dx) )
v[1:-1, 1:-1] = v[1:-1, 1:-1] - dt * (1./rho) * (
(p[1:-1, 2:] - p[1:-1, :-2])/(2*dy) )
return u, v
def calculate_intermediate_velocity(nu, u, v, u_old, v_old, dt, dx, dy):
""" Calculate the intermediate velocities.
"""
# intermediate u
u[1:-1, 1:-1] = u_old[1:-1, 1:-1] - dt * (
# ADVECTION: uu_x + vu_x
u_old[1:-1, 1:-1] *
# see comments in the upwind based solver for advection above
np.select([u_old[1:-1, 1:-1] > 0, u_old[1:-1, 1:-1] < 0],
[np.diff(u_old[:-1, 1:-1], n=1, axis=0)/dx,
np.diff(u_old[1:, 1:-1], n=1, axis=0)/dx]) +
v_old[1:-1, 1:-1] *
np.select([v_old[1:-1, 1:-1] > 0, v_old[1:-1, 1:-1] < 0],
[np.diff(u_old[1:-1, :-1], n=1, axis=1)/dy,
np.diff(u_old[1:-1, 1:], n=1, axis=1)/dy]) ) + (
# DIFFUSION
dt*nu*(np.diff(u_old[:, 1:-1], n=2, axis=0)/(dx**2)
+ np.diff(u_old[1:-1, :], n=2, axis=1)/(dy**2)) )
# intermediate v
v[1:-1, 1:-1] = v_old[1:-1, 1:-1] - dt * (
# ADVECTION: uv_x + vv_x
u_old[1:-1, 1:-1] *
np.select([u_old[1:-1, 1:-1] > 0, u_old[1:-1, 1:-1] < 0],
[np.diff(v_old[:-1, 1:-1], n=1, axis=0)/dx,
np.diff(v_old[1:, 1:-1], n=1, axis=0)/dx]) +
v_old[1:-1, 1:-1] *
np.select([v_old[1:-1, 1:-1] > 0, v_old[1:-1, 1:-1] < 0],
[np.diff(v_old[1:-1, :-1], n=1, axis=1)/dy,
np.diff(v_old[1:-1, 1:], n=1, axis=1)/dy]) ) + (
# DIFFUSION
dt*nu*(np.diff(v_old[:, 1:-1], n=2, axis=0)/(dx**2)
+ np.diff(v_old[1:-1, :], n=2, axis=1)/(dy**2)) )
#apply the velocity boundary condition to the intermediate velocity data
#apply the velocity boundary condition
imax = len(u)
jmax = np.size(u)//imax
u[0, jmax//2:] = u[1,jmax//2:]
u[-1,:] = u[-2,:]
return u, v
def solve_NS(u, v, p, rho, nu, courant, dt_min, t_end, dx, dy, rtol = 1.e-5, logs = False, outint = 100):
""" Solve the incompressible Navier-Stokes equations
using a lot of the numerical choices and approaches we've seen
earlier in this lecture.
"""
t = 0
u_old = u.copy()
v_old = v.copy()
p_RHS = np.zeros_like(X)
time_it = 0
while t < t_end:
time_it+=1
#set dt based on courant number
vel_max = np.max(np.sqrt(u**2+v**2))
if vel_max>0.0:
dt = min(courant*min(dx,dy)/vel_max,dt_min)
else:
dt = dt_min
t += dt
if logs and time_it%outint == 0:
print('\nTime = {:.8f}'.format(t))
# calculate intermediate velocities
u, v = calculate_intermediate_velocity(nu, u, v, u_old, v_old, dt, dx, dy)
# PPM
# calculate RHS for the pressure Poisson problem
p_RHS = calculate_ppm_RHS_central(rho, u, v, p_RHS, dt, dx, dy)
# compute pressure - note that we use the previous p as an initial guess to the solution
p = pressure_poisson_jacobi(p, dx, dy, p_RHS, 1.e-5, logs = (logs and time_it%outint == 0))
# project velocity
u, v = project_velocity(rho, u, v, dt, dx, dy, p)
if logs and time_it%outint == 0:
print('norm(u) = {0:.8f}, norm(v) = {1:.8f}'.format(sl.norm(u),sl.norm(v)))
print('Courant number: {0:.8f}'.format(np.max(np.sqrt(u**2+v**2)) * dt / min(dx,dy)))
#swap pointers without copying data
temp = u_old
u_old = u
u = temp
temp = v_old
v_old = v
v = temp
return u, v, p
# physical parameters
rho = 1
nu = 1./10.
# define spatial mesh
# Size of rectangular domain
Lx = 4.0
Ly = 1.0
P_max = 0.5
# Number of grid points in each direction, including boundary nodes
Nx = 201
Ny = 51
# hence the mesh spacing
dx = Lx/(Nx-1)
dy = Ly/(Ny-1)
# read the docs to see the ordering that mgrid gives us
X, Y = np.mgrid[0:Nx:1, 0:Ny:1]
X = dx*X
Y = dy*Y
# the following is an alternative to the three lines above
#X, Y = np.mgrid[0: Lx + 1e-10: dx, 0: Ly + 1e-10: dy]
# but without the need to add a "small" increement to ensure
# the Lx and Ly end points are included
# initialise independent variables
u = np.zeros_like(X)
v = np.zeros_like(X)
p = np.zeros_like(X)
# Apply Dirichlet BCs to u and v - the code below doesn't touch
# these so we can do this once outside the time loop
u[:, -1] = 0
u[:, 0] = 0
u[0, :] = 0
u[-1, :] = 0
v[:, -1] = 0
v[:, 0] = 0
v[0, :] = 0
v[-1, :] = 0
imax = len(p)
jmax = np.size(p)//imax
p[0, (jmax//2):] = P_max
p[-1, :] = 0.0
#set a Courant number and use dynamic time step
courant = 0.005
dt_min = 1.e-3
t_end = 10.0
import time
start = time.time()
u, v, p = solve_NS(u, v, p, rho, nu, courant, dt_min, t_end, dx, dy, rtol=1.e-6, logs = True, outint = 1000)
end = time.time()
print('Time taken by calculation = ', end - start)
Time = 0.92504142
pressure solve iterations = 9
norm(u) = 5.75588435, norm(v) = 0.82023177
Courant number: 0.00500413
Time = 1.56396603
pressure solve iterations = 4
norm(u) = 7.17385750, norm(v) = 1.01437241
Courant number: 0.00500118
Time = 2.12247992
pressure solve iterations = 2
norm(u) = 7.75132786, norm(v) = 1.09305978
Courant number: 0.00500050
Time = 2.65200717
pressure solve iterations = 2
norm(u) = 8.02592629, norm(v) = 1.13039598
Courant number: 0.00500024
Time = 3.16846933
pressure solve iterations = 1
norm(u) = 8.16407065, norm(v) = 1.14915841
Courant number: 0.00500012
Time = 3.67852842
pressure solve iterations = 1
norm(u) = 8.23535741, norm(v) = 1.15883487
Courant number: 0.00500006
Time = 4.18532702
pressure solve iterations = 1
norm(u) = 8.27260476, norm(v) = 1.16388939
Courant number: 0.00500003
Time = 4.69043372
pressure solve iterations = 1
norm(u) = 8.29219033, norm(v) = 1.16654680
Courant number: 0.00500002
Time = 5.19465397
pressure solve iterations = 1
norm(u) = 8.30252271, norm(v) = 1.16794865
Courant number: 0.00500001
Time = 5.69840745
pressure solve iterations = 1
norm(u) = 8.30798273, norm(v) = 1.16868950
Courant number: 0.00500000
Time = 6.20191444
pressure solve iterations = 1
norm(u) = 8.31087038, norm(v) = 1.16908165
Courant number: 0.00500000
Time = 6.70529101
pressure solve iterations = 1
norm(u) = 8.31239776, norm(v) = 1.16928979
Courant number: 0.00500000
Time = 7.20859847
pressure solve iterations = 1
norm(u) = 8.31320551, norm(v) = 1.16940038
Courant number: 0.00500000
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[5], line 57
55 import time
56 start = time.time()
---> 57 u, v, p = solve_NS(u, v, p, rho, nu, courant, dt_min, t_end, dx, dy, rtol=1.e-6, logs = True, outint = 1000)
58 end = time.time()
59 print('Time taken by calculation = ', end - start)
Cell In[4], line 33, in solve_NS(u, v, p, rho, nu, courant, dt_min, t_end, dx, dy, rtol, logs, outint)
31 p_RHS = calculate_ppm_RHS_central(rho, u, v, p_RHS, dt, dx, dy)
32 # compute pressure - note that we use the previous p as an initial guess to the solution
---> 33 p = pressure_poisson_jacobi(p, dx, dy, p_RHS, 1.e-5, logs = (logs and time_it%outint == 0))
34 # project velocity
35 u, v = project_velocity(rho, u, v, dt, dx, dy, p)
Cell In[2], line 32, in pressure_poisson_jacobi(p, dx, dy, RHS, rtol, logs)
29 p[0, :-(jmax-jmax//2)] = p[1, :-(jmax-jmax//2)]
31 # relative change in pressure
---> 32 tol = sl.norm(p - p_old)/np.maximum(1.0e-10,sl.norm(p))
34 #swap arrays without copying the data
35 temp = p_old
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/scipy/linalg/_misc.py:178, in norm(a, ord, axis, keepdims, check_finite)
175 return lange(*lange_args)
177 # fall back to numpy in every other case
--> 178 return np.linalg.norm(a, ord=ord, axis=axis, keepdims=keepdims)
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/numpy/linalg/linalg.py:2552, in norm(x, ord, axis, keepdims)
2550 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2551 else:
-> 2552 sqnorm = x.dot(x)
2553 ret = sqrt(sqnorm)
2554 if keepdims:
KeyboardInterrupt:
# set up figure
fig = plt.figure(figsize=(21, 7))
ax1 = fig.add_subplot(121)
cont = ax1.contourf(X,Y,p, cmap=cm.coolwarm)
fig.colorbar(cont)
# don't plot at every gird point - every 5th
ax1.quiver(X[::20,::5],Y[::20,::5],u[::20,::5],v[::20,::5])
ax1.set_xlim(-0.1, 4.1)
ax1.set_ylim(-0.1, 1.1)
ax1.set_xlabel('$x$', fontsize=16)
ax1.set_ylabel('$y$', fontsize=16)
ax1.set_title('Pressure driven problem - pressure and velocity vectors', fontsize=16)
ax1 = fig.add_subplot(122)
cont = ax1.contourf(X,Y,np.sqrt(u*u+v*v), cmap=cm.coolwarm)
fig.colorbar(cont)
ax1.set_xlim(-0.1, 4.1)
ax1.set_ylim(-0.1, 1.1)
ax1.set_xlabel('$x$', fontsize=16)
ax1.set_ylabel('$y$', fontsize=16)
ax1.set_title('Pressure driven problem - speed', fontsize=16)
Text(0.5, 1.0, 'Pressure driven problem - speed')
flow_in = np.sum(u[0,:])*Ly
flow_out = np.sum(u[-1,:])*Ly
print('flow in ',flow_in)
print('flow out ',flow_out)
flow in 3.6400883324252575
flow out 3.6920762379752325
# physical parameters
rho = 1
nu = 1./10.
# define spatial mesh
# Size of rectangular domain
Lx = 4.0
Ly = 1.0
P_max = 5.0
# Number of grid points in each direction, including boundary nodes
Nx = 201
Ny = 51
# hence the mesh spacing
dx = Lx/(Nx-1)
dy = Ly/(Ny-1)
# read the docs to see the ordering that mgrid gives us
X, Y = np.mgrid[0:Nx:1, 0:Ny:1]
X = dx*X
Y = dy*Y
# the following is an alternative to the three lines above
#X, Y = np.mgrid[0: Lx + 1e-10: dx, 0: Ly + 1e-10: dy]
# but without the need to add a "small" increement to ensure
# the Lx and Ly end points are included
# initialise independent variables
u = np.zeros_like(X)
v = np.zeros_like(X)
p = np.zeros_like(X)
# Apply Dirichlet BCs to u and v - the code below doesn't touch
# these so we can do this once outside the time loop
u[:, -1] = 0
u[:, 0] = 0
u[0, :] = 0
u[-1, :] = 0
v[:, -1] = 0
v[:, 0] = 0
v[0, :] = 0
v[-1, :] = 0
imax = len(p)
jmax = np.size(p)//imax
p[0, (jmax//2):] = P_max
p[-1, :] = 0.0
#set a Courant number and use dynamic time step
courant = 0.005
dt_min = 1.e-3
t_end = 10.0
import time
start = time.time()
u, v, p = solve_NS(u, v, p, rho, nu, courant, dt_min, t_end, dx, dy, rtol=1.e-6, logs = True, outint = 1000)
end = time.time()
print('Time taken by calculation = ', end - start)
Time = 0.29128222
pressure solve iterations = 4
norm(u) = 26.40120295, norm(v) = 3.82856151
Courant number: 0.00500452
Time = 0.42434206
pressure solve iterations = 2
norm(u) = 35.59942847, norm(v) = 5.05105750
Courant number: 0.00500213
Time = 0.53042189
pressure solve iterations = 3
norm(u) = 42.01154659, norm(v) = 5.87103940
Courant number: 0.00500135
Time = 0.62274982
pressure solve iterations = 2
norm(u) = 47.02031125, norm(v) = 6.49357974
Courant number: 0.00500097
Time = 0.70638201
pressure solve iterations = 1
norm(u) = 51.15035533, norm(v) = 6.99497900
Courant number: 0.00500075
Time = 0.78389566
pressure solve iterations = 2
norm(u) = 54.66746928, norm(v) = 7.41343505
Courant number: 0.00500061
Time = 0.85681524
pressure solve iterations = 1
norm(u) = 57.72801356, norm(v) = 7.77118401
Courant number: 0.00500050
Time = 0.92613075
pressure solve iterations = 1
norm(u) = 60.43299488, norm(v) = 8.08241933
Courant number: 0.00500043
Time = 0.99252778
pressure solve iterations = 1
norm(u) = 62.85199799, norm(v) = 8.35690594
Courant number: 0.00500037
Time = 1.05650399
pressure solve iterations = 2
norm(u) = 65.03529862, norm(v) = 8.60158838
Courant number: 0.00500032
Time = 1.11843438
pressure solve iterations = 2
norm(u) = 67.02059863, norm(v) = 8.82146877
Courant number: 0.00500029
Time = 1.17860935
pressure solve iterations = 1
norm(u) = 68.83695508, norm(v) = 9.02050263
Courant number: 0.00500026
Time = 1.23725804
pressure solve iterations = 1
norm(u) = 70.50726990, norm(v) = 9.20181647
Courant number: 0.00500023
Time = 1.29456085
pressure solve iterations = 1
norm(u) = 72.04986099, norm(v) = 9.36784337
Courant number: 0.00500021
Time = 1.35067099
pressure solve iterations = 2
norm(u) = 73.47988751, norm(v) = 9.52052599
Courant number: 0.00500019
Time = 1.40571684
pressure solve iterations = 2
norm(u) = 74.80994744, norm(v) = 9.66149149
Courant number: 0.00500018
Time = 1.45980673
pressure solve iterations = 1
norm(u) = 76.05061952, norm(v) = 9.79211782
Courant number: 0.00500016
Time = 1.51303302
pressure solve iterations = 2
norm(u) = 77.21088618, norm(v) = 9.91354762
Courant number: 0.00500015
Time = 1.56547517
pressure solve iterations = 1
norm(u) = 78.29844717, norm(v) = 10.02673266
Courant number: 0.00500014
Time = 1.61720209
pressure solve iterations = 2
norm(u) = 79.31995547, norm(v) = 10.13249019
Courant number: 0.00500013
Time = 1.66827399
pressure solve iterations = 1
norm(u) = 80.28120097, norm(v) = 10.23153903
Courant number: 0.00500012
Time = 1.71874387
pressure solve iterations = 2
norm(u) = 81.18725773, norm(v) = 10.32450186
Courant number: 0.00500011
Time = 1.76865862
pressure solve iterations = 1
norm(u) = 82.04260048, norm(v) = 10.41191316
Courant number: 0.00500010
Time = 1.81805994
pressure solve iterations = 2
norm(u) = 82.85119651, norm(v) = 10.49423886
Courant number: 0.00500010
Time = 1.86698512
pressure solve iterations = 2
norm(u) = 83.61657969, norm(v) = 10.57189761
Courant number: 0.00500009
Time = 1.91546766
pressure solve iterations = 2
norm(u) = 84.34191271, norm(v) = 10.64526224
Courant number: 0.00500009
Time = 1.96353776
pressure solve iterations = 2
norm(u) = 85.03003841, norm(v) = 10.71466250
Courant number: 0.00500008
Time = 2.01122277
pressure solve iterations = 2
norm(u) = 85.68352222, norm(v) = 10.78038983
Courant number: 0.00500008
Time = 2.05854752
pressure solve iterations = 1
norm(u) = 86.30468701, norm(v) = 10.84270733
Courant number: 0.00500007
Time = 2.10553465
pressure solve iterations = 2
norm(u) = 86.89564292, norm(v) = 10.90185531
Courant number: 0.00500007
Time = 2.15220487
pressure solve iterations = 2
norm(u) = 87.45831299, norm(v) = 10.95805062
Courant number: 0.00500006
Time = 2.19857713
pressure solve iterations = 1
norm(u) = 87.99445485, norm(v) = 11.01148794
Courant number: 0.00500006
Time = 2.24466887
pressure solve iterations = 2
norm(u) = 88.50567835, norm(v) = 11.06234467
Courant number: 0.00500006
Time = 2.29049616
pressure solve iterations = 2
norm(u) = 88.99346337, norm(v) = 11.11078364
Courant number: 0.00500005
Time = 2.33607384
pressure solve iterations = 1
norm(u) = 89.45917273, norm(v) = 11.15695414
Courant number: 0.00500005
Time = 2.38141564
pressure solve iterations = 1
norm(u) = 89.90406422, norm(v) = 11.20099295
Courant number: 0.00500005
Time = 2.42653428
pressure solve iterations = 1
norm(u) = 90.32930113, norm(v) = 11.24302500
Courant number: 0.00500005
Time = 2.47144161
pressure solve iterations = 1
norm(u) = 90.73596149, norm(v) = 11.28316563
Courant number: 0.00500004
Time = 2.51614863
pressure solve iterations = 2
norm(u) = 91.12504598, norm(v) = 11.32152191
Courant number: 0.00500004
Time = 2.56066561
pressure solve iterations = 2
norm(u) = 91.49748520, norm(v) = 11.35819300
Courant number: 0.00500004
Time = 2.60500213
pressure solve iterations = 1
norm(u) = 91.85414598, norm(v) = 11.39327061
Courant number: 0.00500004
Time = 2.64916715
pressure solve iterations = 2
norm(u) = 92.19583708, norm(v) = 11.42683968
Courant number: 0.00500004
Time = 2.69316908
pressure solve iterations = 1
norm(u) = 92.52331420, norm(v) = 11.45897946
Courant number: 0.00500004
Time = 2.73701579
pressure solve iterations = 2
norm(u) = 92.83728437, norm(v) = 11.48976398
Courant number: 0.00500003
Time = 2.78071467
pressure solve iterations = 2
norm(u) = 93.13841005, norm(v) = 11.51926219
Courant number: 0.00500003
Time = 2.82427266
pressure solve iterations = 2
norm(u) = 93.42731268, norm(v) = 11.54753874
Courant number: 0.00500003
Time = 2.86769632
pressure solve iterations = 1
norm(u) = 93.70457601, norm(v) = 11.57465369
Courant number: 0.00500003
Time = 2.91099179
pressure solve iterations = 2
norm(u) = 93.97074890, norm(v) = 11.60066379
Courant number: 0.00500003
Time = 2.95416491
pressure solve iterations = 2
norm(u) = 94.22634802, norm(v) = 11.62562224
Courant number: 0.00500003
Time = 2.99722115
pressure solve iterations = 1
norm(u) = 94.47186026, norm(v) = 11.64957898
Courant number: 0.00500003
Time = 3.04016571
pressure solve iterations = 2
norm(u) = 94.70774483, norm(v) = 11.67258096
Courant number: 0.00500002
Time = 3.08300349
pressure solve iterations = 2
norm(u) = 94.93443532, norm(v) = 11.69467234
Courant number: 0.00500002
Time = 3.12573913
pressure solve iterations = 2
norm(u) = 95.15234147, norm(v) = 11.71589485
Courant number: 0.00500002
Time = 3.16837704
pressure solve iterations = 1
norm(u) = 95.36185085, norm(v) = 11.73628787
Courant number: 0.00500002
Time = 3.21092138
pressure solve iterations = 2
norm(u) = 95.56333036, norm(v) = 11.75588840
Courant number: 0.00500002
Time = 3.25337613
pressure solve iterations = 2
norm(u) = 95.75712758, norm(v) = 11.77473162
Courant number: 0.00500002
Time = 3.29574503
pressure solve iterations = 1
norm(u) = 95.94357207, norm(v) = 11.79285078
Courant number: 0.00500002
Time = 3.33803167
pressure solve iterations = 1
norm(u) = 96.12297648, norm(v) = 11.81027743
Courant number: 0.00500002
Time = 3.38023944
pressure solve iterations = 1
norm(u) = 96.29563769, norm(v) = 11.82704139
Courant number: 0.00500002
Time = 3.42237157
pressure solve iterations = 2
norm(u) = 96.46183774, norm(v) = 11.84317099
Courant number: 0.00500002
Time = 3.46443115
pressure solve iterations = 2
norm(u) = 96.62184479, norm(v) = 11.85869309
Courant number: 0.00500002
Time = 3.50642109
pressure solve iterations = 1
norm(u) = 96.77591391, norm(v) = 11.87363319
Courant number: 0.00500002
Time = 3.54834421
pressure solve iterations = 2
norm(u) = 96.92428788, norm(v) = 11.88801559
Courant number: 0.00500002
Time = 3.59020316
pressure solve iterations = 2
norm(u) = 97.06719794, norm(v) = 11.90186326
Courant number: 0.00500001
Time = 3.63200049
pressure solve iterations = 1
norm(u) = 97.20486443, norm(v) = 11.91519819
Courant number: 0.00500001
Time = 3.67373862
pressure solve iterations = 2
norm(u) = 97.33749742, norm(v) = 11.92804124
Courant number: 0.00500001
Time = 3.71541987
pressure solve iterations = 2
norm(u) = 97.46529731, norm(v) = 11.94041230
Courant number: 0.00500001
Time = 3.75704645
pressure solve iterations = 1
norm(u) = 97.58845531, norm(v) = 11.95233037
Courant number: 0.00500001
Time = 3.79862049
pressure solve iterations = 2
norm(u) = 97.70715402, norm(v) = 11.96381350
Courant number: 0.00500001
Time = 3.84014400
pressure solve iterations = 2
norm(u) = 97.82156782, norm(v) = 11.97487896
Courant number: 0.00500001
Time = 3.88161893
pressure solve iterations = 2
norm(u) = 97.93186337, norm(v) = 11.98554319
Courant number: 0.00500001
Time = 3.92304712
pressure solve iterations = 2
norm(u) = 98.03820000, norm(v) = 11.99582194
Courant number: 0.00500001
Time = 3.96443035
pressure solve iterations = 1
norm(u) = 98.14073009, norm(v) = 12.00573025
Courant number: 0.00500001
Time = 4.00577032
pressure solve iterations = 1
norm(u) = 98.23959943, norm(v) = 12.01528246
Courant number: 0.00500001
Time = 4.04706866
pressure solve iterations = 2
norm(u) = 98.33494755, norm(v) = 12.02449234
Courant number: 0.00500001
Time = 4.08832693
pressure solve iterations = 1
norm(u) = 98.42690806, norm(v) = 12.03337304
Courant number: 0.00500001
Time = 4.12954661
pressure solve iterations = 1
norm(u) = 98.51560892, norm(v) = 12.04193711
Courant number: 0.00500001
Time = 4.17072915
pressure solve iterations = 2
norm(u) = 98.60117271, norm(v) = 12.05019659
Courant number: 0.00500001
Time = 4.21187592
pressure solve iterations = 1
norm(u) = 98.68371693, norm(v) = 12.05816301
Courant number: 0.00500001
Time = 4.25298824
pressure solve iterations = 1
norm(u) = 98.76335422, norm(v) = 12.06584741
Courant number: 0.00500001
Time = 4.29406738
pressure solve iterations = 2
norm(u) = 98.84019258, norm(v) = 12.07326038
Courant number: 0.00500001
Time = 4.33511454
pressure solve iterations = 1
norm(u) = 98.91433565, norm(v) = 12.08041206
Courant number: 0.00500001
Time = 4.37613091
pressure solve iterations = 2
norm(u) = 98.98588282, norm(v) = 12.08731217
Courant number: 0.00500001
Time = 4.41711760
pressure solve iterations = 2
norm(u) = 99.05492950, norm(v) = 12.09397004
Courant number: 0.00500001
Time = 4.45807568
pressure solve iterations = 1
norm(u) = 99.12156730, norm(v) = 12.10039462
Courant number: 0.00500001
Time = 4.49900619
pressure solve iterations = 1
norm(u) = 99.18588416, norm(v) = 12.10659450
Courant number: 0.00500001
Time = 4.53991012
pressure solve iterations = 2
norm(u) = 99.24796454, norm(v) = 12.11257791
Courant number: 0.00500001
Time = 4.58078844
pressure solve iterations = 2
norm(u) = 99.30788960, norm(v) = 12.11835277
Courant number: 0.00500001
Time = 4.62164206
pressure solve iterations = 2
norm(u) = 99.36573730, norm(v) = 12.12392669
Courant number: 0.00500001
Time = 4.66247186
pressure solve iterations = 2
norm(u) = 99.42158259, norm(v) = 12.12930696
Courant number: 0.00500001
Time = 4.70327869
pressure solve iterations = 1
norm(u) = 99.47549751, norm(v) = 12.13450060
Courant number: 0.00500001
Time = 4.74406338
pressure solve iterations = 1
norm(u) = 99.52755130, norm(v) = 12.13951435
Courant number: 0.00500001
Time = 4.78482670
pressure solve iterations = 1
norm(u) = 99.57781057, norm(v) = 12.14435468
Courant number: 0.00500001
Time = 4.82556942
pressure solve iterations = 1
norm(u) = 99.62633938, norm(v) = 12.14902781
Courant number: 0.00500000
Time = 4.86629226
pressure solve iterations = 2
norm(u) = 99.67319935, norm(v) = 12.15353974
Courant number: 0.00500000
Time = 4.90699593
pressure solve iterations = 2
norm(u) = 99.71844978, norm(v) = 12.15789624
Courant number: 0.00500000
Time = 4.94768110
pressure solve iterations = 2
norm(u) = 99.76214771, norm(v) = 12.16210285
Courant number: 0.00500000
Time = 4.98834841
pressure solve iterations = 2
norm(u) = 99.80434807, norm(v) = 12.16616489
Courant number: 0.00500000
Time = 5.02899851
pressure solve iterations = 1
norm(u) = 99.84510371, norm(v) = 12.17008751
Courant number: 0.00500000
Time = 5.06963198
pressure solve iterations = 1
norm(u) = 99.88446553, norm(v) = 12.17387563
Courant number: 0.00500000
Time = 5.11024941
pressure solve iterations = 2
norm(u) = 99.92248252, norm(v) = 12.17753400
Courant number: 0.00500000
Time = 5.15085136
pressure solve iterations = 2
norm(u) = 99.95920187, norm(v) = 12.18106722
Courant number: 0.00500000
Time = 5.19143838
pressure solve iterations = 1
norm(u) = 99.99466904, norm(v) = 12.18447967
Courant number: 0.00500000
Time = 5.23201097
pressure solve iterations = 1
norm(u) = 100.02892781, norm(v) = 12.18777559
Courant number: 0.00500000
Time = 5.27256964
pressure solve iterations = 2
norm(u) = 100.06202035, norm(v) = 12.19095908
Courant number: 0.00500000
Time = 5.31311487
pressure solve iterations = 1
norm(u) = 100.09398731, norm(v) = 12.19403407
Courant number: 0.00500000
Time = 5.35364714
pressure solve iterations = 1
norm(u) = 100.12486783, norm(v) = 12.19700434
Courant number: 0.00500000
Time = 5.39416687
pressure solve iterations = 2
norm(u) = 100.15469966, norm(v) = 12.19987354
Courant number: 0.00500000
Time = 5.43467452
pressure solve iterations = 1
norm(u) = 100.18351918, norm(v) = 12.20264520
Courant number: 0.00500000
Time = 5.47517049
pressure solve iterations = 2
norm(u) = 100.21136144, norm(v) = 12.20532271
Courant number: 0.00500000
Time = 5.51565519
pressure solve iterations = 1
norm(u) = 100.23826023, norm(v) = 12.20790934
Courant number: 0.00500000
Time = 5.55612901
pressure solve iterations = 1
norm(u) = 100.26424816, norm(v) = 12.21040822
Courant number: 0.00500000
Time = 5.59659231
pressure solve iterations = 2
norm(u) = 100.28935662, norm(v) = 12.21282240
Courant number: 0.00500000
Time = 5.63704547
pressure solve iterations = 1
norm(u) = 100.31361593, norm(v) = 12.21515480
Courant number: 0.00500000
Time = 5.67748882
pressure solve iterations = 2
norm(u) = 100.33705530, norm(v) = 12.21740826
Courant number: 0.00500000
Time = 5.71792271
pressure solve iterations = 1
norm(u) = 100.35970291, norm(v) = 12.21958548
Courant number: 0.00500000
Time = 5.75834745
pressure solve iterations = 1
norm(u) = 100.38158594, norm(v) = 12.22168910
Courant number: 0.00500000
Time = 5.79876337
pressure solve iterations = 2
norm(u) = 100.40273061, norm(v) = 12.22372164
Courant number: 0.00500000
Time = 5.83917075
pressure solve iterations = 1
norm(u) = 100.42316221, norm(v) = 12.22568555
Courant number: 0.00500000
Time = 5.87956990
pressure solve iterations = 1
norm(u) = 100.44290514, norm(v) = 12.22758317
Courant number: 0.00500000
Time = 5.91996109
pressure solve iterations = 1
norm(u) = 100.46198294, norm(v) = 12.22941679
Courant number: 0.00500000
Time = 5.96034460
pressure solve iterations = 2
norm(u) = 100.48041834, norm(v) = 12.23118859
Courant number: 0.00500000
Time = 6.00072068
pressure solve iterations = 2
norm(u) = 100.49823324, norm(v) = 12.23290069
Courant number: 0.00500000
Time = 6.04108959
pressure solve iterations = 2
norm(u) = 100.51544880, norm(v) = 12.23455512
Courant number: 0.00500000
Time = 6.08145157
pressure solve iterations = 2
norm(u) = 100.53208543, norm(v) = 12.23615386
Courant number: 0.00500000
Time = 6.12180686
pressure solve iterations = 1
norm(u) = 100.54816284, norm(v) = 12.23769881
Courant number: 0.00500000
Time = 6.16215568
pressure solve iterations = 1
norm(u) = 100.56370004, norm(v) = 12.23919179
Courant number: 0.00500000
Time = 6.20249826
pressure solve iterations = 1
norm(u) = 100.57871539, norm(v) = 12.24063458
Courant number: 0.00500000
Time = 6.24283480
pressure solve iterations = 2
norm(u) = 100.59322661, norm(v) = 12.24202888
Courant number: 0.00500000
Time = 6.28316550
pressure solve iterations = 2
norm(u) = 100.60725079, norm(v) = 12.24337634
Courant number: 0.00500000
Time = 6.32349058
pressure solve iterations = 1
norm(u) = 100.62080445, norm(v) = 12.24467855
Courant number: 0.00500000
Time = 6.36381021
pressure solve iterations = 2
norm(u) = 100.63390354, norm(v) = 12.24593705
Courant number: 0.00500000
Time = 6.40412458
pressure solve iterations = 1
norm(u) = 100.64656343, norm(v) = 12.24715333
Courant number: 0.00500000
Time = 6.44443388
pressure solve iterations = 1
norm(u) = 100.65879900, norm(v) = 12.24832880
Courant number: 0.00500000
Time = 6.48473826
pressure solve iterations = 1
norm(u) = 100.67062459, norm(v) = 12.24946486
Courant number: 0.00500000
Time = 6.52503789
pressure solve iterations = 2
norm(u) = 100.68205405, norm(v) = 12.25056283
Courant number: 0.00500000
Time = 6.56533294
pressure solve iterations = 1
norm(u) = 100.69310076, norm(v) = 12.25162401
Courant number: 0.00500000
Time = 6.60562357
pressure solve iterations = 2
norm(u) = 100.70377765, norm(v) = 12.25264964
Courant number: 0.00500000
Time = 6.64590991
pressure solve iterations = 1
norm(u) = 100.71409719, norm(v) = 12.25364092
Courant number: 0.00500000
Time = 6.68619211
pressure solve iterations = 1
norm(u) = 100.72407143, norm(v) = 12.25459900
Courant number: 0.00500000
Time = 6.72647031
pressure solve iterations = 2
norm(u) = 100.73371201, norm(v) = 12.25552502
Courant number: 0.00500000
Time = 6.76674466
pressure solve iterations = 1
norm(u) = 100.74303017, norm(v) = 12.25642005
Courant number: 0.00500000
Time = 6.80701527
pressure solve iterations = 1
norm(u) = 100.75203676, norm(v) = 12.25728513
Courant number: 0.00500000
Time = 6.84728227
pressure solve iterations = 1
norm(u) = 100.76074227, norm(v) = 12.25812128
Courant number: 0.00500000
Time = 6.88754578
pressure solve iterations = 2
norm(u) = 100.76915682, norm(v) = 12.25892947
Courant number: 0.00500000
Time = 6.92780593
pressure solve iterations = 2
norm(u) = 100.77729021, norm(v) = 12.25971064
Courant number: 0.00500000
Time = 6.96806282
pressure solve iterations = 2
norm(u) = 100.78515188, norm(v) = 12.26046570
Courant number: 0.00500000
Time = 7.00831656
pressure solve iterations = 1
norm(u) = 100.79275096, norm(v) = 12.26119553
Courant number: 0.00500000
Time = 7.04856726
pressure solve iterations = 1
norm(u) = 100.80009627, norm(v) = 12.26190097
Courant number: 0.00500000
Time = 7.08881502
pressure solve iterations = 1
norm(u) = 100.80719634, norm(v) = 12.26258285
Courant number: 0.00500000
Time = 7.12905994
pressure solve iterations = 1
norm(u) = 100.81405938, norm(v) = 12.26324195
Courant number: 0.00500000
Time = 7.16930212
pressure solve iterations = 1
norm(u) = 100.82069335, norm(v) = 12.26387905
Courant number: 0.00500000
Time = 7.20954164
pressure solve iterations = 1
norm(u) = 100.82710594, norm(v) = 12.26449488
Courant number: 0.00500000
Time = 7.24977859
pressure solve iterations = 1
norm(u) = 100.83330457, norm(v) = 12.26509015
Courant number: 0.00500000
Time = 7.29001307
pressure solve iterations = 1
norm(u) = 100.83929641, norm(v) = 12.26566556
Courant number: 0.00500000
Time = 7.33024515
pressure solve iterations = 1
norm(u) = 100.84508838, norm(v) = 12.26622176
Courant number: 0.00500000
Time = 7.37047491
pressure solve iterations = 2
norm(u) = 100.85068719, norm(v) = 12.26675941
Courant number: 0.00500000
Time = 7.41070244
pressure solve iterations = 2
norm(u) = 100.85609929, norm(v) = 12.26727912
Courant number: 0.00500000
Time = 7.45092780
pressure solve iterations = 2
norm(u) = 100.86133095, norm(v) = 12.26778150
Courant number: 0.00500000
Time = 7.49115107
pressure solve iterations = 1
norm(u) = 100.86638820, norm(v) = 12.26826713
Courant number: 0.00500000
Time = 7.53137232
pressure solve iterations = 1
norm(u) = 100.87127687, norm(v) = 12.26873656
Courant number: 0.00500000
Time = 7.57159162
pressure solve iterations = 2
norm(u) = 100.87600261, norm(v) = 12.26919035
Courant number: 0.00500000
Time = 7.61180903
pressure solve iterations = 2
norm(u) = 100.88057086, norm(v) = 12.26962900
Courant number: 0.00500000
Time = 7.65202461
pressure solve iterations = 1
norm(u) = 100.88498689, norm(v) = 12.27005304
Courant number: 0.00500000
Time = 7.69223843
pressure solve iterations = 1
norm(u) = 100.88925579, norm(v) = 12.27046294
Courant number: 0.00500000
Time = 7.73245054
pressure solve iterations = 1
norm(u) = 100.89338248, norm(v) = 12.27085919
Courant number: 0.00500000
Time = 7.77266100
pressure solve iterations = 2
norm(u) = 100.89737170, norm(v) = 12.27124223
Courant number: 0.00500000
Time = 7.81286988
pressure solve iterations = 2
norm(u) = 100.90122805, norm(v) = 12.27161251
Courant number: 0.00500000
Time = 7.85307721
pressure solve iterations = 1
norm(u) = 100.90495597, norm(v) = 12.27197046
Courant number: 0.00500000
Time = 7.89328305
pressure solve iterations = 1
norm(u) = 100.90855975, norm(v) = 12.27231648
Courant number: 0.00500000
Time = 7.93348745
pressure solve iterations = 2
norm(u) = 100.91204353, norm(v) = 12.27265098
Courant number: 0.00500000
Time = 7.97369046
pressure solve iterations = 2
norm(u) = 100.91541132, norm(v) = 12.27297434
Courant number: 0.00500000
Time = 8.01389213
pressure solve iterations = 1
norm(u) = 100.91866698, norm(v) = 12.27328693
Courant number: 0.00500000
Time = 8.05409250
pressure solve iterations = 2
norm(u) = 100.92181427, norm(v) = 12.27358911
Courant number: 0.00500000
Time = 8.09429161
pressure solve iterations = 1
norm(u) = 100.92485679, norm(v) = 12.27388123
Courant number: 0.00500000
Time = 8.13448950
pressure solve iterations = 1
norm(u) = 100.92779804, norm(v) = 12.27416363
Courant number: 0.00500000
Time = 8.17468622
pressure solve iterations = 2
norm(u) = 100.93064140, norm(v) = 12.27443663
Courant number: 0.00500000
Time = 8.21488181
pressure solve iterations = 1
norm(u) = 100.93339014, norm(v) = 12.27470054
Courant number: 0.00500000
Time = 8.25507630
pressure solve iterations = 2
norm(u) = 100.93604741, norm(v) = 12.27495567
Courant number: 0.00500000
Time = 8.29526973
pressure solve iterations = 1
norm(u) = 100.93861625, norm(v) = 12.27520230
Courant number: 0.00500000
Time = 8.33546213
pressure solve iterations = 2
norm(u) = 100.94109963, norm(v) = 12.27544073
Courant number: 0.00500000
Time = 8.37565355
pressure solve iterations = 1
norm(u) = 100.94350037, norm(v) = 12.27567122
Courant number: 0.00500000
Time = 8.41584400
pressure solve iterations = 2
norm(u) = 100.94582125, norm(v) = 12.27589405
Courant number: 0.00500000
Time = 8.45603353
pressure solve iterations = 1
norm(u) = 100.94806492, norm(v) = 12.27610946
Courant number: 0.00500000
Time = 8.49622216
pressure solve iterations = 2
norm(u) = 100.95023395, norm(v) = 12.27631770
Courant number: 0.00500000
Time = 8.53640993
pressure solve iterations = 1
norm(u) = 100.95233084, norm(v) = 12.27651902
Courant number: 0.00500000
Time = 8.57659686
pressure solve iterations = 2
norm(u) = 100.95435798, norm(v) = 12.27671364
Courant number: 0.00500000
Time = 8.61678298
pressure solve iterations = 1
norm(u) = 100.95631769, norm(v) = 12.27690179
Courant number: 0.00500000
Time = 8.65696832
pressure solve iterations = 2
norm(u) = 100.95821223, norm(v) = 12.27708367
Courant number: 0.00500000
Time = 8.69715291
pressure solve iterations = 1
norm(u) = 100.96004377, norm(v) = 12.27725951
Courant number: 0.00500000
Time = 8.73733676
pressure solve iterations = 2
norm(u) = 100.96181439, norm(v) = 12.27742950
Courant number: 0.00500000
Time = 8.77751991
pressure solve iterations = 1
norm(u) = 100.96352614, norm(v) = 12.27759384
Courant number: 0.00500000
Time = 8.81770237
pressure solve iterations = 1
norm(u) = 100.96518096, norm(v) = 12.27775271
Courant number: 0.00500000
Time = 8.85788418
pressure solve iterations = 2
norm(u) = 100.96678076, norm(v) = 12.27790630
Courant number: 0.00500000
Time = 8.89806535
pressure solve iterations = 1
norm(u) = 100.96832737, norm(v) = 12.27805478
Courant number: 0.00500000
Time = 8.93824589
pressure solve iterations = 2
norm(u) = 100.96982256, norm(v) = 12.27819832
Courant number: 0.00500000
Time = 8.97842585
pressure solve iterations = 1
norm(u) = 100.97126804, norm(v) = 12.27833709
Courant number: 0.00500000
Time = 9.01860522
pressure solve iterations = 1
norm(u) = 100.97266546, norm(v) = 12.27847125
Courant number: 0.00500000
Time = 9.05878404
pressure solve iterations = 2
norm(u) = 100.97401642, norm(v) = 12.27860094
Courant number: 0.00500000
Time = 9.09896232
pressure solve iterations = 1
norm(u) = 100.97532247, norm(v) = 12.27872633
Courant number: 0.00500000
Time = 9.13914008
pressure solve iterations = 1
norm(u) = 100.97658511, norm(v) = 12.27884754
Courant number: 0.00500000
Time = 9.17931734
pressure solve iterations = 2
norm(u) = 100.97780577, norm(v) = 12.27896473
Courant number: 0.00500000
Time = 9.21949411
pressure solve iterations = 1
norm(u) = 100.97898586, norm(v) = 12.27907802
Courant number: 0.00500000
Time = 9.25967041
pressure solve iterations = 1
norm(u) = 100.98012673, norm(v) = 12.27918755
Courant number: 0.00500000
Time = 9.29984625
pressure solve iterations = 2
norm(u) = 100.98122967, norm(v) = 12.27929343
Courant number: 0.00500000
Time = 9.34002165
pressure solve iterations = 2
norm(u) = 100.98229595, norm(v) = 12.27939579
Courant number: 0.00500000
Time = 9.38019663
pressure solve iterations = 1
norm(u) = 100.98332680, norm(v) = 12.27949476
Courant number: 0.00500000
Time = 9.42037120
pressure solve iterations = 2
norm(u) = 100.98432338, norm(v) = 12.27959043
Courant number: 0.00500000
Time = 9.46054537
pressure solve iterations = 2
norm(u) = 100.98528684, norm(v) = 12.27968292
Courant number: 0.00500000
Time = 9.50071916
pressure solve iterations = 1
norm(u) = 100.98621828, norm(v) = 12.27977234
Courant number: 0.00500000
Time = 9.54089257
pressure solve iterations = 1
norm(u) = 100.98711877, norm(v) = 12.27985879
Courant number: 0.00500000
Time = 9.58106563
pressure solve iterations = 2
norm(u) = 100.98798933, norm(v) = 12.27994236
Courant number: 0.00500000
Time = 9.62123834
pressure solve iterations = 2
norm(u) = 100.98883096, norm(v) = 12.28002316
Courant number: 0.00500000
Time = 9.66141071
pressure solve iterations = 2
norm(u) = 100.98964463, norm(v) = 12.28010127
Courant number: 0.00500000
Time = 9.70158276
pressure solve iterations = 1
norm(u) = 100.99043125, norm(v) = 12.28017678
Courant number: 0.00500000
Time = 9.74175449
pressure solve iterations = 1
norm(u) = 100.99119174, norm(v) = 12.28024979
Courant number: 0.00500000
Time = 9.78192592
pressure solve iterations = 2
norm(u) = 100.99192695, norm(v) = 12.28032037
Courant number: 0.00500000
Time = 9.82209706
pressure solve iterations = 2
norm(u) = 100.99263774, norm(v) = 12.28038860
Courant number: 0.00500000
Time = 9.86226792
pressure solve iterations = 2
norm(u) = 100.99332490, norm(v) = 12.28045457
Courant number: 0.00500000
Time = 9.90243850
pressure solve iterations = 1
norm(u) = 100.99398924, norm(v) = 12.28051835
Courant number: 0.00500000
Time = 9.94260881
pressure solve iterations = 1
norm(u) = 100.99463150, norm(v) = 12.28058000
Courant number: 0.00500000
Time = 9.98277887
pressure solve iterations = 1
norm(u) = 100.99525242, norm(v) = 12.28063961
Courant number: 0.00500000
Time taken by calculation = 306.7849705219269
# set up figure
fig = plt.figure(figsize=(21, 7))
ax1 = fig.add_subplot(121)
cont = ax1.contourf(X,Y,p, cmap=cm.coolwarm)
fig.colorbar(cont)
# don't plot at every gird point - every 5th
ax1.quiver(X[::20,::5],Y[::20,::5],u[::20,::5],v[::20,::5])
ax1.set_xlim(-0.1, 4.1)
ax1.set_ylim(-0.1, 1.1)
ax1.set_xlabel('$x$', fontsize=16)
ax1.set_ylabel('$y$', fontsize=16)
ax1.set_title('Pressure driven problem - pressure and velocity vectors', fontsize=16)
ax1 = fig.add_subplot(122)
cont = ax1.contourf(X,Y,np.sqrt(u*u+v*v), cmap=cm.coolwarm)
fig.colorbar(cont)
ax1.set_xlim(-0.1, 4.1)
ax1.set_ylim(-0.1, 1.1)
ax1.set_xlabel('$x$', fontsize=16)
ax1.set_ylabel('$y$', fontsize=16)
ax1.set_title('Pressure driven problem - speed', fontsize=16)
Text(0.5, 1.0, 'Pressure driven problem - speed')
flow_in = np.sum(u[0,:])*Ly
flow_out = np.sum(u[-1,:])*Ly
print('flow in ',flow_in)
print('flow out ',flow_out)
flow in 44.11691653247884
flow out 44.30906374017106