In [26]:
import numpy as np
from IPython.display import Image, display

import matplotlib.pyplot as plt
%matplotlib inline

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Dedalus¶

  • Dedalus is a Python package for solving PDEs using spectral methods

  • Very recently published (Burns et al Phys Rev Research 2019)

  • Idea: specify our PDE in a high-level symbolic language, and then Dedalus will automatically generate optimized code to solve it

  • Requires a few hardware configurations and dependencies

  • Other "high-level" PDE solvers exist, especially for finite-element methods (Ansys, COMSOL, etc), but Dedalus is high quality, actively maintained (difficult to find in academic software), and completely open-source

Installation¶

  • Instructions here
  • Dedalus has a lot of dependencies, so it's best to install it in a conda environment
  • On my Apple Silicon laptop, I had to use the second set of instructions for custom conda installation here to install Dedalus
In [ ]:

The Kuramoto-Sivashinsky equation¶

  • A nonlinear PDE that is a relative of Burgers' equation that describes the dynamics of flame fronts

$$ \frac{\partial u}{\partial t} + \frac{\partial^4 u}{\partial x^4}+ \nu \frac{\partial^2 u}{\partial x^2} + u \frac{\partial u}{\partial x} = 0 $$

In [27]:
import numpy as np
import matplotlib.pyplot as plt
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)


## Specify simulation parameters
Lx = 50
Nx = 1024
nu = 1e0
dealias = 3/2
stop_sim_time = 200
timestepper = d3.SBDF2
timestep = 2e-3
dtype = np.float64


## Set up bases
xcoord = d3.Coordinate('x')
dist = d3.Distributor(xcoord, dtype=dtype)
xbasis = d3.RealFourier(xcoord, size=Nx, bounds=(0, Lx), dealias=dealias) # periodic basis

## Define dynamical Fields
u = dist.Field(name='u', bases=xbasis)

## Define operators
dx = lambda A: d3.Differentiate(A, xcoord)

## State PDE Problem
problem = d3.IVP([u], namespace=locals())
problem.add_equation("dt(u) + dx(dx(dx(dx(u)))) + nu * dx(dx(u)) = -u*dx(u)")


## Make random initial conditions that vanish at the boundaries
x = dist.local_grid(xbasis)
from scipy.ndimage import gaussian_filter1d
ic = gaussian_filter1d(np.random.random(x.shape), 10) * np.sin(np.pi * x / x[-1] )
dx_val = (2 * np.pi)/len(u['g'])
u['g'] = ic / (np.sum(ic) * dx_val)

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

## Main loop
u.change_scales(1)
u_list = [np.copy(u['g'])]
t_list = [solver.sim_time]
while solver.proceed:
    solver.step(timestep)
    if solver.iteration % 500 == 0:
        logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))
    if solver.iteration % 50 == 0:
        u.change_scales(1)
        u_list.append(np.copy(u['g']))
        t_list.append(solver.sim_time)

## Plotting
plt.figure(figsize=(12, 6))
plt.pcolormesh(np.array(t_list), x.ravel(), np.array(u_list).T, cmap='RdBu_r', shading='gouraud', rasterized=True)
plt.xlim(0, stop_sim_time)
plt.ylim(0, Lx)
plt.xlabel('t')
plt.ylabel('x')
plt.title(f'Kuramoto-Sivashinsky, nu={nu}')
plt.tight_layout()

