Discritizing Volume Fraction Transport Equation in 2D¶
The goal is to discritize the volume fraction transport equation using the Crank-Nicolson Method, using the central difference to discritize the spatial derivatives. This is an implicit method in time.
\begin{equation} \frac{\partial \alpha}{\partial t} + \vec{v} \cdot \nabla \alpha = 0 \end{equation}The Unifying Theta Method¶
The unifying theta method is a simplified way of writing the Forward Euler, Backward Euler, and Crank-Nicolson schemes:
\begin{equation} \frac{\partial u}{\partial t} = G(u) \\ \frac{u_i^{n+1}-u_i^n}{\Delta t} = \Theta G(u_i^{n+1}) + (1-\Theta) G(u_i^n) \end{equation}Where $G(u)$ is a spatial discretization method, and $\Theta$ is a constant:
- $\Theta = 0.0$: Forward Euler
- $\Theta = 1.0$: Backward Euler
- $\Theta = 0.5$: Crank-Nicolson
Finite Difference Methods¶
The spatial discretization in the Crank-Nicolson Method is provided by one of the finite difference methods, typically Central Difference for the majority of the matrix and Forward and Backward for the boundaries.
Forward Difference¶
\begin{align} f'(x) = \frac{f(x+h) - f(x)}{h} + O(h^1) \label{eqn:1} \tag{1} \\ f'(x) = \frac{−3f(x) + 4f(x + h) − f(x + 2h)}{2h} + O(h^2) \label{eqn:2} \tag{2} \\\\ f''(x) = \frac{f(x+2h) - 2f(x+h) + f(x)}{h^2} + O(h^1) \label{eqn:3} \tag{3} \\ f''(x) = \frac{2f(x) − 5f(x + h) + 4f(x + 2h) − f(x + 3h)}{h^3} + O(h^2) \label{eqn:4} \tag{4} \end{align}Backward Difference¶
\begin{align} f'(x) = \frac{f(x) - f(x-h)}{h} + O(h^1) \label{eqn:5} \tag{5} \\ f'(x) = \frac{3f(x) − 4f(x − h) + f(x − 2h)}{2h} + O(h^2) \label{eqn:6} \tag{6} \\\\ f''(x) = \frac{f(x+2h)-2f(x+h) + f(x)}{h^2} + O(h^1) \label{eqn:7} \tag{7} \\ f''(x) = \frac{2f(x) − 5f(x − h) + 4f(x − 2h) − f(x − 3h)}{h} + O(h^2) \label{eqn:8} \tag{8} \end{align}Central Difference¶
\begin{align} f'(x) = \frac{f(x + h) − f(x − h)}{2h} + O(h^2) \label{eqn:9} \tag{9} \\ f'(x) = \frac{−f(x + 2h) + 8f(x + h) − 8f(x − h) + f(x − 2h)}{12h} + O(h^4) \label{eqn:10} \tag{10} \\\\ f''(x) = \frac{f(x+h) -2f(x) + f(x-h)}{h^2} + O(h^2) \label{eqn:11} \tag{11} \\ f''(x) = \frac{−f(x + 2h) + 16f(x + h) − 30f(x) + 16f(x − h) − f(x − 2h)}{12h^2} + O(h^4) \label{eqn:12} \tag{12} \end{align}Coefficient Matrix Layout¶
The creation of the coefficient matricies for $\alpha^{n}_{i,j}$ $\alpha^{n+1}_{i,j}$ using the $O(h^2)$ central difference scheme is represented below on a grid of size $Nx = 5$ by $Ny = 4$ elements (note node counts are N+1). The grid is numbered with $(i,j)$ indicies from $0$ to $N_x$, $0$ to $N_y$, respectively. The map from grid points with $(i,j)$ indexing is $j(N_x+1)+i$, so that the matrix index starts a the bottom left node (0,0) and moves across the grid from left to right, then indexes up one row, and starts at the $i = 0$ position again. The grid is annotated with this indicies as $M_i$ for matrix index, and $G_{i,j}$ for grid index.
An arbitrary point $G_{i,j}:(3,2), M_i:15$ is selected to represent the central difference method and the relationship to the coefficient matrix, located by the black dot on both the grid and the matrix. From here, the values of $\alpha^{n+1}_{i_{-1},j}$, $\alpha^{n+1}_{i_{+1},j}$, $\alpha^{n+1}_{i,j_{-1}}$, $\alpha^{n+1}_{i,j_{+1}}$ are computed, shown as the red and blue dots, respectively. Note that to transverse from $\alpha^{n+1}_{i,j}$ (black dot) to either $\alpha^{n+1}_{i,j_{-1}}$ or $\alpha^{n+1}_{i,j_{+1}}$ (blue dot) $N_x+1$ steps must be taken, as represented by the arrows on the spatial grid. This creates an $N_x+1$ offset diagonal on the coefficient matrix.
Note it's shape is diagonal banded Toeplitz Matrix when the grid is unbounded, with bandwidth $N_x + 1$. In a 1D spatial grid, this becomes an Hessenberg Matrix, a special form of the banded matrix with bandwidth = 1. When boundary conditions are applied, the matrix loses its Toeplitz diagonal repetion, and can only be classified as a Band Matrix. This is important when it comes to solving, as there are linear algebra techniques and specific solvers for these special forms which reduce computational cost on these sparce matricies.
In [1]:
from grid_plot import grid_plot
In [2]:
grid_plot()
In [3]:
from sympy import init_printing, symbols, Mul, Add, Matrix, Eq, Float
from sympy import UnevaluatedExpr as ue
import numpy as np
from IPython.display import display, Markdown, HTML
init_printing(use_latex=True)
In [4]:
## Helper Functions
def disp(x,y,mode='display'):
if isinstance(x,str):
x = symbols(x);
if isinstance(y,str):
y = symbols(y);
if isinstance(y,np.ndarray):
y = Matrix(y).T;
if 'disp' in mode:
display(Eq(x,y,evaluate=False));
else:
return Eq(x,y);
In [5]:
## Symbols
# Alpha (n+1) ij Matrix
alpha = Matrix([symbols([('\\alpha^{n+1}_{'+'i{0:+1d}\,j{1:+1d}'.format(i,j)+'}').replace('+0','') for j in range(-1,2)]) for i in range(-1,2)])
# Alpha (n) ij Matrix
alphan = Matrix([symbols([('\\alpha^{n}_{'+'i{0:+1d}\,j{1:+1d}'.format(i,j)+'}').replace('+0','') for j in range(-1,2)]) for i in range(-1,2)])
Pseudo $\alpha$ Grid Matrix¶
In [6]:
display(alpha)
In [7]:
theta,vx,vy,vxn,vyn,dx,dy,dt,Rx,Ry,v = symbols(['Θ','v_x^{n+1}','v_y^{n+1}','v_x^{n}','v_y^{n}','Δx','Δy','Δt','R_x','R_y','vbm'])
In [8]:
def fd(method='',axis='',f=alpha):
method,order = method
if axis == 'x': dir,h = 1,dx
elif axis == 'y': dir,h = -1,dy
if method == 'f':
if order == '1':
return (f[(1+1,1)[::dir]] - f[(1,1)[::dir]])/(h)
elif order == '2':
return (-3*f[(1,1)[::dir]] + 4*f[(1+1,1)[::dir]] -
f[(1+2,1)[::dir]])/(2*h)
if method == 'c':
if order == '2':
return (f[(1+1,1)[::dir]]-f[(1-1,1)[::dir]])/(2*h)
elif order == '4':
return (-f[(1+2,1)[::dir]] + 8*f[(1+1,1)[::dir]] -
8*f[(1-1,1)[::dir]] + f[(1-2,1)[::dir]])/(12*h)
if method == 'b':
if order == '1':
return (f[(1,1)[::dir]] - f[(1-1,1)[::dir]])/(h)
elif order == '2':
return (3*f[(1,1)[::dir]] - 4*f[(1-1,1)[::dir]] +
f[(1-2,1)[::dir]])/(2*h)
Unified Theta Method¶
Crank-Nicolson & 2nd Order Central FDM
In [9]:
# (u^n+1 - u^n)/Δt = Θ*G(U^n+1) + (1-Θ)*G(U^n)
methods = ['f1','c2','b1']
names = [('First','Forward'),('Second','Central'),('First','Backward')]
LHS = []; RHS = []; G = []; Gn = []
Rxs = dt/dx*vx
Rys = dt/dy*vy
## R Coefficients
display(Markdown('### R Coefficients\n' +
'These are the **CLF** numbers for the x and y steps and velocities.'))
display((disp('R_x',Rxs,'eqn'),disp('R_y',Rys,'eqn')))
for n,(method,name) in enumerate(zip(methods,names)):
## Equation Creating
G = theta*(-vx*fd(method,'x') -vy*fd(method,'y'))
Gn = (1-theta)*(-vx*fd(method,'x',alphan)-vy*fd(method,'y',alphan))
# Equations reordered
LHS.append((alpha[1,1] - theta*(-vx*fd(method,'x')*dt*(Rx/Rxs) -vy*fd(method,'y')*dt*(Ry/Rys))).factor(*alpha).simplify())
RHS.append(((1-theta)*(-vx*fd(method,'x',alphan)*dt*(Rx/Rxs) -vy*fd(method,'y',alphan)*dt*(Ry/Rys)) + alphan[1,1]).factor(*alpha).simplify())
## Printing
display(Markdown('### '+' Order '.join(name)+' Finite Difference'))
disp("G(u^{n+1})",G)
disp("G(u^{n})",Gn)
# Pretty Equation
display(Eq((alpha[1,1] - alphan[1,1])/dt,G+Gn))
# Sorted Equation
display(Markdown('Sorted and factored for $\\alpha_{i,j}$ terms:'))
display(Eq(LHS[n],RHS[n]))
print('\n')
In grid boundary areas, central difference cannot be used, as there are no $i-1$ and/or $j-1$ nodes, so forward and backward differences will be required.
Construct Coefficient Vectors¶
These are the coefficients for each diagonal.
In [10]:
def print_vec_table(alpha,vec,title):
fn = lambda x: ['$'+a+'$' for a in x]
display(Markdown('#### '+title))
header = ['|FDM'+(len(alpha)*'| ${:}$ ').format(*alpha.T),(len(alpha)+1)*'|---']
data = [name+(len(alpha)*'| {:~^20s} ').format(*C.astype(str)).replace('~',' ') for C,name in zip(vec,['FWD','CTR','BCK'])]
display(Markdown('|\n'.join(header+data)))
def create_coefficient_vec(alpha,eqns,printing=False):
Cx,Cy,Ca = np.tile(np.array(3*[len(alpha)*[Float(0)]]),(3,1,1,))
for m in range(0,3):
for n,a in enumerate(alpha.T):
for arg in eqns[m].args:
if a in arg.free_symbols:
sep = arg/a
if isinstance(sep,Add):
sep = sep.args
else:
sep = [sep]
for s in sep:
if Rx in s.free_symbols:
Cx[m][n] = (Cx[m][n] + s).factor(theta)
elif Ry in s.free_symbols:
Cy[m][n] = (Cy[m][n] + s).factor(theta)
else:
Ca[m][n] += s
if printing:
sup = str(alpha[0]).split('^')[1].split('_')[0]
[print_vec_table(alpha,C,'$C_'+var+'^'+sup+'$ Coefficient Vector') for C,var in zip([Cx,Cy],['x','y'])];
return Cx,Cy,Ca
Cx,Cy,Ca = create_coefficient_vec(alpha,LHS,printing=True)
Cnx,Cny,Cna = create_coefficient_vec(alphan,RHS,printing=True)
Symbolic Matricies¶
These are shown for a 6,6 grid, i.e.: 5x5 faces.
In [11]:
## Print Partial Matricies
def print_mat_table(mat,title,size=2.5,Nx=6,Ny=6):
nrc = int(size*Nx)
tb_style = 'style="border:1px solid black; font-weight: bold; width:30%;"'
td_style = 'style="text-align: center; width:65px; height: 60px; border:1px solid black;"'
string = ''
for C in mat[0:int(size*Nx)]:
string += ('\t<tr>\n'+(len(C[0:int(size*Nx)])*('\t\t<td '+td_style+'>{:}</td>\n')
).format(*C.astype(str)[0:int(size*Nx)])+'\n\t</tr>\n').replace('>0<','> <').replace('>0.0<','> <')
html = '<table '+tb_style+'>\n'+string+'\n</table>'
display(Markdown('### '+title))
display(HTML(html))
def create_diags(Cx,Cy,Nx,Ny,alpha,return_offsets=False,dtype='O',printing=False):
Nn = Nx*Ny
# Diagonal Offsets
offsets = [list(map(lambda j,i: j*(Nx)+i,*np.where(a==np.array(alpha))))[0]-(Nx+1) for a in alpha]
if return_offsets:
return offsets
# Create Diagonal Coefficient Arrays
xdiags = []; ydiags = [];
for n,(offset,Cxt,Cyt) in enumerate(zip(offsets,Cx.T,Cy.T)):
block = np.zeros(Nx,dtype='O')
diag = np.zeros(Nn,dtype='O')
# Central Difference
block[1:-1] = Cxt[1]
diag[Nx:-Nx] = Cyt[1]
# Forward Difference
block[0] = Cxt[0]
diag[0:Nx] = Cyt[0]
# Backward Difference
block[-1] = Cxt[2]
diag[-Nx::] = Cyt[2]
# Construct Diagonal
#print(np.tile(block,Ny)[0:Nn-abs(offset)])
xdiags.append(np.diag(np.tile(block,Ny)[::np.sign(offset) or 1][0:Nn-abs(offset)][::np.sign(offset) or 1],offset))
# This strips either leading or trailing zeros, making the diagonal the correct length
# np.sign(offset) reverses the array for - offsets, strips from the front, then reverses back
ydiags.append(np.diag(diag[::np.sign(offset) or 1][0:Nn-abs(offset)][::np.sign(offset) or 1],offset))
Cxd = np.sum(xdiags,0,dtype=dtype)
Cyd = np.sum(ydiags,0,dtype=dtype)
Cad = np.eye(Nn,dtype=dtype)
if printing:
sup = str(alpha[0]).split('^')[1].split('_')[0]
[print_mat_table(C,'$C_'+var+'^'+sup+'$ Coefficient Matrix') for C,var in zip([Cxd,Cyd],['x','y'])];
return Cxd,Cyd,Cad
In [12]:
Cxd,Cyd,Cad = create_diags(Cx=Cx,Cy=Cy,Nx=6,Ny=6,alpha=alpha,printing=True)
Cnxd,Cnyd,Cnad = create_diags(Cx=Cnx,Cy=Cny,Nx=6,Ny=6,alpha=alphan,printing=True)
In [13]:
## Write Symbolic Array to an Excel File for viewing
def write_matricies(Cxd,Cyd,Cnxd,Cnyd):
import pandas as pd
from openpyxl import load_workbook
import os
dfx = pd.DataFrame(Cxd).astype(str).replace(['^0$','^0\.0$','1.0\*'],['','',''],regex=True)
dfy = pd.DataFrame(Cyd).astype(str).replace(['^0$','^0\.0$','1.0\*'],['','',''],regex=True)
dfxn = pd.DataFrame(Cnxd).astype(str).replace(['^0$','^0\.0$','1.0\*'],['','',''],regex=True)
dfyn = pd.DataFrame(Cnyd).astype(str).replace(['^0$','^0\.0$','1.0\*'],['','',''],regex=True)
filepath = 'Array Layout.xlsx'
book = load_workbook(filepath)
writer = pd.ExcelWriter(filepath, engine='openpyxl')
writer.book = book
writer.sheets = dict((ws.title, ws) for ws in book.worksheets)
dfx.to_excel(writer,sheet_name='Ax', index=False, header=False, startrow=1, startcol=1)
dfy.to_excel(writer,sheet_name='Ay', index=False, header=False, startrow=1, startcol=1)
dfxn.to_excel(writer,sheet_name='Bx', index=False, header=False, startrow=1, startcol=1)
dfyn.to_excel(writer,sheet_name='By', index=False, header=False, startrow=1, startcol=1)
writer.save()
os.startfile(filepath, 'open')
Required Inputs¶
In [14]:
inputs = [fs for fs in (LHS[1]+RHS[1]).subs([(Rx,Rxs),(Ry,Rys)]).free_symbols if (fs not in alpha) and (fs not in alphan)]
display(inputs)
Numeric Solution¶
In [15]:
from scipy.sparse import spdiags
from scipy.linalg import inv, solve_banded
import numpy as np
from numpy import pi, sin, cos, sinh, cosh, floor
from tqdm import tqdm #tqdm_notebook as tqdm
In [16]:
def array_subs(ary,substitutions,unpack=False):
if not isinstance(ary,list):
ary = [ary]
unpack = True
data = []
for a in ary:
m,n = a.shape
x = np.array([c.subs(substitutions) for c in a.flatten()])
data.append(x.reshape(m,n))
if unpack:
return data[:]
else:
return data
def dbf(A,N,method='indirect'):
# Diagonal Banded Form
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html
temp = np.zeros((1+2*N,N**2))
for i in range(0,len(temp)):
if N-i >= 0:
temp[i,N-i::]=A.diagonal(N-i)
else:
temp[i,0:N-i]=A.diagonal(N-i)
if method == 'direct':
# Returns an array for direct submission to dgbsv
a = np.zeros((3*N + 1, temp.shape[1]), dtype='float64')
a[N:, :] = Adbf
return a
else:
return temp
def normalize(x):
return(x - x.min())/(x - x.min()).max()
In [17]:
## Input Parameters
theta_num,dx_num,dy_num = [0.5,0.04,0.04]
L,H,b,g,amp,N_periods,CFL = [2,1,2,1,0.05,2,0.1]
## Derived Parameters
epsilon = 2*amp/H;
k = 2*pi/L;
omega = np.sqrt(g*k*np.tanh(k*H));
dt_num = round(CFL*dx_num/(epsilon*k*H*g/(2*omega)),2)
x = np.arange(-1,1+dx_num,dx_num)
y = np.arange(-1,1+dx_num,dy_num)
X,Y = np.meshgrid(x,y)
## NOTE THESE ARE N nodes in x and y, not faces!)
Nx = len(x)
Ny = len(y)
Nt = int(floor(2*N_periods*pi/(omega*dt_num)));
Nn = Nx*Ny
## Symbolic Replacements
variables = [(dx,dx_num),(dy,dy_num),(dt,dt_num),(theta,theta_num)]
In [ ]:
In [ ]:
# First Replace Rx, Ry with real variables
# Dividing out by vx, vy as this will be multiplied by the v matricies
Cx,Cy,Cnx,Cny = array_subs([Cx,Cy,Cnx,Cny],[(Rx,Rxs/vx),(Ry,Rys/vy)])
# Now substitute Real Numbers
Cx,Cy,Cnx,Cny = array_subs([Cx,Cy,Cnx,Cny],variables)
# Generate coefficient vectors
Ax,Ay,_ = create_diags(Cx=Cx,Cy=Cy,Nx=Ny,Ny=Nx,alpha=alpha,dtype=float)
Bx,By,_ = create_diags(Cx=Cnx,Cy=Cny,Nx=Nx,Ny=Ny,alpha=alphan,dtype=float)
# Cleanup Symbolic Names
dx,dy,dt,theta = dx_num,dy_num,dt_num,theta_num
#alpha_sym = alpha
In [ ]:
## Construct Sparse Arrays
#from scipy.sparse import spdiags
#Ax = spdiags([x*np.ones(Nn) for x in Cx], offsets, Nn, Nn).toarray()
#Ay = spdiags([x*np.ones(Nn) for x in Cy], offsets, Nn, Nn).toarray()
#Bx = spdiags([x*np.ones(Nn) for x in Cnx], offsets, Nn, Nn).toarray()
#By = spdiags([x*np.ones(Nn) for x in Cny], offsets, Nn, Nn).toarray()
In [ ]:
## Create X,Y Arrays with shape Nx,Ny,Nt+1 ie. X,Y, time -> to n+1
Xt = np.tile(X.T,(Nt+1,1,1)).T
Yt = np.tile(Y.T,(Nt+1,1,1)).T
## Create X,Y Velocity Arrays with shape Nx,Ny,Nt
# NOTE n is 0 to Nt+1 (This is lets us pick the next time step, n+1)
n = np.ones((Nx,Ny,Nt+1))*np.arange(0,Nt+1)
VX = (epsilon*k*H*g/(2*omega))*sin(k*Xt)*cos(omega*n*dt_num)*cosh(k*(Yt+H))/cosh(k*H)
VY = -(epsilon*k*H*g/(2*omega))*cos(k*Xt)*cos(omega*n*dt_num)*sinh(k*(Yt+H))/cosh(k*H)
# Reshape to N^2 X Nt
Vxf = VX.reshape((Nx*Ny,Nt+1))
Vyf = VY.reshape((Nx*Ny,Nt+1))
In [ ]:
## Initialize Alpha
alpha = np.zeros(Yt.shape)
alpha[np.where(Y<=0)] = 1
In [ ]:
In [ ]:
for t in tqdm(range(0,Nt)):
A = (Ax.T*Vxf[:,t+1] + Ay.T*Vyf[:,t+1]).T + np.eye(len(Ax))
B = ((Bx.T*Vxf[:,t+1] + By.T*Vyf[:,t+1]).T + np.eye(len(Ax)))*alpha[:,:,t].flatten()
alpha[:,:,t+1] = normalize(
np.sum(
solve_banded(
(Nx,Ny),
dbf(A,Nx),
B,
check_finite=False,
overwrite_ab=True,
overwrite_b=True,
),
axis=1).reshape(Nx,Ny)
);
Basic Post Processing¶
In [ ]:
from matplotlib import pyplot as plt
import matplotlib
%matplotlib qt
## Animation Settings
cont_args = dict(cmap='Blues',vmin=0,vmax=1,levels=np.linspace(0,1+1/50,50))
quiver_args = dict(pivot='tail', color='k', units='inches',label='quiver', scale=7.5)
fps = 10
## Plotting Update Function
def update_plot(n, data, plot):
ax.collections=[col for col in ax.collections if isinstance(col,matplotlib.quiver.Quiver)]
ax.collections[0].zorder = 2
plot[0] = ax.contourf(X, Y, data[:,:,n],**cont_args);
#plot[0] = ax.pcolor(X, Y, data[:,:,n],**cont_args);
plot[1].set_UVC(VX[:,:,n],-VY[:,:,n]);
plot[2].set_text('Flow Time: {:0>1.3f}s'.format(n*dt));
return plot
## Initialize the Figure and Axes
fig = plt.figure()
cax = fig.add_axes([0.25, 0.95, 0.5, 0.025])
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
# lft bot wth hht
# Custom Color Bar
cax.pcolor(np.tile(np.arange(0,1.1,0.1),(2,1)),
[0,1],
np.tile(np.arange(0,1.1,0.1),(2,1)),
cmap=cont_args['cmap'],edgecolors='k',linewidths=2)
cax.get_yaxis().set_visible(False)
## Initialize the Plot
plot = [
ax.contourf(X,Y,alpha[:,:,0],**cont_args),
#ax.pcolor(X,Y,alpha[:,:,0],**cont_args),
ax.quiver(X,Y,VX[:,:,0],VY[:,:,0], **quiver_args),
ax.text(0.5,0.1, ' ', bbox={'facecolor':'w', 'alpha':0.75, 'pad':5},transform=ax.transAxes, ha="center")
]
## Save Plots as Images Rather than Animation
import os
# Cleanup old images
[os.remove('images/'+filename) for
filename in os.listdir('images') if '.png' in filename];
for n in tqdm(range(0,Nt)):
plot = update_plot(n,alpha,plot)
filename = 'images/vf-{:0>3d}.png'.format(n)
plt.savefig(filename,dpi=300)
while not os.path.exists(filename):
pass
## Write the images to an mp4
from subprocess import Popen, PIPE
if 'vof-transport.mp4' in os.listdir():
os.remove('vof-transport.mp4')
cmd = 'ffmpeg -framerate 31 -i images/vf-%03d.png vof-transport.mp4'
process = Popen(cmd, stdout=PIPE, stderr=PIPE)
stdout, stderr = process.communicate()
Sample Images from Output¶
In [20]:
from IPython.display import Image
Image(filename='max.png')
Out[20]:
In [21]:
Image(filename='vf-020.png')
Out[21]:
