Error in solving PDE with a heterogeneous coefficient

[noirchen]

Hi, I am trying to simulate a PDE on [-50, 50]

$$ \partial_t u = u_{xx} + (a(x)u)_x $$

with BC $u(\pm50)=0$ and IC

$$ u(0,x)=\frac{1}{1+(x-5)^2} $$

The coefficient $a(x)=3v^2-2v$ with $v$ satisfies

$$ \frac{1}{v}+\ln\left\vert\frac{1-v}{v}\right\vert=x $$

My code is

from pde import PDE, CartesianGrid, MemoryStorage, ScalarField, plot_kymograph
import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

# the solution mesh
grid = CartesianGrid([[-50, 50]], 200)


def arrayfun(func):
    """
    Decorator to make a function that takes an array as input and return an array as output.
    """

    def wrapper(a):
        return np.asarray(list(map(func, a))) if isinstance(
            a, np.ndarray) else np.array([func(a)])

    return wrapper


@arrayfun
def degwave(x):
    if x > 2:
        guess = 1 / x
    elif x > -2.5:
        guess = 0.6
    else:
        guess = 1 - np.exp(-2) * np.exp(x)

    return fsolve(lambda v: 1 / v + np.log(abs(
        (1 - v) / v)) - x, guess)[0] if x >= -35 else 1


v = degwave(grid.coordinate_arrays[0])
a = ScalarField(grid, data=3 * v**2 - 2 * v)

# definition of the PDE
eq = PDE({"u": "laplace(u)-a*gradient(u)"}, consts={'a': a}, bc={"value": 0})

ic = ScalarField.from_expression(grid, "1/(1+(x-5)**2)")
storage = MemoryStorage()
result = eq.solve(ic, 10, dt=0.001, tracker=storage.tracker(0.1))
plot_kymograph(storage)

The code complains when initializing the stepper

controller.py:164, in Controller.run(self, state, dt)
162 # initialize the stepper
163 jit_count_init = int(JIT_COUNT)
--> 164 stepper = self.solver.make_stepper(state=state, dt=dt)
165 self.diagnostics["jit_count"] = {
166 "make_stepper": int(JIT_COUNT) - jit_count_init
167 }
168 jit_count_after_init = int(JIT_COUNT)

Any help is appreciated!

[david-zwicker]

The problem is a semantic problem in the definition of the PDE. The way it is written, you're adding the Laplacian of a scalar field to the gradient of the same field, which are fields of different rank and thus cannot work. Instead of taking the gradient, you need to take a derivative in x-direction explicitly. This is as simple as replacing the definition of your PDE by

eq = PDE({"u": "laplace(u) - a*d_dx(u)"}, consts={'a': a}, bc={"value": 0})

Note that a here is not included in the derivative, which you might or might not want.

[noirchen]

Thanks David. One more thing, I found that when solving the PDEs, often I need to explicitly set the value dt <= Δx^2/2, and adaptive=True just causes oscillations. I wonder when should I use adaptive=True.

[david-zwicker]

Yes, this is unfortunately hard to avoid in non-linear PDEs. If you want to use the adaptive stepping, you could try to reduce the tolerance (by passing tolerance=1e-6. However, if you have a clear stability criterion, it might just be best (and more efficient!) to use a fixed time step.