Noise in compiled custom class

[Gregogc]

Hello,

I want to know what would be the best way to implement additive noise in a PDE using the compiled custom class. In particular, I wrote a Ginzburg-Landau model inspired on the example of the Brusselator in the documentation.

import numba as nb
import numpy as np
import os
from pde import FieldCollection, PDEBase, PlotTracker, ScalarField, UnitGrid, MemoryStorage, FileStorage
import subprocess
from pde.grids.operators.cartesian import _make_derivative, _make_derivative2

class GLPDE(PDEBase):
    """GL, whatever"""
    
    explicit_time_dependence = True
    
    def __init__(self, muf=0.1,mui=-1.0,ts=0,bc="natural"):
        super().__init__()
        self.muf = muf
        self.mui = mui
        self.bc = bc  # Boundary condition
        self.ts = ts

    def get_initial_state(self, grid):
         """Prepare a useful initial state"""
         u = ScalarField.random_uniform(grid,vmin=-1,vmax=1, label="Field $u$")
         v = ScalarField.random_uniform(grid,vmin=-1,vmax=1, label="Field $v$")
         return FieldCollection([u, v])

    def evolution_rate(self, state, t=0):
        """Pure Python implementation of the PDE"""
        u, v = state
        rhs = state.copy()
        mui,muf = self.mui,self.muf
        ts=self.ts
        rhs[0] = (muf * float(t>ts+0.0) + mui*float(t<ts+0.0)) * u - u * (u * u + v * v) + u.laplace(self.bc) 
        rhs[1] = (muf * float(t>ts+0.0) + mui*float(t<ts+0.0)) * v - v * (u * u + v * v) + v.laplace(self.bc)

        return rhs

    def _make_pde_rhs_numba(self, state,t=0):
        """Numba-compiled implementation of the PDE"""
        mui,muf = self.mui,self.muf
        ts=self.ts
        laplace = state.grid.make_operator("laplace", bc=self.bc)

        @nb.njit
        def pde_rhs(state_data, t=t):
            u = state_data[0]
            v = state_data[1]
            rate = np.empty_like(state_data)
            rate[0] = (muf * (t>ts+0.0) + mui*(t<ts+0.0)) * u - u * (u * u + v * v) +  laplace(u) 
            rate[1] = (muf * (t>ts+0.0) + mui*(t<ts+0.0)) * v - v * (u * u + v * v) +  laplace(v)
            return rate
        return pde_rhs

# Initialize state
Simsize=64
grid = UnitGrid([Simsize, Simsize])
mui=-1.0
muf=0.5

eqpos = GLPDE(muf= muf, mui=mui ,ts=0.5)
state = FieldCollection([ScalarField.random_uniform(grid,vmin=-0.1 ,vmax=0.1, label="Field $u$"), ScalarField.random_uniform(grid,vmin=-0.1 ,vmax=0.1, label="Field $v$")])

eqpos.solve(state, t_range=100, dt=1e-2, tracker=PlotTracker(interrupts=1),adaptive=True)     

The Ginzburg-Landau model reads \partial_t A=\mu A - AAA+\nabla^2 A with A=u+iv a complex field. I want to implement \partial_t A=\mu A - AAA+\nabla^2 A + \eta, where \eta is a uniform random field (no temporal or spatial correlation) in each real and imaginary part.

When trying to generate a random numpy matrix and adding to the expressions, the compiled and Python results differ, which raises an error. A workaround would be to deactivate the check for the Python and compiled methods, but I would like to know if there is a way to maintain the interpreter-compiler checking and add noise.

Thanks in advance and congratulations for making such a good PDE package!

[david-zwicker]

Noise cannot simply be added in the right hand side of the partial differential equation, since stochastic calculus of partial differential equations is complicated. However, additive noise should be supported by the base class PDEBase. You can simply set the noise variance using the noise argument when you call super().__init__(), e.g., super().__init__(noise=1). You can also give two numbers if you want different variance for the different fields.

[Gregogc]

Thanks a lot. This works perfect.