# Using the transport model classes

## Class hierarchy

`BaseTransportModel1D`

is the abstract base class for all other classes. It includes functions i.e. to retrieve results, set up the model, etc., as far as they are not solver-specific.

Derived model classes are:

`EulerBackwardsModel1D`

, which uses the implicit Euler (first order) time evolution scheme and`CrankNicolsonTransportModel1D`

, which builds on some routines of the former and implements the Crank-Nicolson scheme.

The classes in `ensemble.py`

are wrappers to run multiple submodels of a given transport problem in parallel.
A submodel consists of the transport model for a complete set of radionuclides, that is, a set of radionuclides in which
no radionuclide decays into a radionuclide that is not included in the submodel

## General remarks

The life cycle of transport model objects consists of instantiation/parametrization, running the model, and analyzing the results. Parametrization always happens at instantiation; i.e. parameters should not be changed after instantiation. Similarly, running the same model object twice will usually not work.

## Instantiation/Parametrization

All model classes share, as far as possible, the same interface. For example, the `EulerBackwardsModel1D`

has the following arguments:

```
def __init__(self,
species,
parameters_rock,
parameters_transport_material=None,
geometry=None,
bcs=None,
time_interval=1e6 * u.yr,
options=None,
parameters_source=None
):
```

`species`

is just a list of the radionuclides considered in the model.
The physical parameters of the model domain are provided via
two arguments, `parameters_rock`

and `parameters_transport_material`

. The former is structured as a list of geological units,
where each list entry is a dictionary. An example is shown below

```
import astropy.units as u
parameters_rock = [{'name': "unit1",
'thickness': 1000 * u.m,
'porosity': 0.1,
'head_gradient': 0.0,
'hydraulic_conductivity': 2e-13 * u.m / u.s,
'bulk_density': 2650 * u.kg / u.m**3}]
```

This defines a setup with only one layer, named unit1.
Note the usage of the astropy.units package. We make use of it in many places. `parameters_transport_material`

is another list of dictionaries,
with each entry defining the transport parameters for one geological unit. An example is again shown below.

```
parameters_transport_material = [
{
'I-129':
{'effective_diffusion_coefficient': 2.4e-12 * u.m**2 / u.s,
'sorption_coefficient': 0.0 * u.m**3 / u.kg},
'Ca-41':
{'effective_diffusion_coefficient': 1e-11 * u.m**2 / u.s,
'sorption_coefficient': 6.3e-4 * u.m**3 / u.kg}
}
]
```

In this case, sorption and diffusion coefficients are given for unit1, for the radionuclides I-129 and Ca-41. Both
`parameters_rock`

and `parameters_transport_material`

need to have the same length.

`geometry`

specifies the grid spacing and/or the grid. Internally, this is treated as an GeometryArgs object, but it can also be provided
as a `dict`

object, e.g.

```
geometry = {'dx':1*u.m}
```

This will generate a grid with a grid spacing of 1 m.

`bcs`

will be handled as a `BCArgs`

object, but can be provided a dictionary as well. If None is provided, Dirichlet boundary conditions at 0 mol/m^3 will be assigned per default.

`time_interval`

defines the total time period for the simulation. `options`

will be converted to an `OptionsArgs`

object,
but can again be provided as a dictionary. It can be used to set various options, like flags to include/not include specific
processes, etc. `parameters_source`

will evaluate to a `SourceArgs`

object and can be provided as a dictionary. It controls
how source terms are included in the calculation.

## Running the model

A parametrized model can be run like this:

```
model.run(ics, t_eval)
```

where we have omitted some optional arguments. `ics`

is a data structure that provides the initial conditions, usually a list or a dictionary.
If it is a list, it has to provide the initial concentration vector for each species. If it is a dict, it has to provide the initial concentration vector
for each species as key:value pair. All concentration vectors must have astropy units. Units can be kg/m^3 or mol/m^3.

`t_eval`

is a vector of the times at which outputs of the simulation should be stored, and needs to have units.

Initial conditions need to specify the total mass of radionuclides per rock volume, including the sorbed fraction of radionuclides.

## Analyzing the results

The calculated concentrations are always the total concentration of radionuclides (including any sorbed species), either in mass or amount of substance, per rock volume or per volume with the functions below. Internally, the concentrations are stored as total concentrations in mol/m^3 per fluid volume.

Access to the model results can be obtained by the functions `get_curves_t`

, `get_curves_x`

, and some other derived functions. The
former allows accessing the spatial distribution of concentrations at a fixed time, whereas the latter returns the temporal
evolution of the system at a fixed location. Some more complex routines are already implemented, e.g. `plot_overview_c_x`

in `ensemble.VSGEnsembleModel1D`

.

For example, the following code retrieves the spatial distribution of I-129 (as `c`

), and the location of grid nodes (as `x`

), at the last time stored:

```
x, c = model.get_curves_x("I-129", itime=-1, quantity="molar density", units=u.mol/u.l, per_rock_volume=True)
```