Tutorial Part 1

boxsimu is a simple modelling framework based on the principle of a system of discrete boxes interacting with each other. boxsimu offers a simple interface to define a system, its components and properties. The system definition is done by instantiating several classes of the boxsimu package like Fluid, Variable, Flow, Flux, Process, Reaction, and Box. The temporal evolution of a system defined in boxsimu can then easily be simulated and the result of the simulation can be visualized or further investigated.

After a short primer on what mathematical modelling is and how boxsimu fits in there, this tutorial will lead you through the installation of boxsimu. After the installion, it will briefly explain the basic structure of a system-definition in boxsimu. After that, you should have a basic understanding of the most important classes of boxsimu and you're ready to dive into a first example!

Introduction

Mathematical modelling

In order to understand the dynamics of a system, it is often useful to create a (mathematical) model thereof. A system can be almost anything, e.g. a mountain lake, the Atlantic ocean, the atmosphere, a machine, a company, or even the entire Earth System. A (mathematical) model represents such a system while neglecting a big amount of complexity. Thus a (mathematical) model tries to explain the fundamental system behaviour as simply as possible.

Schwarzenbach et al. (2003) uses the following defintion of a system:

''A model is an imitation of reality that stresses those aspects that are assumed to be important and omits all properties considered to be nonessential''.

Mathematical modelling of a system is a challenging and complex activity. Depending on the system an analytical solution can even be nonexistent. As a consequence one is forced to simulate/solve such systems numerically. This is where boxsimu comes into play:

boxsimu allows the user to easily define systems of simple to intermediate complexity and to simulate their temporal evolution. All that can be done on a high level of abstraction, meaning the user doesn't have to bother about the representation of the system with computer code.

Limitations of boxsimu

At the current status of development of boxsimu, the numerical solver is not very efficient. As a consequence large and complexe systems with a lot of variables, boxes, and processes/reactions will be solved slowly.

Installation

Dependencies

Python: boxsimu was only tested with Python 3.6. Therefore boxsimu may not work as expected with older versions of Python.

The following python packages are required by boxsimu:

Scientifc packages:

  • numpy (1.10 or newer): Fundamental package for scientifc computing, offers powerful and very fast matrix operations.
  • matplotlib (1.5 or newer): Plotting library with publication quality figures.
  • jupyter (4.0 or newer): Offers a browser-based IDE.
  • pandas (0.15 or newer): Allows spreadsheet-like calculation.
  • pint (0.8.1 or newer): Physical units-framework for python.

Other packages:

  • attrdict (1.2 or newer): Extension of the python dict class which allows to access dict items also as instance attributes.
  • svgwrite (1.1 or newer): Generation of SVG graphics.
  • dill (2.7 or newer): Extended pickeling of python objects.

All of the above listed python package are available in the official python repository and thus can be installed using pip. Thus, if a package is missing on your system just use:

pip install <package name>

Additionally an up-to-date C compiler like gcc is needed, since some computational critical parts of boxsimu are written in Cython, which has to be compiled into binary code on your computer.

Installation using pip

boxsimu can easily be installed using pip. On your console type in

pip install boxsimu

this should automatically compile all Cython files and copy all source files into the right directory. If the installation fails check if all of the above mentioned dependencies are met.

Installation from source

Alternatively to the installation via pip, boxsimu can also be installed from source. For this, download the most recent source code from github {{URL}} or clone the repository. Afterwards open a system console and change into the directory where you downloaded boxsimu. Then execute the following commands:

\$ cd boxsimu <br> \$ python setup.py

boxsimu Overview

A system is defined by instantiating the class BoxModelSystem and passing instances of the classes Box, Flow, and Flux to it. The instance of the class BoxModelSystem contains all information on the whole system while its components (Boxes, Flows, and Fluxes) store information on the distribution of Fluids and Variables in the system and how the different compartements (boxes) exchange these quantities. The basic structure of a BoxModelSystem instance is graphically shown below:

