<jupyter_start><jupyter_text>ODYM Example no. 1.  System with two processes, two parameters, one material.A simple MFA system with one material (represented by the indicator element carbon 'C'), a time horizon of 30 years [1980-2010], two processes, and a time-dependent parameter is analysed.The model equations are as follows:+ $a(t) = D(t)$ (exogenous input flow)+ $d(t) = \alpha (t)\cdot b(t)$ (recovery efficiency parameter)+ $a(t) +d(t) = b(t) $ (mass balance process 1)+ $b(t) = c(t) + d(t) $ (mass balance process 2)From these equations the system solution follows:+ $c(t) = a(t) = D(t)$+ $b(t) = \frac{1}{1-\alpha (t)}\cdot D(t)$+ $c(t) = \frac{\alpha}{1-\alpha (t)}\cdot D(t)$We will now programm this solution into ODYM. That is overkill, as ODYM was developed for handling much more complex MFA systems, but instructive. 1) Load ODYM<jupyter_code># Load a local copy of the current ODYM branch:
import sys, os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pickle

# For Ipython Notebook only
%matplotlib inline

sys.path.insert(0, os.path.join('..', 'odym', 'modules')) # add ODYM module directory to system path, relative
sys.path.insert(0, os.path.join(os.getcwd(),'..', 'odym', 'modules')) # add ODYM module directory to system path, absolute
import ODYM_Classes as msc # import the ODYM class file
import ODYM_Functions as msf # import the ODYM function file
import dynamic_stock_model as dsm # import the dynamic stock model library<jupyter_output><empty_output><jupyter_text>2) Define MFA system With the model imported, we can now set up the system definition. The 'classical' elements of the system definition in MFA include: The processes, flows, and stocks, the material, the region, and the time frame studied. Next to these elements, ODYM features/requires the following elements to be specified:+ The list of chemical elements considered+ The classification(s) of the system variables (stocks and flows): Which materials, products, regions, or waste groups are considered?+ An index letter to quickly/directly access a model aspect.+ A dictionary of model parametersFor all these items ODYM has a specific structure, which is used below.First, we define a classification of all objects flowing. In this example we conduct a 'classical' dynamic MFA with just one material/chemical element/indicator element considered, and our classification can therefore be as simple as possible: it contains only one chemical element, in this case, we choose carbon ('C'):<jupyter_code>ModelClassification  = {} # Create dictionary of model classifications

ModelClassification['Time'] = msc.Classification(Name = 'Time', Dimension = 'Time', ID = 1, 
                                                 Items = list(np.arange(2000,2301)))

# Classification for time labelled 'Time' must always be present, 
# with Items containing a list of odered integers representing years, months, or other discrete time intervals

ModelClassification['Element'] = msc.Classification(Name = 'Elements', Dimension = 'Element', ID = 2, 
                                                    Items = ['C'])

# Classification for elements labelled 'Element' must always be present, 
# with Items containing a list of the symbols of the elements covered.

# Get model time start, end, and duration:
Model_Time_Start = int(min(ModelClassification['Time'].Items))
Model_Time_End = int(max(ModelClassification['Time'].Items))
Model_Duration = Model_Time_End - Model_Time_Start<jupyter_output><empty_output><jupyter_text>That dictionary of classifications enters the index table defined for the system. The indext table lists all aspects needed and assigns a classification and index letter to each aspect.<jupyter_code>IndexTable = pd.DataFrame({'Aspect'        : ['Time','Element'], # 'Time' and 'Element' must be present!
                           'Description'   : ['Model aspect "time"', 'Model aspect "Element"'],
                           'Dimension'     : ['Time','Element'], # 'Time' and 'Element' are also dimensions
                           'Classification': [ModelClassification[Aspect] for Aspect in ['Time','Element']],
                           'IndexLetter'   : ['t','e']}) # Unique one letter (upper or lower case) indices to be used later for calculations.

# Default indexing of IndexTable, other indices are produced on the fly
IndexTable.set_index('Aspect', inplace = True) 

