From bab4763336a426b4b9ac930a8c5455a88a6a21d6 Mon Sep 17 00:00:00 2001
From: Samuel Leweke <s.leweke@fz-juelich.de>
Date: Sun, 19 Feb 2023 16:32:20 +0100
Subject: [PATCH] Fix reaction Jacobian in models without bound phase

The binding cell kernel erroneously omitted the reaction model Jacobian
if no bound states are present in the model. The Jacobian is now added
regardless of the number of bound states.
---
 .../model/parts/BindingCellKernel.hpp         | 46 ++++++++++++-------
 1 file changed, 29 insertions(+), 17 deletions(-)

diff --git a/src/libcadet/model/parts/BindingCellKernel.hpp b/src/libcadet/model/parts/BindingCellKernel.hpp
index d22cfe75f..3b3400aa3 100644
--- a/src/libcadet/model/parts/BindingCellKernel.hpp
+++ b/src/libcadet/model/parts/BindingCellKernel.hpp
@@ -176,32 +176,44 @@ void residualKernel(double t, unsigned int secIdx, const ColumnPosition& colPos,
 			}
 		}
 
-		if (wantJac && (params.nTotalBound > 0))
+		if (wantJac)
 		{
-			BufferedArray<double> fluxSolidJacobian = buffer.template array<double>(params.nTotalBound * (params.nTotalBound + params.nComp));
-			linalg::DenseMatrixView dmv(static_cast<double*>(fluxSolidJacobian), nullptr, params.nTotalBound, params.nTotalBound + params.nComp);
-			dmv.setAll(0.0);
-
-			// static_cast should be sufficient here, but this statement is also analyzed when wantJac = false
-			params.dynReaction->analyticJacobianCombinedAdd(t, secIdx, colPos, reinterpret_cast<double const*>(y - params.nComp), reinterpret_cast<double const*>(y), -1.0, jacBase, dmv.row(0, params.nComp), buffer);
-
-			idx = 0;
-			for (unsigned int comp = 0; comp < params.nComp; ++comp)
+			if (params.nTotalBound > 0)
 			{
-				const double invBetaP = (1.0 - static_cast<double>(params.porosity)) / (params.poreAccessFactor ? static_cast<double>(params.poreAccessFactor[comp]) * static_cast<double>(params.porosity) : static_cast<double>(params.porosity));
+				BufferedArray<double> fluxSolidJacobian = buffer.template array<double>(params.nTotalBound * (params.nTotalBound + params.nComp));
+				linalg::DenseMatrixView dmv(static_cast<double*>(fluxSolidJacobian), nullptr, params.nTotalBound, params.nTotalBound + params.nComp);
+				dmv.setAll(0.0);
 
-				for (unsigned int bnd = 0; bnd < params.nBound[comp]; ++bnd, ++idx)
+				// static_cast should be sufficient here, but this statement is also analyzed when wantJac = false
+				params.dynReaction->analyticJacobianCombinedAdd(t, secIdx, colPos, reinterpret_cast<double const*>(y - params.nComp), reinterpret_cast<double const*>(y), -1.0, jacBase, dmv.row(0, params.nComp), buffer);
+
+				idx = 0;
+				for (unsigned int comp = 0; comp < params.nComp; ++comp)
 				{
-					// Add Jacobian row to mobile phase
-					(jacBase + comp).addArray(dmv.rowPtr(idx), -static_cast<int>(comp), dmv.columns(), invBetaP);
+					const double invBetaP = (1.0 - static_cast<double>(params.porosity)) / (params.poreAccessFactor ? static_cast<double>(params.poreAccessFactor[comp]) * static_cast<double>(params.porosity) : static_cast<double>(params.porosity));
 
-					if (!params.qsReaction[idx])
+					for (unsigned int bnd = 0; bnd < params.nBound[comp]; ++bnd, ++idx)
 					{
-						// Add Jacobian row to solid phase
-						(jacBase + params.nComp + idx).addArray(dmv.rowPtr(idx), -static_cast<int>(params.nComp + idx), dmv.columns(), 1.0);
+						// Add Jacobian row to mobile phase
+						(jacBase + comp).addArray(dmv.rowPtr(idx), -static_cast<int>(comp), dmv.columns(), invBetaP);
+
+						if (!params.qsReaction[idx])
+						{
+							// Add Jacobian row to solid phase
+							(jacBase + params.nComp + idx).addArray(dmv.rowPtr(idx), -static_cast<int>(params.nComp + idx), dmv.columns(), 1.0);
+						}
 					}
 				}
 			}
+			else
+			{
+				// We do not have bound states, but still need to obtain the Jacobian for the liquid phase.
+				// So we pass a row iterator for the solid phase that does not point anywhere and hope that the
+				// reaction model does not interact with it.
+
+				// static_cast should be sufficient here, but this statement is also analyzed when wantJac = false
+				params.dynReaction->analyticJacobianCombinedAdd(t, secIdx, colPos, reinterpret_cast<double const*>(y - params.nComp), reinterpret_cast<double const*>(y), -1.0, jacBase, linalg::DenseBandedRowIterator(), buffer);
+			}
 		}
 	}
 }