# Advection-Diffusion#

In this example, we will learn how to perform an advection-diffusion simulation of a given chemical species through a `Cubic`

network. The algorithm can be applied to more complex networks in the same manner as described in this example. For the sake of simplicity, a one layer 3D cubic network is used here. On `OpenPNM`

, 4 different space discretization schemes for the advection-diffusion problem are available and consist of:

Upwind

Hybrid

Powerlaw

Exponential

Depending on the Peclet number characterizing the transport (ratio of advective to diffusive fluxes), the solutions obtained using these schemes may differ. In order to achive a high numerical accuracy, the user should use either the `powerlaw`

or the `exponential`

schemes.

## Generating network#

First, we need to generate a `Cubic`

network. For now, we stick to a one layer 3d network, but you might as well try more complex networks!

```
[1]:
```

```
import numpy as np
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
np.random.seed(10)
%matplotlib inline
ws = op.Workspace()
ws.settings["loglevel"] = 40
np.set_printoptions(precision=5)
shape = [1, 20, 30]
net = op.network.Cubic(shape=shape, spacing=1e-4)
```

## Adding geometry#

Next, we need to add a geometry to the generated network. A geometry contains information about size of the pores/throats in a network. `OpenPNM`

has tons of prebuilt geometries that represent the microstructure of different materials such as Toray090 carbon papers, sand stone, electrospun fibers, etc. For now, we stick to a sample geometry called `SpheresAndCylinders`

that assigns random values to pore/throat diameters.

```
[2]:
```

```
geom = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)
```

## Adding phase#

Next, we need to add a phase to our simulation. A phase object(s) contain(s) thermophysical information about the working fluid(s) in the simulation. `OpenPNM`

has tons of prebuilt phases as well! For this simulation, we use air as our working fluid.

```
[3]:
```

```
air = op.phases.Air(network=net)
```

```
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Input In [3], in <cell line: 1>()
----> 1 air = op.phases.Air(network=net)
AttributeError: module 'openpnm' has no attribute 'phases'
```

## Adding physics#

Finally, we need to add a physics. A physics object contains information about the working fluid in the simulation that depend on the geometry of the network. A good example is diffusive conductance, which not only depends on the thermophysical properties of the working fluid, but also depends on the geometry of pores/throats.

```
[4]:
```

```
phys_air = op.physics.Standard(network=net, phase=air, geometry=geom)
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 phys_air = op.physics.Standard(network=net, phase=air, geometry=geom)
NameError: name 'air' is not defined
```

# Performing Stokes flow#

Note that the advection diffusion algorithm assumes that velocity field is given. Naturally, we solve Stokes flow inside a pore network model to obtain the pressure field, and eventually the velocity field. Therefore, we need to run the `StokesFlow`

algorithm prior to running our advection diffusion. There’s a separate tutorial on how to run `StokesFlow`

in `OpenPNM`

, but here’s a simple code snippet that does the job for us.

```
[5]:
```

```
sf = op.algorithms.StokesFlow(network=net, phase=air)
sf.set_value_BC(pores=net.pores('front'), values=200.0)
sf.set_value_BC(pores=net.pores('back'), values=0.0)
sf.run();
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 sf = op.algorithms.StokesFlow(network=net, phase=air)
2 sf.set_value_BC(pores=net.pores('front'), values=200.0)
3 sf.set_value_BC(pores=net.pores('back'), values=0.0)
NameError: name 'air' is not defined
```

It is essential that you attach the results from `StokesFlow`

(i.e. pressure field) to the corresponding phase, since the results from any algorithm in `OpenPNM`

are by default only attached to the algorithm object (in this case to `sf`

). Here’s how you can update your phase:

```
[6]:
```

```
air.update(sf.results())
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 air.update(sf.results())
NameError: name 'air' is not defined
```

## Performing advection-diffusion#

Now that everything is set up, it’s time to perform our advection-diffusion simulation. For this purpose, we need to add corresponding algorithm to our simulation. As mentioned above, `OpenPNM`

