# Example usage of rdsolver

(c) 2018 Justin Bois. This work is licensed under a [Creative Commons Attribution License CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/). All code contained herein is licensed under an [MIT license](https://opensource.org/licenses/MIT).

`rdsolver` solves the following system of PDEs on a 2D Cartesian domain with periodic boundary conditions. The governing equations are

\begin{align}
\partial_t c_i = D_i(\partial_x^2 + \partial_y^2) c_i + \beta_i + \sum_j\gamma_{ij} c_j + f_i(\mathbf{c}),
\end{align}

where summation over like indices is not implied; all summation is explicit. The last term represents all nonlinear chemical reactions. The $\beta_i + \gamma_{ij}c_j$ terms are linear chemical dynamics. Note that we assume a diagonal constant diffusion tensor.

To specify the problem, the user needs to supply: 

* The physical dimension of the system, which we will call $\mathbf{L} = (L_x, L_y)$
* The number of grid points $\mathbf{n} = (n_x, n_y)$
* The desired time points $t_0, t_1, \ldots$
* The initial concentration profiles of all species, $\mathbf{c}^0(x, y)$
* The values of the parameters $D_i$, $\beta_i$, and $\gamma_{ij}$.
* The function $f_i$ and any parametric arguments that need to be passed to it.

Here, I present an example of how to use `rdsolver`. To learn more about it and installation instructions, see the [README file](https://github.com/justinbois/rdsolver#reaction-diffusion-solver).

## Necessary imports

We work with `rdsolver` in a Jupyter notebook (recommended), you need to import `rdsolver` and `bokeh`, being sure to call `bokeh.io.output_notebook()` to enable interactive plotting in the Jupyter notebook. And it's pretty much automatic to import NumPy!


```python
%load_ext autoreload
%autoreload 2

import numpy as np
import numba

import rdsolver as rd

import bokeh
import bokeh.io
bokeh.io.output_notebook()
```



<div class="bk-root">
    <a href="https://bokeh.org" target="_blank" class="bk-logo bk-logo-small bk-logo-notebook"></a>
    <span id="1001">Loading BokehJS ...</span>
</div>




## The activator-substrate depletion model (ASDM)

The ASDM is a classic system that gives Turing patterns. In its simplest form, the dimensionless equations are

\begin{align}
\partial_t a &= d(\partial_x^2 + \partial_y^2)a + a^2s - a \\[0.5em]
\partial_t s &= (\partial_x^2 + \partial_y^2)s + \mu(1 - a^2s).
\end{align}

Thus, we have

\begin{align}
D &= (d, 1) \\[0.5em]
\beta &= (0, \mu) \\[0.5em]
\gamma &= \begin{pmatrix}
-1 & 0 \\
0 & 0
\end{pmatrix}.
\end{align}

We start by specifying the easy stuff, the physical size of the system, the number of grid points, and the time points we want. Because this particular system does not have very sharp gradients and we are not using a very large physical space, we do not need many grid points at all. Because `rdsolver` uses spectral methods, we can have very accurate calculations even with few grid points. We will use as 32 $\times$ 32 grid here.


```python
# Physical size of system
L = (10, 10)

# Number of grid points in x and y (can often be small with spec. meth.)
n = (32, 32)

# Specify times points we want
t = np.linspace(0, 250, 100)
```

Now, we need to define our parameters. We'll start with $D$, $\beta$, and $\gamma$, choosing $d = 0.05$ and $\mu = 1.4$.


```python
d = 0.05
mu = 1.4

D = (d, 1)
beta = (0, mu)
gamma = np.array([[-1, 0], [0, 0]])
```

Now, we need to define our nonlinear function $f$. This function has call signature `f(u, t, *f_args)`. The first argument, `u`, is an array containing the concentrations that can be unpacked as `a, s = u`. Note, however, that unpacking 3D arrays like this is not yet supported in Numba, so you should use the more verbose version below for Numba'd nonlinear functions.

Next, the function `f` takes the current time as an input. For the ASDM, and indeed for many R-D applications, the nonlinear chemical dynamics do not explicitly depend on time. Finally, `f_args` is a tuple containing any other arguments the function `f` needs.

The function must return an array the same shape as the input array `u` that gives the nonlinear terms in the dynamics. This is most easily accomplished using the `np.stack()` function.


```python
@numba.jit(nopython=True)
def f(u, t, mu):
    """Nonlinear terms for ASDM"""
    a = u[0,:,:]
    s = u[1,:,:]
    return np.stack((a**2 * s, -mu * a**2 * s))
```

We also have to specify the arguments that need to be passed to `f` as a tuple.


```python
# Specify the arguments that need to be passed to f
f_args = (mu, )
```

Now, we need to specify the initial conditions. The initial conditions must be a three-dimensional array of shape $(n_s, n_x, n_y)$, where $n_s$ is the number of chemical species. For convenience, you can specify a homogeneous steady state and use the `rd.initial_condition()` function that will generate an initial condition that is a small perturbation about the specified steady state.


```python
# Homogenous steady state in activator and substrate to perturb
uniform_conc = (1, 1)

# Generate perturbed initial condition
np.random.seed(42)
c0 = rd.initial_condition(uniform_conc=uniform_conc, n=n, L=L)

# Show the shape as a demonstration
c0.shape
```




    (2, 32, 32)



Now, everything is in place. We just have to solve. We do this by calling the `rd.solve()` function. The arguments are obvious from the call below.


```python
# Solve the system
c = rd.solve(c0, t, D=D, beta=beta, gamma=gamma, f=f, f_args=f_args, L=L)

# Look at the shape of the solution
c.shape
```

    /Users/bois/Dropbox/git/rdsolver/rdsolver/rd.py:172: NumbaPerformanceWarning: [1m[1mnp.dot() is faster on contiguous arrays, called on (array(float64, 2d, C), array(float64, 1d, A))[0m[0m
      out[:,i,j] = np.dot(gamma, c[:,i,j])
      0%|          | 0/100 [00:00<?, ?it/s]/Users/bois/Dropbox/git/rdsolver/rdsolver/rd.py:497: NumbaPerformanceWarning: [1m[1mnp.dot() is faster on contiguous arrays, called on (array(complex128, 2d, C), array(complex128, 1d, A))[0m[0m
      + np.dot(A_rhs, c_hat[:,i,j])
    /Users/bois/Dropbox/git/rdsolver/rdsolver/rd.py:497: NumbaPerformanceWarning: [1m[1mnp.dot() is faster on contiguous arrays, called on (array(complex128, 2d, C), array(complex128, 1d, A))[0m[0m
      + np.dot(A_rhs, c_hat[:,i,j])
    100%|██████████| 100/100 [00:14<00:00,  7.00it/s]





    (2, 32, 32, 100)



We note that the solution is of the shape $(n_s, n_x, n_y, n_t)$, where $n_t$ is the number of time points we used. This structure is useful to know for slicing out species and time points of interest.

For plotting purposes, it is useful to interpolate the solution to have smooth concentration profiles. This is purely aesthetic; the solver will give a pixelated, but spectrally accurate, solution. We can use the `rd.viz.interpolate_concs()` function to get the interpolated concentration profiles.


```python
c_interp = rd.viz.interpolate_concs(c)
```

Finally, we are ready to display the solution. We use the `rd.viz.display_notebook()` function that will give a picture of the concentration field with a slider for adjusting the time. By default, for multiple species problems, (like the ASDM), up to three species are shown. The cyan channel is the first, magenta the second, and yellow is the third (though it is not present for the ASDM).


```python
# Use whatever your notebook URL is here
notebook_url = 'localhost:8889'

bokeh.io.show(rd.viz.display_notebook(t, c_interp), notebook_url=notebook_url)
```






If we like, we can look at a single species, in which case the colormap is Viridis.


```python
bokeh.io.show(rd.viz.display_notebook(t, c_interp[0]), notebook_url=notebook_url)
```






## The full ASDM model

We just solved a more simplified ASDM model, but we may consider a more complicated one, as proposed by Koch and Meinhardt (*Rev. Mod. Phys.*, 1994), in which the autocatalysis reactions can saturate.

\begin{align}
\partial_t a &= D_a (\partial_x^2 + \partial_y^2)a + \frac{\rho_a a^2 s}{1 + \kappa_a a^2} - \mu_a a + \sigma_a \\[0.5em]
\partial_t s &= D_s (\partial_x^2 + \partial_y^2)s - \frac{\rho_s a^2 s}{1 + \kappa_a a^2} + \sigma_s
\end{align}

For convenience, `rdsolver` has a growing set of models that you can pre-load. We can use `rd.models.asdm()` to load in the parameters, as well as the homogeneous steady state, for the more complete ASDM model.


```python
D, beta, gamma, f, f_args, homo_ss = rd.models.asdm()
```

The parameters in the above equations are inputted as keyword arguments, with the defaults set to what was used to generate Fig. 2 in the Koch and Meinhardt paper. The outputted function `f` is JITted for performance.

Note that if we wanted to recapitulate the simpler example we already did, we can call the function with the appropriate kwargs.

```python
d = 0.05
mu = 1.4
params = {'D_a': d,
          'D_s': 1,
          'rho_a': 1,
          'rho_s': mu,
          'sigma_a': 0,
          'sigma_s': mu,
          'mu_a': 1,
          'kappa_a': 0}
D, beta, gamma, f, f_args, homo_ss = rd.models.asdm(**params)
```

We already did that one, so we'll do the saturating model now. We'll now define the grid setup and the time points we want.


```python
n = (32, 32)
L = (50, 50)
t = np.linspace(0, 100000, 100)
```

Next, the initial condition, which will be a small perturbation about the steady state.


```python
np.random.seed(42)
c0 = rd.initial_condition(uniform_conc=homo_ss, n=n, L=L)
```

Now we can solve and interpolate....


```python
c = rd.solve(c0, t, D=D, beta=beta, gamma=gamma, f=f, f_args=f_args, L=L)
c_interp = rd.viz.interpolate_concs(c)
```

    100%|██████████| 100/100 [00:09<00:00, 10.46it/s]


And finally, we visualize:


```python
bokeh.io.show(rd.viz.display_notebook(t, c_interp), notebook_url=notebook_url)
```






Finally, if we want to display a single image, we can do so using `rd.viz.display_single_frame()`.


```python
bokeh.io.show(rd.viz.display_single_frame(c_interp, i=-1));
```








<div class="bk-root" id="bba59902-4a07-4c12-9619-3e4c9c7ec8de" data-root-id="4619"></div>





5c3b6e9dbc07863d634681f16e5c6a8f4928926b
ipynb
Jupyter Notebook
notebooks/asdm_example.ipynb
emorisse/rdsolver
89ef35eeadc50bf3618e10fd7e3f1ed0250ead30
MIT
2021-04-27T03:47:17.000Z
2022-01-17T19:30:06.000Z
notebooks/asdm_example.ipynb
emorisse/rdsolver
89ef35eeadc50bf3618e10fd7e3f1ed0250ead30
MIT
2017-07-14T22:52:20.000Z
2017-08-31T22:55:32.000Z
notebooks/asdm_example.ipynb
emorisse/rdsolver
89ef35eeadc50bf3618e10fd7e3f1ed0250ead30
MIT
2021-08-16T14:59:00.000Z
2021-10-14T04:55:48.000Z
Qwen/Qwen-72B
1. YES
2. YES