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 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
4
N=151
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.1 and