IndexTable<jupyter_output><empty_output><jupyter_text>We can now define our MFA system:<jupyter_code>Dyn_MFA_System = msc.MFAsystem(Name = 'NordicBiochar', 
                      Geogr_Scope = 'Nordic regions', 
                      Unit = 'Mt', 
                      ProcessList = [], 
                      FlowDict = {}, 
                      StockDict = {},
                      ParameterDict = {}, 
                      Time_Start = Model_Time_Start, 
                      Time_End = Model_Time_End, 
                      IndexTable = IndexTable, 
                      Elements = IndexTable.loc['Element'].Classification.Items) # Initialize MFA system<jupyter_output><empty_output><jupyter_text>This system has a name, a geographical scope, a system-wide unit, a time frame, an index table with all aspects defined, and a list of chemical elements considered. 3) Inserting data into the MFA systemIt is lacking a list of processes, stocks, flows, and parameters, and these are now defined and inserted into the system:<jupyter_code>Dyn_MFA_System.ProcessList = [] # Start with empty process list, only process numbers (IDs) and names are needed.
Dyn_MFA_System.ProcessList.append(msc.Process(Name = 'Atmosphere', ID   = 0))
Dyn_MFA_System.ProcessList.append(msc.Process(Name = 'Pyrolysis 1'  , ID   = 1))
Dyn_MFA_System.ProcessList.append(msc.Process(Name = 'Decay in soil'  , ID   = 2))

# Print list of processes:
Dyn_MFA_System.ProcessList
ParameterDict = {}

# Define parameter Inflow (D) with indices 'te' (years x element) and matching time series Values (array with size 31 x 1).
# In a more advanced setup the parameters are defined in a data template and then read into the software.
a = np.zeros(301).reshape(301,1)
for i in range(15,35):
    a[i][0] = 10 
ParameterDict['D']  = msc.Parameter(Name = 'Biomass',
                                    ID = 1, P_Res = 1, MetaData = None, Indices = 'te',
                                    Values= a, Unit = 'Mt/yr')

# Define parameter Recovery rate (alpha) with indices 'te' (years x element) and matching time series Values(array with size 31 x 1).
# In a more advanced setup the parameters are defined in a data template and then read into the software.
b = 0.30
ParameterDict['alpha'] = msc.Parameter(Name = 'Biochar yield',
                                       ID = 2, P_Res = 2, MetaData = None, Indices = '',
                                       Values= b, Unit = '1')

# Assign parameter dictionary to MFA system:
Dyn_MFA_System.ParameterDict = ParameterDict
# Define the four flows a,b,c,d of the system, and initialise their values:
Dyn_MFA_System.FlowDict['a'] = msc.Flow(Name = 'Biomass inflow'     , P_Start = 0, P_End = 1, Indices = 't,e', Values=None)
Dyn_MFA_System.FlowDict['b'] = msc.Flow(Name = 'Biochar flow' , P_Start = 1, P_End = 2, Indices = 't,e', Values=None)
Dyn_MFA_System.FlowDict['c'] = msc.Flow(Name = 'Mineralisation'     , P_Start = 2, P_End = 0, Indices = 't,e', Values=None)
Dyn_MFA_System.FlowDict['d'] = msc.Flow(Name = 'Biogenic flow'     , P_Start = 1, P_End = 0, Indices = 't,e', Values=None)


Dyn_MFA_System.StockDict['S']   = msc.Stock(Name = 'Biochar stock', P_Res = 2, Type = 0,
                                              Indices = 't,e', Values=None)
Dyn_MFA_System.StockDict['dS']  = msc.Stock(Name = 'Net biochar stock addition', P_Res = 2, Type = 1,
                                              Indices = 't,e', Values=None)

#Dyn_MFA_System.StockDict['dS_out']  = msc.Stock(Name = 'Biochar stock removal', P_Res = 2, Type = 2,
#                                              Indices = 't,e', Values=None)


