Created
June 1, 2021 13:31
-
-
Save christophernhill/d1f07cec81c19df9903b18461ba38250 to your computer and use it in GitHub Desktop.
Hmm a stretched grid FFT solver - is that possible?
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Setup environment | |
using Pkg | |
Pkg.add("PyPlot") | |
using PyPlot | |
Pkg.add("Printf") | |
using Printf | |
Pkg.add("LinearAlgebra") | |
using LinearAlgebra | |
Pkg.add("FFTW") | |
using FFTW | |
# Simple 1d FFT based solve | |
function mkwaves(N,L) | |
scyc=zeros(N,1); sneu=zeros(N,1); | |
for i in 1:N | |
scyc[i]=(2*sin((i-1)*π/N)/(L/N)).^2 | |
sneu[i]=(2*sin((i-1)*π/(2*(N)))/(L/N)).^2 | |
end | |
return scyc, sneu | |
end | |
## Constant grid spacing, Neumann bc's | |
Nx=10 | |
dxF=ones(Nx) | |
xF=[0 cumsum(dxF)'] | |
xC=(xF[2:end]+xF[1:end-1])*0.5 | |
dxC=zeros(Nx+1) | |
dxC[2:Nx]=(dxF[1:Nx-1]+dxF[2:Nx])*0.5 | |
u=zeros(Nx+1) | |
u[Int(Nx/2)]=1 | |
dudx=(u[2:end]-u[1:end-1])./dxF | |
Lx=Nx | |
sxcyc, sxneu=mkwaves(Nx,Lx) | |
## fF1=FFTW.fft(dudx) | |
fF1=FFTW.r2r(dudx,FFTW.REDFT10) | |
fFi=-fF1 | |
for i=1:Nx | |
if i == 1 | |
fFi[i]=0 | |
else | |
fFi[i]=fFi[i]/sxneu[i] | |
end | |
end | |
## fiF1=FFTW.ifft(fFi) | |
fiF1=FFTW.r2r(fFi,FFTW.REDFT01)/(2*Nx) | |
fiF1=real(fiF1) | |
# | |
# Check the "pressure" gradient corrects the u* flow delta | |
fiF1[2:end]-fiF1[1:end-1] .- u[2:end-1] | |
## Lets try a random u | |
u[2:end-1].=map(x -> rand(), u[2:end-1] ) | |
dudx=(u[2:end]-u[1:end-1])./dxF | |
fF1=FFTW.r2r(dudx,FFTW.REDFT10) | |
fFi=-fF1 | |
for i=1:Nx | |
if i == 1 | |
fFi[i]=0 | |
else | |
fFi[i]=fFi[i]/sxneu[i] | |
end | |
end | |
## fiF1=FFTW.ifft(fFi) | |
fiF1=FFTW.r2r(fFi,FFTW.REDFT01)/(2*Nx) | |
fiF1=real(fiF1) | |
## Now lets try making spacing uneven | |
## As before all we need is to find a field gradient that will balance the | |
## divergence in each cell, and in particular we don't strictly need to know | |
## compute the field in the same physical space as long as we can compute a | |
## gradient that will be valid in the physical space. | |
## | |
## Accordingly if we formulate the divergence/convergence in an integral form | |
## as a volume flux, then we can solve for a field gradient that will ensure | |
## the volume flux integrated over each cell is zero as if in a strictly regular | |
## domain. This same field gradient can then also be used to correct the flux in | |
## an non-uniform domain. I think! | |
imid=Int(Nx/2) | |
dxF=ones(Nx) | |
dxF[imid ]=2 | |
dxF[imid-1]=0.5 | |
dxF[imid+1]=0.5 | |
xF=[0 cumsum(dxF)'] | |
xC=(xF[2:end]+xF[1:end-1])*0.5 | |
dxC=zeros(Nx+1) | |
dxC[2:Nx]=(dxF[1:Nx-1]+dxF[2:Nx])*0.5 | |
u=zeros(Nx+1) | |
u[Int(Nx/2)]=1 | |
### dudx=(u[2:end]-u[1:end-1])./dxF | |
dudx=(u[2:end]-u[1:end-1]) | |
Lx=Nx | |
sxcyc, sxneu=mkwaves(Nx,Lx) | |
## fF1=FFTW.fft(dudx) | |
fF1=FFTW.r2r(dudx,FFTW.REDFT10) | |
fFi=-fF1 | |
for i=1:Nx | |
if i == 1 | |
fFi[i]=0 | |
else | |
fFi[i]=fFi[i]/sxneu[i] | |
end | |
end | |
## fiF1=FFTW.ifft(fFi) | |
fiF1=FFTW.r2r(fFi,FFTW.REDFT01)/(2*Nx) | |
fiF1=real(fiF1) | |
fiF1[2:end]-fiF1[1:end-1] .- u[2:end-1] | |
## Hmmm - seems to work! Thats interesting. | |
# | |
# What about a simple 2d problem e.g. face of a cube | |
# | |
# Lets do a Nx x Ny regular grid first with Neuman bc's | |
Nx=32 | |
Ny=32 | |
dxF=ones(Nx) | |
dyF=ones(Ny) | |
xF=[0 cumsum(dxF)'] | |
xC=(xF[2:end]+xF[1:end-1])*0.5 | |
dxC=zeros(Nx+1) | |
dxC[2:Nx]=(dxF[1:Nx-1]+dxF[2:Nx])*0.5 | |
yF=[0 cumsum(dyF)'] | |
yC=(yF[2:end]+yF[1:end-1])*0.5 | |
dyC=zeros(Ny+1) | |
dyC[2:Ny]=(dyF[1:Ny-1]+dyF[2:Ny])*0.5 | |
u=zeros(Nx+1,Ny ) | |
imid=Int(ceil(Nx/2)) | |
jmid=Int(ceil(Ny/2)) | |
u[imid ,jmid]=1 | |
## u[imid+1,jmid]=-1 | |
v=zeros(Nx ,Ny+1) | |
v[imid,jmid ]=1 | |
## v[imid,jmid+1]=-1 | |
dudi=(u[2:end,:]-u[1:end-1,:]) | |
dvdj=(v[:,2:end]-v[:,1:end-1]) | |
divU=-dudi-dvdj | |
Lx=Nx | |
Ly=Ny | |
sxcyc, sxneu=mkwaves(Nx,Lx) | |
sycyc, syneu=mkwaves(Ny,Ly) | |
fF1=FFTW.r2r(divU,FFTW.REDFT10,1) | |
fF2=FFTW.r2r(fF1 ,FFTW.REDFT10,2) | |
fFi=-fF2 | |
for i=1:Nx | |
for j=1:Ny | |
if i == 1 && j == 1 | |
fFi[i,j]=0 | |
else | |
fFi[i,j]=fFi[i,j]/(sxneu[i]+syneu[j]) | |
end | |
end | |
end | |
fiF2=FFTW.r2r(fFi ,FFTW.REDFT01,2)/(2*Ny) | |
fiF1=FFTW.r2r(fiF2,FFTW.REDFT01,1)/(2*Nx) | |
dpi=zeros(Nx+1,Ny ) | |
dpj=zeros(Nx ,Ny+1) | |
dpi[2:Nx,:]=fiF1[2:Nx,:]-fiF1[1:Nx-1,:] | |
dpj[:,2:Ny]=fiF1[:,2:Ny]-fiF1[:,1:Ny-1] | |
unp1=u+dpi | |
vnp1=v+dpj | |
dunp1di=(unp1[2:end,:]-unp1[1:end-1,:]) | |
dvnp1dj=(vnp1[:,2:end]-vnp1[:,1:end-1]) | |
divUNP1=-dunp1di-dvnp1dj |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment