Skip to content

Commit

Permalink
examples update
Browse files Browse the repository at this point in the history
  • Loading branch information
@GJHSimmons
GJHSimmons committed Sep 13, 2024
1 parent 8872609 commit fef1a0a
Show file tree
Hide file tree
Showing 12 changed files with 1,200 additions and 7 deletions.
3 changes: 0 additions & 3 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,6 @@ share/python-wheels/
*.egg
MANIFEST

# Examples
examples/

# Development
api.py
functions.py
Expand Down
161 changes: 161 additions & 0 deletions docs/examples/mixing_problems/multiple_tanks.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
In the previous example, we considered a mixing problem with a single tank, which had an external flow in and an external flow out. Now we consider a mixing problem consisting of two or more tanks, each of which can have external flows in and out in the same manner as for the single tank, and which furthermore can have flows to and from other tanks in the system.

It will help to first specify the system data. In this case, we will create a system with \( n > 1 \) tanks, such that every tank has a flow in and flow out, plus a flow to and from every other tank. The following data can be freely altered to reduce the external flows and connections, if desired.

```py title="system data"
# Number of tanks
n = 2

# List of tanks created
tanks_id = [i for i in range(1, n + 1)]

# List of pipes between tanks eg. (1,2) is a pipe from tank 1 to tank 2
link_pipes_id = [
(i, j)
for i in tanks_id
for j in tanks_id
if j != i
]

# List of tanks with an external pipe in
pipes_in_id = [1, 2]

# List of tanks with an external pipe out
pipes_out_id = [1, 2]
```

If you followed through the single tank example, the following lists of pipes in and out are created in exactly the same way as there. We create pipes in for every tank in `pipes_in_id` and pipes out for every tank in `pipes_out_id`.

```py
from psymple.build import VariablePortedObject

pipes_in = [
VariablePortedObject(
name=f"in_{i}",
assignments=[
("Q_in", "rate * conc"),
("V_in", "rate"),
],
)
for i in pipes_in_id
]

pipes_out = [
VariablePortedObject(
name=f"out_{i}",
assignments=[
("Q_out", "-rate * Q_out / V_out"),
("V_out", "-rate")
],
)
for i in pipes_out_id
]
```

Next, we need to define connector pipes. These have four variables:

- `Q_in`, the amount of salt entering the pipe from a tank,
- `Q_out`, the amount of salt flowing into the next tank. Note that this has the same concentration
as for `Q_in`,
- `V_in`, the volume of water entering the pipe from the tank,
- `V_out`, the volume of water exiting the pipe into the next tank.

We create a connector pipe for every pair in `link_pipes_id`.

```py
connectors = [
VariablePortedObject(
name=f"pipe_{i}_{j}",
assignments=[
("Q_in", "- rate * Q_in / V_in"),
("Q_out", "rate * Q_in / V_in"),
("V_in", "-rate"),
("V_out", "rate"),
],
)
for i, j in link_pipes_id
]
```

Now we just need to define our system of tanks. Compared to the single tank example, the children, variable ports and input ports are exactly the same (just one for each tank, if it was specified), except for aditionally the specification of the flow rates `rate_i_j` of the connector pipes as additional input ports.

The directed wires simply connect the rates at the input ports to the correct pipes.

The only thing which becomes more complicated is the aggregation of the variables. Each tank `i`, and therefore each variable `Q_i` (and `V_i`), can have:

- External in-flows if `i in pipes_in_id`,
- External out-flows if `i in pipes_out_id`,
- A flow to tank `j` for every `j` such that `(i,j) in link_pipes_id`,
- A flow from tank `j` for every `j` such that `(j,i) in link_pipes_id`.

The variable_wires data simply aggregates all these variables to `Q_i` (and `V_i`).

```py title="tanks model"
from psymple.build import CompositePortedObject

tanks = CompositePortedObject(
name="tanks",
children=pipes_in + pipes_out + connectors,
variable_ports=[f"Q_{i}" for i in tanks_id]
+ [f"V_{i}" for i in tanks_id],
input_ports=[f"rate_{i}_{j}" for i, j in link_pipes_id]
+ [f"conc_in_{i}" for i in pipes_in_id]
+ [f"rate_in_{i}" for i in pipes_in_id]
+ [f"rate_out_{i}" for i in pipes_out_id],
directed_wires=[
(f"rate_{i}_{j}", f"pipe_{i}_{j}.rate")
for i, j in link_pipes_id
]
+ [(f"conc_in_{i}", f"in_{i}.conc") for i in pipes_in_id]
+ [(f"rate_out_{i}", f"out_{i}.rate") for i in pipes_out_id]
+ [(f"rate_in_{i}", f"in_{i}.rate") for i in pipes_in_id],
variable_wires=[
(
([f"in_{i}.Q_in"] if i in pipes_in_id else [])
+ ([f"out_{i}.Q_out"] if i in pipes_out_id else [])
+ [f"pipe_{j}_{i}.Q_out" for j in tanks_id if j != i]
+ [f"pipe_{i}_{j}.Q_in" for j in tanks_id if j != i],
f"Q_{i}",
)
for i in tanks_id
]
+ [
(
([f"in_{i}.V_in"] if i in pipes_in_id else [])
+ ([f"out_{i}.V_out"] if i in pipes_out_id else [])
+ [f"pipe_{j}_{i}.V_out" for j in tanks_id if j != i]
+ [f"pipe_{i}_{j}.V_in" for j in tanks_id if j != i],
f"V_{i}",
)
for i in tanks_id
],
)
```

And that's it!

```py
S = System(tanks)
S.compile()

S.set_parameters(
{
"rate_1_2": 3,
"rate_2_1": 10,
"rate_in_1": "Piecewise((4, T<50), (0, True))",
"rate_in_2": 7,
"conc_in_1": 1 / 2,
"conc_in_2": 0,
"rate_out_1": 11,
#"rate_out_2": 0
}
)

print(S)

sym = S.create_simulation(solver="cts", initial_values={"Q_1": 20, "Q_2": 80, "V_1": 800, "V_2": 1000})

sym.simulate(100, n_steps=100)

sym.plot_solution()
```
130 changes: 130 additions & 0 deletions docs/examples/mixing_problems/single_tank.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
# Mixing problems

In these examples we will consider how to model _mixing problems_ in `psymple`. A mixing problem seeks to understand the evolution of a quantity of solvent in a solute, for example a chemical being mixed in a tank of water. These problems can be highly varied, for example solution can enter a tank at variable concentration or rate, or leave the tank at variable rate. Alternatively, several tanks can be connected together, with their solutions being pumped into each other.

## Assumptions

We will assume that tanks are _well-mixed_, that is, the concentration of the solute is the same throughout the tank at any moment in time.

## Single tank

First, we will consider a single tank with an initial volume of water $V_0\,\mathrm{l}$ and an initial amount of $M_0\,\mathrm{g}$ of salt dissolved in it. A solution with concentration $c(t)\,\mathrm{g}/\mathrm{l}$ of salt flows into the tank at rate $r_0(t)\,\mathrm{l}/\mathrm{s}$ and the mixed solution flows out of the tank at a rate of $r_1(t)\,\mathrm{l}/\mathrm{s}$.

Let $V(t)$ be the volume of solution in the tank at time $t$. Then $V'(t) = r_0(t) - r_1(t)$. Furthermore, let $M(t)$ be the amount of salt in the solution at time $t$. The rate of change of salt in the solution is given by $M'(t) = r_0(t) c(t) - r_1(t) M(t)/V(t)$.

### Modelling flows in `psymple`

In `psymple`, we could model $V(t)$ and $M(t)$ by writing just writing the two differential equations above. In that case, however, our model would become fixed, and if we considered a new situation with more than one in-flow or out-flow, we would have to build something new. Instead, consider what is happening at each point a pipe meets the tank.

For the in-flow pipe, there are fluxes $V'(t) = r_0(t)$ and $M'(t) = r_0 (t) c(t)$. For the out-flow pipe, we similarly get $V'(t) = -r_1(t)$ and $M'(t) = - r_1(t) M(t)/V(t)$. The situation inside the tank is then nothing but the aggregation of these volume and mass variables.

We model this initially with two `VariablePortedObject` instances.

```py
from psymple.build import VariablePortedObject
```

First, the model for the in pipe is:

``` py title="tank in-flow model"
pipe_in = VariablePortedObject(
name="pipe_in",
assignments=[
("V", "r_0"),
("M", "r_0*c"),
],
)
```

and the model for the out pipe is:

```py title="tank out-flow model"
pipe_out = VariablePortedObject(
name="pipe_out",
assignments=[
("V", "-r_1"),
("M", "-r_1 * M/V"),
]
)
```

### Defining the system model

Using a `CompositePortedObject`, we can then aggregate the variables of these two objects together, and expose the parameters `r_0`, `r_1` and `c`.

```py title="single tank model"
from psymple.build import CompositePortedObject

tank = CompositePortedObject(
name="tank",
children=[pipe_in, pipe_out],
input_ports=["r_0", "r_1", "c"],
variable_ports=["V", "M"],
variable_wires=[
(["pipe_in.V", "pipe_out.V"], "V"),
(["pipe_in.M", "pipe_out.M"], "M"),
],
directed_wires=[
("r_0", "pipe_in.r_0"),
("r_1", "pipe_out.r_1"),
("c", "pipe_in.c")
]
)
```

That's it! Let's define and compile a `System` class for `tank`.

### Simulating the model

```py title="tank system"
from psymple.build import System

system = System(tank) # (1)!
system.compile()
```

1. It is also possible to call `System(tank, compile=True)`. In this case, the command `system.compile()` doesn't need to be called.

Before we can simulate, we need to provide:

- Initial values for the variables mass `"V"` and salt amount `"M"`. These are provided using a dictionary passed to the argument `initial_values` when a simulation is created.
- Values for the flow rates and concentration in. We can either do this using the method `system.set_parameters`, or as we do here, by passing a dictionary to the argument `input_parameters` when creating a simulation, allowing us to experiment with multiple scenarios.

We will construct four simulations. For each, we'll set initial values of $V_0 = 1000$ and $M_0 = 20$. The input parameters for each simulation will be:

1. $r_0 = 4 = r_1$ and $c = 0.5$. In this case, the volume of the tank should stay constant,
and the amount of salt should continually increase towards a limit.
2. $r_0 = 2$, $r_1 = 4$ and $c = 0.5$. In this case, the volume of the tank will decrease.
3. Being more creative, we can set $r_0 = 4sin(t) + 4$ and $r_1 = 4$, $c = 0.5$. The
volume of the tank will fluctuate, but stay centred around $1000$.
4. Instead, $r_0 = 4 = r_1$ and $c = 0.5sin(t) + 0.5$.

```py title="Setting up the simulations"
for name, inputs in zip(
["sim_1", "sim_2", "sim_3", "sim_4"], # (1)!
[
{"r_0": 4, "r_1": 4, "c": 0.5}, # (2)!
{"r_0": 2, "r_1": 4, "c": 0.5},
{"r_0": "4*sin(T) + 4", "r_1": 4, "c": 0.5},
{"r_0": 4, "r_1": 4, "c": "0.5*sin(T) + 0.5"},
]
):
sim = system.create_simulation(
name=name,
solver="discrete",
initial_values={"V": 1000, "M": 20}, # (3)!
input_parameters=inputs,
)
sim.simulate(10, n_steps=1000)
```

1. These are the names for each simulation.
2. For each simulation, the set of input parameters is passed as a dictionary.
3. The initial values for each simulation are defined here. They can also be varied in the same way as the input parameters are varied.

Finally, we can visualise each model run by using the `plot_solution` method. Each simulation can be accessed from the dictionary `system.simulations` using the keys `"sim_1"`, `"sim_2"`, `"sim_3"` or `"sim_4"`.

```py title="plotting solutions"
system.simulations["sim_1"].plot_solution({"M"})
system.simulations["sim_1"].plot_solution({"V"})
```
89 changes: 89 additions & 0 deletions docs/examples/population_dynamics/malthusian_population.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,89 @@
#

This example builds up to , \( \frac{dx}{dt} = rx \), which is known as *Malthusian growth*.

In this example, we will:

- Demonstrate how to build a Malthusian growth model in `psymple`,
- Understand how to simulate the system model,
- Understand how to build a reusable interface to the model.



This equation can be formed as a VariablePortedObject with an assignment `{"variable": "x", "expression": "r*x"}`. This entry is equivalent to the shorthand `assg = ("x", "r*x")`.

The value of parameter 'r' could be specified immediately in the assignment, for example we could instead write the assignment ("x", "0.1*x").

Alternatively, psymple allows you to specify a value for 'r' through an input port. This
has two advantages:

1. The variable ported object becomes a general, reusable object,
2. The default value can be overwritten from elsewhere in a larger system.

We will see an example of (2) shortly. For now, to add a default value of 0.1 for the parameter 'r', we add an input port `port = {"name": "r", "default_value": 0.1}`. As for assignments, this is equivalent to the shorthand `port = ("r", 0.1)`.

Putting these together, we define the VariablePortedObject `pop` as follows.

```py
pop = VariablePortedObject(
name="malthusian_pop",
assignments=[("x", "r*x")],
input_ports=[("r", 0.1)],
)
```


The following commands create, run and plot a simulation of our equation with initial condition
'x = 1'.
"""

system_1 = System(pop)
system_1.compile()

sim_1 = system_1.create_simulation("sim", solver="discrete", initial_values={"x": 1})
sim_1.simulate(10, n_steps=10)
sim_1.plot_solution()

"""
Now let's see an example of overriding a default parameter value. Suppose that we individually know
the birth rate 'b' and death rate 'd' for the population, with values 0.4 and 0.2, respectively.
The combined Malthusian rate is 'r = b - d'. We can perform this calculation in a FunctionalPortedObject
with assignment `assg = {"parameter": "r", "expression": "b-d"}`. As before, this is equivalent to the
shorthand `assg = ("r", "b-d")`.

Again, we need to tell psymple what we mean by 'b' and 'd', which we can do as default input ports on
this object.
"""

rate = FunctionalPortedObject(
name="rate",
assignments=[("r", "b-d")],
input_ports=[("b", 0.4), ("d", 0.2)],
)

"""
Next, we need to tell our variable object 'pop' to read its value of 'r' from the value of 'r' produced
by the function 'rate'. We do this in a CompositePortedObject containing both 'pop' and 'rate' as
children, using a directed wire from 'rate.r' to 'malthusian_pop.r'. This is written as
`wire = {"source": "rate.r", "destination": "malthusian_pop.r"}`. There is, of course, an equivalent
shorthand `wire = ("rate.r", "malthusian_pop.r")`.
"""

pop_system = CompositePortedObject(
name="malthusian_pop_system",
children=[pop, rate],
directed_wires=[("rate.r", "malthusian_pop.r")],
)

"""
Notice that we did not have to redefine the variable object 'pop'. Since we defined it generally,
we can reuse its mechanics over and over again. The following commands simulate and plot this
new system.
"""

system_2 = System(pop_system)
system_2.compile()

sim_2 = system_2.create_simulation("sim", solver="discrete", initial_values={"malthusian_pop.x": 1})
sim_2.simulate(10, n_steps=10)
sim_2.plot_solution()
Loading

0 comments on commit fef1a0a

Please sign in to comment.