In this diagram an instance of BoxModelSystem is shown that contains two boxes: 'Box 1' and 'Box 2'. Both boxes contain the same fluid ('Fluid 1') and two instances of the class Variable ('A' and 'B'). Additionally both boxes can contain an arbitrary number of independet processes and reactions of these variables. Finally, the boxes exchange mass of the variables and the fluid via fluxes and flows (the difference between a flow and a flux will be explained further down).

A BoxModelSystem can contain an arbitrary number of boxes, however, the more boxes there are the slower the system's simulation will progress. Similarly the number of fluxes, flows, and variables, processes, reactions within boxes is only limited in regard of computational power. In contrast, every box can only contain zero or one instance of the class Fluid.

The most important class that a user of BoxSimu is interacting with are:

  • BoxModelSystem: Contains the whole system that a user wants to simulate and investigate. Most often this is the last class that is instantiated since instance of the other classes are needed as arguments for the BoxModelSystem-constructor. Once an instance of BoxModelSystem is built, its temporal evolution can be simulated using the instance's 'solve' method.
  • Box: Represents a compartement of the system. Contains information about the mass of variables and fluid within the represented compartement and which processes and reaction take place. Additionally for every box an 'environmental condition' (instance of the class Condition) can be defined. These conditions (e.g. the pH or the temperature within the box) in turn can then be used by user-defined rate functions (dynamic rates of processes/reactions/flows/fluxes/fluid-densities that can be dependent on the condition of the box).
  • Flow: Represent a transport of a fluid between two boxes. This flow can passively transport variables (e.g. a river that passively transport a lot of different chemical substances) but doesn't have to (e.g. evaporation from a lake into the atmosphere where almost pure water is transported).
  • Flux: Represents a transport of a variable that is not associated with a flow of a fluid. Examples of fluxes are for example: sedimentation of particles in a sea/ocean/atmosphere.
  • Variable: Variables are tracer within the system which time evolution is simulated. Thus variables are quantites within the system that are of direct or indirect interst to the user.
  • Process: Represents sink or source mechanisms of a certain variable.
  • Reaction: Represents transformations between different variables. There is no mass conservation constraint. Thus a ''transformation'' of approximately 1kg of phosphate into 114kg of phytoplankton (approx. Redfield ratio) is a valid reaction.

A system is most easily defined following these steps:

  • Define all instances of the classes Variable and Fluid
  • Define all instances of the classes Reaction and Process
  • Define all instances of the class Box
  • Define all instances of the classes Flow and Flux
  • Define the instance of the class BoxModelSystem

Given this first fundamental understanding of boxsimu classes and their sequence of appearence within a system-definition, we jump directly into a first example!

A simple Lake 1-Box Model:

Model description

Our first system consists of a freshwater lake that only has one inflow and one outflow. We want to simulate how the concentration of phosphate in this lake evolves over time. In order to do that we assume the inflow to have a constant concentration of phosphate (PO4) while the outflow has the same concentration of PO4 as the (well-mixed) lake itself. The volume/mass of lake-water is constant over time.

In the following we have a simple depiction of the lake system and important parameters thereof:

Analytical solution

Before we use boxsimu to simulate this system we can solve the govering equations of this system analytically in order to validate the output of boxsimu.

So lets start with the defintion of all needed variables and their physical units:

Variables (symbols are consistent with the figure above):

  • $V$ = Volume of the lake [m^3] ($V$ is constant)
  • $m$ = Mass of phosphate in the lake [kg]
  • $C$ = $\frac{m}{V}$ = Concentration of phosphate in the lake [kg/m^3]
  • $C_0$ = Concentration of phosphate in the lake at the begining of the simulation (t=0) [s]
  • $C_{in}$ = Concentration of phosphate in the inflow [kg/m^3]
  • $Q$ = Volumetric water flow rate of the Inflow and Outflow [m^3/s]
  • $t$ = Time [d]

  • $k_w$ = $\frac{Q}{V}$ = Specific flow rate [1/d]

