Dear Mr Breuer,
Dear Mr Leweke,
thank you very much for your anwers! To add some context this is what my script / the component and reaction system looks like now:
ComponentSystem cs{
{// Format: Component(name).setCharge(...).setMolarMass(...).setTruesdellJonesAlpha(...).setTruesdellJonesBeta(...)
Component("H2O").setMolarMass(18.01528),
Component("H+").setCharge(1).setMolarMass(1.008).setTruesdellJonesParameters(4.78, 0.24),
Component("OH-").setCharge(-1).setMolarMass(17.007).setTruesdellJonesParameters(10.65, 0.21),
Component("Na+").setCharge(1).setMolarMass(22.990).setTruesdellJonesParameters(4.32, 0.06),
Component("Cl-").setCharge(-1).setMolarMass(35.453).setTruesdellJonesParameters(3.71, 0.01),
Component("TrisH+").setCharge(1).setMolarMass(121.14).setTruesdellJonesParameters(4.0, 0.0),
Component("Tris").setMolarMass(121.14),
Component("AcH").setMolarMass(60.052),
Component("Ac-").setCharge(-1).setMolarMass(59.044).setTruesdellJonesParameters(4.5, 0.0),
Component("GlyH2+").setCharge(1).setMolarMass(
75.067), // TODO: No truesdell_jones paramters available -> better to use Debye-Hückel model or 1?
Component("GlyH").setMolarMass(75.067),
Component("Gly-").setCharge(-1).setMolarMass(74.059).setTruesdellJonesParameters(4.0, 0.0),
Component("hIgG").setCharge(1).setMolarMass(150000), // Human Immunoglobulin G (IgG) antibody, roughly 150 kDa
Component("MNP_OH2+").setType(Type::MagneticNanoParticleGroup).setCharge(1).setMolarMass(18.015),
Component("MNP_OH").setType(Type::MagneticNanoParticleGroup).setMolarMass(17.007),
Component("MNP_O-").setType(Type::MagneticNanoParticleGroup).setCharge(-1).setMolarMass(15.999),
Component("MNP_hIgG").setType(Type::MagneticNanoParticleGroup).setMolarMass(150000),
Component("MNP1").setType(Type::MagneticNanoParticle).setRadius(1039e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP2").setType(Type::MagneticNanoParticle).setRadius(1209e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP3").setType(Type::MagneticNanoParticle).setRadius(1406e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP4").setType(Type::MagneticNanoParticle).setRadius(1635e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP5").setType(Type::MagneticNanoParticle).setRadius(1901e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP6").setType(Type::MagneticNanoParticle).setRadius(2211e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP7").setType(Type::MagneticNanoParticle).setRadius(2571e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP8").setType(Type::MagneticNanoParticle).setRadius(2990e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP9").setType(Type::MagneticNanoParticle).setRadius(3477e-9).setDensity(5170).setMagneticSaturation(3.5e5),
Component("MNP10").setType(Type::MagneticNanoParticle).setRadius(4043e-9).setDensity(5170).setMagneticSaturation(3.5e5)}};
// Other physical parameters
const realtype MNP_Ns = 1.08 * 1e-5; // [mol/m²] - number of hydroxyl groups per area
const realtype MNP_specific_surface = 1e5; // m²/kg
// =============================================================
// ========================== REACTIONS ========================
// =============================================================
TruesdellJonesActivityModel truesdellJonesActivityModel(cs);
ReactionSystem reactionSystem(cs, truesdellJonesActivityModel);
realtype K_W = std::pow(10, -14); // [mol²/(m³)²]
realtype c_W = 1000 / cs("H2O").molarMass; // [mol/L] - Concentration of water
realtype Ks_W = K_W / c_W * 1e3; // [mol/m³] - Equilibrium constant for water dissociation in mol/m³
realtype Ks_Tris = std::pow(10, -8.3) * 1e3; // [mol/m³]
realtype Ks_AcH = std::pow(10, -4.75) * 1e3; // [mol/m³]
realtype Ks_GlyH2 = std::pow(10, -2.35) * 1e3; // [mol/m³]
realtype Ks_GlyH = std::pow(10, -9.78) * 1e3; // [mol/m³]
realtype Ks_MNP_OH2_plus = std::pow(10, -(7.9 - 0.5 * 2.5)) * 1e3; // [mol/m³]
realtype Ks_MNP_OH = std::pow(10, -(7.9 + 0.5 * 2.5)) * 1e3; // [mol/m³]
realtype kf_ion = 1e2; // 5.5 (for 10 cells); // [1/s]
realtype kb_W = (kf_ion / Ks_W); // [m³/(mol*s)]
realtype kb_Tris = (kf_ion / Ks_Tris); // [m³/(mol*s)]
realtype kb_AcH = (kf_ion / Ks_AcH); // [m³/(mol*s)]
realtype kb_GlyH2 = (kf_ion / Ks_GlyH2); // [m³/(mol*s)]
realtype kb_GlyH = (kf_ion / Ks_GlyH); // [m³/(mol*s)]
realtype kb_MNP_OH2_plus = (kf_ion / Ks_MNP_OH2_plus); // [m³/(mol*s)]
realtype kb_MNP_OH = (kf_ion / Ks_MNP_OH); // [m³/(mol*s)]
reactionSystem.add(massActionLaw(cs, "H2O <=> H+ + OH-", kf_ion,
kb_W)); // H2O concentration is included in the equilibrium constant
reactionSystem.add(massActionLaw(cs, "TrisH+ <=> Tris + H+", kf_ion, kb_Tris));
reactionSystem.add(massActionLaw(cs, "AcH <=> H+ + Ac-", kf_ion, kb_AcH));
reactionSystem.add(massActionLaw(cs, "GlyH2+ <=> H+ + GlyH", kf_ion, kb_GlyH2));
reactionSystem.add(massActionLaw(cs, "GlyH <=> H+ + Gly-", kf_ion, kb_GlyH));
Plus two additional ion reactions between the hydroxl groups of the Magnetic Nano Particles (MNP) which include an extra correlation factor derived after the Gouy-Chapman model (more details on that: Modeling magnetic separation for hIgG purification: A theoretical and experimental study - ScienceDirect ).
Most importantly, I want to model the adsorption from hIgG to the MNP for which we use the Langmuir model with pH depented K_b and q_max. The other reactions are only needed to determine the pH value for this adsorption reaction.
So far, I added normal ODE terms as such:
Reaction massActionLaw(std::vector<std::pair<size_t, int>> reactant_terms,
std::vector<std::pair<size_t, int>> product_terms,
realtype k_forward,
realtype k_backward) {
return Reaction([reactant_terms, product_terms, k_forward, k_backward](realtype t, const ConstArrayMap& activities,
ArrayMap& dc_dt) {
const auto n_cells = activities.rows();
const auto n_components = activities.cols();
// Calculate forward and backward rates for all cells simultaneously
ColVector rate_fwd = ColVector::Constant(n_cells, k_forward);
ColVector rate_bwd = ColVector::Constant(n_cells, k_backward);
// Calculate forward reaction rates (vectorized)
for (const auto& [idx, power] : reactant_terms) {
for (int i = 0; i < power; ++i) {
rate_fwd *= activities.col(idx);
}
}
// Calculate backward reaction rates (vectorized)
for (const auto& [idx, power] : product_terms) {
for (int i = 0; i < power; ++i) {
rate_bwd *= activities.col(idx);
}
}
// Calculate net reaction rates for all cells
ColVector net_rate = rate_fwd - rate_bwd;
// Apply stoichiometric changes (vectorized)
for (const auto& [idx, power] : reactant_terms) {
dc_dt.col(idx) -= power * net_rate;
}
for (const auto& [idx, power] : product_terms) {
dc_dt.col(idx) += power * net_rate;
}
});
}
I used the jacobians via automatic differentiation, so I have no implementations for them since the derivation should be doable but quite difficult trivial because of the activity model, surface complexiation and Langmuir parameter interpolation.
Also, I read the chapter 2.4 of your PhD thesis and find it a very interesting and elegant approach because it solves both problems together in one solver. On the other hand, I am unsure if it is the right approach for me since:
- I surely would have some difficulties to deeply understand and apply this quite advanced approach on my system
- Since I noticed that the other parts of the system are not to stiff, (I can solve the 6500 s of realtime experiment in about 30 s with ERKStep). I am thinking about solving the reactions and the rest in an alternating manner which would not be possible with IDA.
Specifically, I am thinking about modelling my reaction system algebraicly (equilibrium constraints + mass conservation + charge neutrality) which should be quite easy. And after each (or some) ERK step(s) I could compute a residual e.g. r_H+ = H+_current - H+_approx_eq = H+_current - K (HAc / Ac–) (equivalently for the other equations) which estimates the error of the H+ and hIgG-MNP concentrations. If this error surpasses a tolerance I could solve the system with KINSOL for the respective cells which would bring the jacobians down to 27*27 (instead of 2500*2500).
What do you think about this and do you have other ideas / would you recommend using only one solver and e.g. try analytic or sparse jacobians?
And would you recommend to model H2O as a component which (very slightly) varying concentrations or as a constant, which it basically is since we deal we relatively small concentrations, although this means the reactions generate and destroy mass (H+ and OH-) out of nothing?
Thank you!