Hmm a stretched grid FFT solver - is that possible?

stretchy-solver.jl

# 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