Example
#First import numpy, since all arrays are based on numpy
import numpy as np
#Import matplotlib, for plotting
%matplotlib inline
import matplotlib.pyplot as plt
#Import the main software library
import simex_lib as simex
All parameters must be adopted in the same units of space/time/mass/volume here we adopt:
#Experiment basename (for filenames and figures)
name = "magd"
The first step is to set the domain, that is, the position of the interfaces and boundaries.
#Domain definition (position of interfaces)
# x0 x1 x2 x3
# |----|----|---|
#x=np.array([initial_position, interface_1, interface_2, final_position])
x=np.array([0, 0.3, 0.305, 0.355]) # cm
Now we set the names of the compartments.
#Names of compartments
#xnames = ["", "", ""]
xnames = ["sample", "headspace", "extract"]
Set initial concentrations, constant in each compartment.
#Initial concentrations in each space
# Assumed constant for each space
# C0 C1 C2
# |----|----|----|
#C=np.array([C0, C1, C2]) # (mug/mL)
C=np.array([2, 0.0, 0.0]) # (mug/mL)
Diffusion coefficients for each compartment, constant at each one.
#Diffusion coefficients for each space
# D0 D1 D2
# |----|----|----|
#D=np.array([D0, D1, D2]) # (cm^2/s)
D=np.array([219e-6, 0.11, 219e-6]) # (cm^2/s)
Set the partition coefficients. Always enforce zero at external boundaries for no flow condition.
#Partition coefficients K for each interface
# Set 0.0 for beggining and endpoints
# K0 K1 K2 K3
# |---|---|----|
#K=np.array([K0, K1, K2, K3]) # 2 boundaries and 2 interface coefficient, non-dimensional
K=np.array([0.0, 13.2, 1/13.2, 0.0])
Now configure runtime options for this simulation.
#Max time definition (how long to run the simulation)
maxtime = 6000 # seconds
#Time step size
dt = 0.1 # seconds - depending on your diffusion coefficients, this may be required to be small
#Plotting time spots (snapshots of the concentrations are saved in these time instants)
iplot_time=np.array([0.0, 30, 60, 120, 300, 600]) # seconds
#iplot_time=np.linspace(0, 50, 21, endpoint=True) #Equally spaced plots, for animations
#Space discretization (number of grid points) - depending on your device, this may need to be incresed
N = 500
We now build the device with these parameters (we have given the same names of the variables as the device function arguments, but we could have called the variables differently).
p = simex.device(
D = D,
K = K,
C = C,
xspace = x,
xnames = xnames,
name = name,
N = N, dt = dt, maxtime = maxtime, iplot_time = iplot_time)
-------------------------------------------------------------- Simex - Simulation of Extraction Processes -------------------------------------------------------------- You defined a device with 3 compartment(s). Mechanism layout/interfaces (x): [0. 0.3 0.305 0.355] Initial concentrations: [2. 0. 0.] Diffusion coefficients: [0.000219 0.11 0.000219] Interface coefficients: [ 0. 13.2 0.07575758 0. ] Output basename: output/magd/magd Compartment 0 setup Name: sample Local Domain: [0. 0.3] Difusion (neigbours): [0. 0.000219 0.11 ] Border/Interfaces Coef: [ 0. 13.2] Compartment 1 setup Name: headspace Local Domain: [0.3 0.305] Difusion (neigbours): [0.000219 0.11 0.000219] Border/Interfaces Coef: [13.2 0.07575758] Compartment 2 setup Name: extract Local Domain: [0.305 0.355] Difusion (neigbours): [0.11 0.000219 0. ] Border/Interfaces Coef: [0.07575758 0. ] Proposed number of control volumes (grid points): 500 Adjusted number of grid points: 500 Number of degrees of freedom: 496 Time-space info (dx, dt, Nt, maxD, dx/maxD): 0.0007099999999999999 0.1 60000 0.11 0.006454545454545454 Equilibrium concentrations: [1.7124324324324325, 0.12972972972972974, 1.7124324324324325] ------------------------------------------------
Now that the device structure is built, we can actually use in a few ways.
Lets loop in time the model to see how the concentrations evolve in time. The main function here is run(), which calculates all time steps of an implicit solver for the model of the device. It will save and return snapshots at all time set in iplot_time.
concentrations = p.run()
It: 0 Time: 0.0 Mass: 0.59928 %Dif Eq: 100.0000 % It: 300 Time: 30.0 Mass: 0.60162 %Dif Eq: 23.1960 % It: 600 Time: 60.0 Mass: 0.60186 %Dif Eq: 13.4087 % It: 1200 Time: 120.0 Mass: 0.60199 %Dif Eq: 5.7680 % It: 3000 Time: 300.0 Mass: 0.60205 %Dif Eq: 1.1937 % It: 6000 Time: 600.0 Mass: 0.60206 %Dif Eq: 0.4138 % Saving concentrations as output/magd/magd_data.csv %Eq Time 50.0 % 9.00 60.0 % 13.30 70.0 % 20.60 80.0 % 34.60 90.0 % 70.00 95.0 % 112.20 99.0 % 334.90 99.5 % 529.90
The snapshots were saved in the csv files and also returned to the list of snapshots given at "concentrations". We can use this output to plot the concentrations, as follows.
fig, axes = plt.subplots(1, 1, constrained_layout=True, figsize=(15,10))
plt.xlabel("x-distance (cm)")
plt.ylabel("concentration ($\mu g/mL$)")
plt.title("membrane-aided gas-diffusion")
for i, c in enumerate(concentrations):
#Plot
t = iplot_time[i]
axes.plot(p.x, c, label=str(t)+"s")
for i, name in enumerate(xnames):
xlabel=0.5*p.xspace[i]+0.5*p.xspace[i+1]
plt.text(xlabel, -0.3, name)
axes.plot(p.x, p.u_equi_ext, 'k--', label="Theoretical \n Equilibrium")
axes.legend()
plt.savefig(p.basename+"_final.png")