From 56975a7e0c6ee92cec5342a9e036fab72914f75f Mon Sep 17 00:00:00 2001
From: David Parkhurst <dlparkhurst@gmail.com>
Date: Sat, 17 Jul 2021 08:49:17 -0600
Subject: [PATCH] Saved surface for numerical derivatives

---
 model.cpp  |  9 +++++++++
 pitzer.cpp | 27 ++++++++++++++++++---------
 sit.cpp    | 10 +++++++++-
 3 files changed, 36 insertions(+), 10 deletions(-)

diff --git a/model.cpp b/model.cpp
index 88d33609..39955f26 100644
--- a/model.cpp
+++ b/model.cpp
@@ -5560,6 +5560,7 @@ numerical_jacobian(void)
 	std::vector<class phase> base_phases;
 	double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
 	cxxGasPhase base_gas_phase;
+	cxxSurface base_surface;
 
 	if (!
 		(numerical_deriv ||
@@ -5571,6 +5572,10 @@ numerical_jacobian(void)
 
 	calculating_deriv = TRUE;
 	//jacobian_sums();
+	if (use.Get_surface_ptr() != NULL)
+	{
+		base_surface = *use.Get_surface_ptr();
+	}
 	if (use.Get_gas_phase_ptr() != NULL)
 	{
 		//cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
@@ -5761,6 +5766,10 @@ numerical_jacobian(void)
 			x[i]->moles -= d2;
 			break;
 		}
+		if (use.Get_surface_ptr() != NULL)
+		{
+			base_surface = *use.Get_surface_ptr();
+		}
 		if (use.Get_gas_phase_ptr() != NULL)
 		{
 			*use.Get_gas_phase_ptr() = base_gas_phase;
diff --git a/pitzer.cpp b/pitzer.cpp
index 21d667aa..b84aa850 100644
--- a/pitzer.cpp
+++ b/pitzer.cpp
@@ -1951,6 +1951,7 @@ jacobian_pz(void)
 	std::vector<class phase*> phase_ptrs;
 	std::vector<class phase> base_phases;
 	cxxGasPhase base_gas_phase;
+	cxxSurface base_surface;
 	LDBLE d, d1, d2;
 	int i, j;
 Restart:
@@ -1962,6 +1963,10 @@ Restart:
 	}
 	mb_sums();
 	residuals();
+	if (use.Get_surface_ptr() != NULL)
+	{
+		base_surface = *use.Get_surface_ptr();
+	}
 	if (use.Get_gas_phase_ptr() != NULL)
 	{
 		cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
@@ -2141,6 +2146,10 @@ Restart:
 			reset();
 			break;
 		}
+		if (use.Get_surface_ptr() != NULL)
+		{
+			base_surface = *use.Get_surface_ptr();
+		}
 		if (use.Get_gas_phase_ptr() != NULL)
 		{
 			*use.Get_gas_phase_ptr() = base_gas_phase;
@@ -2155,15 +2164,15 @@ Restart:
 		mb_sums();
 		residuals();
 	}
-	//for (i = 0; i < count_unknowns; i++)
-	//{
-	//	//Debugging
-	//	if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
-	//		fabs(residual[i]) + fabs(base[i]) > 1e-8)
-	//	{
-	//		std::cerr << i << ": " << x[i]->description << "  " << residual[i] << "  " << base[i] << std::endl;
-	//	}
-	//}
+	for (i = 0; i < count_unknowns; i++)
+	{
+		//Debugging
+		if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
+			fabs(residual[i]) + fabs(base[i]) > 1e-6)
+		{
+			std::cerr << i << ": " << x[i]->description << "  " << residual[i] << "  " << base[i] << std::endl;
+		}
+	}
 	base.clear();
 	calculating_deriv = 0;
 	return OK;
diff --git a/sit.cpp b/sit.cpp
index a1b140b6..ec6d328f 100644
--- a/sit.cpp
+++ b/sit.cpp
@@ -967,7 +967,12 @@ jacobian_sit(void)
 	std::vector<class phase> base_phases;
 	double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
 	cxxGasPhase base_gas_phase;
+	cxxSurface base_surface;
 Restart:
+	if (use.Get_surface_ptr() != NULL)
+	{
+		base_surface = *use.Get_surface_ptr();
+	}
 	if (use.Get_gas_phase_ptr() != NULL)
 	{
 		cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
@@ -1120,7 +1125,10 @@ Restart:
 			x[i]->moles -= d2;
 			break;
 		}
-
+		if (use.Get_surface_ptr() != NULL)
+		{
+			base_surface = *use.Get_surface_ptr();
+		}
 		if (use.Get_gas_phase_ptr() != NULL)
 		{
 			*use.Get_gas_phase_ptr() = base_gas_phase;