diff --git a/.gitignore b/.gitignore index 381e0b6..e789b5e 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ *.jl.*.cov *.jl.mem deps/deps.jl +.ipynb_checkpoints/ +*.swp diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..dd71e2e --- /dev/null +++ b/Dockerfile @@ -0,0 +1,41 @@ +FROM jupyter/minimal-notebook + +LABEL maintainer="Pierre Navaro " + +USER root + +# Julia dependencies +# install Julia packages in /opt/julia instead of $HOME +ENV JULIA_PKGDIR=/opt/julia +ENV JULIA_VERSION=1.0.2 + +RUN mkdir /opt/julia-${JULIA_VERSION} && \ + cd /tmp && \ + wget -q https://julialang-s3.julialang.org/bin/linux/x64/`echo ${JULIA_VERSION} | cut -d. -f 1,2`/julia-${JULIA_VERSION}-linux-x86_64.tar.gz && \ + tar xzf julia-${JULIA_VERSION}-linux-x86_64.tar.gz -C /opt/julia-${JULIA_VERSION} --strip-components=1 && \ + rm /tmp/julia-${JULIA_VERSION}-linux-x86_64.tar.gz +RUN ln -fs /opt/julia-*/bin/julia /usr/local/bin/julia + +# Show Julia where conda libraries are \ +RUN mkdir /etc/julia && \ + echo "push!(Libdl.DL_LOAD_PATH, \"$CONDA_DIR/lib\")" >> /etc/julia/juliarc.jl && \ + # Create JULIA_PKGDIR \ + mkdir $JULIA_PKGDIR && \ + chown $NB_USER $JULIA_PKGDIR && \ + fix-permissions $JULIA_PKGDIR + +COPY . ${HOME} +RUN chown -R ${NB_UID} ${HOME} +USER $NB_UID + +RUN julia -e 'import Pkg; Pkg.update()' && \ + julia -e 'import Pkg; Pkg.add("IJulia")' && \ + # Precompile Julia packages \ + julia -e 'using IJulia' && \ + # move kernelspec out of home \ + mv $HOME/.local/share/jupyter/kernels/julia* $CONDA_DIR/share/jupyter/kernels/ && \ + chmod -R go+rx $CONDA_DIR/share/jupyter && \ + rm -rf $HOME/.local && \ + fix-permissions $JULIA_PKGDIR $CONDA_DIR/share/jupyter + +RUN cd ${HOME} && julia -q ${HOME}/deps_install.jl diff --git a/README.md b/README.md index 6ddb16c..78c0bdc 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,5 @@ # VlasovNotebooks + Notebooks with some examples + +[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/juliavlasov/VlasovNotebooks/master) diff --git a/REQUIRE b/REQUIRE new file mode 100644 index 0000000..b4d971e --- /dev/null +++ b/REQUIRE @@ -0,0 +1,6 @@ +BenchmarkTools +FFTW +Plots +ProgressMeter +QuadGK +FINUFFT diff --git a/Vlasov.Ampere.Fourier.ipynb b/Vlasov.Ampere.Fourier.ipynb new file mode 100644 index 0000000..c25cc7a --- /dev/null +++ b/Vlasov.Ampere.Fourier.ipynb @@ -0,0 +1,1025 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1D Vlasov–Ampere system\n", + "\n", + "$$\n", + "\\frac{\\partial f}{\\partial t} + \\upsilon \\frac{\\partial f}{\\partial x}\n", + "- E(t,x) \\frac{\\partial f}{\\partial \\upsilon} = 0\n", + "$$\n", + "\n", + "$$\n", + "\\frac{\\partial E}{\\partial t} = - J = \\int f\\upsilon \\; d\\upsilon\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Algorithm \n", + "\n", + "- For each $j$ compute discrete Fourier transform in $x$ of $f^n(x_i,\\upsilon_j)$ yielding $f_k^n(\\upsilon_j)$, \n", + "\n", + "- For $ k \\neq 0 $\n", + "\n", + " - Compute \n", + " \n", + " $$f^{n+1}_k(\\upsilon_j) = e^{−2i\\pi k \\upsilon\n", + " \\Delta t/L} f_n^k(\\upsilon_j),$$\n", + " \n", + " - Compute \n", + " \n", + " $$\\rho_k^{n+1} = \\Delta \\upsilon \\sum_j􏰄 f^{n+1}_k(\\upsilon_j),$$\n", + " \n", + " - Compute\n", + " \n", + " $$E^{n+1}_k = \\rho^{n+1}_k L/(2i\\pi k \\epsilon_0),$$\n", + " \n", + "- For $k = 0$ do nothing: \n", + "\n", + "$$f_{n+1}(\\upsilon_j) = f^n_k(\\upsilon_j), E^{n+1}_k = E^n_k$$.\n", + "\n", + "- Perform inverse discrete Fourier transform of $E^{n+1}_k$ and for each $j$ of $f^{n+1}_k (\\upsilon_j)$." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "using ProgressMeter, FFTW, Plots, LinearAlgebra\n", + "using BenchmarkTools, Statistics" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "UniformMesh" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + "\n", + " UniformMesh(start, stop, length)\n", + "\n", + " 1D uniform mesh data.\n", + "\n", + " length : number of points\n", + " length-1 : number of cells\n", + "\n", + " To remove the last point, set endpoint=false\n", + "\n", + "\"\"\"\n", + "struct UniformMesh\n", + "\n", + " start :: Float64\n", + " stop :: Float64\n", + " length :: Int\n", + " step :: Float64\n", + " points :: Vector{Float64}\n", + " endpoint :: Bool\n", + "\n", + " function UniformMesh(start, stop, length::Int; endpoint=true)\n", + "\n", + " if (endpoint)\n", + " points = range(start, stop=stop, length=length)\n", + " else\n", + " points = range(start, stop=stop, length=length+1)[1:end-1]\n", + " end\n", + " step = points[2]-points[1]\n", + "\n", + " new( start, stop, length, step, points, endpoint)\n", + "\n", + " end\n", + "\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "using FFTW, LinearAlgebra\n", + "\n", + "abstract type AbstractAdvection end\n", + "abstract type AdvectionType end" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [], + "source": [ + "import FFTW: fft!, ifft!\n", + "import LinearAlgebra: transpose\n", + "\n", + "export Ampere\n", + "\n", + "struct Ampere end\n", + "\n", + "\"\"\"\n", + "\n", + " advection!( fᵀ, mesh1, mesh2, E, dt, type, axis ) \n", + "\n", + " ∂f/∂t − υ∂f/∂x = 0\n", + " ∂E/∂t = −J = ∫ fυ dυ\n", + " ∂f/∂t − E(x) ∂f/∂υ = 0\n", + "\n", + "\"\"\"\n", + "struct AmpereAdvection <: AbstractAdvection\n", + " \n", + " mesh :: UniformMesh\n", + " kx :: Vector{Float64}\n", + "\n", + " function AmpereAdvection( mesh )\n", + " \n", + " nx = mesh.length\n", + " dx = mesh.step\n", + " Lx = mesh.stop - mesh.start\n", + " kx = zeros(Float64, nx)\n", + " kx .= 2π/Lx * [0:nx÷2-1;-nx÷2:-1]\n", + " new( mesh, kx)\n", + " \n", + " end\n", + "\n", + "end\n", + "\n", + "function (adv :: AmpereAdvection)( fᵗ :: Array{ComplexF64,2}, \n", + "\t\t e :: Vector{ComplexF64}, \n", + "\t\t dt :: Float64 )\n", + " fft!(fᵗ, 1)\n", + " fᵗ .= fᵗ .* exp.(-1im * dt * adv.kx * transpose(e))\n", + " ifft!(fᵗ, 1)\n", + "\n", + "end\n", + "\n", + "function (adv :: AmpereAdvection)( f :: Array{ComplexF64,2}, \n", + "\t\t e :: Vector{ComplexF64}, \n", + "\t\t v :: Vector{Float64}, \n", + "\t\t dt :: Float64 )\n", + " \n", + " ev = exp.(-1im*dt * adv.kx * transpose(v)) \n", + " \n", + " fft!(f,1)\n", + " f .= f .* ev\n", + " dv = v[2]-v[1]\n", + " ρ = dv * vec(sum(f,dims=2)) \n", + " for i in 2:length(e)\n", + " e[i] = -1im * ρ[i] ./ adv.kx[i]\n", + " end\n", + " e[1] = 0.0\n", + " ifft!(f,1)\n", + " ifft!(e)\n", + " e .= real(e)\n", + "end\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "BSplineAdvection" + ] + }, + "execution_count": 42, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using FFTW, LinearAlgebra\n", + "\n", + "\"\"\"\n", + " bspline(p, j, x)\n", + "\n", + "Return the value at x in [0,1[ of the B-spline with integer nodes of degree p with support starting at j.\n", + "Implemented recursively using the [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)\n", + "\n", + "```math\n", + "B_{i,0}(x) := \\\\left\\\\{\n", + "\\\\begin{matrix}\n", + "1 & \\\\mathrm{if} \\\\quad t_i ≤ x < t_{i+1} \\\\\\\\\n", + "0 & \\\\mathrm{otherwise} \n", + "\\\\end{matrix}\n", + "\\\\right.\n", + "```\n", + "\n", + "```math\n", + "B_{i,p}(x) := \\\\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x) \n", + "+ \\\\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).\n", + "```\n", + "\"\"\"\n", + "function bspline(p::Int, j::Int, x::Float64)\n", + " if p == 0\n", + " if j == 0\n", + " return 1.0\n", + " else\n", + " return 0.0\n", + " end\n", + " else\n", + " w = (x - j) / p\n", + " w1 = (x - j - 1) / p\n", + " end\n", + " return (w * bspline(p - 1, j, x) + (1 - w1) * bspline(p - 1, j + 1, x))\n", + "end\n", + "\n", + "export BSpline\n", + "\n", + "struct BSpline <: AdvectionType\n", + " p :: Int64\n", + "end\n", + "\n", + "\n", + "\"\"\"\n", + " BsplineAdvection(p, mesh)\n", + "\n", + "Advection type\n", + "\n", + "\"\"\"\n", + "mutable struct BSplineAdvection <: AbstractAdvection\n", + " \n", + " p :: Int64 \n", + " mesh :: UniformMesh\n", + " modes :: Vector{Float64}\n", + " eig_bspl :: Vector{Float64}\n", + " eigalpha :: Vector{Complex{Float64}}\n", + " \n", + " function BSplineAdvection( p, mesh )\n", + "\n", + " println(\" Create BSL advection with bspline of degree $p \")\n", + " @show nx = mesh.length\n", + "\n", + " modes = zeros(Float64, nx)\n", + " modes .= [2π * i / nx for i in 0:nx-1]\n", + " eig_bspl = zeros(Float64, nx)\n", + " eig_bspl = zeros(Float64, nx)\n", + " eig_bspl .= bspline(p, -div(p+1,2), 0.0)\n", + " for i in 1:div(p+1,2)-1\n", + " eig_bspl .+= bspline(p, i-(p+1)÷2, 0.0) * 2 .* cos.(i * modes)\n", + " end\n", + " eigalpha = zeros(Complex{Float64}, nx)\n", + " new( p, mesh, modes, eig_bspl, eigalpha )\n", + " end\n", + " \n", + "end\n", + "\n", + "\"\"\"\n", + "\n", + " advection( f, v, dt)\n", + "\n", + "Instantiate an advection type \n", + "\n", + "```julia\n", + "advection! = Advection( p, mesh)\n", + "\n", + "advection!(f, v, dt)\n", + "```\n", + "\n", + "\"\"\"\n", + "function (adv :: BSplineAdvection)(f :: Array{Complex{Float64},2}, \n", + " v :: Vector{Float64}, \n", + " dt :: Float64)\n", + " \n", + " nx = adv.mesh.length\n", + " nv = length(v)\n", + " dx = adv.mesh.step\n", + " \n", + " fft!(f,1)\n", + " \n", + " for j in 1:nv\n", + "\n", + " @inbounds alpha = dt * v[j] / dx\n", + " # compute eigenvalues of cubic splines evaluated at displaced points\n", + " ishift = floor(-alpha)\n", + " beta = -ishift - alpha\n", + " fill!(adv.eigalpha,0.0im)\n", + " for i in -div(adv.p-1,2):div(adv.p+1,2)\n", + " adv.eigalpha .+= (bspline(adv.p, i-div(adv.p+1,2), beta) \n", + " .* exp.((ishift+i) * 1im .* adv.modes))\n", + " end\n", + " \n", + " # compute interpolating spline using fft and properties of circulant matrices\n", + " \n", + " @inbounds f[:,j] .*= adv.eigalpha ./ adv.eig_bspl\n", + " \n", + " end\n", + " \n", + " ifft!(f,1)\n", + " \n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Advection (generic function with 2 methods)" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function Advection( b::BSpline, mesh:: UniformMesh )\n", + "\n", + " BSplineAdvection( b.p, mesh)\n", + "\n", + "end\n", + "\n", + "function Advection( b::Ampere, mesh:: UniformMesh )\n", + "\n", + " AmpereAdvection(mesh)\n", + "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "compute_rho" + ] + }, + "execution_count": 44, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + "\n", + " compute_rho( mesh, f)\n", + "\n", + " Compute charge density\n", + "\n", + " ρ(x,t) = ∫ f(x,v,t) delta2\n", + "\n", + " return ρ - ρ̄ if neutralized=true\n", + "\n", + "\"\"\"\n", + "function compute_rho(meshv::UniformMesh, f, neutralized=true)\n", + "\n", + " local dv = meshv.step\n", + " ρ = dv * vec(sum(real(f), dims=2))\n", + " if (neutralized)\n", + " ρ .- mean(ρ)\n", + " else\n", + " ρ\n", + " end\n", + "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "compute_e (generic function with 1 method)" + ] + }, + "execution_count": 45, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function compute_e(mesh::UniformMesh, ρ)\n", + "\n", + " local n = mesh.length\n", + " local k = 2π / (mesh.stop - mesh.start)\n", + " local modes = zeros(Float64, n)\n", + " modes .= k * vcat(0:n÷2-1,-n÷2:-1)\n", + " modes[1] = 1.0\n", + " ρ̂ = fft(ρ)./modes\n", + " vec(real(ifft(-1im*ρ̂)))\n", + "\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "landau (generic function with 1 method)" + ] + }, + "execution_count": 46, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function landau( ϵ, kx, x, v )\n", + " \n", + " (1.0.+ϵ*cos.(kx*x))/sqrt(2π) .* transpose(exp.(-0.5*v.*v))\n", + " \n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "vlasov_ampere (generic function with 1 method)" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function vlasov_ampere( nx, nv, xmin, xmax, vmin, vmax , tf, nt)\n", + "\n", + " meshx = UniformMesh(xmin, xmax, nx)\n", + " meshv = UniformMesh(vmin, vmax, nv)\n", + " \n", + " # Initialize distribution function\n", + " x = meshx.points\n", + " v = meshv.points\n", + " ϵ, kx = 0.001, 0.5\n", + " \n", + " f = zeros(Complex{Float64},(nx,nv))\n", + " fᵀ= zeros(Complex{Float64},(nv,nx))\n", + " \n", + " f .= landau( ϵ, kx, x, v)\n", + " \n", + " transpose!(fᵀ,f)\n", + " \n", + " ρ = compute_rho(meshv, f)\n", + " e = zeros(ComplexF64, nx)\n", + " e .= compute_e(meshx, ρ)\n", + " \n", + " nrj = Float64[]\n", + " \n", + " dt = tf / nt\n", + " \n", + " advection_x! = Advection( Ampere(), meshx )\n", + " advection_v! = Advection( Ampere(), meshv )\n", + " \n", + " for i in 1:nt\n", + " advection_v!(fᵀ, e, 0.5*dt)\n", + " transpose!(f,fᵀ)\n", + " advection_x!( f, e, v, dt)\n", + " push!(nrj, log(sqrt((sum(e.^2))*meshx.step)))\n", + " transpose!(fᵀ,f)\n", + " advection_v!(fᵀ, e, 0.5*dt)\n", + " end\n", + " nrj\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "vlasov_poisson (generic function with 1 method)" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function vlasov_poisson( nx, nv, xmin, xmax, vmin, vmax , tf, nt)\n", + "\n", + " meshx = UniformMesh(xmin, xmax, nx)\n", + " meshv = UniformMesh(vmin, vmax, nv)\n", + " \n", + " # Initialize distribution function\n", + " x = meshx.points\n", + " v = meshv.points\n", + " ϵ, kx = 0.001, 0.5\n", + " \n", + " e = zeros(Float64, nx)\n", + " ρ = zeros(Float64, nx)\n", + " f = zeros(Complex{Float64},(nx,nv))\n", + " fᵀ= zeros(Complex{Float64},(nv,nx))\n", + " \n", + " f .= landau( ϵ, kx, x, v)\n", + "\n", + " nrj = Float64[]\n", + " \n", + " dt = tf / nt\n", + " \n", + " advection_x! = Advection( BSpline(5), meshx )\n", + " advection_v! = Advection( BSpline(5), meshv )\n", + " \n", + " for i in 1:nt\n", + " advection_x!(f, v, 0.5dt)\n", + " ρ .= compute_rho(meshv, f)\n", + " e .= compute_e(meshx, ρ)\n", + " push!(nrj, log(sqrt((sum(e.^2))*meshx.step)))\n", + " transpose!(fᵀ, f)\n", + " advection_v!(fᵀ, e, dt)\n", + " transpose!(f, fᵀ)\n", + " advection_x!( f, v, 0.5dt)\n", + " end\n", + " nrj\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "600" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nx, nv = 256, 256\n", + "xmin, xmax = 0., 4*π\n", + "vmin, vmax = -6., 6.\n", + "tf = 60\n", + "nt = 600" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Create BSL advection with bspline of degree 5 \n", + "nx = mesh.length = 256\n", + " Create BSL advection with bspline of degree 5 \n", + "nx = mesh.length = 256\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "0\n", + "\n", + "\n", + "10\n", + "\n", + "\n", + "20\n", + "\n", + "\n", + "30\n", + "\n", + "\n", + "40\n", + "\n", + "\n", + "50\n", + "\n", + "\n", + "60\n", + "\n", + "\n", + "-18\n", + "\n", + "\n", + "-16\n", + "\n", + "\n", + "-14\n", + "\n", + "\n", + "-12\n", + "\n", + "\n", + "-10\n", + "\n", + "\n", + "-8\n", + "\n", + "\n", + "-6\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "y1\n", + "\n", + "\n", + "\n", + "ampere\n", + "\n", + "\n", + "\n", + "poisson\n", + "\n", + "\n" + ] + }, + "execution_count": 50, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "t = range(0,stop=tf,length=nt)\n", + "plot(t, -0.1533*t.-5.48)\n", + "plot!(t, vlasov_ampere(nx, nv, xmin, xmax, vmin, vmax, tf, nt), label=:ampere )\n", + "plot!(t, vlasov_poisson(nx, nv, xmin, xmax, vmin, vmax, tf, nt), label=:poisson )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.1.0", + "language": "julia", + "name": "julia-1.1" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.1.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Vlasov.Poisson.BSL.ipynb b/Vlasov.Poisson.BSL.ipynb new file mode 100644 index 0000000..392007f --- /dev/null +++ b/Vlasov.Poisson.BSL.ipynb @@ -0,0 +1,2026 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Semi-Lagrangian method\n", + "\n", + "Let us consider an abstract scalar advection equation of the form\n", + "$$\n", + "\\frac{\\partial f}{\\partial t}+ a(x, t) \\cdot \\nabla f = 0. \n", + "$$\n", + "The characteristic curves associated to this equation are the solutions of the ordinary differential equations\n", + "$$\n", + "\\frac{dX}{dt} = a(X(t), t)\n", + "$$\n", + "We shall denote by $X(t, x, s)$ the unique solution of this equation associated to the initial condition $X(s) = x$.\n", + "\n", + "The classical semi-Lagrangian method is based on a backtracking of characteristics. Two steps are needed to update the distribution function $f^{n+1}$ at $t^{n+1}$ from its value $f^n$ at time $t^n$ :\n", + "1. For each grid point $x_i$ compute $X(t^n; x_i, t^{n+1})$ the value of the characteristic at $t^n$ which takes the value $x_i$ at $t^{n+1}$.\n", + "2. As the distribution solution of first equation verifies \n", + "$$f^{n+1}(x_i) = f^n(X(t^n; x_i, t^{n+1})),$$\n", + "we obtain the desired value of $f^{n+1}(x_i)$ by computing $f^n(X(t^n;x_i,t^{n+1})$ by interpolation as $X(t^n; x_i, t^{n+1})$ is in general not a grid point.\n", + "\n", + "*[Eric Sonnendrücker - Numerical methods for the Vlasov equations](http://www-m16.ma.tum.de/foswiki/pub/M16/Allgemeines/NumMethVlasov/Num-Meth-Vlasov-Notes.pdf)*" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "using Plots, FFTW, LinearAlgebra, Statistics, BenchmarkTools" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Bspline\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "bspline" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + " bspline(p, j, x)\n", + "\n", + "Return the value at x in [0,1[ of the B-spline with integer nodes of degree p with support starting at j.\n", + "Implemented recursively using the [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)\n", + "\n", + "```math\n", + "B_{i,0}(x) := \\\\left\\\\{\n", + "\\\\begin{matrix}\n", + "1 & \\\\mathrm{if} \\\\quad t_i ≤ x < t_{i+1} \\\\\\\\\n", + "0 & \\\\mathrm{otherwise} \n", + "\\\\end{matrix}\n", + "\\\\right.\n", + "```\n", + "\n", + "```math\n", + "B_{i,p}(x) := \\\\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x) \n", + "+ \\\\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).\n", + "```\n", + "\"\"\"\n", + "function bspline(p::Int, j::Int, x::Float64)\n", + " if p == 0\n", + " if j == 0\n", + " return 1.0\n", + " else\n", + " return 0.0\n", + " end\n", + " else\n", + " w = (x - j) / p\n", + " w1 = (x - j - 1) / p\n", + " end\n", + " return (w * bspline(p - 1, j, x) + (1 - w1) * bspline(p - 1, j + 1, x))\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "search: \u001b[0m\u001b[1mb\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1mi\u001b[22m\u001b[0m\u001b[1mn\u001b[22m\u001b[0m\u001b[1me\u001b[22m\n", + "\n" + ] + }, + { + "data": { + "text/latex": [ + "\\begin{verbatim}\n", + "bspline(p, j, x)\n", + "\\end{verbatim}\n", + "Return the value at x in [0,1[ of the B-spline with integer nodes of degree p with support starting at j. Implemented recursively using the \\href{https://en.wikipedia.org/wiki/De_Boor%27s_algorithm}{De Boor's Algorithm}\n", + "\n", + "$$B_{i,0}(x) := \\left\\{\n", + "\\begin{matrix}\n", + "1 & \\mathrm{if} \\quad t_i ≤ x < t_{i+1} \\\\\n", + "0 & \\mathrm{otherwise} \n", + "\\end{matrix}\n", + "\\right.$$\n", + "$$B_{i,p}(x) := \\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x) \n", + "+ \\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).$$\n" + ], + "text/markdown": [ + "```\n", + "bspline(p, j, x)\n", + "```\n", + "\n", + "Return the value at x in [0,1[ of the B-spline with integer nodes of degree p with support starting at j. Implemented recursively using the [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)\n", + "\n", + "$$\n", + "B_{i,0}(x) := \\left\\{\n", + "\\begin{matrix}\n", + "1 & \\mathrm{if} \\quad t_i ≤ x < t_{i+1} \\\\\n", + "0 & \\mathrm{otherwise} \n", + "\\end{matrix}\n", + "\\right.\n", + "$$\n", + "\n", + "$$\n", + "B_{i,p}(x) := \\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x) \n", + "+ \\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).\n", + "$$\n" + ], + "text/plain": [ + "\u001b[36m bspline(p, j, x)\u001b[39m\n", + "\n", + " Return the value at x in [0,1[ of the B-spline with integer nodes of degree\n", + " p with support starting at j. Implemented recursively using the De Boor's\n", + " Algorithm (https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)\n", + "\n", + "\u001b[35mB_{i,0}(x) := \\left\\{\u001b[39m\n", + "\u001b[35m\\begin{matrix}\u001b[39m\n", + "\u001b[35m1 & \\mathrm{if} \\quad t_i ≤ x < t_{i+1} \\\\\u001b[39m\n", + "\u001b[35m0 & \\mathrm{otherwise} \u001b[39m\n", + "\u001b[35m\\end{matrix}\u001b[39m\n", + "\u001b[35m\\right.\u001b[39m\n", + "\n", + "\u001b[35mB_{i,p}(x) := \\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x) \u001b[39m\n", + "\u001b[35m+ \\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).\u001b[39m" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "?bspline" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "UniformMesh" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + " UniformMesh(xmin, xmax, nx)\n", + "\n", + "1D uniform mesh data\n", + "\"\"\"\n", + "struct UniformMesh\n", + " xmin :: Float64\n", + " xmax :: Float64\n", + " nx :: Int\n", + " dx :: Float64\n", + " x :: Vector{Float64}\n", + " function UniformMesh(xmin, xmax, nx)\n", + " dx = (xmax - xmin) / nx\n", + " x = range(xmin, stop=xmax, length=nx+1)[1:end-1] \n", + " new( xmin, xmax, nx, dx, x)\n", + " end\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.6283185307179586, 10)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mesh = UniformMesh(-π, π, 10)\n", + "mesh.dx, mesh.nx" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "advection!" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + " advection!(f, p, mesh, v, dt)\n", + "\n", + "function to advect the distribution function `f` with velocity `v`\n", + "during a time step `dt`. Interpolation method uses bspline periodic.\n", + "\"\"\"\n", + "function advection!(f :: Array{Complex{Float64},2}, \n", + " p :: Int64, \n", + " mesh :: UniformMesh, \n", + " v :: Vector{Float64}, \n", + " dt :: Float64)\n", + " \n", + " nx = mesh.nx\n", + " nv = length(v)\n", + " dx = mesh.dx\n", + " modes = [2π * i / nx for i in 0:nx-1]\n", + " \n", + " # compute eigenvalues of degree p b-spline matrix\n", + " eig_bspl = zeros(Float64, nx)\n", + " eig_bspl .= bspline(p, -div(p+1,2), 0.0)\n", + " for i in 1:div(p+1,2)-1\n", + " eig_bspl .+= bspline(p, i - (p+1)÷2, 0.0) * 2 .* cos.(i * modes)\n", + " end\n", + " eigalpha = zeros(Complex{Float64}, nx)\n", + " \n", + " fft!(f,1)\n", + " \n", + " for j in 1:nv\n", + " @inbounds alpha = dt * v[j] / dx\n", + " # compute eigenvalues of cubic splines evaluated at displaced points\n", + " ishift = floor(-alpha)\n", + " beta = -ishift - alpha\n", + " fill!(eigalpha,0.0im)\n", + " for i in -div(p-1,2):div(p+1,2)\n", + " eigalpha .+= (bspline(p, i-div(p+1,2), beta) \n", + " .* exp.((ishift+i) * 1im .* modes))\n", + " end\n", + " \n", + " # compute interpolating spline using fft and properties of circulant matrices\n", + " \n", + " @inbounds f[:,j] .*= eigalpha ./ eig_bspl\n", + " \n", + " end\n", + " \n", + " ifft!(f,1)\n", + " \n", + "end " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Vlasov-Poisson equation\n", + "We consider the dimensionless Vlasov-Poisson equation for one species\n", + "with a neutralizing background.\n", + "\n", + "$$ \n", + "\\frac{\\partial f}{\\partial t}+ v\\cdot \\nabla_x f + E(t,x) \\cdot \\nabla_v f = 0, \\\\\n", + "- \\Delta \\phi = 1 - \\rho, E = - \\nabla \\phi \\\\\n", + "\\rho(t,x) = \\int f(t,x,v)dv.\n", + "$$\n", + "\n", + "- [Vlasov Equation - Wikipedia](https://en.wikipedia.org/wiki/Vlasov_equation)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "compute_rho" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + " compute_rho(meshv, f)\n", + "\n", + "Compute charge density\n", + "ρ(x,t) = ∫ f(x,v,t) dv\n", + "\"\"\"\n", + "function compute_rho(meshv::UniformMesh, \n", + " f::Array{Complex{Float64},2})\n", + " \n", + " dv = meshv.dx\n", + " rho = dv * sum(real(f), dims=2)\n", + " vec(rho .- mean(rho)) # vec squeezes the 2d array returned by sum function\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "compute_e" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + " compute_e(meshx, rho)\n", + "compute Ex using that -ik*Ex = rho \n", + "\"\"\"\n", + "function compute_e(meshx::UniformMesh, rho::Vector{Float64})\n", + " nx = meshx.nx\n", + " k = 2π / (meshx.xmax - meshx.xmin)\n", + " modes = zeros(Float64, nx)\n", + " modes .= k * vcat(0:div(nx,2)-1,-div(nx,2):-1)\n", + " modes[1] = 1.0\n", + " rhok = fft(rho)./modes\n", + " rhok .*= -1im\n", + " ifft!(rhok)\n", + " real(rhok)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Landau Damping\n", + "\n", + "[Landau damping - Wikipedia](https://en.wikipedia.org/wiki/Landau_damping)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "landau (generic function with 2 methods)" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function landau( ϵ, kx, meshx, meshv)\n", + " nx = meshx.nx\n", + " nv = meshv.nx\n", + " x = meshx.x\n", + " v = meshv.x\n", + " f = zeros(Complex{Float64},(nx,nv))\n", + " f .= (1.0.+ϵ*cos.(kx*x))/sqrt(2π) .* transpose(exp.(-0.5*v.^2))\n", + " f\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "advection_v! (generic function with 1 method)" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function advection_x!( f, p, meshx, v, dt) \n", + " advection!(f, p, meshx, v, dt)\n", + "end\n", + "\n", + "function advection_v!( f, fᵗ, p, meshx, meshv, ℰ, dt)\n", + " dx = meshx.dx\n", + " nx = meshx.nx\n", + " ρ = compute_rho(meshv, f)\n", + " e = compute_e(meshx, ρ)\n", + " push!(ℰ, 0.5*log(sum(e.*e)*dx))\n", + " transpose!(fᵗ, f)\n", + " advection!(fᵗ, p, meshv, e, dt)\n", + " transpose!(f, fᵗ)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "landau_damping (generic function with 1 method)" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function landau_damping(tf::Float64, nt::Int64)\n", + " \n", + " # Set grid\n", + " p = 3\n", + " nx, nv = 128, 256\n", + " xmin, xmax = 0.0, 4π\n", + " vmin, vmax = -6., 6.\n", + " meshx = UniformMesh(xmin, xmax, nx)\n", + " meshv = UniformMesh(vmin, vmax, nv)\n", + " x, v = meshx.x, meshv.x\n", + " \n", + " # Set distribution function for Landau damping\n", + " ϵ, kx = 0.001, 0.5\n", + " f = landau( ϵ, kx, meshx, meshv)\n", + " \n", + " fᵗ = zeros(Complex{Float64},(nv,nx))\n", + " \n", + " # Set time domain\n", + " dt = tf / nt\n", + " \n", + " # Run simulation\n", + " ℰ = Float64[]\n", + " \n", + " for it in 1:nt\n", + " \n", + " advection_x!( f, p, meshx, v, 0.5dt)\n", + " advection_v!( f, fᵗ, p, meshx, meshv, ℰ, dt)\n", + " advection_x!( f, p, meshx, v, 0.5dt)\n", + " \n", + " end\n", + " \n", + " ℰ\n", + "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 9.649130 seconds (1.64 M allocations: 1.872 GiB, 3.24% gc time)\n" + ] + } + ], + "source": [ + "nt = 1000\n", + "tf = 100.0\n", + "t = range(0.0, stop=tf, length=nt)\n", + "@time nrj = landau_damping(tf, nt);" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "0\n", + "\n", + "\n", + "25\n", + "\n", + "\n", + "50\n", + "\n", + "\n", + "75\n", + "\n", + "\n", + "100\n", + "\n", + "\n", + "-20\n", + "\n", + "\n", + "-15\n", + "\n", + "\n", + "-10\n", + "\n", + "\n", + "-5\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "E\n", + "\n", + "\n", + "\n", + "-0.1533t.-5.5\n", + "\n", + "\n" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "plot( t, nrj; label = \"E\")\n", + "plot!(t, -0.1533*t.-5.50; label=\"-0.1533t.-5.5\")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Callable type" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Advection" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + " Advection(f, p, mesh, v, nv, dt)\n", + "\n", + "Advection type\n", + "\n", + "\"\"\"\n", + "mutable struct Advection\n", + " \n", + " p :: Int64 \n", + " mesh :: UniformMesh\n", + " modes :: Vector{Float64}\n", + " eig_bspl :: Vector{Float64}\n", + " eigalpha :: Vector{Complex{Float64}}\n", + " \n", + " function Advection( p, mesh )\n", + " nx = mesh.nx\n", + " modes = zeros(Float64, nx)\n", + " modes .= [2π * i / nx for i in 0:nx-1]\n", + " eig_bspl = zeros(Float64, nx)\n", + " eig_bspl = zeros(Float64, nx)\n", + " eig_bspl .= bspline(p, -div(p+1,2), 0.0)\n", + " for i in 1:div(p+1,2)-1\n", + " eig_bspl .+= bspline(p, i-(p+1)÷2, 0.0) * 2 .* cos.(i * modes)\n", + " end\n", + " eigalpha = zeros(Complex{Float64}, nx)\n", + " new( p, mesh, modes, eig_bspl, eigalpha )\n", + " end\n", + " \n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [], + "source": [ + "function (adv :: Advection)(f :: Array{Complex{Float64},2}, \n", + " v :: Vector{Float64}, \n", + " dt :: Float64)\n", + " \n", + " nx = adv.mesh.nx\n", + " nv = length(v)\n", + " dx = adv.mesh.dx\n", + " \n", + " fft!(f,1)\n", + " \n", + " @simd for j in 1:nv\n", + " @inbounds alpha = dt * v[j] / dx\n", + " # compute eigenvalues of cubic splines evaluated at displaced points\n", + " ishift = floor(-alpha)\n", + " beta = -ishift - alpha\n", + " fill!(adv.eigalpha,0.0im)\n", + " for i in -div(adv.p-1,2):div(adv.p+1,2)\n", + " adv.eigalpha .+= (bspline(adv.p, i-div(adv.p+1,2), beta) \n", + " .* exp.((ishift+i) * 1im .* adv.modes))\n", + " end\n", + " \n", + " # compute interpolating spline using fft and properties of circulant matrices\n", + " \n", + " @inbounds f[:,j] .*= adv.eigalpha ./ adv.eig_bspl\n", + " \n", + " end\n", + " \n", + " ifft!(f,1)\n", + " \n", + "end " + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "landau_damping_hl (generic function with 1 method)" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function landau_damping_hl(tf::Float64, nt::Int64)\n", + " \n", + " # Set grid\n", + " p = 3\n", + " nx, nv = 128, 256\n", + " xmin, xmax = 0.0, 4π\n", + " vmin, vmax = -6., 6.\n", + " meshx = UniformMesh(xmin, xmax, nx)\n", + " meshv = UniformMesh(vmin, vmax, nv)\n", + " x, v = meshx.x, meshv.x \n", + " dx = meshx.dx\n", + " \n", + " # Set distribution function for Landau damping\n", + " ϵ, kx = 0.001, 0.5\n", + " f = landau( ϵ, kx, meshx, meshv)\n", + " fᵗ = zeros(Complex{Float64},(nv,nx))\n", + " \n", + " # Instantiate advection functions\n", + " advection_x! = Advection(p, meshx)\n", + " advection_v! = Advection(p, meshv)\n", + " \n", + " # Set time step\n", + " dt = tf / nt\n", + " \n", + " # Run simulation\n", + " ℰ = Float64[]\n", + " \n", + " for it in 1:nt\n", + " \n", + " advection_x!(f, v, 0.5dt)\n", + "\n", + " ρ = compute_rho(meshv, f)\n", + " e = compute_e(meshx, ρ)\n", + " \n", + " push!(ℰ, 0.5*log(sum(e.*e)*dx))\n", + " \n", + " transpose!(fᵗ, f)\n", + " advection_v!(fᵗ, e, dt)\n", + " transpose!(f, fᵗ)\n", + " \n", + " advection_x!(f, v, 0.5dt)\n", + " \n", + " end\n", + " \n", + " ℰ\n", + "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 9.974339 seconds (1.43 M allocations: 1.843 GiB, 2.86% gc time)\n" + ] + } + ], + "source": [ + "nt = 1000\n", + "tf = 100.0\n", + "t = range(0.0, stop=tf, length=nt)\n", + "@time nrj = landau_damping_hl(tf, nt);" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "0\n", + "\n", + "\n", + "25\n", + "\n", + "\n", + "50\n", + "\n", + "\n", + "75\n", + "\n", + "\n", + "100\n", + "\n", + "\n", + "-20\n", + "\n", + "\n", + "-15\n", + "\n", + "\n", + "-10\n", + "\n", + "\n", + "-5\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "E\n", + "\n", + "\n", + "\n", + "-0.1533t.-5.5\n", + "\n", + "\n" + ] + }, + "execution_count": 40, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "plot( t, nrj; label = \"E\")\n", + "plot!(t, -0.1533*t.-5.50; label=\"-0.1533t.-5.5\")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Metaprogramming" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "@Strang" + ] + }, + "execution_count": 41, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\n", + "\n", + " @Strang( push_t, push_v )\n", + "\n", + " Apply the second order Strang splitting\n", + "\n", + " push_t and push_v are two function calls with\n", + " `dt` as argument.\n", + "\n", + "\"\"\"\n", + "macro Strang(push_t, push_v)\n", + " return esc(quote\n", + " local full_dt = dt\n", + " dt = 0.5full_dt\n", + " $push_t\n", + " dt = full_dt\n", + " $push_v\n", + " dt = 0.5full_dt\n", + " $push_t\n", + " dt = full_dt\n", + " end)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "landau_with_macro (generic function with 1 method)" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function landau_with_macro(tf::Float64, nt::Int64)\n", + " \n", + " # Set grid\n", + " p = 3\n", + " nx, nv = 128, 256\n", + " xmin, xmax = 0.0, 4π\n", + " vmin, vmax = -6., 6.\n", + " meshx = UniformMesh(xmin, xmax, nx)\n", + " meshv = UniformMesh(vmin, vmax, nv)\n", + " x, v = meshx.x, meshv.x \n", + " dx = meshx.dx\n", + " \n", + " # Set distribution function for Landau damping\n", + " ϵ, kx = 0.001, 0.5\n", + " f = landau( ϵ, kx, meshx, meshv)\n", + " fᵗ = zeros(Complex{Float64},(nv,nx))\n", + " \n", + " # Instantiate distribution functions\n", + " advection_x! = Advection(p, meshx)\n", + " advection_v! = Advection(p, meshv)\n", + " \n", + " push_t!( f, v, dt ) = advection_x!(f, v, dt)\n", + "\n", + " function push_v!( f, fᵗ, meshx, meshv, dt, ℰ)\n", + " \n", + " ρ = compute_rho(meshv, f)\n", + " e = compute_e(meshx, ρ)\n", + " \n", + " push!(ℰ, 0.5*log(sum(e.*e)*dx))\n", + " \n", + " transpose!(fᵗ, f)\n", + " advection_v!(fᵗ, e, dt)\n", + " transpose!(f, fᵗ)\n", + " \n", + " end\n", + " \n", + " # Set time step\n", + " dt = tf / nt\n", + " \n", + " # Run simulation\n", + " ℰ = Float64[]\n", + " \n", + " for it in 1:nt\n", + " \n", + " @Strang( push_t!( f, v, dt ),\n", + " push_v!( f, fᵗ, meshx, meshv, dt, ℰ))\n", + " \n", + " end\n", + " \n", + " ℰ\n", + "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 9.808677 seconds (1.14 M allocations: 1.830 GiB, 2.92% gc time)\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "0\n", + "\n", + "\n", + "25\n", + "\n", + "\n", + "50\n", + "\n", + "\n", + "75\n", + "\n", + "\n", + "100\n", + "\n", + "\n", + "-20\n", + "\n", + "\n", + "-15\n", + "\n", + "\n", + "-10\n", + "\n", + "\n", + "-5\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "E\n", + "\n", + "\n", + "\n", + "-0.1533t.-5.5\n", + "\n", + "\n" + ] + }, + "execution_count": 44, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nt = 1000\n", + "tf = 100.0\n", + "t = range(0.0, stop=tf, length=nt)\n", + "@time nrj = landau_with_macro(tf, nt);\n", + "plot( t, nrj; label = \"E\")\n", + "plot!(t, -0.1533*t.-5.50; label=\"-0.1533t.-5.5\")" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Julia 1.1.0", + "language": "julia", + "name": "julia-1.1" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.1.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/deps_install.jl b/deps_install.jl new file mode 100644 index 0000000..432ead9 --- /dev/null +++ b/deps_install.jl @@ -0,0 +1,14 @@ +using Pkg + +function install_deps() + io = open("REQUIRE", "r") + deps = read(io, String) + + for dep in split(deps) + Pkg.add(String(dep)) + end +end + +install_deps() + +using Plots, FFTW, LinearAlgebra, Statistics, BenchmarkTools