# Assign empty arrays to flows according to dimensions.
Dyn_MFA_System.Initialize_FlowValues()
Dyn_MFA_System.Initialize_StockValues()
# Check whether flow value arrays match their indices, etc. See method documentation.
Dyn_MFA_System.Consistency_Check()<jupyter_output><empty_output><jupyter_text>4) Programming a solution of the MFA systemNow the system definition is complete, and we can program the model solution:<jupyter_code>Dyn_MFA_System.FlowDict['a'].Values = Dyn_MFA_System.ParameterDict['D'].Values

Dyn_MFA_System.FlowDict['b'].Values = Dyn_MFA_System.ParameterDict['alpha'].Values * Dyn_MFA_System.ParameterDict['D'].Values
Dyn_MFA_System.FlowDict['d'].Values = (1-Dyn_MFA_System.ParameterDict['alpha'].Values) * Dyn_MFA_System.ParameterDict['D'].Values

S_calc = np.diag(Dyn_MFA_System.FlowDict['b'].Values.reshape(301))
S_calc.shape
i = 0
C_stable = 0.95
C_labile = 1-C_stable
k1 = 3.3945
k2 = 0.0018

for col in S_calc.T:
    if col[i] != 0:
        for t in range(i+1,301):
            S_calc[i][t] = col[i]*( C_labile*np.exp(-k1*(t-i)) + C_stable*np.exp(-k2*(t-i)) )
    i=i+1
    
e_sum = np.ones(301).reshape(301,1)
e_sum.shape
S = S_calc.T.dot(e_sum)
S.shape

O_C = np.zeros(301).reshape(301,1)
O_C[1:300] = S[1:300]-S[0:299]
O_C = Dyn_MFA_System.FlowDict['b'].Values - O_C
O_C

dS = np.zeros(301).reshape(301,1)
dS[1:300] = S[1:300]-S[0:299]
dS.shape
Dyn_MFA_System.StockDict['S'].Values = S
Dyn_MFA_System.StockDict['dS'].Values = dS
Dyn_MFA_System.FlowDict['c'].Values = O_C # annual outflow
#Dyn_MFA_System.StockDict['dS_out'].Values = dS<jupyter_output><empty_output><jupyter_text>5) Mass-balance-check, analyse, and store the model solutionOne major advantage of the ODYM system structure is that mass balance checks can be performed automatically using unit-tested routines without further programming need:<jupyter_code>Bal = Dyn_MFA_System.MassBalance()
print(Bal.shape) # dimensions of balance are: time step x process x chemical element
print(np.abs(Bal).sum()) # reports the sum of all absolute balancing errors.<jupyter_output>(301, 3, 1)
120.0<jupyter_text>The ODYM mass balance array reports the balance for each chemical element, each year, and each process, including the system balance (process 0).<jupyter_code>fig, ax = plt.subplots()
ax.plot(Dyn_MFA_System.IndexTable['Classification']['Time'].Items, Dyn_MFA_System.FlowDict['a'].Values)
ax.set_ylabel('Flow a in Mt/yr', fontsize =16)
fig, ax = plt.subplots()
ax.plot(Dyn_MFA_System.IndexTable['Classification']['Time'].Items, Dyn_MFA_System.FlowDict['c'].Values)
ax.set_ylabel('Flow c in Mt/yr', fontsize =16)
fig, ax = plt.subplots()
ax.plot(Dyn_MFA_System.IndexTable['Classification']['Time'].Items, Dyn_MFA_System.StockDict['S'].Values)
ax.set_ylabel('Stock S in Mt/yr', fontsize =16)<jupyter_output><empty_output><jupyter_text>Save entire system:<jupyter_code>pickle.dump({'BiocharCseq': Dyn_MFA_System}, open("BiocharCseq.p", "wb") )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1135125317835089, 'fraction_numerical': 0.015346894297130403, 'mean_word_length': 3.8218038528896674, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '217119', 'license': 'MIT', 'n_tokens_mistral': 3619, 'n_tokens_neox': 3371, 'n_words': 1074, 'path': 'docs/.ipynb_checkpoints/ODYM Attempt No. 1. Biochar C stocks-checkpoint.ipynb', 'repo_name': 'ntropy-esa/ODYM_biochar', 'size': '11012'}