Energy term in wave equation simulation

General Discussion
1

Dear users,

I have to implement an energy term defined by

E(t)= integral2d (u_t^2+c^2*(u_x^2+u_y^2));

The domain of the integral is the square [0,3]X[0,3]:

I defined:
real x0=0.0,y0=0.0;
real xF=3.0,yF=3.0;
int n=7;
int N =2^n;
real h=(xF-x0)/N;
real dt=0.163*h;// step in time
td=floor(3/dt);

macro grad(u) [dx(u), dy(u)]//EOM
mesh Th = square(N,N,[x0+(xF-x0)*x,y0+(yF-y0)y]);
real[int] Energy(td+1);
Energy[T]=
int2d(Th)(((u-uold)/dt)^2+c^2grad(uold)’*grad(uold));

u, uold I obtain for each time of the simulation from solving the variational formulation of the wave equation in this domain.

I would like to ask please, is there an error in the calculation of the energy term in this way?

Thanks in advance,
Mordechai.
Forward problem all boundaries Dirichlet.edp (2.0 KB)

2

Hmm the index of your energy T in Energy[T] is real. i guess it should be an integer: Energy[s]… isn’t it?

3

Hi Julien,
Thank you! yes, you are right it should be Energy[s] I changed it now.
Still, the result of the energy term in the forward run should be constant and I obviously don’t obtain it with my calculation, I don’t know where is my mistake?

4

In your code you evaluate the energy after updating uold. Hence the kinetic energy is alway zero:

wave;
uvold=uold;
uold =u;
Energy[s]=int2d(Th)(((u-uold)/dt)^2+c^2*grad(uold)'*grad(uold));

Instead evaluate the energy before updating uold

wave;
Energy[s]=int2d(Th)(((u-uold)/dt)^2+c^2*grad(uold)'*grad(uold));
uvold=uold;
uold =u;

then the energy is pretty much conserved (also maybe you should us u instead of uold with the gradients)

5

Thank you! Julien, I will changed it.

Best regards,
Mordechai.

6

Hi Julien,

I changed in the code of the forward problem the corrections you suggested, thanks a lot!
now the energy is constant for each time of the forward simulation.
I would like to ask please, if you can please check my backward code for the wave equation, because the energy in the backward simulation should also be constant, however it gives me constants results for certain regions.
Backward problem all boundaries Dirichlet without noise.edp (2.4 KB) Forward problem all boundaries Dirichlet without noise.edp (2.0 KB)

Thanks in advance,
Mordechai.