Assumptions:

  • $C_{in}$ = const.
  • $Q$ = const.
  • $V$ = const.

Based on the definition of our system given above we can set up the following differential equation:

$\frac{dm}{dt} = \sum sources - \sum sinks = Q \cdot C_{in} - Q \cdot C$.

The rate of change of the phosphate concentration in the lake is given by the sum of all source terms minus the sum of all sink terms. In this case the only source term is the inflow of phosphate from the river upstream, while the only sink term is the outflow of phosphate through the river downstream.

Where $Q \cdot C_{in}$ and $Q \cdot C$ represents the phosphate-gain and phosphate-loss of the lake per time, respectively.

Next, we devide both sides by the volume of the lake ($V$) and use the specific flow rate $k_w$ on the right hand side (r.h.s.):

$\frac{1}{V}\frac{dm}{dt} = k_w \cdot C_{in} - k_w \cdot C = k_w (C_{in} - C)$

Now, since the volume of the box is constant, we can incorporate the volume into the time-derivative and end up with:

$\frac{1}{V}\frac{dm}{dt} = \frac{d(m/V)}{dt} = \frac{dC}{dt} = k_w (C_{in} - C)$

The solution of this linear, inhomogene, ordinary differential equation is:

$C(t) = (C_0 - C_{in}) e^{-k_wt} + C_{in}$

where we also used the initial condition $C(t=0) = C_0$.

The solution is plotted below for a system with the following parameters:

  • $V = 10^7m^3$
  • $Q = 10^5\frac{m^3}{d}$
  • $C_0 = 10^{-2}\frac{kg}{m^3}$
  • $C_{in} = 3 \cdot 10^{-1}\frac{kg}{m^3}$

We define a function that calculates and returns the concentration of phosphate at a time $t$. We also vectorize this function using numpy in order to be able to apply it on arrays. The resulting array is then plotted as a function of time:


In [5]:
import matplotlib.pyplot as plt
import numpy as np

@np.vectorize
def C(t):
    V = 1e7
    Q = 1e5
    C0 = 1e-2
    Cin = 3e-1
    kw = Q/V
    return (C0-Cin)*np.exp(-kw*t) + Cin

t = np.linspace(0, 8e2, 1000)
c_phosphate = C(t)
plt.plot(t, c_phosphate)
plt.xlabel('Time [days]')
plt.ylabel('PO4 concentration [kg/m^3]')


Out[5]:
<matplotlib.text.Text at 0x7f175c2a6cc0>

We can see that the system reaches a steady-state after about $400$ days.

Now we want to use boxsimu to simulate this system...

boxsimu

Since boxsimu accepts some quantities only in certain units (dimensionalities) we first have to calculate the system parameters in the right dimensionalities/units:

  • The amount of fluid inside a box has to be given in mass units: $m_{water} = V \cdot \rho = 10^7\,m^3 \cdot 10^3\,\frac{kg}{m^3} = 10^{10}kg$
  • The inital amount of phosphate also has to be given in mass units: $m(t=0) = V \cdot C_0 = 10^7\,m^3 \cdot 10^{-2}\,\frac{kg}{m^3} = 10^5\,kg$
  • The flow rate has to be given in units of mass per time: $J = Q \cdot \rho = 10^5\,\frac{m^3}{d} \cdot 10^3\,\frac{kg}{m^3} = 10^8\,\frac{kg}{d}$

Now we define these parameters as python variables. But before we can do that we have to import boxsimu (if not already happened) and the instance of pint.UnitRegistry that is used by boxsimu (boxsimu.ur).

If you use boxsimu you have to use the UnitRegistry from boxsimu, since different UnitRegistry instances are incompatible. Thus use

from boxsimu import ur

instead of importing the UnitRegistry directly from the pint libary.


