Introduction
CADET is a fast and accurate solver to model and simulate process engineering applications. One of these applications is chromatography in which the separation of different components in a mixture takes place. The separation of components of a mixture is driven by adsorption behaviour in a chromatography column. There are almost 20+ adsorption models implemented in CADET to model this behaviour, but there are other few models that are not implemented. The purpose of this post is to aid the users of CADET to implement their own or some other adsorption model in CADET.
The implementation of a new adsorption model occurs in two steps:
- Registering the model in CADET binding model factory
- Implementation of adsorption model equations, and
- Testing
Example-Case
In this guide, the native form of the multi-component Langmuir adsorption model is used as an example. Mathematically this model is given as:
\frac{dq_i}{dt}=k_{a,i} c_{p,i} q_{max,i} \left(1- \sum_{j=0}^{Ncomp-1}q_j/q_{max,j} \right)-k_{d,i} q_i ~~ ~~ i=0,ā¦,N_{comp}-1 ~~ ~~ ~~ (1)
where, k_a,k_d and q_{max} are the adsorption rate, desorption rate and maximum binding capacity respectively. These parameters are generally known as mechanistic parameters and they are provided by the user. The above-mentioned equation is the analytical form of the Langmuir model, but numerically this equation this solved by ensuring the following condition:
\frac{dq_i}{dt} + k_{d,i} q_i -k_{a,i} c_{p,i} q_{max,i} \left(1- \sum_{j=0}^{Ncomp-1}q_j/q_{max,j} \right) = 0 ~~ ~~ ~~ (2)
The combined second and third term of the above equation is termed as the residual, and this is what will be implemented in the binding model module. This can be written as:
F_i=res_i= k_{d,i} q_i- k_{a,i} c_{p,i} q_{max,i} \left(1- \sum_{j=0}^{N_{comp}-1} \frac{q_j}{q_{max,j}} \right) ~~ ~~ ~~ (3)
For solving a set of nonlinear differential-algebraic equations (DAE), solutions of linear systems require the Jacobian (differentiation of output with input) of the system and one of which comes from the binding model modules. In CADET, the jacobian of the binding model can either be calculated via automatic differentiation (AD) or via solving an analytical Jacobian. To implement the analytical Jacobian of the Langmuir model, the residual term should be differentiated w.r.t all the variables i.e., q_i,q_j,c_{p,i} and c_{p,j} and they are given as:
\frac{\partial F_i}{\partial c_{p,i}}= -k_{a,i}*q_{max,i}*\left(1- \sum_{j=0}^{N_{comp}-1}\frac{q_j}{q_{max,j}} \right) ~~ ~~ ~~ (4)
\frac{\partial F_i}{\partial c_{p,j}}=0.0 ~~ ~~ ~~ (5)
\frac{\partial F_i}{\partial q_{i}}=k_{d,i} ~~ ~~ ~~ (6)
\frac{\partial F_i}{\partial q_{j}}=k_{a,i}*c_{p,i}*\frac{q_{max,i}}{q_{max,j}} ~~ ~~ ~~ (7)
Equations 3,4, 6, and 7 are everything that is needed to implement the binding model in CADET.
1. Registring a new binding model
To register a new binding model we first need to create a new binding model file in the CADET binding model source directory. The working directory for this purpose is:
\src\libcadet\model\binding\
Under this binding directory create a new .cpp file and name it as follows:
YourModelNameBinding.cpp
For example, if the name of the new binding model is Example then the name of the new binding model file will be ExampleBinding.cpp. In case of Langmuir isotherm, the name of the binding model will be LangmuirBinding.cpp.
Even though itās an empty file (contents of this file will be added in the later section) but letās first register this binding model in the Binding Model Factory of CADET in your working brach or forked repository. Open the following file:
\src\libcadet\BindingModelFactory.cpp
Add the following lines after the last statement in the scope of namespace binding {ā¦}
void registerExampleModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings);
Then go to the last statement in BindingModelFactory::BindingModelFactory() definition and add following lines:
model::binding::registerExampleModel(_bindingModels);
The next step in registering the binding model is to make changes in CMakeLists.txt file which is located at the following location:
Open this file and navigate to Line 79 or line with the following statement:
set(LIBCADET_BINDINGMODEL_SOURCES).
This section holds all the source files of the binding models registered in CADET. Within this scope go to the last line and add the following line of code:
${CMAKE_SOURCE_DIR} /src/libcadet/model/binding/ExampleBinding.cpp
After this point, the binding model is registered in the CADET binding model factory.
2. Implementation
To aid the users in implementation, a sample binding model file has been created. Users are advised to use this file as a starting point.
Implementation of any adsorption model occurs in two parts:
- Binding model parameters configuration, where the relavant mechanistic parameters reading interface is defined.
- Implementation of adsorption flux and Jacobian
In a multi-component Langmuir isotherm, three parameters need to be configured and they are, adsorption rate constant k_a , desorption rate constant k_d and maximum binding capacity q_{max}. To set up the configuration of isotherm parameters a macro (.json script) has been included in the code, which generates the relevant code section when the user defines the parameters in the scope of this script. To modify the script go to Line 30 and you will see the following piece of code:
/*<codegen>
{
"name": "LangmuirParamHandler",
"externalName": "ExtLangmuirParamHandler",
"parameters":
[
{ "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "MCL_KA"},
{ "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "MCL_KD"},
{ "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCL_QMAX"}
]
}
</codegen>*/
The above script has three scopes namely, ānameā (defining the name of binding model parameter handler class, in this case, we are implementing Langmuir isotherm so as per CADET naming convention the name is defined as "LangmiurParamHandler" ), āexternalNameā (defining the name of the class for external function support, as per the CADET naming convention it is defined as āExtLangmuirParamHandlerā) and āparametersā. To understand the scope field under parameters consider the following line of code
{ "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "MCL_KA"},
If you have ever tried to define a variable in C++ or C then the above line will be easy to understand. For example, to define a variable in the code the user has to write the type of the variable followed by the name of the variable like double varName; Similarly, in the āparametersā scope, we first defined the " type" of the parameter, which followed the name of the variable (āvarNameā) and in addition, a configuration name of the parameter is given in the scope āconfNameā. The string defined in āconfNameā is the field that the user will define in the CADET input file to assign the value to the corresponding variable. As per CADET naming convention, the āconfNameā should follow following semantics:
Abbreviation of isotherm name_varName
In case of multi-component Langmuir, the abbreviation is MCL so the confName will take the form:
MCL_KA
Note: āconfNameā scope entry should be defined in CAPITAL LETTERS. Similarly, all other parameters can be defined by adding more fields in āparametersā scope.
To configure the name of the binding model in CADET, go to Line 56, which is a class method definition, and it is given as follows:
inline const char* LangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR"; }
The user needs to make two changes in the above line of code when implementing their binding model. The first change is to modify the name of the class to the one that has been defined in ānameā section of .json script. For example, in this example we have used the class name as
āLangmuirParamHandlerā.
The second change is the configuration name of the binding model which should be defined as return value of this function. The return value should be given as string in CAPITAL LETTERS and in CADET input file the same string should be used to call this binding model.
Now we move on to configuration validation methodologies that will check whether the user input a correct number of entries in the parameter list. For example, for two-component analysis, the user should input two values each for MCL_KA, MCL_KD and MCL_QMAX, and in case if the user input fewer or more entries in these fields then CADET should throw an exception and inform the user about this error. To see how this is coded and how you can make changes related to your binding model consider following the LangmuirParamHandler class method defined in between Lines 58-64.
inline bool LangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
{
if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp))
throw InvalidParameterException("MCL_KA, MCL_KD, and MCL_QMAX have to have the same size");
return true;
}
The if condition in this function checks for all possible conditions that checks whether the defined parameter size are compatible with each other or not. One should notice here, that the name of the parameters used are nothing but the ones that has been defined in the .json script field āvarNameā i.e., _kA, _kD and _qMax (_ should be added in front of these variables as these varibles are defined like that in LangmuirParamHandler class). For external function support of this binding model a similar methodology is adopted, as can be seen in Lines 66-74.
inline const char* ExtLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR"; }
inline bool ExtLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
{
if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp))
throw InvalidParameterException("EXT_MCL_KA, EXT_MCL_KD, and EXT_MCL_QMAX have to have the same size");
return true;
}
Now, we move on to the part where the flux and Jacobian of the Langmuir implemented will be implemented. To get familiarized with the CADET source code, the physical quantities in an adsorption isotherm is represented by the following variables:
c_{p,i} (liquid phase concentration of ith component) \rightarrow yCp[i]
q_i (solid phase conentration of ith component) \rightarrow y[i]
q_{max,i} (maximum binding capacity of ith component) \rightarrow p->qMax[i]
The adsorption flux related to Langmuir model, Equation (3), should be implemented in the following function:
int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y,
CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const
For the ease of implementation, the summation term \left(1 - \sum_{j=0}^{N_{comp}-1} q_j/q_{max,j}\right) in Equation (3) is defined and stored in a variable defined in Line 117 as:
ResidualType qsum = 1
where ResidualType is the predefined data type of the variable. The following lines of code defines the summation term:
ResidualType qSum = 1.0;
unsigned int bndIdx = 0;
for (int i = 0; i < _nComp; ++i)
{
// Skip components without bound states (bound state index bndIdx is not advanced)
if (_nBoundStates[i] == 0)
continue;
qSum -= y[bndIdx] / static_cast<ParamType>(p->qMax[i]);
// Next bound component
++bndIdx;
}
And the complete residual term (Equation 3) is defined in Line 139 as:
res[bndIdx] = static_cast<ParamType>(p->kD[i]) * y[bndIdx] - static_cast<ParamType>(p->kA[i]) * yCp[i] * static_cast<ParamType>(p->qMax[i]) * qSum;
Note: One might notice that the index of the residual variable is not loop iterating variable āiā but ābndIdxā to take into account the fact that not every component is a bounding component.
So, thatās all about implmenting the flux of binding model and you should only make changes in the above discussed code sections.
The next step in the implementation of the binding model is the implementation of the Jacobian of the Langmuir adsorption model. For this purpose, Equations 4, 6 and 7 are will be implemented in the void jacobianImpl() function. Similar to flux implementation, the summation term in Eq. (4) is defined in Lines 156-168, such as:
double qSum = 1.0;
int bndIdx = 0;
for (int i = 0; i < _nComp; ++i)
{
// Skip components without bound states (bound state index bndIdx is not advanced)
if (_nBoundStates[i] == 0)
continue;
qSum -= y[bndIdx] / static_cast<double>(p->qMax[i]);
// Next bound component
++bndIdx;
}
After that, the Eq. (4) \left(\frac{\partial F_i}{\partial c_{p,i}}\right) is defined in the Jacobian matrix (jac[]) in Line 181 in the following manner:
// dres_i / dc_{p,i}
jac[i - bndIdx - offsetCp] = -static_cast<double>(p->kA[i]) * static_cast<double>(p->qMax[i]) * qSum;
where āiā is the iteration index over the number of components, ābndIdxā is the index of the bounding component and āoffsetCpā is the offset in the Jacobian matrix corresponding to the pore-phase concentration. To reach the Jacobian matrix entry corresponding to the c_{p,i}, -bndIdx takes us to q_0, another offset with āoffsetCp takes us to c_{p,0} and then finally an increment with +i leads to c_{p,i}. In short jac[i ā bndIdx - offsetCp] corresponds to c_{p,i}.
Similarly, Eq. (7) \left(\frac{\partial F_i}{\partial q_{j}}\right) is implemented in the Line 194 such that:
// dres_i / dq_j
jac[bndIdx2 - bndIdx] = static_cast<double>(p->kA[i]) * yCp[i] * static_cast<double>(p->qMax[i]) / static_cast<double>(p->qMax[j]);
To reach the Jacobian matrix entry corresponding to the q_{j}, -bndIdx takes us to q_0, another offset with +bndIdx2 to q_j. In short jac[bndIdx2 - bndIdx] corresponds to q_{j}.
Finally, Eq. (6) \left(\frac{\partial F_i}{\partial q_{i}}\right) is implemented in Line 201 in the following manner:
// Add to dres_i / dq_i
jac[0] += kd;
where k_d is a local variable defined in Line 178.
When implementing any other multi-component binding model where components are interacting with each other, it is expected that the changes are made in only the places that are discussed above.
Last but not the least, do not forget to make modification in Lines 211-221 which are given as:
typedef LangmuirBindingBase<LangmuirParamHandler> LangmuirBinding;
typedef LangmuirBindingBase<ExtLangmuirParamHandler> ExternalLangmuirBinding;
namespace binding
{
void registerLangmuirModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings)
{
bindings[LangmuirBinding::identifier()] = []() { return new LangmuirBinding(); };
bindings[ExternalLangmuirBinding::identifier()] = []() { return new ExternalLangmuirBinding(); };
}
}
In the above code, replace the string Langmuir with your own binding model name. For example, if your binding model is Example then the above piece of code would like this:
typedef LangmuirBindingBase<LangmuirParamHandler> LangmuirBinding;
typedef LangmuirBindingBase<ExtLangmuirParamHandler> ExternalLangmuirBinding;
namespace binding
{
void registerLangmuirModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings)
{
bindings[LangmuirBinding::identifier()] = []() { return new LangmuirBinding(); };
bindings[ExternalLangmuirBinding::identifier()] = []() { return new ExternalLangmuirBinding(); };
}
} // namespace binding
One should take note, that this example is focussed on implementing multi-component Langmuir isotherm and thatās why naming convention relating to Langmuir is adopted in the implementation part of this tutorial. In case of implemeting any arbitrary binding model (for example: Example) one should give great care to change the name Langmuir with your binding model name. Any mistake in this regard can cause an error when you try to compile your code.