2022-10-14 00:00:43,425 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 4.2e+01/s
2022-10-14 00:00:43,556 __main__ 0/1 INFO :: Iteration=500, Time=1.000000e+00, dt=2.000000e-03
2022-10-14 00:00:43,686 __main__ 0/1 INFO :: Iteration=1000, Time=2.000000e+00, dt=2.000000e-03
2022-10-14 00:00:43,820 __main__ 0/1 INFO :: Iteration=1500, Time=3.000000e+00, dt=2.000000e-03
2022-10-14 00:00:43,952 __main__ 0/1 INFO :: Iteration=2000, Time=4.000000e+00, dt=2.000000e-03
2022-10-14 00:00:44,078 __main__ 0/1 INFO :: Iteration=2500, Time=5.000000e+00, dt=2.000000e-03
2022-10-14 00:00:44,205 __main__ 0/1 INFO :: Iteration=3000, Time=6.000000e+00, dt=2.000000e-03
2022-10-14 00:00:44,333 __main__ 0/1 INFO :: Iteration=3500, Time=7.000000e+00, dt=2.000000e-03
2022-10-14 00:00:44,467 __main__ 0/1 INFO :: Iteration=4000, Time=8.000000e+00, dt=2.000000e-03
2022-10-14 00:00:44,601 __main__ 0/1 INFO :: Iteration=4500, Time=9.000000e+00, dt=2.000000e-03
2022-10-14 00:00:44,736 __main__ 0/1 INFO :: Iteration=5000, Time=1.000000e+01, dt=2.000000e-03
2022-10-14 00:00:44,863 __main__ 0/1 INFO :: Iteration=5500, Time=1.100000e+01, dt=2.000000e-03
2022-10-14 00:00:44,989 __main__ 0/1 INFO :: Iteration=6000, Time=1.200000e+01, dt=2.000000e-03
2022-10-14 00:00:45,114 __main__ 0/1 INFO :: Iteration=6500, Time=1.300000e+01, dt=2.000000e-03
2022-10-14 00:00:45,239 __main__ 0/1 INFO :: Iteration=7000, Time=1.400000e+01, dt=2.000000e-03
2022-10-14 00:00:45,365 __main__ 0/1 INFO :: Iteration=7500, Time=1.500000e+01, dt=2.000000e-03
2022-10-14 00:00:45,488 __main__ 0/1 INFO :: Iteration=8000, Time=1.600000e+01, dt=2.000000e-03
2022-10-14 00:00:45,615 __main__ 0/1 INFO :: Iteration=8500, Time=1.700000e+01, dt=2.000000e-03
2022-10-14 00:00:45,742 __main__ 0/1 INFO :: Iteration=9000, Time=1.800000e+01, dt=2.000000e-03
2022-10-14 00:00:45,869 __main__ 0/1 INFO :: Iteration=9500, Time=1.900000e+01, dt=2.000000e-03
2022-10-14 00:00:45,999 __main__ 0/1 INFO :: Iteration=10000, Time=2.000000e+01, dt=2.000000e-03
2022-10-14 00:00:46,129 __main__ 0/1 INFO :: Iteration=10500, Time=2.100000e+01, dt=2.000000e-03
2022-10-14 00:00:46,255 __main__ 0/1 INFO :: Iteration=11000, Time=2.200000e+01, dt=2.000000e-03
2022-10-14 00:00:46,383 __main__ 0/1 INFO :: Iteration=11500, Time=2.300000e+01, dt=2.000000e-03
2022-10-14 00:00:46,512 __main__ 0/1 INFO :: Iteration=12000, Time=2.400000e+01, dt=2.000000e-03
2022-10-14 00:00:46,639 __main__ 0/1 INFO :: Iteration=12500, Time=2.500000e+01, dt=2.000000e-03
2022-10-14 00:00:46,766 __main__ 0/1 INFO :: Iteration=13000, Time=2.600000e+01, dt=2.000000e-03
2022-10-14 00:00:46,893 __main__ 0/1 INFO :: Iteration=13500, Time=2.700000e+01, dt=2.000000e-03
2022-10-14 00:00:47,016 __main__ 0/1 INFO :: Iteration=14000, Time=2.800000e+01, dt=2.000000e-03
2022-10-14 00:00:47,139 __main__ 0/1 INFO :: Iteration=14500, Time=2.900000e+01, dt=2.000000e-03
2022-10-14 00:00:47,264 __main__ 0/1 INFO :: Iteration=15000, Time=3.000000e+01, dt=2.000000e-03
2022-10-14 00:00:47,387 __main__ 0/1 INFO :: Iteration=15500, Time=3.100000e+01, dt=2.000000e-03
2022-10-14 00:00:47,510 __main__ 0/1 INFO :: Iteration=16000, Time=3.200000e+01, dt=2.000000e-03
2022-10-14 00:00:47,632 __main__ 0/1 INFO :: Iteration=16500, Time=3.300000e+01, dt=2.000000e-03
2022-10-14 00:00:47,756 __main__ 0/1 INFO :: Iteration=17000, Time=3.400000e+01, dt=2.000000e-03
2022-10-14 00:00:47,879 __main__ 0/1 INFO :: Iteration=17500, Time=3.500000e+01, dt=2.000000e-03
2022-10-14 00:00:48,003 __main__ 0/1 INFO :: Iteration=18000, Time=3.600000e+01, dt=2.000000e-03
2022-10-14 00:00:48,133 __main__ 0/1 INFO :: Iteration=18500, Time=3.700000e+01, dt=2.000000e-03
2022-10-14 00:00:48,261 __main__ 0/1 INFO :: Iteration=19000, Time=3.800000e+01, dt=2.000000e-03
2022-10-14 00:00:48,389 __main__ 0/1 INFO :: Iteration=19500, Time=3.900000e+01, dt=2.000000e-03
2022-10-14 00:00:48,518 __main__ 0/1 INFO :: Iteration=20000, Time=4.000000e+01, dt=2.000000e-03
2022-10-14 00:00:48,647 __main__ 0/1 INFO :: Iteration=20500, Time=4.100000e+01, dt=2.000000e-03
2022-10-14 00:00:48,774 __main__ 0/1 INFO :: Iteration=21000, Time=4.200000e+01, dt=2.000000e-03
2022-10-14 00:00:48,900 __main__ 0/1 INFO :: Iteration=21500, Time=4.300000e+01, dt=2.000000e-03
2022-10-14 00:00:49,027 __main__ 0/1 INFO :: Iteration=22000, Time=4.400000e+01, dt=2.000000e-03
2022-10-14 00:00:49,161 __main__ 0/1 INFO :: Iteration=22500, Time=4.500000e+01, dt=2.000000e-03
2022-10-14 00:00:49,291 __main__ 0/1 INFO :: Iteration=23000, Time=4.600000e+01, dt=2.000000e-03
2022-10-14 00:00:49,417 __main__ 0/1 INFO :: Iteration=23500, Time=4.700000e+01, dt=2.000000e-03
2022-10-14 00:00:49,543 __main__ 0/1 INFO :: Iteration=24000, Time=4.800000e+01, dt=2.000000e-03
2022-10-14 00:00:49,669 __main__ 0/1 INFO :: Iteration=24500, Time=4.900000e+01, dt=2.000000e-03
2022-10-14 00:00:49,795 __main__ 0/1 INFO :: Iteration=25000, Time=5.000000e+01, dt=2.000000e-03
2022-10-14 00:00:49,922 __main__ 0/1 INFO :: Iteration=25500, Time=5.100000e+01, dt=2.000000e-03
2022-10-14 00:00:50,044 __main__ 0/1 INFO :: Iteration=26000, Time=5.200000e+01, dt=2.000000e-03
2022-10-14 00:00:50,164 __main__ 0/1 INFO :: Iteration=26500, Time=5.300000e+01, dt=2.000000e-03
2022-10-14 00:00:50,286 __main__ 0/1 INFO :: Iteration=27000, Time=5.400000e+01, dt=2.000000e-03
2022-10-14 00:00:50,407 __main__ 0/1 INFO :: Iteration=27500, Time=5.500000e+01, dt=2.000000e-03
2022-10-14 00:00:50,526 __main__ 0/1 INFO :: Iteration=28000, Time=5.600000e+01, dt=2.000000e-03
2022-10-14 00:00:50,645 __main__ 0/1 INFO :: Iteration=28500, Time=5.700000e+01, dt=2.000000e-03
2022-10-14 00:00:50,772 __main__ 0/1 INFO :: Iteration=29000, Time=5.800000e+01, dt=2.000000e-03
2022-10-14 00:00:50,894 __main__ 0/1 INFO :: Iteration=29500, Time=5.900000e+01, dt=2.000000e-03
2022-10-14 00:00:51,014 __main__ 0/1 INFO :: Iteration=30000, Time=6.000000e+01, dt=2.000000e-03
2022-10-14 00:00:51,134 __main__ 0/1 INFO :: Iteration=30500, Time=6.100000e+01, dt=2.000000e-03
2022-10-14 00:00:51,254 __main__ 0/1 INFO :: Iteration=31000, Time=6.200000e+01, dt=2.000000e-03
2022-10-14 00:00:51,374 __main__ 0/1 INFO :: Iteration=31500, Time=6.300000e+01, dt=2.000000e-03
2022-10-14 00:00:51,495 __main__ 0/1 INFO :: Iteration=32000, Time=6.400000e+01, dt=2.000000e-03
2022-10-14 00:00:51,615 __main__ 0/1 INFO :: Iteration=32500, Time=6.500000e+01, dt=2.000000e-03
2022-10-14 00:00:51,735 __main__ 0/1 INFO :: Iteration=33000, Time=6.600000e+01, dt=2.000000e-03
2022-10-14 00:00:51,854 __main__ 0/1 INFO :: Iteration=33500, Time=6.700000e+01, dt=2.000000e-03
2022-10-14 00:00:51,974 __main__ 0/1 INFO :: Iteration=34000, Time=6.800000e+01, dt=2.000000e-03
2022-10-14 00:00:52,094 __main__ 0/1 INFO :: Iteration=34500, Time=6.900000e+01, dt=2.000000e-03
2022-10-14 00:00:52,214 __main__ 0/1 INFO :: Iteration=35000, Time=7.000000e+01, dt=2.000000e-03
2022-10-14 00:00:52,335 __main__ 0/1 INFO :: Iteration=35500, Time=7.100000e+01, dt=2.000000e-03
2022-10-14 00:00:52,453 __main__ 0/1 INFO :: Iteration=36000, Time=7.200000e+01, dt=2.000000e-03
2022-10-14 00:00:52,573 __main__ 0/1 INFO :: Iteration=36500, Time=7.300000e+01, dt=2.000000e-03
2022-10-14 00:00:52,692 __main__ 0/1 INFO :: Iteration=37000, Time=7.400000e+01, dt=2.000000e-03
2022-10-14 00:00:52,815 __main__ 0/1 INFO :: Iteration=37500, Time=7.500000e+01, dt=2.000000e-03
2022-10-14 00:00:52,936 __main__ 0/1 INFO :: Iteration=38000, Time=7.600000e+01, dt=2.000000e-03
2022-10-14 00:00:53,055 __main__ 0/1 INFO :: Iteration=38500, Time=7.700000e+01, dt=2.000000e-03
2022-10-14 00:00:53,173 __main__ 0/1 INFO :: Iteration=39000, Time=7.800000e+01, dt=2.000000e-03
2022-10-14 00:00:53,293 __main__ 0/1 INFO :: Iteration=39500, Time=7.900000e+01, dt=2.000000e-03
2022-10-14 00:00:53,412 __main__ 0/1 INFO :: Iteration=40000, Time=8.000000e+01, dt=2.000000e-03
2022-10-14 00:00:53,531 __main__ 0/1 INFO :: Iteration=40500, Time=8.100000e+01, dt=2.000000e-03
2022-10-14 00:00:53,650 __main__ 0/1 INFO :: Iteration=41000, Time=8.200000e+01, dt=2.000000e-03
2022-10-14 00:00:53,769 __main__ 0/1 INFO :: Iteration=41500, Time=8.300000e+01, dt=2.000000e-03
2022-10-14 00:00:53,888 __main__ 0/1 INFO :: Iteration=42000, Time=8.400000e+01, dt=2.000000e-03
2022-10-14 00:00:54,008 __main__ 0/1 INFO :: Iteration=42500, Time=8.500000e+01, dt=2.000000e-03
2022-10-14 00:00:54,126 __main__ 0/1 INFO :: Iteration=43000, Time=8.600000e+01, dt=2.000000e-03
2022-10-14 00:00:54,245 __main__ 0/1 INFO :: Iteration=43500, Time=8.700000e+01, dt=2.000000e-03
2022-10-14 00:00:54,363 __main__ 0/1 INFO :: Iteration=44000, Time=8.800000e+01, dt=2.000000e-03
2022-10-14 00:00:54,484 __main__ 0/1 INFO :: Iteration=44500, Time=8.900000e+01, dt=2.000000e-03
2022-10-14 00:00:54,602 __main__ 0/1 INFO :: Iteration=45000, Time=9.000000e+01, dt=2.000000e-03
2022-10-14 00:00:54,721 __main__ 0/1 INFO :: Iteration=45500, Time=9.100000e+01, dt=2.000000e-03
2022-10-14 00:00:54,839 __main__ 0/1 INFO :: Iteration=46000, Time=9.200000e+01, dt=2.000000e-03
2022-10-14 00:00:54,959 __main__ 0/1 INFO :: Iteration=46500, Time=9.300000e+01, dt=2.000000e-03
2022-10-14 00:00:55,080 __main__ 0/1 INFO :: Iteration=47000, Time=9.400000e+01, dt=2.000000e-03
2022-10-14 00:00:55,198 __main__ 0/1 INFO :: Iteration=47500, Time=9.500000e+01, dt=2.000000e-03
2022-10-14 00:00:55,317 __main__ 0/1 INFO :: Iteration=48000, Time=9.600000e+01, dt=2.000000e-03
2022-10-14 00:00:55,436 __main__ 0/1 INFO :: Iteration=48500, Time=9.700000e+01, dt=2.000000e-03
2022-10-14 00:00:55,554 __main__ 0/1 INFO :: Iteration=49000, Time=9.800000e+01, dt=2.000000e-03
2022-10-14 00:00:55,674 __main__ 0/1 INFO :: Iteration=49500, Time=9.900000e+01, dt=2.000000e-03
2022-10-14 00:00:55,794 __main__ 0/1 INFO :: Iteration=50000, Time=1.000000e+02, dt=2.000000e-03
2022-10-14 00:00:55,913 __main__ 0/1 INFO :: Iteration=50500, Time=1.010000e+02, dt=2.000000e-03
2022-10-14 00:00:56,030 __main__ 0/1 INFO :: Iteration=51000, Time=1.020000e+02, dt=2.000000e-03
2022-10-14 00:00:56,147 __main__ 0/1 INFO :: Iteration=51500, Time=1.030000e+02, dt=2.000000e-03
2022-10-14 00:00:56,266 __main__ 0/1 INFO :: Iteration=52000, Time=1.040000e+02, dt=2.000000e-03
2022-10-14 00:00:56,385 __main__ 0/1 INFO :: Iteration=52500, Time=1.050000e+02, dt=2.000000e-03
2022-10-14 00:00:56,503 __main__ 0/1 INFO :: Iteration=53000, Time=1.060000e+02, dt=2.000000e-03
2022-10-14 00:00:56,621 __main__ 0/1 INFO :: Iteration=53500, Time=1.070000e+02, dt=2.000000e-03
2022-10-14 00:00:56,740 __main__ 0/1 INFO :: Iteration=54000, Time=1.080000e+02, dt=2.000000e-03
2022-10-14 00:00:56,858 __main__ 0/1 INFO :: Iteration=54500, Time=1.090000e+02, dt=2.000000e-03
2022-10-14 00:00:56,976 __main__ 0/1 INFO :: Iteration=55000, Time=1.100000e+02, dt=2.000000e-03
2022-10-14 00:00:57,096 __main__ 0/1 INFO :: Iteration=55500, Time=1.110000e+02, dt=2.000000e-03
2022-10-14 00:00:57,215 __main__ 0/1 INFO :: Iteration=56000, Time=1.120000e+02, dt=2.000000e-03
2022-10-14 00:00:57,336 __main__ 0/1 INFO :: Iteration=56500, Time=1.130000e+02, dt=2.000000e-03
2022-10-14 00:00:57,456 __main__ 0/1 INFO :: Iteration=57000, Time=1.140000e+02, dt=2.000000e-03
2022-10-14 00:00:57,577 __main__ 0/1 INFO :: Iteration=57500, Time=1.150000e+02, dt=2.000000e-03
2022-10-14 00:00:57,697 __main__ 0/1 INFO :: Iteration=58000, Time=1.160000e+02, dt=2.000000e-03
2022-10-14 00:00:57,816 __main__ 0/1 INFO :: Iteration=58500, Time=1.170000e+02, dt=2.000000e-03
2022-10-14 00:00:57,936 __main__ 0/1 INFO :: Iteration=59000, Time=1.180000e+02, dt=2.000000e-03
2022-10-14 00:00:58,055 __main__ 0/1 INFO :: Iteration=59500, Time=1.190000e+02, dt=2.000000e-03
2022-10-14 00:00:58,175 __main__ 0/1 INFO :: Iteration=60000, Time=1.200000e+02, dt=2.000000e-03
2022-10-14 00:00:58,294 __main__ 0/1 INFO :: Iteration=60500, Time=1.210000e+02, dt=2.000000e-03
2022-10-14 00:00:58,415 __main__ 0/1 INFO :: Iteration=61000, Time=1.220000e+02, dt=2.000000e-03
2022-10-14 00:00:58,535 __main__ 0/1 INFO :: Iteration=61500, Time=1.230000e+02, dt=2.000000e-03
2022-10-14 00:00:58,654 __main__ 0/1 INFO :: Iteration=62000, Time=1.240000e+02, dt=2.000000e-03
2022-10-14 00:00:58,773 __main__ 0/1 INFO :: Iteration=62500, Time=1.250000e+02, dt=2.000000e-03
2022-10-14 00:00:58,893 __main__ 0/1 INFO :: Iteration=63000, Time=1.260000e+02, dt=2.000000e-03
2022-10-14 00:00:59,012 __main__ 0/1 INFO :: Iteration=63500, Time=1.270000e+02, dt=2.000000e-03
2022-10-14 00:00:59,130 __main__ 0/1 INFO :: Iteration=64000, Time=1.280000e+02, dt=2.000000e-03
2022-10-14 00:00:59,247 __main__ 0/1 INFO :: Iteration=64500, Time=1.290000e+02, dt=2.000000e-03
2022-10-14 00:00:59,366 __main__ 0/1 INFO :: Iteration=65000, Time=1.300000e+02, dt=2.000000e-03
2022-10-14 00:00:59,484 __main__ 0/1 INFO :: Iteration=65500, Time=1.310000e+02, dt=2.000000e-03
2022-10-14 00:00:59,601 __main__ 0/1 INFO :: Iteration=66000, Time=1.320000e+02, dt=2.000000e-03
2022-10-14 00:00:59,719 __main__ 0/1 INFO :: Iteration=66500, Time=1.330000e+02, dt=2.000000e-03
2022-10-14 00:00:59,836 __main__ 0/1 INFO :: Iteration=67000, Time=1.340000e+02, dt=2.000000e-03
2022-10-14 00:00:59,956 __main__ 0/1 INFO :: Iteration=67500, Time=1.350000e+02, dt=2.000000e-03
2022-10-14 00:01:00,077 __main__ 0/1 INFO :: Iteration=68000, Time=1.360000e+02, dt=2.000000e-03
2022-10-14 00:01:00,197 __main__ 0/1 INFO :: Iteration=68500, Time=1.370000e+02, dt=2.000000e-03
2022-10-14 00:01:00,317 __main__ 0/1 INFO :: Iteration=69000, Time=1.380000e+02, dt=2.000000e-03
2022-10-14 00:01:00,436 __main__ 0/1 INFO :: Iteration=69500, Time=1.390000e+02, dt=2.000000e-03
2022-10-14 00:01:00,554 __main__ 0/1 INFO :: Iteration=70000, Time=1.400000e+02, dt=2.000000e-03
2022-10-14 00:01:00,672 __main__ 0/1 INFO :: Iteration=70500, Time=1.410000e+02, dt=2.000000e-03
2022-10-14 00:01:00,789 __main__ 0/1 INFO :: Iteration=71000, Time=1.420000e+02, dt=2.000000e-03
2022-10-14 00:01:00,907 __main__ 0/1 INFO :: Iteration=71500, Time=1.430000e+02, dt=2.000000e-03
2022-10-14 00:01:01,027 __main__ 0/1 INFO :: Iteration=72000, Time=1.440000e+02, dt=2.000000e-03
2022-10-14 00:01:01,147 __main__ 0/1 INFO :: Iteration=72500, Time=1.450000e+02, dt=2.000000e-03
2022-10-14 00:01:01,266 __main__ 0/1 INFO :: Iteration=73000, Time=1.460000e+02, dt=2.000000e-03
2022-10-14 00:01:01,384 __main__ 0/1 INFO :: Iteration=73500, Time=1.470000e+02, dt=2.000000e-03
2022-10-14 00:01:01,502 __main__ 0/1 INFO :: Iteration=74000, Time=1.480000e+02, dt=2.000000e-03
2022-10-14 00:01:01,621 __main__ 0/1 INFO :: Iteration=74500, Time=1.490000e+02, dt=2.000000e-03
2022-10-14 00:01:01,739 __main__ 0/1 INFO :: Iteration=75000, Time=1.500000e+02, dt=2.000000e-03
2022-10-14 00:01:01,856 __main__ 0/1 INFO :: Iteration=75500, Time=1.510000e+02, dt=2.000000e-03
2022-10-14 00:01:01,974 __main__ 0/1 INFO :: Iteration=76000, Time=1.520000e+02, dt=2.000000e-03
2022-10-14 00:01:02,093 __main__ 0/1 INFO :: Iteration=76500, Time=1.530000e+02, dt=2.000000e-03
2022-10-14 00:01:02,211 __main__ 0/1 INFO :: Iteration=77000, Time=1.540000e+02, dt=2.000000e-03
2022-10-14 00:01:02,330 __main__ 0/1 INFO :: Iteration=77500, Time=1.550000e+02, dt=2.000000e-03
2022-10-14 00:01:02,448 __main__ 0/1 INFO :: Iteration=78000, Time=1.560000e+02, dt=2.000000e-03
2022-10-14 00:01:02,567 __main__ 0/1 INFO :: Iteration=78500, Time=1.570000e+02, dt=2.000000e-03
2022-10-14 00:01:02,684 __main__ 0/1 INFO :: Iteration=79000, Time=1.580000e+02, dt=2.000000e-03
2022-10-14 00:01:02,802 __main__ 0/1 INFO :: Iteration=79500, Time=1.590000e+02, dt=2.000000e-03
2022-10-14 00:01:02,919 __main__ 0/1 INFO :: Iteration=80000, Time=1.600000e+02, dt=2.000000e-03
2022-10-14 00:01:03,037 __main__ 0/1 INFO :: Iteration=80500, Time=1.610000e+02, dt=2.000000e-03
2022-10-14 00:01:03,155 __main__ 0/1 INFO :: Iteration=81000, Time=1.620000e+02, dt=2.000000e-03
2022-10-14 00:01:03,273 __main__ 0/1 INFO :: Iteration=81500, Time=1.630000e+02, dt=2.000000e-03
2022-10-14 00:01:03,392 __main__ 0/1 INFO :: Iteration=82000, Time=1.640000e+02, dt=2.000000e-03
2022-10-14 00:01:03,510 __main__ 0/1 INFO :: Iteration=82500, Time=1.650000e+02, dt=2.000000e-03
2022-10-14 00:01:03,628 __main__ 0/1 INFO :: Iteration=83000, Time=1.660000e+02, dt=2.000000e-03
2022-10-14 00:01:03,745 __main__ 0/1 INFO :: Iteration=83500, Time=1.670000e+02, dt=2.000000e-03
2022-10-14 00:01:03,862 __main__ 0/1 INFO :: Iteration=84000, Time=1.680000e+02, dt=2.000000e-03
2022-10-14 00:01:03,980 __main__ 0/1 INFO :: Iteration=84500, Time=1.690000e+02, dt=2.000000e-03
2022-10-14 00:01:04,098 __main__ 0/1 INFO :: Iteration=85000, Time=1.700000e+02, dt=2.000000e-03
2022-10-14 00:01:04,215 __main__ 0/1 INFO :: Iteration=85500, Time=1.710000e+02, dt=2.000000e-03
2022-10-14 00:01:04,333 __main__ 0/1 INFO :: Iteration=86000, Time=1.720000e+02, dt=2.000000e-03
2022-10-14 00:01:04,450 __main__ 0/1 INFO :: Iteration=86500, Time=1.730000e+02, dt=2.000000e-03
2022-10-14 00:01:04,568 __main__ 0/1 INFO :: Iteration=87000, Time=1.740000e+02, dt=2.000000e-03
2022-10-14 00:01:04,685 __main__ 0/1 INFO :: Iteration=87500, Time=1.750000e+02, dt=2.000000e-03
2022-10-14 00:01:04,802 __main__ 0/1 INFO :: Iteration=88000, Time=1.760000e+02, dt=2.000000e-03
2022-10-14 00:01:04,922 __main__ 0/1 INFO :: Iteration=88500, Time=1.770000e+02, dt=2.000000e-03
2022-10-14 00:01:05,041 __main__ 0/1 INFO :: Iteration=89000, Time=1.780000e+02, dt=2.000000e-03
2022-10-14 00:01:05,159 __main__ 0/1 INFO :: Iteration=89500, Time=1.790000e+02, dt=2.000000e-03
2022-10-14 00:01:05,277 __main__ 0/1 INFO :: Iteration=90000, Time=1.800000e+02, dt=2.000000e-03
2022-10-14 00:01:05,396 __main__ 0/1 INFO :: Iteration=90500, Time=1.810000e+02, dt=2.000000e-03
2022-10-14 00:01:05,515 __main__ 0/1 INFO :: Iteration=91000, Time=1.820000e+02, dt=2.000000e-03
2022-10-14 00:01:05,635 __main__ 0/1 INFO :: Iteration=91500, Time=1.830000e+02, dt=2.000000e-03
2022-10-14 00:01:05,755 __main__ 0/1 INFO :: Iteration=92000, Time=1.840000e+02, dt=2.000000e-03
2022-10-14 00:01:05,873 __main__ 0/1 INFO :: Iteration=92500, Time=1.850000e+02, dt=2.000000e-03
2022-10-14 00:01:05,992 __main__ 0/1 INFO :: Iteration=93000, Time=1.860000e+02, dt=2.000000e-03
2022-10-14 00:01:06,113 __main__ 0/1 INFO :: Iteration=93500, Time=1.870000e+02, dt=2.000000e-03
2022-10-14 00:01:06,233 __main__ 0/1 INFO :: Iteration=94000, Time=1.880000e+02, dt=2.000000e-03
2022-10-14 00:01:06,351 __main__ 0/1 INFO :: Iteration=94500, Time=1.890000e+02, dt=2.000000e-03
2022-10-14 00:01:06,471 __main__ 0/1 INFO :: Iteration=95000, Time=1.900000e+02, dt=2.000000e-03
2022-10-14 00:01:06,590 __main__ 0/1 INFO :: Iteration=95500, Time=1.910000e+02, dt=2.000000e-03
2022-10-14 00:01:06,708 __main__ 0/1 INFO :: Iteration=96000, Time=1.920000e+02, dt=2.000000e-03
2022-10-14 00:01:06,827 __main__ 0/1 INFO :: Iteration=96500, Time=1.930000e+02, dt=2.000000e-03
2022-10-14 00:01:06,946 __main__ 0/1 INFO :: Iteration=97000, Time=1.940000e+02, dt=2.000000e-03
2022-10-14 00:01:07,066 __main__ 0/1 INFO :: Iteration=97500, Time=1.950000e+02, dt=2.000000e-03
2022-10-14 00:01:07,186 __main__ 0/1 INFO :: Iteration=98000, Time=1.960000e+02, dt=2.000000e-03
2022-10-14 00:01:07,305 __main__ 0/1 INFO :: Iteration=98500, Time=1.970000e+02, dt=2.000000e-03
2022-10-14 00:01:07,422 __main__ 0/1 INFO :: Iteration=99000, Time=1.980000e+02, dt=2.000000e-03
2022-10-14 00:01:07,540 __main__ 0/1 INFO :: Iteration=99500, Time=1.990000e+02, dt=2.000000e-03
2022-10-14 00:01:07,659 __main__ 0/1 INFO :: Iteration=100000, Time=2.000000e+02, dt=2.000000e-03
2022-10-14 00:01:07,659 solvers 0/1 INFO :: Simulation stop time reached.
No description has been provided for this image