In [7]:
import boxsimu
from boxsimu import ur

Now you are able to specify the units of a python variable by simply multiplying by pint units:


In [8]:
length = 3 * ur.meter

Pint units can easily be transformed to other units of the same physical dimensionality:


In [9]:
length.to(ur.km) # transforms length to kilometer


Out[9]:
0.003 kilometer

Now lets define the system parameters calculated above:


In [10]:
m_water = 1e10 * ur.kg
m_0 = 1e5 * ur.kg
flow_rate = 1e8 * ur.kg / ur.day

Next we define the boxsimu system! In the following you see the complete model defintion in boxsimu:


In [11]:
# FLUIDS
freshwater = boxsimu.Fluid('freshwater', rho=1000*ur.kg/ur.meter**3)

# VARIABLES
po4 = boxsimu.Variable('po4')

# PROCESSES
# No processes in this system

# REACTIONS
# No reactions in this system

# BOXES
lake = boxsimu.Box(
    name='lake',
    description='Little Lake',
    fluid=freshwater.q(m_water),
    variables=[po4.q(m_0)],
)

# FLOWS
inflow = boxsimu.Flow(
    name='Inflow', 
    source_box=None,
    target_box=lake,
    rate=flow_rate,
    tracer_transport=True,
    concentrations={po4: 3e-1 * ur.gram / ur.kg}, 
)

outflow = boxsimu.Flow(
    name='Outflow',
    source_box=lake,
    target_box=None,
    rate=flow_rate,
    tracer_transport=True,
)

# FLUXES
# No fluxes in this system

# BOXMODELSYSTEM
system = boxsimu.BoxModelSystem(
    name='lake_system', 
    description='Simple Lake Box Model',
    boxes=[lake,], 
    flows=[inflow, outflow,],
)

Code overview

We will go through this code line by line. But before we dive into the details, have a quick look at the sequence of definitions of boxsimu classes:

In a first step, instances of the classes Fluid and Variable are instantiated. A fluid represents the solvent of a box in which variables can be dissolved (e.g. the fluid of a lake is water, in which a myrade of chemicals can be dissolved). A box can, but doesn't have to be filled with a fluid. This allows to also simulate quantites that are not typically dissolved (e.g. to simulate a population of animals). Variables are quantities of the system that are of interest to the user and that are simulated. A variable defined for one box, flow, flux, process, or reaction will automatically be created for all boxes.

At the time of instantiation (creation of the instance) instances of the classes Fluid and Variable are not quantified, that means the instances contain no information about the mass of fluid or variable they represent. However, before they are passed to an instance of the class Box they need to be quantified. This is done by calling the method ''q(mass)'' for a fluid or variable.

In our simple lake system the fluids and variables are:

  • Fluid: The freshwater that is present in the lake (see line 2 in the code).
  • Variable: Phosphate whose concentration we want to analyse as a function of time (see line 5 in the code).

Next, if needed, the classes Process and Reaction are instantiated (created). However, in this simple lake system there are no processes and no reactions.

The created instances of Fluid, Variable, Process, and Reaction are then assigned to instances of the class Box that represent a compartement of the system (see lines 14-19 in the code).

