A coupled PDE/ODE IVP system

[James Taylor]

Dear Dedalus Users,

I'd like to solve a 1D IVP problem formulated in the attached PDF file. The problem is described by a coupled PDE/ODE system and involves a time-dependent boundary condition. An analytical solution exists. I'd like to recover this solution numerically.

I'm trying to implement the problem in Dedalus. My attempt is in the attached script. I've defined 's' as a field with only an ξ-basis (since it doesn't depend on t). The script produces the following error message:

ValueError: Pencil () has 3 constant equations for 1 constant variables plus 1 differential equations / tau terms.

Can you provide guidance or suggest modifications to my script that would make it work?

Thanks in advance for all your help!  

James

[Keaton Burns]

Hi James,

In Dedalus v2, the differential order of the system is inferred from the number of derivatives on the LHS. You can overwrite this, and mark equations as differential, with the tau=True keyword.  I’d suggest doing this the dt(u) equation. However, I’d also suggest switching to Dedalus v3, since we’re phasing out support for v2 going forward. 

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/d2dcb12b-a6b5-41bd-ad78-0ec807cc8d36n%40googlegroups.com.

[James Taylor]

Hi Keaton,

Thanks for your feedback. I've switched my script to the Dedalus 3 format (see the script below). Specifically, I've introduced the tau fields 'tau_1' and 'tau_2' for the variable u(ξ,t) (governed by the 2nd-order PDE in ξ). I've also introduced the field 's' for the variable s(t) (governed by the ODE in t).

Running the script, I'm getting the error "ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()". This might suggest that the script is mixing up dedalus fields and numpy arrays. However, the script only uses the Dedalus field objects, and not the numpy arrays of the local data.

Any advice on what to do to make it work would be appreciated!

James

===

import numpy as np
import dedalus.public as d3
import matplotlib.pyplot as plt
import logging
logger = logging.getLogger(__name__)

# Parameters
nξ = 128
timestepper = d3.RK443
timestep = 1e-3
stop_sim_time = 20
dealias = 3/2
dtype = np.float64

# Bases
ξcoord = d3.Coordinate('ξ')
dist = d3.Distributor(ξcoord, dtype=dtype)
ξbasis = d3.Chebyshev(ξcoord, size=nξ, bounds=(0,1), dealias=dealias)
ξ = dist.local_grid(ξbasis)

# Fields
u = dist.Field(name='u', bases=(ξbasis))
s = dist.Field(name='s')
tau_1 = dist.Field(name='tau_1')
tau_2 = dist.Field(name='tau_2')

# Substitutions
dξ = lambda A: d3.Differentiate(A, ξcoord)
lift_basis = ξbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
uξ = dξ(u) + lift(tau_1) # substitution for uξ
uξξ = dξ(uξ) + lift(tau_2) # substitution for uξξ

# Problem
problem = d3.IVP([u, tau_1, tau_2], namespace=locals())
problem.add_equation("dt(u) = ξ*( -(1/s)*uξ(ξ=1) )*(1/s)*uξ + (1/s**2)*uξξ") # PDE for u(ξ,t) with dt(s) substituted
problem.add_equation("dt(s) = -(1/s)*uξ(ξ=1)") # ODE for s(t)
problem.add_equation("dξ(u)(ξ=0) = -s*exp(t)") # left BC
problem.add_equation("u(ξ=1) = 0") # right BC

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time