Fix Jacobians and equations for LDF binding models

Fixes several problems with the model equations and Jacobians of the LDF binding models. Also fixes the tests, which all pass now.
modsim · Dec 16, 2022 · 160b4f2 · 160b4f2
1 parent 9f5dca9
commit 160b4f2
Show file tree

Hide file tree

Showing 5 changed files with 75 additions and 84 deletions.
diff --git a/src/libcadet/model/binding/BiLangmuirLDFBinding.cpp b/src/libcadet/model/binding/BiLangmuirLDFBinding.cpp
@@ -1,7 +1,7 @@
 // =============================================================================
 //  CADET
 //  
-//  Copyright © 2008-2021: The CADET Authors
+//  Copyright © 2008-2022: The CADET Authors
 //            Please see the AUTHORS and CONTRIBUTORS file.
 //  
 //  All rights reserved. This program and the accompanying materials
@@ -245,7 +245,7 @@ namespace cadet
 					active const* const localKkin = p->kkin[site];
 					active const* const localQmax = p->qMax[site];
 
-					ResidualType qSum = 1.0;
+					ResidualType cpSum = 1.0;
 					unsigned int bndIdx = 0;
 
 					// bndIdx is used as a counter inside one binding site type
@@ -257,7 +257,7 @@ namespace cadet
 						if (_nBoundStates[i] == 0)
 							continue;
 
-						qSum += yCp[i]* static_cast<ParamType>(localKeq[i]);
+						cpSum += yCp[i] * static_cast<ParamType>(localKeq[i]);
 
 						// Next bound component
 						++bndIdx;
@@ -271,7 +271,7 @@ namespace cadet
 							continue;
 
 						// Residual
-						res[bndIdx * nSites] = static_cast<ParamType>(localKkin[i]) * y[bndIdx * nSites] - (static_cast<ParamType>(localKkin[i])*static_cast<ParamType>(localKeq[i]) * yCp[i] * static_cast<ParamType>(localQmax[i]) / qSum);
+						res[bndIdx * nSites] = static_cast<ParamType>(localKkin[i]) * (y[bndIdx * nSites] - static_cast<ParamType>(localKeq[i]) * yCp[i] * static_cast<ParamType>(localQmax[i]) / cpSum);
 
 						// Next bound component
 						++bndIdx;
@@ -304,15 +304,15 @@ namespace cadet
 					active const* const localKkin = p->kkin[site];
 					active const* const localQmax = p->qMax[site];
 
-					double qSum = 1.0;
+					double cpSum = 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 += yCp[i] * static_cast<double>(localKeq[i]);
+						cpSum += yCp[i] * static_cast<double>(localKeq[i]);
 
 						// Next bound component
 						++bndIdx;
@@ -329,12 +329,13 @@ namespace cadet
 						const double kkin = static_cast<double>(localKkin[i]);
 
 						// dres_i / dc_{p,i}
-						jac[i - site - offsetCp - nSites * bndIdx] = -kkin*keq * static_cast<double>(localQmax[i]) / qSum;
+						jac[i - site - offsetCp - nSites * bndIdx] = -kkin * keq * static_cast<double>(localQmax[i]) / cpSum;
 						// Getting to c_{p,i}: -nSites * bndIdx takes us to q_{0,site}, another -site to q_{0,0}. From there, we
 						//                     take a -offsetCp to reach c_{p,0} and a +i to arrive at c_{p,i}.
 						//                     This means jac[i - site - offsetCp - nSites * bndIdx] corresponds to c_{p,i}.
 
 						// Fill dres_i / dc_j
