Add tests for flux Jacobians of column models

Adds tests that check the flux part of the Jacobians of column models
against finite differences.

test/ColumnTests.cpp

@@ -525,6 +525,49 @@ namespace column
 		destroyModelBuilder(mb);
 	}
 
+	void testArrowHeadJacobianFD(const std::string& uoType, double h, double absTol, double relTol)
+	{
+		cadet::JsonParameterProvider jpp = createColumnWithTwoCompLinearBinding(uoType);
+
+		cadet::IModelBuilder* const mb = cadet::createModelBuilder();
+		REQUIRE(nullptr != mb);
+
+		cadet::IUnitOperation* const unitAna = createAndConfigureUnit(uoType, *mb, jpp, cadet::Weno::maxOrder());
+		cadet::IUnitOperation* const unitFD = createAndConfigureUnit(uoType, *mb, jpp, cadet::Weno::maxOrder());
+
+		// Obtain offset to fluxes
+		const unsigned int fluxOffset = fluxOffsetOfColumnUnitOp(unitFD);
+
+		// Setup matrices
+		unitAna->notifyDiscontinuousSectionTransition(0.0, 0u, nullptr, nullptr, 0u);
+		unitFD->notifyDiscontinuousSectionTransition(0.0, 0u, nullptr, nullptr, 0u);
+
+		// Obtain memory for state, Jacobian multiply direction, Jacobian column
+		std::vector<double> y(unitFD->numDofs(), 0.0);
+		std::vector<double> jacDir(unitFD->numDofs(), 0.0);
+		std::vector<double> jacCol1(unitFD->numDofs(), 0.0);
+		std::vector<double> jacCol2(unitFD->numDofs(), 0.0);
+
+		// Fill state vector with some values
+		util::populate(y.data(), [](unsigned int idx) { return std::abs(std::sin(idx * 0.13)) + 1e-4; }, unitAna->numDofs());
+//		util::populate(y.data(), [](unsigned int idx) { return 1.0; }, unitAna->numDofs());
+
+		// Compute state Jacobian
+		unitAna->residualWithJacobian(0.0, 0u, 1.0, y.data(), nullptr, jacDir.data(), nullptr, nullptr, 0u);
+		unitFD->residualWithJacobian(0.0, 0u, 1.0, y.data(), nullptr, jacDir.data(), nullptr, nullptr, 0u);
+		std::fill(jacDir.begin(), jacDir.end(), 0.0);
+
+		// Compare Jacobians
+		cadet::test::compareJacobianArrowHeadFD(
+			[=](double const* lDir, double* res) -> void { unitFD->residual(0.0, 0u, 1.0, lDir, nullptr, res); }, 
+			[&](double const* lDir, double* res) -> void { unitAna->multiplyWithJacobian(0.0, 0u, 1.0, y.data(), nullptr, lDir, 1.0, 0.0, res); }, 
+			y.data(), jacDir.data(), jacCol1.data(), jacCol2.data(), unitFD->numDofs(), fluxOffset, h, absTol, relTol);
+
+		mb->destroyUnitOperation(unitAna);
+		mb->destroyUnitOperation(unitFD);
+		destroyModelBuilder(mb);
+	}
+
 	void testFwdSensJacobians(const std::string& uoType, double h, double absTol, double relTol)
 	{
 		cadet::IModelBuilder* const mb = cadet::createModelBuilder();

test/ColumnTests.hpp

@@ -128,6 +128,16 @@ namespace column
 	 */
 	void testTimeDerivativeJacobianFD(const std::string& uoType, double h = 1e-6, double absTol = 0.0, double relTol = std::numeric_limits<float>::epsilon() * 100.0);
 
+	/**
+	 * @brief Checks the bottom macro row and right macro column of the Jacobian against FD
+	 * @details Uses centered finite differences to check the flux part of the Jacobian.
+	 * @param [in] uoType Unit operation type
+	 * @param [in] h Step size of centered finite differences
+	 * @param [in] absTol Absolute error tolerance
+	 * @param [in] relTol Relative error tolerance
+	 */
+	void testArrowHeadJacobianFD(const std::string& uoType, double h = 1e-6, double absTol = 0.0, double relTol = std::numeric_limits<float>::epsilon() * 100.0);
+
 	/**
 	 * @brief Checks the forward sensitivity residual using analytic Jacobians
 	 * @details Uses centered finite differences.

test/GeneralRateModel.cpp

@@ -50,6 +50,11 @@ TEST_CASE("GRM time derivative Jacobian vs FD", "[GRM],[UnitOp],[Residual],[Jaco
 	cadet::test::column::testTimeDerivativeJacobianFD("GENERAL_RATE_MODEL", 1e-6, 0.0, 9e-4);
 }
 
+TEST_CASE("GRM flux Jacobian vs FD", "[GRM],[UnitOp],[Residual],[Jacobian]")
+{
+	cadet::test::column::testArrowHeadJacobianFD("GENERAL_RATE_MODEL", 1e-6, 2e-9);
+}
+
 TEST_CASE("GRM sensitivity Jacobians", "[GRM],[UnitOp],[Sensitivity]")
 {
 	cadet::test::column::testFwdSensJacobians("GENERAL_RATE_MODEL", 1e-4, 6e-7);

test/LumpedRateModelWithPores.cpp

@@ -50,6 +50,11 @@ TEST_CASE("LRMP time derivative Jacobian vs FD", "[LRMP],[UnitOp],[Residual],[Ja
 	cadet::test::column::testTimeDerivativeJacobianFD("LUMPED_RATE_MODEL_WITH_PORES");
 }
 
+TEST_CASE("LRMP flux Jacobian vs FD", "[LRMP],[UnitOp],[Residual],[Jacobian]")
+{
+	cadet::test::column::testArrowHeadJacobianFD("LUMPED_RATE_MODEL_WITH_PORES");
+}
+
 TEST_CASE("LRMP sensitivity Jacobians", "[LRMP],[UnitOp],[Sensitivity]")
 {
 	cadet::test::column::testFwdSensJacobians("LUMPED_RATE_MODEL_WITH_PORES", 1e-4, 6e-7);