OrdinaryDiffEq.jl 1D wave equation

begin
    import Pkg as _Pkg
    haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
    using Test
    using Revise
    using Printf
    using VoronoiFVM
    using OrdinaryDiffEqBDF
    using OrdinaryDiffEqRosenbrock
    using OrdinaryDiffEqSDIRK
    using LinearAlgebra
    using PlutoUI
    using ExtendableGrids
    using DataStructures
    using GridVisualize, CairoMakie
    CairoMakie.activate!(type = "png")
end

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 nothing
end;

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)", Plotter = CairoMakie, resolution = (600, 150))
end

let
    vis = GridVisualizer(Plotter = CairoMakie)
    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

Built with Julia 1.11.2 and