Created
December 2, 2018 19:31
-
-
Save christophernhill/5d9c9aafe3fab413ded728da034cb01e to your computer and use it in GitHub Desktop.
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
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 47, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m registry at `~/.julia/registries/General`\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m git-repo `https://github.com/JuliaRegistries/General.git`\n", | |
"\u001b[2K\u001b[?25h\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Project.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Manifest.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n", | |
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Project.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Manifest.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n", | |
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Project.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Manifest.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n", | |
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Project.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n", | |
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.0/Manifest.toml`\n", | |
"\u001b[90m [no changes]\u001b[39m\n" | |
] | |
} | |
], | |
"source": [ | |
"# Setup environment\n", | |
"using Pkg\n", | |
"Pkg.add(\"PyPlot\")\n", | |
"using PyPlot\n", | |
"Pkg.add(\"Printf\")\n", | |
"using Printf\n", | |
"Pkg.add(\"LinearAlgebra\")\n", | |
"using LinearAlgebra\n", | |
"Pkg.add(\"FFTW\")\n", | |
"using FFTW\n", | |
"tol=1e-12;" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 34, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Make 3-d PPP operator \n", | |
"function mkA_PPP(Nx,Ny,Nz)\n", | |
" NN=Nx*Ny*Nz;\n", | |
" A=zeros(NN,NN)\n", | |
" # Modulo and offset for 1 based index\n", | |
" MOD(i,n)=mod(i-1,n)+1\n", | |
" OFF(i,j,k,ni,nj,nk)= (k-1)*ni*nj + (j-1)*ni + (i-1) + 1 \n", | |
" for k=1:Nz\n", | |
" for j=1:Ny\n", | |
" for i=1:Nx\n", | |
" ic=i ; iw=MOD(i-1,Nx); ie=MOD(i+1,Nx)\n", | |
" jc=j ; js=MOD(j-1,Ny); jn=MOD(j+1,Ny)\n", | |
" kc=k ; ku=MOD(k-1,Nz); kd=MOD(k+1,Nz)\n", | |
" offc=OFF(i , j, k, Nx, Ny, Nz)\n", | |
" offw=OFF(iw, j, k, Nx, Ny, Nz)\n", | |
" offe=OFF(ie, j, k, Nx, Ny, Nz)\n", | |
" offs=OFF( i,js, k, Nx, Ny, Nz)\n", | |
" offn=OFF( i,jn, k, Nx, Ny, Nz)\n", | |
" offu=OFF( i, j,ku, Nx, Ny, Nz)\n", | |
" offd=OFF( i, j,kd, Nx, Ny, Nz)\n", | |
" A[offc,offc]=-6\n", | |
" A[offc,offw]= A[offc,offw]+1\n", | |
" A[offc,offe]= A[offc,offe]+1\n", | |
" A[offc,offs]= A[offc,offs]+1\n", | |
" A[offc,offn]= A[offc,offn]+1\n", | |
" A[offc,offu]= A[offc,offu]+1\n", | |
" A[offc,offd]= A[offc,offd]+1 \n", | |
" end\n", | |
" end\n", | |
" end\n", | |
"\n", | |
" # show(IOContext(stdout), \"text/plain\", Matrix(A))\n", | |
" return A, Nx, Ny, Nz, NN\n", | |
"end\n", | |
"APPP,Nx,Ny,Nz,NN=mkA_PPP(3,3,3);\n", | |
"\n", | |
"A=APPP;\n", | |
"show(IOContext(stdout), \"text/plain\", Int.(Matrix(A)))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 87, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"2×2 Array{Int64,2}:\n", | |
" -1 1\n", | |
" 1 -1" | |
] | |
} | |
], | |
"source": [ | |
"# Make 3-d PPP operator \n", | |
"function mkA_PPN(Nx,Ny,Nz)\n", | |
" NN=Nx*Ny*Nz;\n", | |
" A=zeros(NN,NN)\n", | |
" # Modulo and offset for 1 based index\n", | |
" MOD(i,n)=mod(i-1,n)+1\n", | |
" OFF(i,j,k,ni,nj,nk)= (k-1)*ni*nj + (j-1)*ni + (i-1) + 1 \n", | |
" for k=1:Nz\n", | |
" for j=1:Ny\n", | |
" for i=1:Nx\n", | |
" ic=i ; iw=MOD(i-1,Nx); ie=MOD(i+1,Nx)\n", | |
" jc=j ; js=MOD(j-1,Ny); jn=MOD(j+1,Ny)\n", | |
" kc=k ; ku=MOD(k-1,Nz); kd=MOD(k+1,Nz)\n", | |
" offc=OFF(i , j, k, Nx, Ny, Nz)\n", | |
" offw=OFF(iw, j, k, Nx, Ny, Nz)\n", | |
" offe=OFF(ie, j, k, Nx, Ny, Nz)\n", | |
" offs=OFF( i,js, k, Nx, Ny, Nz)\n", | |
" offn=OFF( i,jn, k, Nx, Ny, Nz)\n", | |
" offu=OFF( i, j,ku, Nx, Ny, Nz)\n", | |
" offd=OFF( i, j,kd, Nx, Ny, Nz)\n", | |
" A[offc,offc]=-6\n", | |
" A[offc,offw]= A[offc,offw]+1\n", | |
" A[offc,offe]= A[offc,offe]+1\n", | |
" A[offc,offs]= A[offc,offs]+1\n", | |
" A[offc,offn]= A[offc,offn]+1\n", | |
" if k == 1\n", | |
" A[offc,offu]= A[offc,offu]+0\n", | |
" A[offc,offc]= A[offc,offc]+1\n", | |
" A[offc,offd]= A[offc,offd]+1 \n", | |
" elseif k == Nz\n", | |
" A[offc,offu]= A[offc,offu]+1\n", | |
" A[offc,offc]= A[offc,offc]+1\n", | |
" A[offc,offd]= A[offc,offd]+0 \n", | |
" else\n", | |
" A[offc,offu]= A[offc,offu]+1\n", | |
" A[offc,offd]= A[offc,offd]+1 \n", | |
" end\n", | |
" end\n", | |
" end\n", | |
" end\n", | |
"\n", | |
" # show(IOContext(stdout), \"text/plain\", Matrix(A))\n", | |
" return A, Nx, Ny, Nz, NN\n", | |
"end\n", | |
"APPN,Nx,Ny,Nz,NN=mkA_PPN(1,1,2);\n", | |
"\n", | |
"A=APPN;\n", | |
"show(IOContext(stdout), \"text/plain\", Int.(Matrix(A)))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 88, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"eigen values \n", | |
"2-element Array{Float64,1}:\n", | |
" -2.0\n", | |
" 0.0\n", | |
"\n", | |
"eigen vectors \n", | |
"2×2 Array{Float64,2}:\n", | |
" 0.707107 0.707107\n", | |
" -0.707107 0.707107\n", | |
"\n" | |
] | |
} | |
], | |
"source": [ | |
"# Show eigen vectors\n", | |
"E=eigen(A);El=E.values;Ev=E.vectors\n", | |
"println(\"eigen values \")\n", | |
"show(IOContext(stdout), \"text/plain\", El )\n", | |
"println(\"\\n\")\n", | |
"println(\"eigen vectors \")\n", | |
"show(IOContext(stdout), \"text/plain\", Matrix(Ev) )\n", | |
"println(\"\\n\")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 90, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"1×1×2 Array{Float64,3}:\n", | |
"[:, :, 1] =\n", | |
" 0.7071067811865475\n", | |
"\n", | |
"[:, :, 2] =\n", | |
" -0.7071067811865475" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
}, | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
" f = [0.707107, -0.707107]\n", | |
" L2(f)0.9999999999999998\n" | |
] | |
} | |
], | |
"source": [ | |
"# Create a RHS\n", | |
"N=size(A)[1]\n", | |
"f=rand(N,1);\n", | |
"# f=ones(N,1);\n", | |
"# f=[1/3 1/3 1/3]';\n", | |
"f=E.vectors[:,1];\n", | |
"display(reshape(f,Nx,Ny,Nz))\n", | |
"\n", | |
"f=f.-sum(f)/N;\n", | |
"println(\" f = \",f)\n", | |
"println(\" L2(f)\", sum(f.*f))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 91, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"1×1×2 Array{Float64,3}:\n", | |
"[:, :, 1] =\n", | |
" 0.0\n", | |
"\n", | |
"[:, :, 2] =\n", | |
" 1.9999999999999996" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"# Show FFT wave number amplitudes\n", | |
"fND=reshape(f,Nx,Ny,Nz)\n", | |
"ff=FFTW.fft(fND)\n", | |
"# println( ff )\n", | |
"display( real.(ff .* conj.(ff)) )" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 92, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"f projected onto eigen vectors[1.0, 0.0]\n", | |
"ϕ projected onto eigen vectors[-0.5, 0.0]\n", | |
"ϕ (eigen vectors projected onto Φ)[-0.353553, 0.353553]\n", | |
"[0.707107, -0.707107]\n", | |
"[0.707107, -0.707107]\n", | |
"[-1.11022e-16, 1.11022e-16]\n" | |
] | |
} | |
], | |
"source": [ | |
"# Solve Eigen style\n", | |
"F=(E.vectors)'*f;\n", | |
"println(\"f projected onto eigen vectors\",F)\n", | |
"λ=E.values;\n", | |
"rλ=map(x -> if (abs(x)>tol) 1.0/x; else 0. ; end , λ);\n", | |
"Φ=F.*rλ;\n", | |
"println(\"ϕ projected onto eigen vectors\",Φ)\n", | |
"ϕ=(E.vectors)*Φ\n", | |
"println(\"ϕ (eigen vectors projected onto Φ)\",ϕ)\n", | |
"println(A*ϕ)\n", | |
"println(f)\n", | |
"println(A*ϕ-f)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 93, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"[0.707107, -0.707107]\n", | |
"[0.707107, -0.707107]\n", | |
"[1.0, 1.0]\n" | |
] | |
} | |
], | |
"source": [ | |
"# Solve FFT style\n", | |
"function mkwaves(N,L)\n", | |
" scyc=zeros(N,1); sneu=zeros(N,1);\n", | |
" for i in 1:N\n", | |
" scyc[i]=(2*sin((i-1)*π/N)/(L/N)).^2\n", | |
" sneu[i]=(2*sin((i-1)*π/(2*(N)))/(L/N)).^2 \n", | |
" end \n", | |
" return scyc, sneu\n", | |
"end\n", | |
"\n", | |
"Lx=Nx;\n", | |
"Ly=Ny;\n", | |
"Lz=Nz;\n", | |
"\n", | |
"# fz = FFTW.r2r(f,FFTW.REDFT10,3)\n", | |
"# fxyz= FFTW.fft(fz,[1,2])\n", | |
"\n", | |
"# fF1=FFTW.fft(reshape(f,Nx,Ny,Nz),[1,2]);\n", | |
"# fF2= FFTW.r2r(real.(fF1),FFTW.REDFT10)\n", | |
"\n", | |
"fF1=FFTW.r2r(reshape(f,Nx,Ny,Nz),FFTW.REDFT10,3)\n", | |
"fF2= FFTW.fft(fF1,[1,2])\n", | |
"\n", | |
"sxcyc, sxneu=mkwaves(Nx,Lx)\n", | |
"sycyc, syneu=mkwaves(Ny,Ly)\n", | |
"szcyc, szneu=mkwaves(Nz,Lz)\n", | |
"sx=sxcyc;\n", | |
"sy=sycyc;\n", | |
"sz=szneu;\n", | |
"\n", | |
"# sx=sxneu;\n", | |
"# sy=syneu;\n", | |
"\n", | |
"fFi=-fF2;\n", | |
"\n", | |
"for i=1:Nx\n", | |
" for j=1:Ny\n", | |
" for k=1:Nz\n", | |
" if i == 1 && j == 1 && k == 1\n", | |
" fFi[i,j,k]=0;\n", | |
" else\n", | |
" fFi[i,j,k]=fFi[i,j,k]./(sx[i]+sy[j]+sz[k])\n", | |
" end\n", | |
" end\n", | |
" end\n", | |
"end\n", | |
"\n", | |
"fiF1=FFTW.ifft(fFi,[1,2])\n", | |
"fiF2=FFTW.r2r(real.(fiF1),FFTW.REDFT01,3)/(2*Nz) \n", | |
" \n", | |
"fiF=reshape(fiF2,Nx*Ny*Nz);\n", | |
"println(f)\n", | |
"println(A*fiF)\n", | |
"println(f./(A*fiF))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Julia 1.0.1", | |
"language": "julia", | |
"name": "julia-1.0" | |
}, | |
"language_info": { | |
"file_extension": ".jl", | |
"mimetype": "application/julia", | |
"name": "julia", | |
"version": "1.0.1" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 2 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment