Wave equation with Dirichlet bc blows up

36 views
Skip to first unread message

Ju Dith

unread,
Jun 12, 2024, 6:27:17 AMJun 12
to Dedalus Users
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

Keaton Burns

unread,
Jun 12, 2024, 9:11:21 AMJun 12
to dedalu...@googlegroups.com
Hi Judith,

This is the right approach for reducing the system to first-order in time using “wt”, but that equation actually doesn’t need a tau term. Both tau terms should go in the equation with the second-order spatial derivatives, which you could do either in:

  • first order form, making the substitution wx = dx(w)+lift(tauw,-1) and using dx(wx) in the equation.

or

  • second order form, adding lift(tauw,-2) to the equation, and making lift_basis = xbasis.derivative_basis(2).

Best,
-Keaton


--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/a5d9c535-63c8-42ea-ac92-6f56344c38dcn%40googlegroups.com.

Ju Dith

unread,
Jun 12, 2024, 10:29:43 AMJun 12
to dedalu...@googlegroups.com
Hi Keaton,

Thanks a lot, this works!

Judith

You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/bgSseJJA9pM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBhidU3C%2B10P9fHRszQ%3DToCPPgLkHAnr6pSZziznLssWCw%40mail.gmail.com.
Reply all
Reply to author
Forward
0 new messages