In a next step we have to define how the boxes in our system interact with each other and with the environement (outside the system). That means how the boxes are exchaning mass of fluids and variables. Therefore we instantiate the classes Flow and Flux (see lines 22-40 in the code). The difference between a flow and a flux is that flows represent a transport of fluid that again can (but doesn't have to) passively transport dissolved variables, whereas fluxes represent a transport of variables that is not associated with any fluid transport.

Finally we create an instance of the class BoxModelSystem that glues together all elements of the system (see lines 43-46 in the following code).

Code Line by Line

Lets go through this code line by line:

Line 2:

freshwater = boxsimu.Fluid('freshwater', rho=1000*ur.kg/ur.meter**3)

An instance of the class Fluid, representing the freshwater of the lake, is created. The class Fluid has the following signature:

Fluid(name, rho, description=None)
  • name: String which also has to be a valid python variable name since it is internally used to identify specific fluid instances. Thus, 'fluid123' or 'fluid_123' would be valid, while 'fluid 123' or '123fluid' would not be valid name arguments.
  • rho: Pint quantity or a user-defined function that returns a pint quantity. You find more information about user-defined functions in the second part of the tutorial.
  • description: String which gives a short description of the fluid.

(For further information on the class Fluid confront the documentation: DOCUMENTATION_URL)

Line 5:

po4 = boxsimu.Variable('po4')

An instance of the class Variable is instantiated (created). This variable represents the substance we are interested in: phosphate. Every system needs to define at least one variable since its the temporal evolution of these variables that is simulated with the boxsimu-solver.

The class Variable has the following signature:

Variable(name, molar_mass=None, mobility=True, description=None)
  • name: String which also has to be a valid python variable name since it is internally used to identify specific fluid instances. Thus, 'variable123' or 'variable_123' would be valid, while 'variable 123' or '123variable' would not be valid name arguments.
  • molar_mass: Pint quantity that specifies the molar mass of the substance that is represented by the variable instance. Must be of dimensionality [M/N] (e.g. mol/kg).
  • mobility: Boolean or function that returns boolean. Specifies if the variable is dissolved in a fluid (if one is present) and therefore it is transported by fluid flow.
  • description: String which gives a short description of the variable.

Line 14-19:

lake = boxsimu.Box(
    name='lake',
    name_long='Little Lake',
    fluid=freshwater.q(m_water),
    variables=[po4.q(m_0)],
)

An instance of the class Box is instantiated (created). A box represents a compartement (region/part of the system which is somehow separated or at least distinguishable from the rest of the system.

The class Box has the following signature:

Box(name, description, fluid=None, condition=None, 
    variables=None, processes=None, reactions=None)
  • name: String which also has to be a valid python variable name since it is internally used to identify specific box instances. Thus, 'box123' or 'box_123' would be valid, while 'box 123' or '123box' would not be valid name arguments.
  • description: String which gives a short description of the box.
  • fluid: Instance of the class Fluid or None. Specifies the fluid that is present within a box. If None is given, the box is assumed to not contain a fluid at all.
  • condition: Instance of the class boxsimu.Condition. Specifies the condition within the box (e.g. the pH, the temperature, and more many more properties of a box that remain constant).
  • variables: List of quantified (!) instances of the class Variable. Specifies the variables that are present within the box. The given variables must be quantified - that means a variable's quantified (alias q) must have been called.
  • processes: List of instances of the class Process. Specifies the processes that are running within the box.
  • reactions: List of instances of the class Reaction. Specifies the reactions that are running within the box.

Line 22-37:

inflow = boxsimu.Flow(
    name='Inflow',
    source_box=None,
    target_box=lake,
    rate=flow_rate,
    tracer_transport=True,
    concentrations={po4: 3e-4},
)

outflow = boxsimu.Flow(
    name='Outflow',
    source_box=lake,
    target_box=None,
    rate=flow_rate,
    tracer_transport=True,
)

Two instances of the class Flow are instantiated. A flow represents a transport of a fluid from one box to another, or an exchange of fluid-mass with the environment of the system (the ''the outside''). The first flow (inflow) represents a river that is flowing into the lake, and has an associated phosphate concentratin of $3\cdot10^{-4}\frac{kg}{kg}$. The second flow represent the outflow of the lake.

The class Flow has the following signature:

Flow(name, source_box, target_box, rate, 
     tracer_transport=True, concentrations={}):
  • name: String which also has to be a valid python variable name since it is internally used to identify specific flow instances. Thus, 'flow123' or 'flow_123' would be valid, while 'flow 123' or '123flow' would not be valid name arguments.
  • source_box: Instance of the class Box or None. Specifies where the flow originates. If source_box is set to None, the flow is assumed to come from outside the system.
  • target_box: Instance of the class Box or None. Specifies where the flow ends. If target_box is set to None, the flow is assumed to go outside the system.
  • rate: Pint Quantity or function that returns Pint Quantity of dimensionality [M/T]. Defines the mass-transport per time of the flow.
  • tracer_transport: Boolean. Specifies whether variables are passively transported by this flow.
  • concentrations: Dict of Variables and associated concentrations (Pint Quantity of dimensionality [M/M]). Specifies the concentration of variables within the flow. This attribute can only be set for flows that originate form outside the system (flows that originate from a box within the system always have concentrations equal to the concentration of the box itself).

Line 43-48:

system = boxsimu.BoxModelSystem(
    name='lake_system', 
    description='Simple Lake Box Model',
    boxes=[lake,], 
    flows=[inflow, outflow,],
)

An instance of the class BoxModelSystem is instantiated. A BoxModelSystem represents the system that will be simulated with boxsimu.

The class BoxModelSystem has the following signature:

BoxModelSystem

Solving/Simulating the system

The system defined in boxsimu can now easily be solved using the solve method of the BoxModelSystem class:


In [13]:
sol = system.solve(800*ur.day, 1*ur.day)


Start solving the BoxModelSystem...
- total integration time: 800 day
- dt (time step): 1 day
- number of time steps: 800
10.0%
20.0%
30.0%
40.0%
50.0%
60.0%
70.0%
80.0%
90.0%
Function "solve(...)" used 5.345s

The simulation output is accessible as a pandas dataframe: Solution.df. In our example we can retrieve the model output as follows:


In [18]:
sol.df[:10]


Out[18]:
Box lake
Quantity mass volume po4
Timestep
0 1.000000e+10 10000000.0 129000.000000
1 1.000000e+10 10000000.0 157710.000000
2 1.000000e+10 10000000.0 186132.900000
3 1.000000e+10 10000000.0 214271.571000
4 1.000000e+10 10000000.0 242128.855290
5 1.000000e+10 10000000.0 269707.566737
6 1.000000e+10 10000000.0 297010.491070
7 1.000000e+10 10000000.0 324040.386159
8 1.000000e+10 10000000.0 350799.982297
9 1.000000e+10 10000000.0 377291.982474

Example: The first 10 timesteps of the po4 mass in the box lake can be accessed by:


In [27]:
sol.df[('lake', 'po4')].head(10)


Out[27]:
Timestep
0    129000.000000
1    157710.000000
2    186132.900000
3    214271.571000
4    242128.855290
5    269707.566737
6    297010.491070
7    324040.386159
8    350799.982297
9    377291.982474
Name: (lake, po4), dtype: float64

The output can now be visualized using several methods of the class Solution. In the following just a few examples are given.


In [29]:
sol.plot_variable_concentration(po4)


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-29-bddec9a88ff5> in <module>()
----> 1 sol.plot_variable_concentration(po4)

~/Documents/MyPrivateRepo/boxsimu/boxsimu/solution.py in plot_variable_concentration(self, variable, boxes, figsize, yaxis_log, volumetric, units)
    192         if not self.time_units:
    193             self.time_units = self.time[0].units
--> 194         if not self.time_magnitude:
    195             self.time_magnitude = [t.magnitude for t in self.time]
    196         if not figsize:

AttributeError: 'Solution' object has no attribute 'time_magnitude'

In [ ]:


In [ ]:

References

  • Schwarzenbach, Rene P., Philip M. Gschwend, and Dieter M. Imboden. "Environmental Organic Chemistry. 2003, Hoboken." 808-811.

Add custom css


In [2]:
from IPython.core.display import HTML
#import urllib2
#HTML( urllib2.urlopen('http://bit.ly/1Bf5Hft').read() )
HTML(open('../costum_jupyter_look.css', 'r').read())


Out[2]:

In [ ]: