# Demo

This notebook is a demonstration of **```ngen-kutral```** framework [1,2]. This work is based on model described in [3] also detailed in [4,5] among others.


```python
import wildfire # Framework
from demo_extras import * # Extra functions for experiments
import ipywidgets as widgets # Interactive 
```

## Mathematical Model

Let $u(\mathbf{x},t)$ be the numerical value of the temperature and $\beta(\mathbf{x}, t)$ the amount of fuel,
both at spatial coordinate $\mathbf{x}=(x,y)$ and at time $t$, where $\mathbf{x}\in\Omega=[x_{min},x_{max}]\times[y_{min},y_{max}]\subseteq\mathbb{R}^2$ 
and $t\in [0,t_{\text{max}}]$. 
Let $\mathbf{v}(\mathbf{x},t) = \mathbf{w}(\mathbf{x}, t)+\nabla T(\mathbf{x})$ the vector field that models the effect of wind and topography, and $f(u, \beta)$ a nonlinear heat source, then the model is defined as,

\begin{equation}
    \begin{array}{rcll}
        u_t &=& \kappa\,\Delta u - \mathbf{v} \cdot \nabla u + f(u, \beta), 
            & \textrm{in} ~ \Omega\,\times\, ]0, t_{\text{max}}], \\
        \beta_t &=& g(u, \beta), & \textrm{in} ~ \Omega \, \times \, ]0, t_{\text{max}}], \\
        u(\mathbf{x}, t) &=& h_1(\mathbf{x},t), & \textrm{on} ~ \Gamma \, \times \, ]0, t_{\text{max}}], \\
        \beta(\mathbf{x}, t) &=& h_2(\mathbf{x},t), & \textrm{on} ~ \Gamma \, \times \, ]0, t_{\text{max}}], \\
        u(\mathbf{x}, 0) &=& u_0(\mathbf{x}), & \textrm{in} ~ \Omega, \hfill \\
        \beta(\mathbf{x}, 0) &=& \beta_0(\mathbf{x}), & \textrm{in} ~ \Omega,
    \end{array}
    \tag{1}
\end{equation}

where $\Gamma$ represents the boundary of $\Omega$ and,

\begin{align}
    f(u, \beta) &= H_{pc}(u)\beta \exp\left(\frac{u}{1 + \varepsilon u}\right) - \alpha u \nonumber, \\ 
    g(u, \beta) &= -H_{pc}(u)\frac{\varepsilon}{q}\beta\exp\left(\frac{u}{1 + \varepsilon u}\right) \nonumber,
\end{align}

with

\begin{equation}
    H_{pc}(u) = 
    \begin{cases}
        1, & \text{if} ~ u \geq u_{pc}, \\
        0, & \text{otherwise}.
    \end{cases}
    \nonumber
\end{equation}

In this case the gradient is defined as 
$\nabla = \left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}\right)$
and the Laplace operator $\Delta = \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$.
This model is presented in a nondimensional form using the Frank-Kamenetskii change of variable 
detailed in [1], so the nondimensional parameters of the model $(1)$ are: diffusion coefficient $\kappa$, the inverse of activation energy of fuel $\varepsilon$,
natural convection $\alpha$, reaction heat $q$ and the phase change threshold $u_{pc}$. 

## Difference with Asension-Ferragut 2002

Diffusion term used by [3] is $\operatorname{div}(K(u)\,\nabla u)=\nabla \cdot (K(u)\,\nabla u)$ instead of $\kappa\,\Delta u$, where $K(u)=\bar{\kappa}\,(1 + \varepsilon\, u)^3 + 1$. The diffusion constant coefficient is assumed in works such as [4] and [5].
Expanding $\operatorname{div}(K(u)\,\nabla u)$ we have,

