eigvfftneuz.ipynb

{
 "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
}