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 OrdinaryDiffEq
    using LinearAlgebra
    using PlutoUI
    using ExtendableGrids
    using DataStructures
    using GridVisualize,CairoMakie
    CairoMakie.activate!(type="svg")
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 oder 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
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}, Matrix{Float64}}(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.10.2 and

CairoMakie 0.11.10
DataStructures 0.18.20
ExtendableGrids 1.5.0
GridVisualize 1.5.1
OrdinaryDiffEq 6.74.1
Pkg 1.10.0
PlutoUI 0.7.59
Revise 3.5.14
VoronoiFVM 1.20.0