+						const double commonFactor = kkin * keq * static_cast<double>(localQmax[i]) * yCp[i] / (cpSum * cpSum);
 						int bndIdx2 = 0;
 						for (int j = 0; j < _nComp; ++j)
 						{
@@ -343,7 +344,7 @@ namespace cadet
 								continue;
 
 							// dres_i / dc_j
-							jac[j - site - offsetCp - nSites * bndIdx] += static_cast<double>(localKeq[j])*kkin * keq * static_cast<double>(localQmax[i]) *yCp[i]/ (qSum*qSum);
+							jac[j - site - offsetCp - nSites * bndIdx] += static_cast<double>(localKeq[j]) * commonFactor;
 							// Getting to c_{p,j}: -nSites * bndIdx takes us to q_{0,site}, another -site to q_{0,0}. From there, we
 							//                     take a -offsetCp to reach c_{p,0} and a +j to arrive at c_{p,j}.
 							//                     This means jac[j - site - offsetCp - nSites * bndIdx] corresponds to c_{p,j}.

diff --git a/src/libcadet/model/binding/FreundlichLDFBinding.cpp b/src/libcadet/model/binding/FreundlichLDFBinding.cpp
@@ -1,7 +1,7 @@
 // =============================================================================
 //  CADET
 //  
-//  Copyright © 2008-2021: The CADET Authors
+//  Copyright © 2008-2022: The CADET Authors
 //            Please see the AUTHORS and CONTRIBUTORS file.
 //  
 //  All rights reserved. This program and the accompanying materials
@@ -97,46 +97,42 @@ namespace cadet
 			using ParamHandlerBindingModelBase<ParamHandler_t>::_nComp;
 			using ParamHandlerBindingModelBase<ParamHandler_t>::_nBoundStates;
 
-			const double threshold = 1e-14;
+			const double _threshold = 1e-14;
 
 			template <typename StateType, typename CpStateType, typename ResidualType, typename ParamType>
 			int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y,
 				CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const
 			{
 				using std::abs;
+				using std::pow;
+
 				typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
 
-				/* Protein Flux: dq/dt = k_kin * (q* – q) with
-				q* = k_F * c^(1/n)
-				*/
+				// Protein Flux: dq/dt = k_kin * (q* - q) with q* = k_F * c^(1/n)
 				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;
+
 					const ParamType n_param = static_cast<ParamType>(p->n[i]);
 					const ParamType kF = static_cast<ParamType>(p->kF[i]);
 					const ParamType kkin = static_cast<ParamType>(p->kkin[i]);
-
-					const ParamType alpha_1 = (((2 * n_param - 1) / n_param) * kF * pow(threshold, (1 - n_param) / n_param));
-					const ParamType alpha_2 = (((1 - n_param) / n_param) * kF * pow(threshold, (1 - 2 * n_param) / n_param));
 
 					// Residual
-					if (n_param > 1)
+					if ((n_param > 1) && (abs(yCp[i]) < _threshold))
 					{
-						if (abs(yCp[i]) < threshold)
-							res[bndIdx] = kkin * (y[bndIdx] - kkin* alpha_1 * yCp[i] - alpha_2 * yCp[i] * yCp[i]);
-						else
-							res[bndIdx] = kkin * (y[bndIdx] -  kF * pow(abs(yCp[i]), 1.0 / n_param));
+						const ParamType alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param);
+						const ParamType alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param);
+						res[bndIdx] = kkin * (y[bndIdx] - yCp[i] * (alpha_1 + alpha_2 * yCp[i]));
 					}
 					else
 					{
-						res[bndIdx] = kkin * (y[bndIdx] -  kF * pow(abs(yCp[i]), 1.0 / n_param));
+						res[bndIdx] = kkin * (y[bndIdx] - kF * pow(abs(yCp[i]), 1.0 / n_param));
 					}
 
-
 					// Next bound component
 					++bndIdx;
 				}
@@ -147,43 +143,37 @@ namespace cadet
 			template <typename RowIterator>
 			void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const
 			{
-				typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
 				using std::abs;
+				using std::pow;
+
+				typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
 
 				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;
-					// dres / dq_i
+
 					const double n_param = static_cast<double>(p->n[i]);
 					const double kF = static_cast<double>(p->kF[i]);
 					const double kkin = static_cast<double>(p->kkin[i]);
 
-					double const alpha_1 = ((2 * n_param - 1) / n_param) * kF * pow(threshold, (1 - n_param) / (n_param));
-					double const alpha_2 = ((1 - n_param) / n_param) * kF * pow(threshold, (1 - 2 * n_param) / (n_param));
+					// dres / dq_i
+					jac[0] = kkin;
 
-					jac[0] = static_cast<double>(p->kkin[i]);
 					// dres / dc_{p,i}
-					// This isotherm is non-differentiable at yCp = 0 so following segment of code deals with mitigating this issue.
-					if (n_param > 1)
+					if ((n_param > 1) && (abs(yCp[i]) < _threshold))
 					{
-						if (abs(yCp[i]) < threshold)
-							jac[i - bndIdx - offsetCp] = -alpha_1 - 2 * alpha_2 * yCp[i];
-						else
-							jac[i - bndIdx - offsetCp] = -(1.0 / n_param) * kkin * kF * pow(abs(yCp[i]), (1.0 - n_param) / n_param);
+						double const alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param);
+						double const alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param);
+						jac[i - bndIdx - offsetCp] = -kkin * (alpha_1 + 2.0 * alpha_2 * yCp[i]);
 					}
 					else
 					{
 						jac[i - bndIdx - offsetCp] = -(1.0 / n_param) * kkin * kF * pow(abs(yCp[i]), (1.0 - n_param) / n_param);
 					}
 