\begin{equation}
    \begin{split}
        \nabla \cdot (K(u)\,\nabla u) &= \left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}\right)
            \cdot \left(K(u)\,\left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}\right)\right)\\
            &= \left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}\right)
            \cdot \left(K(u)\frac{\partial u}{\partial x}, K(u)\frac{\partial u}{\partial y}\right)\\
            &= \frac{\partial}{\partial x}\left(K(u)\frac{\partial u}{\partial x}\right) + 
               \frac{\partial}{\partial y}\left(K(u)\frac{\partial u}{\partial y}\right) \\
            &=\frac{\partial K(u)}{\partial x}\frac{\partial u}{\partial x} + K(u)\frac{\partial^2 u}{\partial x^2} +
              \frac{\partial K(u)}{\partial y}\frac{\partial u}{\partial y} + K(u)\frac{\partial^2 u}{\partial y^2} \\
            &= \frac{\partial K(u)}{\partial u}\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}
                + \frac{\partial K(u)}{\partial u}\frac{\partial u}{\partial y}\frac{\partial u}{\partial y}
                + K(u)\,\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \\
            &= \frac{\partial K(u)}{\partial u}\left(\frac{\partial u}{\partial x}\right)^2 
                + \frac{\partial K(u)}{\partial u}\left(\frac{\partial u}{\partial y}\right)^2
                + K(u)\,\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \\
            &= \frac{\partial K(u)}{\partial u}\left(\left(\frac{\partial u}{\partial x}\right)^2  
                + \left(\frac{\partial u}{\partial y}\right)^2 \right) + K(u)\,\Delta u
    \end{split}
\end{equation}

and $\dfrac{\partial K(u)}{\partial u} = 3\bar{\kappa}\varepsilon\,(1+\varepsilon\,u)^2$. To use this model, set the physical parameter ```complete = True```.

## Physical and nondimensional relation

### Variables

* Space and time:

\begin{equation}
    x = \frac{\bar{x}}{l_0}, \quad y = \frac{\bar{y}}{l_0}, \quad t=\frac{\bar{t}}{t_0}
\end{equation}

* Model components:

\begin{equation}
    u = \frac{\bar{u}-u_{\infty}}{\varepsilon u_{\infty}}, \quad \beta = \frac{\bar{\beta}}{\beta_0}, 
        \quad \mathbf{v} = \frac{t_0}{l_0}\bar{\mathbf{v}}
\end{equation}

$l_0$ is a characteristic length, $t_0$ is characteristic time, $u_{\infty}$ is a reference temperature and $\beta_0$ is a reference fuel. The bar over variables indicates physical magnitude. According to [3,5] $\beta\in[0,1]$.

### Units

|Unit of | Symbol   | Basic SI|
|:------:|:--------:|:-------------------------------------------:|
|Amount of substance|$\text{mol}$||
|Distance|$\text{m}$||
|Energy|$\text{J}$|$\dfrac{\text{kg}\cdot\text{m}^2}{\text{s}^{2}}$|
|Mass|$\text{kg}$||
|Temperature|$\text{K}$||
|Time|$\text{s}$||
| Power  |$\text{W}$|$\dfrac{\text{kg}\cdot\text{m}^2}{\text{s}^{3}}$|



### Constants
* **Universal gas constant**: $R=8.31446261815324\,[\text{J}\cdot\text{K}^{-1} \cdot\text{mol}^{-1}]$.
* **Stefan-Boltzmann constant**: $\sigma = 5.670374419... \times 10^{-8} \,[\text{W}\cdot\text{m}^{-2}\cdot\text{K}^{-4}]$.

### Definitions
* **Thermal conductivity coefficient** $k$, units $[\text{W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}]$.
    - *Wood:* $0.059\,[\text{W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}]$ for maple or oak to $1.17\,[\text{W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}]$ for sawdust.
    - *Air:* $0.024\,[\text{W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}]$
* **Optical path length for radiation through the substance** $\delta$, units $[\text{m}]$.
* **Density** $\rho$, units $[\text{kg}\cdot\text{m}^{-3}]$.
    - *Wood:* Density varies from $420\,[\text{kg}\cdot\text{m}^{-3}]$ for fir to $640\,[\text{kg}\cdot\text{m}^{-3}]$ for yellow pine. 
    - *Air:* At atmospheric pressure and ambient temperature $\rho=1.1774\,[\text{kg}\cdot\text{m}^{-3}]$.
