findiff icon indicating copy to clipboard operation
findiff copied to clipboard

Time-dependency

Open engsbk opened this issue 1 year ago • 2 comments

Hi! Thank you for the fantastic library!

I'm trying to simulate the wave equation using findiff, but I'm not getting the expected results even after changing the accuracy multiple times. The equation I'm trying to simulate is: Screen Shot 2022-08-17 at 12 33 26 AM

The results I'm getting are: output

I'm expecting a wave that starts from the left and propagates all the way to the right of the domain, but this is not the case. Also, I noticed that the middle-time instances have unreasonable values for u(x,t). Normally, it should not be more than 1.0 or less than -1.0.

Please advise. Below is my code:


import numpy as np
import matplotlib.pyplot as plt

from findiff import Gradient, Divergence, Laplacian, Curl
from findiff import FinDiff, PDE, Coef, BoundaryConditions, Identity
from findiff.stencils import Stencil
from itertools import product


shape = (100, 100)
x, t = np.linspace(0, 1, shape[0]), np.linspace(0, 1, shape[1])
dx, dt = x[1]-x[0], t[1]-t[0]
X, T = np.meshgrid(x, t, indexing='ij')


freq = 2
c = 1
a = 10 # accuracy
L = FinDiff(1, dt, 2, acc=a) - Coef(c**2) * FinDiff(0, dx, 2, acc=a)
f = np.zeros(shape)

bc = BoundaryConditions(shape)

#IC
bc[:, 0] = 0.   # Dirichlet BC

#BC
bc[0, :] =  np.sin(2*np.pi*freq*t)  # Dirichlet BC, Left
bc[-1, :] = 0.   # Dirichlet BC, Right

pde = PDE(L, f, bc)
u = pde.solve()


fig = plt.figure(figsize=(8, 15), dpi=80)
gs = fig.add_gridspec(5, hspace=0.5)
ax = gs.subplots()
plt.xlim(x[0],x[-1])

ax[0].plot(x, u[:,0])
ax[0].set_title('t=0.0')

ax[1].plot(x, u[:,25])
ax[1].set_title('t=0.25')

ax[2].plot(x, u[:,50])
ax[2].set_title('t=0.5')

ax[3].plot(x, u[:,75])
ax[3].set_title('t=0.75')

ax[4].plot(x, u[:,99])
ax[4].set_title('t=1.0')

engsbk avatar Aug 17 '22 04:08 engsbk

Hi,

I think you are actually solving a different problem. The PDE class is suitable for pure boundary value problems. However, the wave equation poses an initial value + boundary value problem. In addition to the the boundary conditions for x, you need two initial conditions. The first is for the u(t=0, x) and the second for \partial_t u(t=0, x). I wouldn't use the PDE class for this kind of problem, but instead derive an iterative scheme using the Stencil class (there will be an example on the docs pages, soon).

One other thing I noticed: the left boundary condition is a time dependent Dirichlet BC and the right one is a zero Dirichlet BC. This doesn't have the effect you want. You said you want to model a wave passing by. But in your case, the left BC is continuously pumping energy into the system sending out waves in both directions and the waves traveling along the positive x-axis will be reflected at the right bounday (u=0 there).

maroba avatar Aug 17 '22 13:08 maroba

think you are actually solving a different problem. The PDE class is suitable for pure boundary value problems. However, the wave equation poses an initial value + boundary value problem. In addition to the the boundary conditions for x, you need two initial conditions. The first is for the u(t=0, x) and the second for \partial_t u(t=0, x). I wouldn't use the PDE class for this kind of problem, but instead derive an iterative scheme using the Stencil class (there will be an example on the docs pages, soon).

Looking forward for that!

One other thing I noticed: the left boundary condition is a time dependent Dirichlet BC and the right one is a zero Dirichlet BC. This doesn't have the effect you want. You said you want to model a wave passing by. But in your case, the left BC is continuously pumping energy into the system sending out waves in both directions and the waves traveling along the positive x-axis will be reflected at the right bounday (u=0 there).

Yes, my apologies, I forgot to mention that it should be a continuous wave starting from the left. However, it should be absorbed at the right boundary instead of getting reflected.

engsbk avatar Aug 17 '22 14:08 engsbk