Solving wave equation with Leap-frog scheme#

Keywords: PDE:Wave Method:FEM FE:Lagrange Condition:Dirichlet Solver:Direct

So far, only the harmonic problems were considered. Time problem may also be solved using XLiFE++. But there is no specific tools dedicated to. Users have to implement the time loop related to the finite difference time scheme they choose.

As an example, consider the wave equation:

\begin{array}{r} {\begin{cases} \frac{\partial^{2} u}{\partial t^{2}} - c^{2} Δ u = f & in Ω \times] 0, T [ \\ \frac{\partial u}{\partial n} = 0 & in \partial Ω \times] 0, T [ \\ u (x, 0) = \frac{\partial u}{\partial t} (x, 0) = 0 & in Ω \end{cases} \end{array}

Using classical leap-frog scheme with time discretization $t^{n} = n Δ t$ , leads to ( $u^{n}$ approximates $u (x, t^{n})$ ):

\begin{array}{r} {\begin{cases} u^{n + 1} = 2 u^{n} - u^{n - 1} + (c Δ t)^{2} Δ u^{n} + (Δ t)^{2} f^{n} & in Ω, \forall n > 1 \\ \frac{\partial u^{n}}{\partial n} = 0 & in \partial Ω, \forall n > 1 \\ u^{0} = u^{1} = 0 & in Ω \end{cases} \end{array}

or, in variational form, $\forall v \in V = H^{1} (Ω)$ :

\begin{array}{r} {\begin{cases} \int_{Ω} u^{n + 1} v = 2 \int_{Ω} u^{n} v - \int_{Ω} u^{n - 1} v - (c Δ t)^{2} \int_{Ω} \nabla u^{n} . \nabla v + (Δ t)^{2} \int_{Ω} f^{n} v \forall n > 1 \\ u^{0} = u^{1} = 0 in Ω \end{cases} \end{array}

When approximating space $V$ by a finite dimension space $V_{h}$ with basis $(w_{i})_{i = 1, p}$ , the variational formulation is reinterpreted in terms of matrices and vectors as follows:

\begin{array}{r} {\begin{cases} U^{n + 1} = 2 U^{n} - U^{n - 1} - M^{- 1} ((c Δ t)^{2} K U^{n} - (Δ t)^{2} F^{n}) \forall n > 1 \\ U^{0} = U^{1} = 0 in Ω \end{cases} \end{array}

where

M_{i j} = \int_{Ω} w_{i} w_{j}, K_{i j} = \int_{Ω} \nabla w_{i} . \nabla w_{j}, (F^{n})_{i} = \int_{Ω} f^{n} w_{i} .

The XLiFE++ implementation of this scheme on the unity square when using P1 Lagrange interpolation looks like ( $f (x, t) = h (t) g (x)$ ):

#include "xlife++.h"
using namespace xlifepp;

Real g(const Point& P, Parameters& pa = defaultParameters)
{
  Real d=P.distance(Point(0.5, 0.5));
  Real R= 0.02;        // source radius
  Real amp= 1./(pi_*R*R); // source amplitude (constant power)
  if (d<=0.02) return amp; else return 0.;
}

Real h(const Real& t)
{
  Real a=10000, t0=0.04 ; //gaussian slope and center
  return exp(-a*(t-t0)*(t-t0));
}

int main(int argc, char** argv)
{
  init(argc, argv, _lang=en);
  // create a mesh and domain omega
  SquareGeo sq(_origin=Point(0., 0.), _length=1, _nnodes=70);
  Mesh mesh2d(sq, _shape=triangle, _generator=structured);
  Domain omega=mesh2d.domain("Omega");
  // create interpolation
  Space V(_domain=omega, _interpolation=P1, _name="V");
  Unknown u(V, _name="u");
  TestFunction v(u, _name="v");
  // define FE terms
  TermMatrix A(intg(omega, grad(u)|grad(v)), _name="A"), M(intg(omega, u*v), _name="M");
  TermVector G(intg(omega, g*v), _name="G");
  TermMatrix L; ldltFactorize(M, L);
  // leap-frog scheme
  Real c=1, dt=0.004, dt2=dt*dt, cdt2=c*c*dt2;
  Number nbt=200;
  TermVectors U(nbt);  // to store solution at t=ndt
  TermVector zeros(u, omega, 0.);  U(1)=zeros;  U(2)=zeros;
  Real t=dt;
  for (Number n=2; n<nbt; n++, t+=dt)
  {
    U(n+1)=2.*U(n)-U(n-1)-factSolve(L, cdt2*(A*U(n))-dt2*h(t)*G);
  }
  saveToFile("U", U, vtu);
  return 0;
}

Note the very simple syntax taken into account the leap-frog scheme. The following figure represents the solution at different instants for a constant source localized in disk with center $(0.5, 0.5)$ , radius $R = 0.02$ and time excitation that is a Gaussian function. For chosen parameter $d t = 0.04$ , the leap-frog scheme is stable (it satisfies the CFL condition) but dispersion effects obviously appear.

../_images/wave2d.png — Fig. 12 Solution of the wave equation at different instants for a constant source localized in disk with center $(0.5, 0.5)$ #