The Poisson equation with mixed boundary conditions¶

  • A 2D problem that is a relative of the Laplace equation, but with mixed boundary conditions

$$ \nabla^2 u(\mathbf{r}) = f(\mathbf{r}) $$

  • The top and bottom boundaries are Dirichlet, and the left and right boundaries are periodic

  • The current version of Dedalus requires boundary conditions to be specified via the tau method

  • This example is adapted from the Dedalus tutorials

In [25]:
import numpy as np
import matplotlib.pyplot as plt
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)


# Parameters
Lx, Ly = 2*np.pi, np.pi
Nx, Ny = 256, 128
dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))
ybasis = d3.Chebyshev(coords['y'], size=Ny, bounds=(0, Ly))

## Fields
u = dist.Field(name='u', bases=(xbasis, ybasis))
tau_1 = dist.Field(name='tau_1', bases=xbasis)
tau_2 = dist.Field(name='tau_2', bases=xbasis)

## Forcing
x, y = dist.local_grids(xbasis, ybasis)
f = dist.Field(bases=(xbasis, ybasis))
g = dist.Field(bases=xbasis)
h = dist.Field(bases=xbasis)
f.fill_random('g', seed=40)
f.low_pass_filter(shape=(64, 32))
g['g'] = np.sin(8*x) * 0.025
h['g'] = 0

