{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Plots" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 1d1v vlasov equation for electrons is given by:\n", "\n", "\\begin{align} &\\frac{\\partial f}{\\partial t} + v \\frac{\\partial f}{\\partial x} - E \\frac{\\partial f}{\\partial v} = 0 \\\\\n", "&- \\frac{\\partial^2 \\phi}{\\partial x^2} = 1 - \\int f \\ \\rm{d}v, \\qquad E = - \\frac{\\partial \\phi}{\\partial x}, \n", "\\end{align} \n", "\n", "with $\\ f > 0, t \\in [0, T]. x \\in [0, L_x], v \\in [v_{min}, v_{max}]$. Where $f$ is $L_x$ and $L_v$ periodic and $\\phi$ also $L_x$ periodic. \n", "\n", "We introduce the splitting given by: \n", "\n", "\\begin{align} \\rm{(A): } &\\frac{\\partial f^*}{\\partial t} - E \\frac{\\partial f^*}{\\partial v} = 0, \\\\\n", "& f^*(t_n) = f(t_n). \n", "\\end{align} \n", "\n", "\n", "\\begin{align} \\rm{(B): } &\\frac{\\partial f}{\\partial t} - v \\frac{\\partial f}{\\partial x} = 0, \\\\\n", "& f(t_n) = f^*(t_{n+1}). \n", "\\end{align} \n", "\n", "One can show that (A) is a constant-coefficient advection problem. So in every timestep we solve the poisson equation for $E$ and solve both advection equations successively." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "function solve_damped(h, λ, T)\n", "#DEFINE PARAMETERS\n", "Lₓ = 12\n", "vₘₐₓ = 6\n", "vₘᵢₙ = -6\n", "Lᵥ = vₘₐₓ - vₘᵢₙ \n", "f(x,v) = (1+0.01*cos.(2*pi*x/Lₓ)).*(1/(2*pi)^(1/2)).*exp.(-v.^2 / 2);\n", "M(v) = (1/(2*pi)^(1/2))*exp.(-v.^2 / 2);\n", "\n", "# for simplicity same space and velocity discr.\n", "dv = h\n", "τ = λ*h; # space-time coupling\n", "x = collect(0:h:Lₓ)\n", "v = collect(vₘᵢₙ:dv:vₘₐₓ)\n", " \n", "Nₕ = length(x)\n", "Nᵥ = length(v)\n", "\n", "# make sure that N is even for the spectral solvers \n", "if(rem(Nₕ,2) != 0)\n", "x = collect(0:h:Lₓ+h)\n", "Nₕ = length(x)\n", "end\n", "\n", "if(rem(Nᵥ,2) != 0)\n", " v = collect(vₘᵢₙ:dv:vₘₐₓ+dv)\n", " Nᵥ = length(v)\n", "end\n", "\n", "#INITIALIZE VECTORS\n", "FieldEnergy = zeros(Int(round(T/τ)))\n", "\n", "sol = zeros(Nᵥ, Nₕ)\n", "for i in 1:Nₕ\n", " for j in 1:Nᵥ \n", " sol[j,i] = f(x[i], v[j])\n", " end\n", "end\n", " \n", "maxw = M(v); \n", "\n", "# weights for our trapezoidal rule integrals\n", "weights = [1,1]/2\n", "λ = zeros(Nᵥ)\n", "for k=1:Nᵥ-1\n", " λ[k:k+1] += weights \n", "end\n", "\n", "weights = [1,1]/2\n", "λ2 = zeros(Nₕ)\n", "for k=1:Nₕ-1\n", " λ2[k:k+1] += weights \n", "end\n", "\n", "\n", "#anim = @animate #infront of the loop if you want to create a gif\n", "for t=1:Int(round(T/τ))\n", " \n", " #CALC. INTEGRAL w/ Trapezoidal Rule\n", " ∫f = zeros(Nₕ)\n", " for k=1:Nₕ\n", " ∫f[k] = dv*dot(λ, sol[:,k])\n", " end\n", "\n", " #SOLVE POISSON w/ spectral solver\n", " u = fftshift(fft(1 - ∫f))\n", " l = -Nₕ/2:1:Nₕ/2-1\n", " i = findin(l,0)[1]\n", " u[1:i-1] = u[1:i-1]./((2pi.*l[1:i-1]./Lₓ).^2);\n", " u[i+1:end] = u[i+1:end]./((2pi.*l[i+1:end]./Lₓ).^2);\n", " u[1] = u[end];\n", " ∫E = real(ifft(fftshift(u)));\n", " \n", " #DIFFERENTIATE w/ FD (there may be a better way to do that, e.g. directly in the fourier space)\n", " E = zeros(Nₕ)\n", " for j = 2:(Nₕ-1)\n", " E[j] = -1/(2h)*(∫E[j+1] - ∫E[j-1])\n", " end\n", " E[end] = -1/(2h)*(∫E[1] - ∫E[end-1]);\n", " E[1] = -1/(2h)*(∫E[2] - ∫E[end]);\n", " \n", " #CALCULATE FIELD ENERGY w/ Trapezoidal Rule\n", " FieldEnergy[t] = h*dot(λ2, (abs.(E).^2)/2);\n", "\n", " #SOLVE ADVECTION IN V SPACE w/ spectral solver\n", " u = fftshift(fft(sol,1),1)\n", " l = -Nᵥ/2:1:Nᵥ/2-1\n", " i = findin(l,0)[1]\n", " u[1:i-1,:] = exp.(-2pi*im./(2Lᵥ).*l[1:i-1].*E[1:i-1].*τ).*u[1:i-1, :]\n", " u[i+1:end,:] = exp.(-2pi*im./(2Lᵥ).*l[i+1:end].*E[i+1:end].*τ).*u[i+1:end, :]\n", " u[1,:] = u[end,:]\n", " sol = real(ifft(fftshift(u,1),1));\n", " \n", " #SOLVE ADVECTION IN X SPACE w/ spectral solver\n", " u = fftshift(fft(sol,2),2)\n", " l = -Nₕ/2:1:Nₕ/2-1\n", " i = findin(l,0)[1]\n", " u[:, 1:i-1] = exp.(-2pi*im./(2Lₓ).*l'[1:i-1].*v'[1:i-1].*τ).*u[:, 1:i-1]\n", " u[:, i+1:end] = exp.(-2pi*im./(2Lₓ).*l'[i+1:end].*v'[i+1:end].*τ).*u[:, i+1:end]\n", " u[:,1] = u[:,end]\n", " sol = real(ifft(fftshift(u,2),2));\n", " \n", " \n", " #PREPARE PLOT (IF NEEDED)\n", " # we are only interessted in how the damping behaves over time\n", " p1 = plot(x,v, sol.-maxw, xlim=(0,Lₓ), ylim=(vₘᵢₙ,vₘₐₓ), xlab = \"space\", ylab = \"velocity\",\n", " title = \"Damping of the 1D1V Vlasov\")\n", " p2 = plot(τ.*(1:t), FieldEnergy[1:t], xlim=(0,T), xlab = \"time\", ylab = \"field energy\",title = \"Field Energy over Time\", label=\"Field Energy\")\n", " \n", " #DISPLAY WHILE IN LOOP (IF NEEDED)\n", " if rem(t,10) == 0\n", " IJulia.clear_output(true)\n", " Plots.display(plot(p1,p2, layout=(2,1), size = (900,700)))\n", " #plot(p1,p2, layout=(2,1), size = (900,700)) #FOR GIF OPTION\n", " end\n", " \n", "end \n", "#gif(anim, \"/tmp/anim_fps30.gif\", fps = 30) #gif output\n", " \n", " return x,v,sol,maxw, FieldEnergy\n", " end;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x,v,sol,maxw, F = solve_damped(0.1, 1, 300); # h = dv = 0.1, λ = 1, T= 300 play with input parameters" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.2", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }