From ccff9328e6356c120146a1dabc22cfe6dd971949 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 23 Aug 2022 08:26:10 +0200 Subject: [PATCH] KdV example in the docs (#144) * KdV example in the docs * use animation from the comment * set version to 0.5.21 --- Project.toml | 2 +- docs/make.jl | 1 + docs/src/tutorials/advection_diffusion.md | 13 ++ .../tutorials/constant_linear_advection.md | 14 ++ docs/src/tutorials/kdv.md | 158 ++++++++++++++++++ .../tutorials/variable_linear_advection.md | 14 ++ docs/src/tutorials/wave_equation.md | 14 ++ 7 files changed, 215 insertions(+), 1 deletion(-) create mode 100644 docs/src/tutorials/kdv.md diff --git a/Project.toml b/Project.toml index 1f6872f8..1eb16e9f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SummationByPartsOperators" uuid = "9f78cca6-572e-554e-b819-917d2f1cf240" author = ["Hendrik Ranocha"] -version = "0.5.20" +version = "0.5.21" [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" diff --git a/docs/make.jl b/docs/make.jl index 49937b5a..dbeca512 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -57,6 +57,7 @@ makedocs( "tutorials/advection_diffusion.md", "tutorials/variable_linear_advection.md", "tutorials/wave_equation.md", + "tutorials/kdv.md", ], "Applications & references" => "applications.md", "Benchmarks" => "benchmarks.md", diff --git a/docs/src/tutorials/advection_diffusion.md b/docs/src/tutorials/advection_diffusion.md index b107a787..3fe07197 100644 --- a/docs/src/tutorials/advection_diffusion.md +++ b/docs/src/tutorials/advection_diffusion.md @@ -96,3 +96,16 @@ end ![example_advection_diffusion_animation](https://user-images.githubusercontent.com/12693098/119226459-7b242c00-bb09-11eb-848b-d09590aa1c31.gif) + +## Package versions + +These results were obtained using the following versions. + +```@example advection_diffusion +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["SummationByPartsOperators", "OrdinaryDiffEq"], + mode=PKGMODE_MANIFEST) +``` diff --git a/docs/src/tutorials/constant_linear_advection.md b/docs/src/tutorials/constant_linear_advection.md index 53c93bac..d149a5be 100644 --- a/docs/src/tutorials/constant_linear_advection.md +++ b/docs/src/tutorials/constant_linear_advection.md @@ -146,3 +146,17 @@ savefig("example_linear_advection_Galerkin.png"); ``` ![](example_linear_advection_Galerkin.png) + + +## Package versions + +These results were obtained using the following versions. + +```@example linear_advection +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["SummationByPartsOperators", "OrdinaryDiffEq"], + mode=PKGMODE_MANIFEST) +``` diff --git a/docs/src/tutorials/kdv.md b/docs/src/tutorials/kdv.md new file mode 100644 index 00000000..649ab29b --- /dev/null +++ b/docs/src/tutorials/kdv.md @@ -0,0 +1,158 @@ +# Linear advection diffusion equation with periodic boundary conditions + +Let's consider the Korteweg-de Vries (KdV) equation + +```math +\begin{aligned} + \partial_t u(t,x) + \partial_x \frac{u(t,x)^2}{2} + \partial_x^3 u(t,x) &= 0, && t \in (0,T), x \in (x_{min}, x_{max}), \\ + u(0,x) &= u_0(x), && x \in (x_{min}, x_{max}), \\ +\end{aligned} +``` + +with periodic boundary conditions. The KdV equation has the quadratic invariant + +```math + J = \frac{1}{2} \int u(t,x)^2 \mathrm{d}x. +``` + +A classical trick to conserve this invariant is to use following split form + +```math + u_t + \frac{1}{3} (u^2)_x + \frac{1}{3} u u_x = 0. +``` + +Indeed, integration by parts with periodic boundary conditions yields + +```math +\begin{aligned} + \partial_t J + &= + \int u u_t + = + -\frac{1}{3} \int u (u^2)_x - \frac{1}{3} \int u^2 u_x + \\ + &= + 0 + \frac{1}{3} \int u_x u^2 - \frac{1}{3} \int u^2 u_x + = + 0. +\end{aligned} +``` + +## Basic example using finite difference SBP operators + +Let's create an appropriate discretization of this equation step by step. At first, +we load packages that we will use in this example. + +```@example kdv +using SummationByPartsOperators, OrdinaryDiffEq +using LaTeXStrings; using Plots: Plots, plot, plot!, savefig +``` + +Next, we specify the initial data as Julia function as well as the +spatial domain and the time span. Here, we use an analytic soliton solution +of the KdV equation for the initial data. + +```@example kdv +# traveling wave solution (soliton) +get_xmin() = 0.0 # left boundary of the domain +get_xmax() = 80.0 # right boundary of the domain +get_c() = 2 / 3 # wave speed +function usol(t, x) + xmin = get_xmin() + xmax = get_xmax() + μ = (xmax - xmin) / 2 + c = get_c() + A = 3 * c + K = sqrt(1/c - 1) + x_t = mod(x - c*t - xmin, xmax - xmin) + xmin - μ + + A / cosh(sqrt(3*A) / 6 * x_t)^2 +end + +tspan = (0.0, (get_xmax() - get_xmin()) / (3 * get_c()) + 10 * (get_xmax() - get_xmin()) / get_c()) +``` + +Next, we implement the semidiscretization using the interface of +[OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) +which is part of [DifferentialEquations.jl](https://diffeq.sciml.ai/latest/). +For simplicity, we just use the out-of-place version here since we do not have +to worry about appropriate temporary buffers when using automatic differentiation +in implicit time integration methods. + +```@example kdv +function kdv(u, parameters, t) + D1 = parameters.D1 + D3 = parameters.D3 + + # conservative semidiscretization using a split form + return (-1 / 3) * (u .* (D1 * u) + D1 * (u.^2)) - D3 * u +end +``` + +Next, we choose first- and third-derivative SBP operators `D1, D3`, evaluate +the initial data on the grid, and set up the semidiscretization as an ODE problem. + +```@example kdv +N = 128 # number of grid points +D1 = periodic_derivative_operator(derivative_order=1, accuracy_order=8, + xmin=get_xmin(), xmax=get_xmax(), N=N) +D3 = periodic_derivative_operator(derivative_order=3, accuracy_order=8, + xmin=get_xmin(), xmax=get_xmax(), N=N) +u0 = usol.(first(tspan), grid(D1)) +parameters = (; D1, D3) + +ode = ODEProblem(kdv, u0, tspan, parameters); +``` + +Finally, we can solve the ODE using a Rosenbrock method with adaptive time stepping. +We use such a linearly implicit time integration method since the third-order +derivative makes the system stiff. + +```@example kdv +sol = solve(ode, Rodas5(), saveat=range(first(tspan), stop=last(tspan), length=200)) + +plot(xguide=L"x", yguide=L"u") +plot!(evaluate_coefficients(sol[1], D1), label=L"u_0") +plot!(evaluate_coefficients(sol[end], D1), label=L"u_\mathrm{numerical}") +savefig("example_kdv.png"); +``` + +![](example_kdv.png) + + +## Advanced visualization + +Let's create an animation of the numerical solution. + +```julia +using Printf; using Plots: Animation, frame, gif + +let anim = Animation() + idx = 1 + x, u = evaluate_coefficients(sol[idx], D1) + fig = plot(x, u, xguide=L"x", yguide=L"u", xlim=extrema(x), ylim=(-0.05, 2.05), + label="", title=@sprintf("\$t = %6.2f \$", sol.t[idx])) + for idx in 1:length(sol.t) + fig[1] = x, sol.u[idx] + plot!(title=@sprintf("\$t = %6.2f \$", sol.t[idx])) + frame(anim) + end + gif(anim, "example_kdv.gif") +end +``` + +![example_kdv_animation](https://user-images.githubusercontent.com/12693098/186075685-a4f12cf0-df6d-486a-ba5c-04d7a7466743.gif) + + +## Package versions + +These results were obtained using the following versions. + +```@example kdv +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["SummationByPartsOperators", "OrdinaryDiffEq"], + mode=PKGMODE_MANIFEST) +``` diff --git a/docs/src/tutorials/variable_linear_advection.md b/docs/src/tutorials/variable_linear_advection.md index 43adbc3f..76006f85 100644 --- a/docs/src/tutorials/variable_linear_advection.md +++ b/docs/src/tutorials/variable_linear_advection.md @@ -61,3 +61,17 @@ savefig("example_linear_advection.png"); ``` ![](example_linear_advection.png) + + +## Package versions + +These results were obtained using the following versions. + +```@example variable_linear_advection +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["SummationByPartsOperators", "OrdinaryDiffEq"], + mode=PKGMODE_MANIFEST) +``` diff --git a/docs/src/tutorials/wave_equation.md b/docs/src/tutorials/wave_equation.md index cecac42d..3b6d000f 100644 --- a/docs/src/tutorials/wave_equation.md +++ b/docs/src/tutorials/wave_equation.md @@ -105,3 +105,17 @@ create_gif(Val(:HomogeneousNeumann), Val(:NonReflecting)) ``` ![wave_equation_HomogeneousNeumann_NonReflecting](https://user-images.githubusercontent.com/12693098/119228041-5633b700-bb11-11eb-9c17-bc56c906dae3.gif) + + +## Package versions + +These results were obtained using the following versions. + +```@example wave_equation +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["SummationByPartsOperators", "OrdinaryDiffEq"], + mode=PKGMODE_MANIFEST) +```