# Substitutions
dy = lambda A: d3.Differentiate(A, coords['y'])
lift_basis = ybasis.derivative_basis(2)
lift = lambda A, n: d3.Lift(A, lift_basis, n)

## Specify PDE Problem. Notice that we need to add explicit boundary conditions
## along the y direction, since we have Dirichlet BCs there.
problem = d3.LBVP([u, tau_1, tau_2], namespace=locals())
# problem.add_equation("lap(u) = f")
problem.add_equation("lap(u) + lift(tau_1,-1) + lift(tau_2,-2) = f")
problem.add_equation("u(y=0) = g")
problem.add_equation("dy(u)(y=Ly) = h")

## Build and runSolver
solver = problem.build_solver()
solver.solve()

# Gather global data
x = xbasis.global_grid()
y = ybasis.global_grid()
ug = u.allgather_data('g')

# Plot
if dist.comm.rank == 0:
    plt.figure(figsize=(6, 4))
    plt.pcolormesh(x.ravel(), y.ravel(), ug.T, cmap='viridis', shading='gouraud', rasterized=True)
    plt.gca().set_aspect('equal')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title("Randomly forced Poisson equation")
    plt.tight_layout()
    plt.savefig('poisson.pdf')
    plt.savefig('poisson.png', dpi=200)

