2D Stokes problem with periodic BC¶
In [1]:
import numpy as np
from numpy.linalg import multi_dot
import matplotlib
import matplotlib.pyplot as plt
In [17]:
N = 2**8
#k = np.append(np.arange(0,N//2+1).tolist(),np.arange(-N//2+1,0).tolist())
#k = [kk + 0j for kk in k]
k = N*np.fft.fftfreq(N)
k1,k2 = np.meshgrid(k,k)
mu = 1.
# Source term f
#x = np.linspace(0,2*np.pi,N)
x = np.linspace(-np.pi,np.pi,N+1)[:-1]
x1, x2 = np.meshgrid(x, x, sparse=False)
# Manufactured solution (validation)
a = 1.; c1 = 1.; c2 = 1.; #velocity coeffs
b = 1.; l1 = 1.; l2 = 1.; #pressure coeffs
um1 = a*np.sin(c1*x1)*np.cos(c2*x2) # vel. x1-component
um2 = -a*c1/c2*np.cos(c1*x1)*np.sin(c2*x2) # vel. x2-component
pm = b*np.sin(l1*x1)*np.cos(l2*x2) # pressure
d1pm = b*l1*np.cos(l1*x1)*np.cos(l2*x2)
d2pm = -b*l2*np.sin(l1*x1)*np.sin(l2*x2)
# forcing
#f1 = mu*(c1**2+c2**2)*um1+d1pm
#f2 = -mu*c1/c2*(c1**2+c2**2)*um2+d2pm
#f1 = np.sin(2*x1)**2*np.cos(x2)**2 + 0*np.random.randn(x1.shape[0],x2.shape[0])
#f2 = 0*np.sin(2*x1)**2*np.cos(x2)**2 + 0*np.random.randn(x1.shape[0],x2.shape[0])
# funny terms (wierd shape)
#f1 = (np.sin((x1-1)**2+x2**2))**2#+1*np.sin(x1*x2)
#f2 = (np.sin((x1+1)**2+x2**2))**2
#f1 = np.random.randn(x1.shape[0],x2.shape[0])+0.2*np.cos(x1**2+x2**2)#-0.075*x2/(x1**2+x2**2+1e-8)**0.2#+0.005*(np.sin((x1-1)**2+x2**2))**2
#f2 = np.random.randn(x1.shape[0],x2.shape[0])#+0.075*x1/(x1**2+x2**2+1e-8)**0.7#+0.01*(np.sin((x1+1)**2+x2**2))**2
f1 = np.random.randn(x1.shape[0],x2.shape[0])
f2 = np.random.randn(x1.shape[0],x2.shape[0])
#f1 = np.random.randn(x1.shape[0],x2.shape[0])-.5*np.cos((x1**2+x2**2))**2
#f2 = np.random.randn(x1.shape[0],x2.shape[0])#0/(x1**2+x2**2+1)
#f1 = np.sin(x1)*np.sin(x2)
#f2 = np.cos(x1)*np.cos(x2)
f1_ = np.fft.fft2(f1)
f2_ = np.fft.fft2(f2)
u_ = np.zeros((2,N,N),dtype=complex)
# for i in range(N):
# for j in range(N):
# k1 = k[i]; k2 = k[j]
# if k1 == 0 and k2 == 0:
# continue
# kn = k1**2+k2**2
# u_[0][i][j] = 1./mu/kn*((1./(2*np.pi)-k1**2/kn)*f1_[i][j]-(k1*k2/kn)*f2_[i][j])
# u_[1][i][j] = 1./mu/kn*((1./(2*np.pi)-k2**2/kn)*f1_[i][j]-(k2*k1/kn)*f2_[i][j])
kk = k1**2+k2**2; print(kk[0,0])
kk[0,0] = 1;
c = 2*np.pi
#u_[0] = 1./mu/kk*((1.-k1*k1/kk)*f1_-k1*k2/kk*f2_)
#u_[1] = 1./mu/kk*((1.-k2*k2/kk)*f1_-k2*k1/kk*f2_)
#u_[0] = 1./mu/kk*((1./c-k1*k1/kk)*f1_-k1*k2/kk*f2_)
#u_[1] = 1./mu/kk*((1./c-k2*k2/kk)*f2_-k2*k1/kk*f1_)
u_[0] = 1./mu/kk**1*(f1_-k1/kk*(f1_*k1+f2_*k2))
u_[1] = 1./mu/kk**1*(f2_-k2/kk*(f1_*k1+f2_*k2))
u_[:][0][0] = 0
u1 = np.fft.ifft2(u_[0])
u2 = np.fft.ifft2(u_[1])
0.0
In [18]:
# Plotting
from matplotlib import cm
#norm = cm.colors.Normalize(vmax=abs(Z).max(), vmin=-abs(Z).max())
cmap = cm.magma; size = 8
fig, ax = plt.subplots(1,2,figsize=(2*size,size+2));
im=ax[0].contourf(x,x,(f1**2+f2**2)**0.5,40,cmap=cm.get_cmap(cmap,20))
fig.colorbar(im,ax=ax[0],orientation="horizontal")
lw = (f1.real**2+f2.real**2)**0.5/np.max((f1**2+f2**2)**0.5)
ax[0].streamplot(x,x,f1.real,f2.real,color='purple')
ax[0].set_xlabel('x1'); ax[0].set_ylabel('x2'); ax[0].set_title('Source term $f$')
#im=ax[1].contourf(x,x,(abs(u1)**2+abs(u2)**2)**0.5-(abs(um1)**2+abs(um2)**2)**0.5,40, cmap=cm.get_cmap(cmap,))
im=ax[1].contourf(x,x,(u1.real**2+u2.real**2)**0.5,40, cmap=cm.get_cmap(cmap,))
# im = ax[1].imshow((u1.real**2+u2.real**2)**0.5, cmap=cm.get_cmap(cmap,),interpolation='bilinear')
lw = (abs(u1)**2+abs(u2)**2)**0.5/np.max((abs(u1)**2+abs(u2)**2)**0.5)
#ax[1].streamplot(x,x,u1.real-um1.real,u2.real-um2.real,color='purple')
#ax[1].streamplot(x,x,u1.real,u2.real,color='cyan',density=2)
fig.colorbar(im,ax=ax[1],orientation="horizontal")
ax[1].set_xlabel('x1'); ax[1].set_ylabel('x2'); ax[1].set_title('Stokes solution (or not)')
Out[18]:
Text(0.5, 1.0, 'Stokes solution (or not)')
In [ ]:
In [ ]: