Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Time-dependency #63

Open
engsbk opened this issue Aug 17, 2022 · 2 comments
Open

Time-dependency #63

engsbk opened this issue Aug 17, 2022 · 2 comments

Comments

@engsbk
Copy link

engsbk commented Aug 17, 2022

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')

@maroba
Copy link
Owner

maroba commented Aug 17, 2022

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).

@engsbk
Copy link
Author

engsbk commented Aug 17, 2022

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants