OrdinaryDiffEq.jl 1D wave equation

# Convert this cell to markdown in order to enable Pluto's inbuilt package managerif isdefined(Main, :PlutoRunner)    using Pkg    docsdir = joinpath(@__DIR__, "..", "docs")    if isdir(docsdir)        Pkg.activate(docsdir)    endend

begin    using Test    using Revise    using Printf    using VoronoiFVM    using OrdinaryDiffEqBDF    using OrdinaryDiffEqRosenbrock    using OrdinaryDiffEqSDIRK    using LinearAlgebra    using PlutoUI    using ExtendableGrids    using DataStructures    using GridVisualize    if isdefined(Main, :PlutoRunner)        using CairoMakie        default_plotter!(CairoMakie)        CairoMakie.activate!(type = "png")    endend

The wave equation as system of equations

This is the n-dimensional wave equation:

$$u_{tt}- c^2 \Delta u = 0$$

We can create a system of first order in time PDEs out of this:

$$ \begin{aligned} u_t - v&=0\\ v_t -c^2\Delta u&=0 \end{aligned}$$

This allows for a quick implementation in VoronoiFVM (which may be not the optimal way, in particular with respect to time discretization).

const iu = 1; const iv = 2;

storage(y, u, node, data) = y .= u;

reaction(y, u, node, data) = y[iu] = -u[iv];

flux(y, u, edge, data) = y[iv] = data.c^2 * (u[iu, 1] - u[iu, 2]);

Implementation of transparent or mirror bc

function brea(y, u, node, data)    if node.region == 2        if data.bctype == :transparent            y[iu] = data.c * u[iu]        elseif data.bctype == :mirror            boundary_dirichlet!(y, u, node, species = 1, region = 2, value = 0)        end    end    return nothingend;

Wave velocity:

const c = 0.1

0.1

Domain length:

L = 4

N = 151

dt = 1.0e-2; tend = 100.0;

grid = simplexgrid(range(-L, L, length = N));

sys = VoronoiFVM.System(grid, storage = storage, flux = flux, breaction = brea, reaction = reaction, data = (c = c, bctype = Symbol(bc2type)), species = [1, 2])

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=151, ncells=150,
  nbfaces=2),  
  physics = Physics(data=@NamedTuple{c::Float64, bctype::Symbol}, flux=flux,
  storage=storage, reaction=reaction, breaction=brea, ),  
  num_species = 2)

Perturbation in the center of the domain:

begin    inival = unknowns(sys, inival = 0)    inival[1, :] .= map(x -> cos(κ * π * x) * exp(-x^2 / 0.1), grid)end;

problem = ODEProblem(sys, inival, (0.0, tend));

tsol = solve(    problem, diffeqmethods[method]();    force_dtmin = true,    adaptive = true,    reltol = 1.0e-4,    abstol = 1.0e-5,    dtmin = dt,    progress = true,    progress_steps = 1,    dt = dt);

Boundary condition at x=L:

Reflection (Neumann) bc $\partial_x u|_{x=L}=0$

Package wave number κ: method:

t=49.95

let    u = tsol1(t)    scalarplot(grid, u[1, :], flimits = (-1, 1), clear = true, show = true, title = "t=$(t)", resolution = (600, 150))end

let    vis = GridVisualizer()    scalarplot!(vis, sys, tsol1, colormap = :bwr, limits = (-1, 1), levels = (-0.9:0.2:0.9))    reveal(vis)end

tsol1 = reshape(tsol, sys);

diffeqmethods = OrderedDict(    "QNDF2" => QNDF2,    "FBDF" => FBDF,    "Rosenbrock23 (Rosenbrock)" => Rosenbrock23,    "Implicit Euler" => ImplicitEuler,    "Implicit Midpoint" => ImplicitMidpoint,)

OrderedDict{String, UnionAll} with 5 entries:
  "QNDF2"                     => QNDF2
  "FBDF"                      => FBDF
  "Rosenbrock23 (Rosenbrock)" => Rosenbrock23
  "Implicit Euler"            => ImplicitEuler
  "Implicit Midpoint"         => ImplicitMidpoint