⚡ Improved convergence

- Improves convergence of all solvers assuming a value of 1.0 for
  std::pow(0, 0) instead of NaN in df and dg calculations.

include/chemmisol/chemical/solver.h

@@ -542,22 +542,25 @@ namespace chemmisol {
 							T activity = activities[index];
 							// d x^a/dx = a * x^(a-1)
 							if(species_coefficient_in_reactions < 0) {
-								if(activity != T(0)) {
 								d_f[index] = -species_coefficient_in_reactions
-									* std::pow(
+									* d_products;
+								if(species_coefficient_in_reactions < -1) {
+									// Prevents std::pow(0, 0) issues
+									d_f[index] *= std::pow(
 											activity,
 											-species_coefficient_in_reactions-1
-											) * d_products;
-								} else {
-									d_f[index] = -species_coefficient_in_reactions;
+											);
 								}
-
 							} else {
 								d_f[index] = species_coefficient_in_reactions
-									* std::pow(
+									* d_reactives;
+								if(species_coefficient_in_reactions > 1) {
+									// Prevents std::pow(0, 0) issues
+									d_f[index] *= std::pow(
 											activity,
 											species_coefficient_in_reactions-1
-											) * d_reactives;
+											);
+								}
 							}
 						}
 					}

src/chemmisol/chemical/solver.cpp

@@ -57,8 +57,9 @@ namespace chemmisol {
 				dg[i].resize(reduced_system.fxSize());
 				// All other derivatives are equal to 0.0
 				if(reduced_activities[i] != CX::value_type(0)) {
-					dg[i][i] = a[i] * degrees[i]
-						* std::pow(reduced_activities[i], degrees[i]-1);
+					dg[i][i] = a[i] * degrees[i];
+					if(degrees[i] > 1)
+						dg[i][i] *= std::pow(reduced_activities[i], degrees[i]-1);
 				}
 			}
 			return dg;

tests/chemmisol/tests/chemical/system/agcl_test.cpp

@@ -380,26 +380,27 @@ TEST_F(AgClChemicalSystemTest, solve_equilibrium_absolute_newton) {
 	checkPh(default_ph);
 }
 