-					// The distance from liquid phase to solid phase is reduced for each non-binding component
-					// since a bound state is neglected. The number of neglected bound states so far is i - bndIdx.
-					// Thus, by going back offsetCp - (i - bndIdx) = -[ i - bndIdx - offsetCp ] we get to the corresponding
-					// liquid phase component.
-
 					++bndIdx;
 					++jac;
 				}

diff --git a/src/libcadet/model/binding/LangmuirLDFBinding.cpp b/src/libcadet/model/binding/LangmuirLDFBinding.cpp
@@ -1,7 +1,7 @@
 // =============================================================================
 //  CADET
 //  
-//  Copyright © 2008-2021: The CADET Authors
+//  Copyright © 2008-2022: The CADET Authors
 //            Please see the AUTHORS and CONTRIBUTORS file.
 //  
 //  All rights reserved. This program and the accompanying materials
@@ -107,8 +107,8 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 		std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
 
 		// Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j}
-		double qSum = 1.0;
-		double qSumT = 0.0;
+		double cpSum = 1.0;
+		double cpSumT = 0.0;
 		unsigned int bndIdx = 0;
 		for (int i = 0; i < _nComp; ++i)
 		{
@@ -117,8 +117,8 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 				continue;
 
 			const double summand = y[bndIdx] / static_cast<double>(p->qMax[i]);
-			qSum -= summand;
-			qSumT += summand / static_cast<double>(p->qMax[i]) * static_cast<double>(dpDt->qMax[i]);
+			cpSum -= summand;
+			cpSumT += summand / static_cast<double>(p->qMax[i]) * static_cast<double>(dpDt->qMax[i]);
 
 			// Next bound component
 			++bndIdx;
@@ -133,9 +133,9 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 
 			// Residual
 			dResDt[bndIdx] = static_cast<double>(dpDt->kD[i]) * y[bndIdx] 
-				- yCp[i] * (static_cast<double>(dpDt->kA[i]) * static_cast<double>(p->qMax[i]) * qSum
-				           + static_cast<double>(p->kA[i]) * static_cast<double>(dpDt->qMax[i]) * qSum
-				           + static_cast<double>(p->kA[i]) * static_cast<double>(p->qMax[i]) * qSumT);
+				- yCp[i] * (static_cast<double>(dpDt->kA[i]) * static_cast<double>(p->qMax[i]) * cpSum
+				           + static_cast<double>(p->kA[i]) * static_cast<double>(dpDt->qMax[i]) * cpSum
+				           + static_cast<double>(p->kA[i]) * static_cast<double>(p->qMax[i]) * cpSumT);
 
 			// Next bound component
 			++bndIdx;
@@ -157,15 +157,15 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 		typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
 
 		// Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j}
-		ResidualType qSum = 1.0;
+		ResidualType cpSum = 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 += yCp[i]* static_cast<ParamType>(p->keq[i]); //\sum b_j*c_j
+			cpSum += yCp[i] * static_cast<ParamType>(p->keq[i]);
 
 			// Next bound component
 			++bndIdx;
@@ -179,7 +179,7 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 				continue;
 
 			// Residual
-			res[bndIdx] = static_cast<ParamType>(p->kkin[i]) * (y[bndIdx] - (static_cast<ParamType>(p->keq[i]) * yCp[i] * static_cast<ParamType>(p->qMax[i]) / qSum));
+			res[bndIdx] = static_cast<ParamType>(p->kkin[i]) * (y[bndIdx] - static_cast<ParamType>(p->keq[i]) * yCp[i] * static_cast<ParamType>(p->qMax[i]) / cpSum);
 
 			// Next bound component
 			++bndIdx;
@@ -194,15 +194,15 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 		typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
 
 		// Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j}
-		double qSum = 1.0;
+		double cpSum = 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 += yCp[i] * static_cast<double>(p->keq[i]);
+			cpSum += yCp[i] * static_cast<double>(p->keq[i]);
 
 			// Next bound component
 			++bndIdx;
@@ -217,13 +217,14 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 
 			const double keq = static_cast<double>(p->keq[i]);
 			const double kkin = static_cast<double>(p->kkin[i]);
-			//(keq *keq *static_cast<double>(p->qMax[i]) * yCp[i] / (qSum*qSum))
+			//(keq *keq *static_cast<double>(p->qMax[i]) * yCp[i] / (cpSum*cpSum))
 			// dres_i / dc_{p,i}
-			jac[i - bndIdx - offsetCp] = -kkin*( (keq* static_cast<double>(p->qMax[i]) /qSum)  );
+			jac[i - bndIdx - offsetCp] = -kkin * keq * static_cast<double>(p->qMax[i]) / cpSum;
 			// Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}.
 			//                     This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}.
 
 			// Fill dres_i / dc_j
+			const double commonFactor = keq * kkin * static_cast<double>(p->qMax[i]) * yCp[i] / (cpSum * cpSum);
 			int bndIdx2 = 0;
 			for (int j = 0; j < _nComp; ++j)
 			{
@@ -232,7 +233,7 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase<ParamHandler_
 					continue;
 
 				// dres_i / dc_j
-				jac[j - bndIdx - offsetCp] += static_cast<double>(p->keq[j])*keq * kkin * static_cast<double>(p->qMax[i])* yCp[i] / (qSum*qSum);
+				jac[j - bndIdx - offsetCp] += static_cast<double>(p->keq[j]) * commonFactor;
 				// Getting to c_{p,j}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +j to c_{p,j}.
 				//                     This means jac[j - bndIdx - offsetCp] corresponds to c_{p,j}.
 

diff --git a/src/libcadet/model/binding/LangmuirLDFCBinding.cpp b/src/libcadet/model/binding/LangmuirLDFCBinding.cpp
@@ -1,7 +1,7 @@
 // =============================================================================
 //  CADET
 //  
-//  Copyright © 2008-2021: The CADET Authors
+//  Copyright © 2008-2022: The CADET Authors
 //            Please see the AUTHORS and CONTRIBUTORS file.
 //  
 //  All rights reserved. This program and the accompanying materials
@@ -165,7 +165,7 @@ class LangmuirLDFCBindingBase : public ParamHandlerBindingModelBase<ParamHandler
 			if (_nBoundStates[i] == 0)
 				continue;
 
-			qSum -= y[bndIdx]/ static_cast<ParamType>(p->qMax[i]); 
+			qSum -= y[bndIdx] / static_cast<ParamType>(p->qMax[i]); 
 
 			// Next bound component
 			++bndIdx;
@@ -179,7 +179,7 @@ class LangmuirLDFCBindingBase : public ParamHandlerBindingModelBase<ParamHandler
 				continue;
 
 			// Residual
-			res[bndIdx] = -static_cast<ParamType>(p->kkin[i]) * (yCp[bndIdx] - ( y[bndIdx] / (static_cast<ParamType>(p->keq[i]) * static_cast<ParamType>(p->qMax[i]) - qSum) ));
+			res[bndIdx] = -static_cast<ParamType>(p->kkin[i]) * (yCp[i] - y[bndIdx] / (static_cast<ParamType>(p->keq[i]) * static_cast<ParamType>(p->qMax[i]) * qSum));
 
 			// Next bound component
 			++bndIdx;
@@ -218,13 +218,13 @@ class LangmuirLDFCBindingBase : public ParamHandlerBindingModelBase<ParamHandler
 			const double keq = static_cast<double>(p->keq[i]);
 			const double kkin = static_cast<double>(p->kkin[i]);
 			const double qMax = static_cast<double>(p->qMax[i]);
-			//(keq *keq *static_cast<double>(p->qMax[i]) * yCp[i] / (qSum*qSum))
+
 			// dres_i / dc_{p,i}
-			jac[i - bndIdx - offsetCp] = -kkin;//*( (keq* static_cast<double>(p->qMax[i]) /qSum)  );
+			jac[i - bndIdx - offsetCp] = -kkin;
 			// Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}.
 			//                     This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}.
 
-			// Fill dres_i / dc_j
+			// Fill dres_i / dq_j
 			int bndIdx2 = 0;
 			for (int j = 0; j < _nComp; ++j)
 			{
@@ -233,16 +233,15 @@ class LangmuirLDFCBindingBase : public ParamHandlerBindingModelBase<ParamHandler
 					continue;
 
 				// dres_i / dq_j 
-				jac[bndIdx2 - bndIdx] += -y[i] * kkin / (keq * qMax * static_cast<double>(p->qMax[j]) * (1 - qSum) * (1 - qSum));
+				jac[bndIdx2 - bndIdx] += y[bndIdx] * kkin / (keq * qMax * static_cast<double>(p->qMax[j]) * qSum * qSum);
 				// Getting to q_j: -bndIdx takes us to q_{0}, another +bndIdx2 to q_{j}.
-					// This means jac[bndIdx2 - bndIdx] corresponds to q_{j}.
-
+				// This means jac[bndIdx2 - bndIdx] corresponds to q_{j}.
 
 				++bndIdx2;
 			}
 
 			// Add to dres_i / dq_i
-			jac[0] += kkin / (keq * qMax  * (1 - qSum));
+			jac[0] += kkin / (keq * qMax * qSum);
 
 			// Advance to next flux and Jacobian row
 			++bndIdx;