Gray-Scott Equation¶

Gray-Scott Equation: is a reaction-diffusion system that describes the interaction between two chemical substances, denoted as $\phi_0$ and $\phi_1$. The equation is given by:

$$ \begin{aligned} \frac{\partial \phi_0}{\partial t} &= \nu_0 \Delta \phi_0 - \phi_0 \phi_1^2 + f (1 - \phi_0) \\ \frac{\partial \phi_1}{\partial t} &= \nu_1 \Delta \phi_1 + \phi_0 \phi_1^2 - (f + k) \phi_1 \end{aligned} $$

where $\nu_0$ and $\nu_1$ are the diffusion coefficients, $f$ is the feed rate, and $k$ is the kill rate.

Solution in 2D¶

In [1]:

from torchfsm.operator import ChannelWisedDiffusion, GrayScottSource, Operator

def GrayScott(nu_0: float, nu_1:float, feed_rate:float, kill_rate:float) -> Operator:
    op = ChannelWisedDiffusion([nu_0,nu_1])+ GrayScottSource(feed_rate, kill_rate)
    op.set_de_aliasing_rate(0.5) # GrayScottSource requires lower de-aliasing rate due to the cubic nonlinearity
    return op

gray_scott= GrayScott(
    nu_0=2e-5, 
    nu_1=1e-5, 
    feed_rate=0.04, 
    kill_rate=0.06
)

from torchfsm.operator import ChannelWisedDiffusion, GrayScottSource, Operator def GrayScott(nu_0: float, nu_1:float, feed_rate:float, kill_rate:float) -> Operator: op = ChannelWisedDiffusion([nu_0,nu_1])+ GrayScottSource(feed_rate, kill_rate) op.set_de_aliasing_rate(0.5) # GrayScottSource requires lower de-aliasing rate due to the cubic nonlinearity return op gray_scott= GrayScott( nu_0=2e-5, nu_1=1e-5, feed_rate=0.04, kill_rate=0.06 )

In [2]:

import torch
from torchfsm.mesh import MeshGrid
from torchfsm.plot import plot_traj
from torchfsm.traj_recorder import AutoRecorder, IntervalController
from torchfsm.field import random_gaussian_blobs
device='cuda' if torch.cuda.is_available() else 'cpu'
L=1.0; N=128; 

import torch from torchfsm.mesh import MeshGrid from torchfsm.plot import plot_traj from torchfsm.traj_recorder import AutoRecorder, IntervalController from torchfsm.field import random_gaussian_blobs device='cuda' if torch.cuda.is_available() else 'cpu' L=1.0; N=128;

In [7]:

mesh=MeshGrid([(0,L,N)]*2,device=device)
x=mesh.bc_mesh_grid()
u_0=torch.cat(
    [1-random_gaussian_blobs(mesh),random_gaussian_blobs(mesh)]
    ,dim=1
)
traj=gray_scott.integrate(
    u_0=u_0,
    mesh=mesh,
    dt=1,
    step=1200,
    trajectory_recorder=AutoRecorder(control_func=IntervalController(30),include_initial_state=True),
)

mesh=MeshGrid([(0,L,N)]*2,device=device) x=mesh.bc_mesh_grid() u_0=torch.cat( [1-random_gaussian_blobs(mesh),random_gaussian_blobs(mesh)] ,dim=1 ) traj=gray_scott.integrate( u_0=u_0, mesh=mesh, dt=1, step=1200, trajectory_recorder=AutoRecorder(control_func=IntervalController(30),include_initial_state=True), )

In [9]:

plot_traj(traj)

plot_traj(traj)

Out[9]:

Solution in 3D¶

In [5]:

mesh=MeshGrid([(0,L,N)]*3,device=device)
x=mesh.bc_mesh_grid()
u_0=torch.cat(
    [1-random_gaussian_blobs(mesh),random_gaussian_blobs(mesh)]
    ,dim=1
)
traj=gray_scott.integrate(
    u_0=u_0,
    mesh=mesh,
    dt=1,
    step=1200,
    trajectory_recorder=AutoRecorder(control_func=IntervalController(30),include_initial_state=True),
)

mesh=MeshGrid([(0,L,N)]*3,device=device) x=mesh.bc_mesh_grid() u_0=torch.cat( [1-random_gaussian_blobs(mesh),random_gaussian_blobs(mesh)] ,dim=1 ) traj=gray_scott.integrate( u_0=u_0, mesh=mesh, dt=1, step=1200, trajectory_recorder=AutoRecorder(control_func=IntervalController(30),include_initial_state=True), )

In [6]:

plot_traj(traj)

plot_traj(traj)

Out[6]: