mountie

An R package providing an efficient C++ implementation of the Reactive Multi-Particle Collisions (RMPC) algorithm.

Installation

You can install the development version from GitHub with:

devtools::install_github("dbarrows/mountie")

Examples

library(mountie)

Diffusion-only system (MPC)

Setup:

(diffusion <- rdsys_examples("diffusion"))
#> $volume
#> [1] "S1"
#> # dims: 100 x 5 x 5
#> # h: 0.01
#> # states: 
#> # # A tibble: 2,500 x 4
#>       x     y     z    S1
#>   <int> <int> <int> <dbl>
#> 1     1     1     1     0
#> 2     2     1     1     2
#> 3     3     1     1     1
#> 4     4     1     1     0
#> 5     5     1     1     0
#> # … with 2,495 more rows
#> 
#> $D
#> S1 
#>  1 
#> 
#> $T
#> [1] 0.1
#> 
#> $tau
#> [1] 1e-04
#> 
#> $boundaries
#>       x            y          z         
#> lower "reflective" "periodic" "periodic"
#> upper "reflective" "periodic" "periodic"

Solve:

sol <- rmpc(diffusion)
#> Starting RMPC simulation with the following parameters:
#>  - D:  
#>     S1 = 1
#>  - tau: 0.0001
#>  - T:   0.1
#> ....................................................................................................

names(sol)
#> [1] "sys" "t"   "u"
sol$u[[1]]
#> # A tibble: 27,000 x 7
#>    Type          P.x      P.y     P.z     V.x     V.y     V.z
#>    <chr>       <dbl>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
#>  1 S1-bath 0.00310   0.00363  0.00879  -0.943  -39.7    41.2 
#>  2 S1-bath 0.000319  0.00813  0.00839  88.6    -12.4  -120.  
#>  3 S1-bath 0.00952   0.000979 0.00554 -43.8     24.4    -4.04
#>  4 S1-bath 0.00429   0.000438 0.00252 -16.0     92.1   213.  
#>  5 S1-bath 0.00296   0.00184  0.00425 -71.9     -2.96 -185.  
#>  6 S1-bath 0.00406   0.00369  0.00242  20.8   -209.    -59.9 
#>  7 S1-bath 0.00552   0.00422  0.00692 -91.5   -130.   -109.  
#>  8 S1-bath 0.000410  0.00678  0.00830  58.9     -6.90    2.60
#>  9 S1-bath 0.0000600 0.00699  0.00921 106.      64.0  -186.  
#> 10 S1-bath 0.00696   0.00727  0.00487 147.      36.2  -105.  
#> # … with 26,990 more rows

Solution plot:

library(ggplot2)
library(wplot)
library(patchwork)
theme_set(theme_wc())

px_plots <- lapply(1:length(sol$u), function(i) {
    position_plot(sol, "S1", frame = i) +
        ylim(0, 2.5)
})
plot_indicies <- round(seq(1, length(px_plots), length.out = 3))
wrap_plots(px_plots[plot_indicies], ncol = 1)

Schnakenberg system

Setup:

(schnakenberg <- rdsys_examples("schnakenberg"))
#> $network
#> # Reaction network: 4 reactions x 2 species
#>     Reactants    Products  Rate
#> 1           0 -> U            2
#> 2           U -> 0            6
#> 3           0 -> V            8
#> 4      2U + V -> 3U           6
#> 
#> $volume
#> [1] "U" "V"
#> # dims: 200 x 1 x 1
#> # h: 0.005
#> # states: 
#> # # A tibble: 200 x 5
#>       x     y     z     U     V
#>   <int> <int> <int> <dbl> <dbl>
#> 1     1     1     1     0     0
#> 2     2     1     1     0     0
#> 3     3     1     1     0     0
#> 4     4     1     1     0     0
#> 5     5     1     1     0     0
#> # … with 195 more rows
#> 
#> $D
#>   U   V 
#> 1.0 0.1 
#> 
#> $T
#> [1] 0.01
#> 
#> $tau
#> [1] 1e-04
#> 
#> $boundaries
#>       x            y            z           
#> lower "reflective" "reflective" "reflective"
#> upper "reflective" "reflective" "reflective"

Solve:

sol <- rmpc(schnakenberg)
#> Starting RMPC simulation with the following parameters:
#>  - D:  
#>     U = 1
#>     V = 0.1
#>  - tau: 0.0001
#>  - T:   0.01
#> ....................................................................................................

Solution plot:

library(ggplot2)
library(wplot)
library(patchwork)
theme_set(theme_wc())

px_plots <- function(species) {
    plot_indicies <- round(seq(3, length(sol$u)))
    plots <- lapply(plot_indicies, function(i) {
        position_plot(sol, species, frame = i) +
            ylab(species) +
            ylim(0, 600)
    })
}
wrap_plots(c(px_plots("U"), px_plots("V")), ncol = 2, byrow = FALSE)