* **Specific heat** $C$, units $[\text{m}^2\cdot\text{K}^{-1}\cdot\text{s}^{-2}]$. 
    - *Wood:* $2400\,[\text{m}^2\cdot\text{K}^{-1}\cdot\text{s}^{-2}]$ for maple to $2800\,[\text{m}^2\cdot\text{K}^{-1}\cdot\text{s}^{-2}]$ for yellow pine.
    - *Air:* $1005.7\,[\text{m}^2\cdot\text{K}^{-1}\cdot\text{s}^{-2}]$
* **Natural convection coefficient** $h$, units $[\text{W}\cdot\text{m}^{-2}\cdot\text{K}^{-1}]$.
* **Activation Energy** $E_A$, units $[\text{J}\cdot\text{mol}^{-1}]$. 
    - Typical values $E_A\approx 83680\,[\text{J}\cdot\text{mol}^{-1}]$.
* **Heat of combustion** $H$, units $[\text{k}\text{J}\cdot\text{kg}^{-1}]$.
    - For cellulose $H=66525.6\,[\text{k}\text{J}\cdot\text{kg}^{-1}]$.
    - Realistic cases may consider smaller values because there are other components in fuel with lower heats of combustion.
* **Pre-exponential factor** $A$, units $[\text{s}^{-1}]$. 
    - Admissible value $A=10^9\,[\text{s}^{-1}]$
* **Reference temperature** $u_{\infty}$, units $[\text{K}]$. 
    - Typical value $u_{\infty}=300\,[\text{K}]$
* **Phase change temperature**, units $[\text{K}]$.
    - Estimated value $\bar{u}_{pc}=550\,[\text{K}]$.
    
### Parameters
* **Inverse of conductivity** $\kappa=\dfrac{4\sigma\delta u_{\infty}^3}{k}$.
* **Inverse of activation energy** $\varepsilon = \dfrac{R~u_{\infty}}{E_A}$.
* **Reaction Heat** $q=\dfrac{H\beta_0}{C\,u_{\infty}}$.
* **Natural convection** $\alpha=\dfrac{t_0h}{\rho C}$.
* **Phase change temperature** $u_{pc}=\dfrac{\bar{u}_{pc}-u_{\infty}}{\varepsilon u_{\infty}}$

### Characteristics
* **Time** $t_0=\dfrac{\varepsilon}{qA}\exp\left(\dfrac{1}{\varepsilon}\right)\,[\text{s}]$.
* **Length** $l_0=\sqrt{\dfrac{t_0 k}{\rho C}}\,[\text{m}]$.

#### Used in [3]:

##### Physical
* $\rho=10^2\,[\text{kg}\cdot\text{m}^{-3}]$
* $C=10^3\,[\text{m}^2\cdot\text{K}^{-1}\cdot\text{s}^{-2}]$
* $k=1\,[\text{W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}]$
* $H=66525.6\,[\text{k}\text{J}\cdot\text{kg}^{-1}]$
* $u_{\infty}=300\,[\text{K}]$
* $\bar{u}_{pc}=550\,[\text{K}]$
* $E_A\approx 83680\,[\text{J}\cdot\text{mol}^{-1}]$
* $t_0=8987\,[\text{s}]$
* $l_0=0.3\,[\text{m}]$

##### Non-dimensional
* $\bar{\kappa}=0.1$
* $\varepsilon \approx 0.03$
* $q=1$
* $\alpha=10^{-3}$

Using this values, we may estimate a constant $\kappa$ to replace $K(u)$, $K(u)=\bar{\kappa}\,(1 + \varepsilon\, u)^3 + 1=0.1\,(1+0.03\cdot 300)^3+1\approx 101$. 