supports 4 different discretizations that may be used with the `AdvectionDiffusion`

and `Dispersion`

algorithms. Setting the discretization scheme can be performed when defining the physics model as follows:

```
[7]:
```

```
mod = op.models.physics.ad_dif_conductance.ad_dif
phys_air.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw')
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [7], in <cell line: 2>()
1 mod = op.models.physics.ad_dif_conductance.ad_dif
----> 2 phys_air.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw')
NameError: name 'phys_air' is not defined
```

Then, the advection-diffusion algorithm is defined by:

```
[8]:
```

```
ad = op.algorithms.AdvectionDiffusion(network=net, phase=air)
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 ad = op.algorithms.AdvectionDiffusion(network=net, phase=air)
NameError: name 'air' is not defined
```

Note that `network`

and `phase`

are required parameters for pretty much every algorithm we add, since we need to specify on which network and for which phase do we want to run the algorithm.

Note that you can also specify the discretization scheme by modifying the `settings`

of our `AdvectionDiffusion`

algorithm. You can choose between `upwind`

, `hybrid`

, `powerlaw`

, and `exponential`

. It is important to note that the scheme specified within the algorithm’s settings is only used when calling the `rate`

method for post processing.

## Adding boundary conditions#

Next, we need to add some boundary conditions to the simulation. By default, `OpenPNM`

assumes zero flux for the boundary pores.

```
[9]:
```

```
inlet = net.pores('front')
outlet = net.pores(['back', 'top', 'bottom'])
ad.set_value_BC(pores=inlet, values=100.0)
ad.set_value_BC(pores=outlet, values=0.0)
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [9], in <cell line: 3>()
1 inlet = net.pores('front')
2 outlet = net.pores(['back', 'top', 'bottom'])
----> 3 ad.set_value_BC(pores=inlet, values=100.0)
4 ad.set_value_BC(pores=outlet, values=0.0)
NameError: name 'ad' is not defined
```

`set_value_BC`

applies the so-called “Dirichlet” boundary condition to the specified pores. Note that unless you want to apply a single value to all of the specified pores (like we just did), you must pass a list (or `ndarray`

) as the `values`

parameter.

## Running the algorithm#

Now, it’s time to run the algorithm. This is done by calling the `run`

method attached to the algorithm object.

```
[10]:
```

```
ad.run();
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 ad.run()
NameError: name 'ad' is not defined
```

# Post processing#

When an algorithm is successfully run, the results are attached to the same object. To access the results, you need to know the quantity for which the algorithm was solving. For instance, `AdvectionDiffusion`

solves for the quantity `pore.concentration`

, which is somewhat intuitive. However, if you ever forget it, or wanted to manually check the quantity, you can take a look at the algorithm `settings`

:

```
[11]:
```

```
print(ad.settings)
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 print(ad.settings)
NameError: name 'ad' is not defined
```

Now that we know the quantity for which `AdvectionDiffusion`

was solved, let’s take a look at the results:

```
[12]:
```

```
c = ad['pore.concentration']
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 c = ad['pore.concentration']
NameError: name 'ad' is not defined
```

## Heatmap#

Since the network is 2d, we can simply reshape the results in form of a 2d array similar to the shape of the network and plot the heatmap of it using `matplotlib`

.

```
[13]:
```

```
print('Network shape:', shape)
c2d = c.reshape(shape)
```

```
Network shape: [1, 20, 30]
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [13], in <cell line: 2>()
1 print('Network shape:', shape)
----> 2 c2d = c.reshape(shape)
NameError: name 'c' is not defined
```

```
[14]:
```

```
import matplotlib.pyplot as plt
plt.imshow(c2d[0,:,:]);
plt.title('Concentration (mol/m$^3$)');
plt.colorbar();
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [14], in <cell line: 2>()
1 import matplotlib.pyplot as plt
----> 2 plt.imshow(c2d[0,:,:]);
3 plt.title('Concentration (mol/m$^3$)');
4 plt.colorbar()
NameError: name 'c2d' is not defined
```