-TEST_F(AgClChemicalSystemTest, solve_equilibrium_homotopy) {
-	solver::HomotopyContinuation<std::minstd_rand> homotopy(
-			std::minstd_rand {}, 1000, 100
-			);
+class AgClChemicalSystemHomotopyTest :
+	public AgClChemicalSystemTest, public WithParamInterface<int> {
+		protected:
+			// seeds are defined in utils.h
+			solver::HomotopyContinuation<std::minstd_rand> homotopy {
+					std::minstd_rand {seeds[GetParam()]}, 500, 100
+					};
+};
 
-	// Run the solver several times to ensure the algorithm converges from any
-	// random a/b selected to build G.
-	for(std::size_t i = 0; i < 10; i++) {
-		chemical_system.solveEquilibrium(homotopy);
-		checkEquilibrium();
-		checkPh(default_ph);
-	}
+TEST_P(AgClChemicalSystemHomotopyTest, solve_equilibrium_homotopy) {
+	chemical_system.solveEquilibrium(homotopy);
+	checkEquilibrium();
+	checkPh(default_ph);
 }
+INSTANTIATE_TEST_SUITE_P(HomotopyTest, AgClChemicalSystemHomotopyTest, Range(0, 10));
 
 TEST_F(AgClChemicalSystemTest, solve_equilibrium_brutal_ph_homotopy) {
 	solver::HomotopyContinuation<std::minstd_rand> homotopy(
-			std::minstd_rand {}, 1000, 100
+			std::minstd_rand {}, 500, 100
 			);
 
-	chemical_system.setMaxIteration(1000);
 	chemical_system.fixPH(10);
 	chemical_system.solveEquilibrium(homotopy);
 	checkEquilibrium();

tests/chemmisol/tests/chemical/system/nacl_test.cpp

@@ -368,19 +368,21 @@ TEST_F(NaClChemicalSystemTest, homotopy_G0) {
 	}
 }
 
-TEST_F(NaClChemicalSystemTest, solve_equilibrium_homotopy) {
-	solver::HomotopyContinuation<std::minstd_rand> homotopy(
-			std::minstd_rand {}, 1000, 100
-			);
+class NaClChemicalSystemHomotopyTest :
+	public NaClChemicalSystemTest, public WithParamInterface<int> {
+		protected:
+			// seeds are defined in utils.h
+			solver::HomotopyContinuation<std::minstd_rand> homotopy {
+					std::minstd_rand {seeds[GetParam()]}, 100, 100
+					};
+};
 
-	// Run the solver several times to ensure the algorithm converges from any
-	// random a/b selected to build G.
-	for(std::size_t i = 0; i < 10; i++) {
-		chemical_system.solveEquilibrium(homotopy);
-		checkEquilibrium();
-		checkPh(default_ph);
-	}
+TEST_P(NaClChemicalSystemHomotopyTest, solve_equilibrium_homotopy) {
+	chemical_system.solveEquilibrium(homotopy);
+	checkEquilibrium();
+	checkPh(default_ph);
 }
+INSTANTIATE_TEST_SUITE_P(HomotopyTest, NaClChemicalSystemHomotopyTest, Range(0, 10));
 
 TEST_F(NaClChemicalSystemTest, solve_equilibrium_brutal_ph_homotopy) {
 	solver::HomotopyContinuation<std::minstd_rand> homotopy(

tests/chemmisol/tests/chemical/system/po4_test.cpp

@@ -286,15 +286,14 @@ TEST_F(PO4ChemicalSystemTest, solve_equilibrium_absolute_newton) {
 class PO4ChemicalSystemHomotopyTest :
 	public PO4ChemicalSystemTest, public WithParamInterface<int> {
 		protected:
+			// seeds are defined in utils.h
 			solver::HomotopyContinuation<std::minstd_rand> homotopy {
 					std::minstd_rand {seeds[GetParam()]}, 100, 100
 					};
 };
 
 
 TEST_P(PO4ChemicalSystemHomotopyTest, solve_equilibrium_homotopy) {
-	// Run the solver several times to ensure the algorithm converges from any
-	// random a/b selected to build G.
 	chemical_system.solveEquilibrium(homotopy);
 	checkEquilibrium();
 	checkPh(default_ph);
@@ -303,7 +302,7 @@ INSTANTIATE_TEST_SUITE_P(HomotopyTest, PO4ChemicalSystemHomotopyTest, Range(0, 1
 
 TEST_F(PO4ChemicalSystemTest, solve_equilibrium_brutal_ph_homotopy) {
 	solver::HomotopyContinuation<std::minstd_rand> homotopy(
-			std::minstd_rand {}, 1000, 100
+			std::minstd_rand {}, 100, 100
 			);
 
 	chemical_system.fixPH(10);

tests/chemmisol/tests/chemical/system/soh_test.cpp

@@ -288,11 +288,18 @@ TEST_F(SohTest, solve_equilibrium_absolute_newton) {
 	checkEquilibrium();
 }
 
-TEST_F(SohTest, solve_equilibrium_homotopy) {
-	solver::HomotopyContinuation<std::minstd_rand> homotopy(
-			std::minstd_rand {}, 10, 10
-			);
+class SohHomotopyTest :
+	public SohTest, public WithParamInterface<int> {
+		protected:
+			// seeds are defined in utils.h
+			solver::HomotopyContinuation<std::minstd_rand> homotopy {
+					std::minstd_rand {seeds[GetParam()]}, 100, 100
+					};
+};
+
+TEST_P(SohHomotopyTest, solve_equilibrium_homotopy) {
 	chemical_system.solveEquilibrium(homotopy);
 
 	checkEquilibrium();
 }
+INSTANTIATE_TEST_SUITE_P(HomotopyTest, SohHomotopyTest, Range(0, 10));