```python
widget_params_ = widgets.interact_manual(checkParams, k=k_w, delta=delta_w, rho=rho_w, 
                                        C=C_w, h=h_w, E_A=E_A_w, H=H_w, A=A_w, U_inf=U_inf_w, B_0=B_0_w)
```


    interactive(children=(FloatText(value=1.0, description='$k$'), FloatText(value=0.01632918494453172, descriptio…


## Physical parameters

Parameters used in our works.


```python
# Physical parameters for the model
physical_parameters = {
    'jcc2018': {
        'kap': 1e-1, # diffusion coefficient
        'eps': 3e-1, # inverse of activation energy
        'upc': 3, # u phase change
        'q': 1, # reaction heat
        'alp': 1e-3, # natural convection,
        # Domain
        'x_lim': (0, 90), # x-axis domain 
        'y_lim': (0, 90), # y-axis domain
        't_lim': (0, 30), # time domain
    },
    'jcc2019': {
        'kap': 1e-1, # diffusion coefficient
        'eps': 3e-1, # inverse of activation energy
        'upc': 3, # u phase change
        'q': 1, # reaction heat
        'alp': 1e-3, # natural convection,
        # Domain
        'x_lim': (0, 90), # x-axis domain 
        'y_lim': (0, 90), # y-axis domain
        't_lim': (0, 20), # time domain
    },
    'af2002': {    
        'kap': 1e-1, # diffusion coefficient
        'eps': 3e-2, # inverse of activation energy
        'upc': 0, # u phase change
        'q': 1.0, # reaction heat
        'alp': 1e-3, # natural convection,
        # Domain
        'x_lim': (0, 300), # x-axis domain 
        'y_lim': (0, 300), # y-axis domain
        't_lim': (0, 0.2), # time domain
        'complete': True
    }
}
```

## Numerical Parameters

### Asension & Ferragut 2002 [3]
* Temperature initial condition: 

\begin{equation}
    u_0(x,y)=3.4\exp\left(-\dfrac{(x-50)^2 + (y-50)^2}{100}\right)
\end{equation}

* Fuel initial condition: 

\begin{equation}
    \beta_0(x,y) \sim \text{Uniform}((0,1)^2)
\end{equation}

* Vector field: 

\begin{equation}
    \mathbf{v}(x, y, t) = (300, 300)
\end{equation}

### San Martin & Torres 2018 [1]
* Temperature initial condition: 

\begin{equation}
    u_0(x,y)=6\exp\left(-\dfrac{(x-20)^2 + (y-20)^2}{20}\right)
\end{equation}

* Fuel initial condition: 

\begin{equation}
    \beta_0(x,y) \sim \text{Uniform}((0,1)^2)
\end{equation}

* Vector field: 

\begin{equation}
    \mathbf{v}(x, y, t) = \left(\cos\left(\frac{\pi}{4}\right), \sin\left(\frac{\pi}{4}\right)\right)
\end{equation}

### San Martín & Torres 2019 [2]
* Temperature initial condition: 

\begin{equation}
    u_0(x,y)=7\exp\left(-\dfrac{(x-20)^2 + (y-20)^2}{20}\right) + 4\exp\left(-\dfrac{(x-80)^2 + (y-70)^2}{20}\right) + 4\exp\left(-\dfrac{(x-20)^2 + (y-35)^2}{70}\right)
\end{equation}

* Fuel initial condition:

\begin{equation}
    \beta_0(x,y) = 1
\end{equation}

* Vector field:

\begin{equation}
    \begin{split}
        \mathbf{w}(x, y, t) &= \left(\cos\left(\dfrac{\pi}{4}+0.025 t\right), 
            \sin\left(\dfrac{\pi}{4} + 0.025 t\right)\right) \\ 
        T(x,y) &= 1.5\,\bigg[3\exp\left(\dfrac{(x-45)^2+ (y-45)^2}{40}\right) 
            + 2\exp\left(\dfrac{(x-30)^2+(y-30)^2}{60}\right) \\
              &\quad + 3\exp\left(\dfrac{(x-70)^2 + (y-70)^2}{60}\right) 
                + 2\exp\left(\dfrac{(x-20)^2 + (y-70)^2}{70}\right)\bigg] \\
        \mathbf{v}(x, y, t) &= \mathbf{w}(x, y, t) 
            + \left(\dfrac{\partial T(x,y)}{\partial x}, \dfrac{\partial T(x,y)}{\partial y}\right)
    \end{split}
\end{equation}


```python
solve_fun = {
    'af2002': {
        'u0': u0_af2002,
        'b0': b0_af2002,
        'V' : V_af2002
    },
    'jcc2018': {
        'u0': u0_jcc2018,
        'b0': b0_af2002,
        'V' : V_jcc2018
    },
    'jcc2019': {
        'u0': u0_jcc2019,
        'b0': b0_jcc2019,
        'V' : V_jcc2019
    },
}
```


```python
def simulation(experiment, space_method, Nx, Ny, acc, time_method, Nt):
    wildfire_ = wildfire.Fire(**physical_parameters[experiment])
    u0 = solve_fun[experiment]['u0']
    b0 = solve_fun[experiment]['b0']
    V  = solve_fun[experiment]['V']
    t, X, Y, U, B = wildfire_.solvePDE(Nx, Ny, Nt, u0, b0, V, space_method, time_method, last=False, acc=acc)
    time_steps_sl = widgets.IntSlider(value=0, min=0, max=Nt, step=1, continuous_update=False)
    # PHYSICAL VALUES
    t0 = 8987#11067.041552426204
    l0 = 0.3#0.33267163318242515
    U = wildfire_.getTemperature(U)
    B = wildfire_.getFuel(B)
    V = wildfire_.getV(V, l0, t0)
    X, Y = wildfire_.getSpace(X, Y, l0)
    t = wildfire_.getTime(t, t0)
    widgets.interact(showSimulation, t=widgets.fixed(t), X=widgets.fixed(X), Y=widgets.fixed(Y),
        U=widgets.fixed(U), B=widgets.fixed(B), V=widgets.fixed(V), k=time_steps_sl)
```


```python
widget_ = widgets.interact_manual(simulation, experiment=experiment_dp, space_method=space_method_dp, 
    Nx=space_nx_text, Ny=space_ny_text, acc=space_acc_dp, 
    time_method=time_method_dp_2, Nt=time_nt_text)
```


    interactive(children=(Dropdown(description='Experiment', index=1, options=(('Asensio & Ferragut 2002', 'af2002…


# References

* [1] San Martín, D., & Torres, C. E. (2018). Ngen-Kütral: Toward an Open Source Framework for Chilean Wildfire Spreading. In 2018 37th International Conference of the Chilean Computer Science Society (SCCC) (pp. 1–8). https://doi.org/10.1109/SCCC.2018.8705159
* [2] San Martín, D., & Torres, C. E. (2019). Exploring a Spectral Numerical Algorithm for Solving a Wildfire Mathematical Model. In 2019 38th International Conference of the Chilean Computer Science Society (SCCC) (pp. 1–7). https://doi.org/10.1109/SCCC49216.2019.8966412
* [3] Asensio, M. I., & Ferragut, L. (2002). On a wildland fire model with radiation. International Journal for Numerical Methods in Engineering, 54(1), 137-157. https://doi.org/10.1002/nme.420
* [4] Weber, R. O. (1991). Toward a comprehensive wildfire spread model. International Journal of Wildland Fire, 1(4), 245–248. https://doi.org/10.1071/WF9910245
* [5] Eberle, S., Freeden, W., & Matthes, U. (2015). Forest fire spreading. In Freeden Willi, M. Z. and Nashed, & and Sonar Thomas (Eds.), Handbook of Geomathematics: Second Edition (pp. 1349–1385). Berlin, Heidelberg: Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-54551-1_70

3237c2b75bfb5e33c28258c768be35f7a7313288
ipynb
Jupyter Notebook
examples/demo.ipynb
dsanmartin/ngen-kutral
4e724b8108a698dec25aa122b831f42c9b91ba9a
BSD-3-Clause
examples/demo.ipynb
dsanmartin/ngen-kutral
4e724b8108a698dec25aa122b831f42c9b91ba9a
BSD-3-Clause
examples/demo.ipynb
dsanmartin/ngen-kutral
4e724b8108a698dec25aa122b831f42c9b91ba9a
BSD-3-Clause
2019-06-12T18:46:06.000Z
2019-06-12T18:46:06.000Z
Qwen/Qwen-72B
1. YES
2. YES