Simple change in boundary conditions breaks code

[HennesHajduk]

Hi, I am new to dedalus and spectral methods in general but I do have experience with PDEs and numerics. I was able to set up the 2D system I want to solve using RealFourier bases in x and y. Now I am trying to port my script to cases involving non-periodic boundaries.

To be honest, I only have a rough idea of what the tau-method is, to me it feels like adding a penalization term in the PDE to weakly enforce the BC. Anyway to get a better understanding, I wanted to implement progressively more complicated BCs to finally be able to setup the correct BCs for the model that I am actually interested in.

As a first step, I simply wanted to change the provided Poisson equation script to be periodic in x instead of y-direction, which works in serial but in parallel I get the error: ValueError: Problem is coupled along distributed dimensions: (0,)

The poisson script in particular seems to be prone to fail in parallel upon minor changes that do not make it fail in serial.

Here's my code, thanks in advance!

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.Chebyshev(coords['x'], size=Nx, bounds=(0, Lx))
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(0, Ly))

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

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

# Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
lift_basis = xbasis.derivative_basis(2)
lift = lambda A, n: d3.Lift(A, lift_basis, n)

# Problem
problem = d3.LBVP([u, tau_1, tau_2], namespace=locals())
problem.add_equation("lap(u) + lift(tau_1,-1) + lift(tau_2,-2) = f")
problem.add_equation("u(x=0) = g")
problem.add_equation("dx(u)(x=Lx) = h")

# Solver
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.show()

[kburns]

Chebyshev dimensions need to be on the trailing coordinates (in the order the coordinates are entered in the coordinate system) to run in parallel. You can still call that coordinate x, but you should do d3.CartesianCoordinates('y', 'x'). Please give that a try, and if there are still issues please follow up on the user mailing list, which is a better spot for getting help with model design and debugging.