2022-10-13 23:41:55,261 subsystems 0/1 INFO :: Building subproblem matrices 1/128 (~1%) Elapsed: 0s, Remaining: 4s, Rate: 3.3e+01/s
2022-10-13 23:41:55,411 subsystems 0/1 INFO :: Building subproblem matrices 13/128 (~10%) Elapsed: 0s, Remaining: 2s, Rate: 7.2e+01/s
2022-10-13 23:41:55,568 subsystems 0/1 INFO :: Building subproblem matrices 26/128 (~20%) Elapsed: 0s, Remaining: 1s, Rate: 7.7e+01/s
2022-10-13 23:41:55,731 subsystems 0/1 INFO :: Building subproblem matrices 39/128 (~30%) Elapsed: 1s, Remaining: 1s, Rate: 7.8e+01/s
2022-10-13 23:41:55,912 subsystems 0/1 INFO :: Building subproblem matrices 52/128 (~41%) Elapsed: 1s, Remaining: 1s, Rate: 7.6e+01/s
2022-10-13 23:41:56,090 subsystems 0/1 INFO :: Building subproblem matrices 65/128 (~51%) Elapsed: 1s, Remaining: 1s, Rate: 7.6e+01/s
2022-10-13 23:41:56,268 subsystems 0/1 INFO :: Building subproblem matrices 78/128 (~61%) Elapsed: 1s, Remaining: 1s, Rate: 7.5e+01/s
2022-10-13 23:41:56,436 subsystems 0/1 INFO :: Building subproblem matrices 91/128 (~71%) Elapsed: 1s, Remaining: 0s, Rate: 7.5e+01/s
2022-10-13 23:41:56,597 subsystems 0/1 INFO :: Building subproblem matrices 104/128 (~81%) Elapsed: 1s, Remaining: 0s, Rate: 7.6e+01/s
2022-10-13 23:41:56,768 subsystems 0/1 INFO :: Building subproblem matrices 117/128 (~91%) Elapsed: 2s, Remaining: 0s, Rate: 7.6e+01/s
2022-10-13 23:41:56,905 subsystems 0/1 INFO :: Building subproblem matrices 128/128 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 7.6e+01/s
No description has been provided for this image
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: