diff --git a/notebooks/Vlasov-Poisson-1D1V.ipynb b/notebooks/Vlasov-Poisson-1D1V.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..1c513af8d0567db3ed11aaafdd52021f35179dcf --- /dev/null +++ b/notebooks/Vlasov-Poisson-1D1V.ipynb @@ -0,0 +1,574 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "using Plots, LinearAlgebra, Statistics, BenchmarkTools\n", + "using HermiteGF\n", + "using SparseArrays" + ] + }, + { + "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": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "SparseMatrixCSC{Float64,Int64}" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nx = 64\n", + "matrix = spdiagm(-1 => ones(Float64,nx-2),\n", + " 0 => -2*ones(Float64,nx),\n", + " 1 => ones(Float64,nx-2))\n", + "\n", + "typeof(matrix)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "struct UniformMesh\n", + " nx :: Int64\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( nx, dx, x)\n", + " end\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "struct ChebyshevMesh\n", + " nx :: Int64\n", + " x :: Vector{Float64}\n", + " dx :: Vector{Float64}\n", + " function ChebyshevMesh( nodes::HermiteGF.NodesType)\n", + " dx =[ x1 - x0 for (x1,x0) in zip(nodes.xk[2:end],nodes[1:end-1]) ]\n", + " x = nodes.xk\n", + " new( nx, dx, x)\n", + " end\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "f_i'' \\approx \\frac{2 \\Big[ f_{i+1}\n", + "- \\big(1+\\frac{h_i}{h_{i-1}}\\big) f_i\n", + "+\\frac{h_i}{h_{i-1}}f_{i-1} \n", + " \\Big]}\n", + "{ h_i h_{i-1} (1+\\frac{h_i}{h_{i-1}})}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "## Julia implementation of solving Poisson on a square grid using Chebyshev collocation\n", + "# julia> load(\"logg.jl\")\n", + "# julia> loggProfile(10000)\n", + "# error: 3.310154100518143e-7\n", + "# Average time (us): 132.39421844482422\n", + "\n", + "function logg(m)\n", + " (D,xb) = cheb(m+1) # Chebyshev differentiation matrix and nodes\n", + " D = 2*D; xb = 0.5 .+ 0.5*xb # Remap to interval [0,1]\n", + " D2 = -(D*D)[2:end-1,2:end-1] # Negative 1D Laplacian, with Dirichlet BCs\n", + " xd = xb[2:end-1] # Remove Dirichlet BCs\n", + " sinx = sin.(pi*xd'')\n", + " xexact = kron(sinx', sinx) # sin(pi*x)*sin(pi*y)\n", + " b = 2*pi^2*xexact\n", + " # Compute action of inverse of L using sum factorization\n", + " # See Lynch, Rice, and Thomas (1964)\n", + " (Lambda, R) = eig(D2)\n", + " e = ones(m)\n", + " Diag = kron(Lambda,e) + kron(e,Lambda)\n", + " iDiag = 1/Diag\n", + " Rinv = inv(R)\n", + " # L = kron(D2,eye(m)) + kron(eye(m),D2)\n", + " # Linv = kron(R,R) * diagm(reshape(iDiag,m^2,1)) * kron(Rinv,Rinv)\n", + " x = R * (iDiag .* (Rinv*b*Rinv')) * R' # x = Linv * b (with reshaping)\n", + " norm(reshape(x - xexact,m^2,1)) # L2 norm of the vector\n", + "end\n", + "\n", + "function cheb(N)\n", + " x = cos.(pi*(0:N)/N)\n", + " c = ones(N+1)\n", + " c[1] = 2\n", + " c[end] = 2\n", + " for i in 2:2:N+1\n", + " c[i] = - c[i]\n", + " end\n", + " X = repeat(x,1,N+1)\n", + " dX = X .- X'\n", + " D = (c*(1/c')) ./ (dX+eye(N+1)) # off-diagonal entries\n", + " D = D - diagm(sum(D,dims=2)) # diagonal entries\n", + " (D, x)\n", + "end\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [ + { + "ename": "MethodError", + "evalue": "MethodError: no method matching /(::Int64, ::Adjoint{Float64,Array{Float64,1}})\nClosest candidates are:\n /(::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}, !Matched::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}) at int.jl:59\n /(::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}, !Matched::BigInt) at gmp.jl:466\n /(::T<:Integer, !Matched::T<:Integer) where T<:Integer at int.jl:57\n ...", + "output_type": "error", + "traceback": [ + "MethodError: no method matching /(::Int64, ::Adjoint{Float64,Array{Float64,1}})\nClosest candidates are:\n /(::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}, !Matched::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}) at int.jl:59\n /(::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}, !Matched::BigInt) at gmp.jl:466\n /(::T<:Integer, !Matched::T<:Integer) where T<:Integer at int.jl:57\n ...", + "", + "Stacktrace:", + " [1] cheb(::Int64) at ./In[34]:33", + " [2] logg(::Int64) at ./In[34]:8", + " [3] top-level scope at util.jl:156", + " [4] top-level scope at In[36]:1" + ] + } + ], + "source": [ + "@time logg(7)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "struct Poisson\n", + " \n", + " nx :: Int64\n", + " dx :: Float64\n", + " Φ :: Array{Float64,1}\n", + " matrix :: SparseMatrixCSC{Float64,Int64}\n", + " \n", + " function Poisson( meshx )\n", + " nx = meshx.nx\n", + " dx = meshx.dx\n", + " Φ = zeros(Float64,nx)\n", + " matrix = spdiagm(-1 => -ones(Float64,nx-2),\n", + " 0 => +2*ones(Float64,nx),\n", + " 1 => -ones(Float64,nx-2))\n", + " new( nx, dx, Φ, matrix)\n", + "\n", + " end\n", + " \n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "function (p::Poisson)( rho::Array{Float64,1} )\n", + " p.Φ .= p.matrix \\ rho\n", + " (circshift(p.Φ, 1) - circshift(p.Φ, -1)) ./ (2*p.dx)\n", + "end " + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "compute_rho" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\"\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 * vec(sum(real(f), dims=2))\n", + " rho .- mean(rho)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Landau Damping\n", + "\n", + "[Landau damping - Wikipedia](https://en.wikipedia.org/wiki/Landau_damping)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "landau (generic function with 1 method)" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function landau(tf::Float64, nt::Int64)\n", + " \n", + " nx, nv = 128, 256\n", + " xmin, xmax = 0.0, 4*pi\n", + " vmin, vmax = -6., 6.\n", + " meshx = UniformMesh(xmin, xmax, nx)\n", + " meshv = UniformMesh(vmin, vmax, nv)\n", + " x = meshx.x\n", + " v = meshv.x\n", + " dx = meshx.dx\n", + " \n", + " # Create Vlasov-Poisson simulation\n", + " poisson = Poisson(meshx)\n", + " \n", + " eps, kx = 0.001, 0.5\n", + " f = zeros(Complex{Float64},(nx,nv))\n", + " f .= (1.0.+eps*cos.(kx*x))/sqrt(2Ï€) * transpose(exp.(-0.5*v.^2))\n", + "\n", + " Ï = compute_rho( meshv, f)\n", + " \n", + " e = poisson( Ï )\n", + " \n", + " ## Set time domain\n", + " #dt = tf / nt\n", + " #\n", + " ## Run simulation\n", + " #â„° = Float64[]\n", + " #\n", + " #for it in 1:nt\n", + " # advection!(f, p, meshx, v, nv, 0.5*dt)\n", + " # rho = compute_rho(meshv, f)\n", + " # e = compute_e(meshx, rho)\n", + " # push!(â„°, 0.5*log(sum(e.*e)*dx))\n", + " # transpose!(fáµ—, f)\n", + " # advection!(fáµ—, p, meshv, e, nx, dt)\n", + " # transpose!(f, fáµ—)\n", + " # advection!(f, p, meshx, v, nv, 0.5*dt)\n", + " #end\n", + " # \n", + " #â„°\n", + " meshx.x, e\n", + "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n", + "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n", + "<defs>\n", + " <clipPath id=\"clip9500\">\n", + " <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n", + " </clipPath>\n", + "</defs>\n", + "<defs>\n", + " <clipPath id=\"clip9501\">\n", + " <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n", + " </clipPath>\n", + "</defs>\n", + "<polygon clip-path=\"url(#clip9501)\" points=\"\n", + "0,1600 2400,1600 2400,0 0,0 \n", + " \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n", + "<defs>\n", + " <clipPath id=\"clip9502\">\n", + " <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n", + " </clipPath>\n", + "</defs>\n", + "<polygon clip-path=\"url(#clip9501)\" points=\"\n", + "243.029,1503.47 2321.26,1503.47 2321.26,47.2441 243.029,47.2441 \n", + " \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n", + "<defs>\n", + " <clipPath id=\"clip9503\">\n", + " <rect x=\"243\" y=\"47\" width=\"2079\" height=\"1457\"/>\n", + " </clipPath>\n", + "</defs>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 301.847,1503.47 301.847,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 616.342,1503.47 616.342,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 930.838,1503.47 930.838,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 1245.33,1503.47 1245.33,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 1559.83,1503.47 1559.83,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 1874.32,1503.47 1874.32,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 2188.82,1503.47 2188.82,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 243.029,1462.4 2321.26,1462.4 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 243.029,1118.88 2321.26,1118.88 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 243.029,775.359 2321.26,775.359 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 243.029,431.84 2321.26,431.84 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n", + " 243.029,88.3202 2321.26,88.3202 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 243.029,1503.47 2321.26,1503.47 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 243.029,1503.47 243.029,47.2441 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 301.847,1503.47 301.847,1481.63 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 616.342,1503.47 616.342,1481.63 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 930.838,1503.47 930.838,1481.63 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 1245.33,1503.47 1245.33,1481.63 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 1559.83,1503.47 1559.83,1481.63 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 1874.32,1503.47 1874.32,1481.63 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 2188.82,1503.47 2188.82,1481.63 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 243.029,1462.4 274.202,1462.4 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 243.029,1118.88 274.202,1118.88 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 243.029,775.359 274.202,775.359 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 243.029,431.84 274.202,431.84 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 243.029,88.3202 274.202,88.3202 \n", + " \"/>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 301.847, 1557.47)\" x=\"301.847\" y=\"1557.47\">0</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 616.342, 1557.47)\" x=\"616.342\" y=\"1557.47\">2</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 930.838, 1557.47)\" x=\"930.838\" y=\"1557.47\">4</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1245.33, 1557.47)\" x=\"1245.33\" y=\"1557.47\">6</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1559.83, 1557.47)\" x=\"1559.83\" y=\"1557.47\">8</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1874.32, 1557.47)\" x=\"1874.32\" y=\"1557.47\">10</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2188.82, 1557.47)\" x=\"2188.82\" y=\"1557.47\">12</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 219.029, 1479.9)\" x=\"219.029\" y=\"1479.9\">-0.002</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 219.029, 1136.38)\" x=\"219.029\" y=\"1136.38\">-0.001</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 219.029, 792.859)\" x=\"219.029\" y=\"792.859\">0.000</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 219.029, 449.34)\" x=\"219.029\" y=\"449.34\">0.001</text>\n", + "</g>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 219.029, 105.82)\" x=\"219.029\" y=\"105.82\">0.002</text>\n", + "</g>\n", + "<polyline clip-path=\"url(#clip9503)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 301.847,766.938 317.284,741.655 332.722,708.031 348.16,674.57 363.598,641.352 379.035,608.456 394.473,575.963 409.911,543.949 425.349,512.494 440.786,481.671 \n", + " 456.224,451.556 471.662,422.222 487.1,393.738 502.538,366.173 517.975,339.594 533.413,314.065 548.851,289.647 564.289,266.399 579.726,244.378 595.164,223.635 \n", + " 610.602,204.222 626.04,186.185 641.477,169.567 656.915,154.408 672.353,140.745 687.791,128.612 703.228,118.036 718.666,109.044 734.104,101.657 749.542,95.8928 \n", + " 764.979,91.7658 780.417,89.2856 795.855,88.4582 811.293,89.2856 826.73,91.7658 842.168,95.8928 857.606,101.657 873.044,109.044 888.481,118.036 903.919,128.612 \n", + " 919.357,140.745 934.795,154.408 950.233,169.567 965.67,186.185 981.108,204.222 996.546,223.635 1011.98,244.378 1027.42,266.399 1042.86,289.647 1058.3,314.065 \n", + " 1073.73,339.594 1089.17,366.173 1104.61,393.738 1120.05,422.222 1135.49,451.556 1150.92,481.671 1166.36,512.494 1181.8,543.949 1197.24,575.963 1212.67,608.456 \n", + " 1228.11,641.352 1243.55,674.57 1258.99,708.031 1274.43,741.655 1289.86,775.359 1305.3,809.064 1320.74,842.688 1336.18,876.149 1351.61,909.367 1367.05,942.263 \n", + " 1382.49,974.756 1397.93,1006.77 1413.37,1038.23 1428.8,1069.05 1444.24,1099.16 1459.68,1128.5 1475.12,1156.98 1490.55,1184.55 1505.99,1211.12 1521.43,1236.65 \n", + " 1536.87,1261.07 1552.31,1284.32 1567.74,1306.34 1583.18,1327.08 1598.62,1346.5 1614.06,1364.53 1629.49,1381.15 1644.93,1396.31 1660.37,1409.97 1675.81,1422.11 \n", + " 1691.24,1432.68 1706.68,1441.68 1722.12,1449.06 1737.56,1454.83 1753,1458.95 1768.43,1461.43 1783.87,1462.26 1799.31,1461.43 1814.75,1458.95 1830.18,1454.83 \n", + " 1845.62,1449.06 1861.06,1441.68 1876.5,1432.68 1891.94,1422.11 1907.37,1409.97 1922.81,1396.31 1938.25,1381.15 1953.69,1364.53 1969.12,1346.5 1984.56,1327.08 \n", + " 2000,1306.34 2015.44,1284.32 2030.88,1261.07 2046.31,1236.65 2061.75,1211.12 2077.19,1184.55 2092.63,1156.98 2108.06,1128.5 2123.5,1099.16 2138.94,1069.05 \n", + " 2154.38,1038.23 2169.82,1006.77 2185.25,974.756 2200.69,942.263 2216.13,909.367 2231.57,876.149 2247,851.109 2262.44,809.064 \n", + " \"/>\n", + "<polygon clip-path=\"url(#clip9501)\" points=\"\n", + "1958.43,251.724 2249.26,251.724 2249.26,130.764 1958.43,130.764 \n", + " \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 1958.43,251.724 2249.26,251.724 2249.26,130.764 1958.43,130.764 1958.43,251.724 \n", + " \"/>\n", + "<polyline clip-path=\"url(#clip9501)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n", + " 1982.43,191.244 2126.43,191.244 \n", + " \"/>\n", + "<g clip-path=\"url(#clip9501)\">\n", + "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2150.43, 208.744)\" x=\"2150.43\" y=\"208.744\">y1</text>\n", + "</g>\n", + "</svg>\n" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nt = 1000\n", + "tf = 100.0\n", + "t = range(0.0, stop=tf, length=nt)\n", + "x, e = landau(tf, nt)\n", + "plot( x, e)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "ename": "MethodError", + "evalue": "MethodError: no method matching adjoint(::Type{Poisson})\nClosest candidates are:\n adjoint(!Matched::Missing) at missing.jl:79\n adjoint(!Matched::Number) at number.jl:193\n adjoint(!Matched::Adjoint{#s177,#s176} where #s176<:Union{StaticArray{Tuple{N},T,1} where T where N, StaticArray{Tuple{N,M},T,2} where T where M where N} where #s177) at /Users/navaro/.julia/packages/StaticArrays/WmJnA/src/linalg.jl:78\n ...", + "output_type": "error", + "traceback": [ + "MethodError: no method matching adjoint(::Type{Poisson})\nClosest candidates are:\n adjoint(!Matched::Missing) at missing.jl:79\n adjoint(!Matched::Number) at number.jl:193\n adjoint(!Matched::Adjoint{#s177,#s176} where #s176<:Union{StaticArray{Tuple{N},T,1} where T where N, StaticArray{Tuple{N,M},T,2} where T where M where N} where #s177) at /Users/navaro/.julia/packages/StaticArrays/WmJnA/src/linalg.jl:78\n ...", + "", + "Stacktrace:", + " [1] \\(::Type, ::Array{Float64,1}) at ./operators.jl:536", + " [2] (::Poisson)(::Array{Float64,1}) at ./In[5]:2", + " [3] landau(::Float64, ::Int64) at ./In[7]:21", + " [4] top-level scope at util.jl:156", + " [5] top-level scope at In[9]:6" + ] + } + ], + "source": [ + "using Profile\n", + "\n", + "nt = 1000\n", + "tf = 100.0\n", + "t = range(0.0, stop=tf, length=nt)\n", + "@time nrj = landau(tf, nt)\n", + "plot( t, nrj; label = \"E\")\n", + "plot!(t, -0.1533*t.-5.50; label=\"-0.1533t.-5.5\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using ProfileView\n", + "ProfileView.v" + ] + } + ], + "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 +}