Hi Dedalus experts,
I tried implementing a nonlinear wave equation
w_{tt} = w_{xx} + alpha*w^2
with homogeneous Dirichlet boundary conditions.
The relevant code lines are these:
a = 0 # even unstable without nonlinearity
# Bases and domain
xcoord = d3.Coordinate('x')
dist = d3.Distributor(xcoord, comm=MPI.COMM_SELF, dtype=float)
xbasis = d3.Chebyshev(xcoord, size=Nx, bounds=(xmin, xmax), dealias=dealias)
domain = Domain(dist, bases=[xbasis])
# Fields
w = dist.Field(bases=xbasis, name='w')
wt = dist.Field(bases=xbasis, name='wt')
tauw = dist.Field(name='tauw')
tauwt = dist.Field(name='tauw1')
t = dist.Field()
tempul = dist.Field()
lift_basis = xbasis.derivative_basis(1)
x = dist.local_grid(xbasis)
def lift(A, n): return d3.Lift(A, lift_basis, n)
def dx(p): return d3.Differentiate(p, xcoord)
# Problems
nonlinProblem = d3.IVP([w, wt, tauw, tauwt], time=t, namespace=locals())
nonlinProblem.add_equation("dt(w) - wt + tauw = 0")
nonlinProblem.add_equation("dt(wt) + lift(tauwt, -1) - dx(dx(w)) = a*w**2")
nonlinProblem.add_equation("w(x='left') = 0")
nonlinProblem.add_equation("w(x='right') = 0")
nonlinSolver = nonlinProblem.build_solver(d3.RK443)
w['g'] = np.exp(-100*x**2)
wt['g'] = 0
This blows up after a short time. I guess something's messed up with the tau terms, but I don't know what. A finer time step (like 1e-5) does not help a lot.
Any help would be very much appreciated!
Best,
Judith