Merge pull request #78 from dlparkhurst/mix_error

Mix error

mytest/CMakeLists.txt

@@ -257,6 +257,7 @@ set(TESTS
   minteqv4
   missing_surf_related_equi
   missing_surf_related_kin
+  mix_unequal
   mixes
   mmb
   mmb2

mytest/mix_unequal

@@ -0,0 +1,37 @@
+database ../database/phreeqc.dat
+SELECTED_OUTPUT 101
+    -file                 mix_unequal_101.sel
+USER_PUNCH 101
+    -headings Mu SC TC
+    -start
+10 PUNCH STR_F$(MU, 20, 12)
+20 PUNCH STR_F$(SC, 20, 10)
+30 PUNCH STR_F$(TC, 20, 10)
+    -end
+SOLUTION 1
+-temp 10
+Na 1
+Cl 1
+-water 0.4
+END
+SOLUTION 2
+-temp 20
+K  10
+Br 10
+-water 0.6
+END
+MIX 1
+# T = (0.4*0.4*10 + 0.6*0.6*20)/(0.16+0.36) = 16.92307692, mass = 0.52
+# Cl = (0.4*0.4*1)/0.52 = 0.3077 mmol/kgw
+# K = (0.6*0.6*10)/0.52 = 6.92 mmol/kgw
+1 0.4
+2 0.6
+END
+
+# T = (0.4*0.2*10 + 0.6*0.7*20)/(0.08+0.34) = 18.4, mass = (0.4*0.2+0.6*0.7) = 0.5
+# Cl = (0.4*0.2*1)/0.5 = 0.3077 mmol/kgw = 0.16
+# K = (0.6*0.7*10)/0.5 = 6.92 mmol/kgw = 8.4
+MIX
+1 0.2
+2 0.7
+END

mytest/mix_unequal_101.sel

@@ -0,0 +1,5 @@
+          Mu	          SC	          TC	
+      0.001000066900	       86.7123676711	       10.0000000000	
+      0.010000092689	     1291.7967478544	       20.0000000000	
+      0.007230851816	      873.0591329232	       16.9230769231	
+      0.008560087209	     1068.6105947361	       18.4000000000	

src/Solution.cxx

@@ -1409,7 +1409,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
 	this->ph = f1 * this->ph + f2 * addee.ph;
 	this->pe = f1 * this->pe + f2 * addee.pe;
 	this->mu = f1 * this->mu + f2 * addee.mu;
-	this->ah2o = f1 * this->mu + f2 * addee.ah2o;
+	this->ah2o = f1 * this->ah2o + f2 * addee.ah2o;
 	this->total_h += addee.total_h * extensive;
 	this->total_o += addee.total_o * extensive;
 	this->cb += addee.cb * extensive;

src/step.cpp

@@ -620,6 +620,7 @@ add_mix(cxxMix *mix_ptr)
  *   and other variables.
  */
 	LDBLE sum_fractions, intensive, extensive;
+	LDBLE sum_fractions_water=0, sum_positive_water=0, intensive_water=0, extensive_water=0;
 	cxxSolution *solution_ptr;
 	int count_positive;
 	LDBLE sum_positive;
@@ -634,10 +635,21 @@ add_mix(cxxMix *mix_ptr)
 	std::map<int, LDBLE>::const_iterator it;
 	for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
 	{
+		solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
+		if (solution_ptr == NULL)
+		{
+			error_string = sformatf("Mix solution not found, %d.",
+				it->first);
+			error_msg(error_string, CONTINUE);
+			input_error++;
+			continue;
+		}
 		sum_fractions += it->second;
+		sum_fractions_water += it->second * solution_ptr->Get_mass_water();
 		if (it->second > 0)
 		{
 			sum_positive += it->second;
+			sum_positive_water += it->second * solution_ptr->Get_mass_water();
 			count_positive++;
 		}
 	}
@@ -653,8 +665,72 @@ add_mix(cxxMix *mix_ptr)
 			continue;
 		}
 		extensive = it->second;
+		extensive_water = it->second * solution_ptr->Get_mass_water();
 		intensive = extensive / sum_fractions;
+		intensive_water = extensive_water / sum_fractions_water;
+
 		if (count_positive < (int) mix_ptr->Get_mixComps().size())
+		{
+			if (it->second > 0)
+			{
+				intensive = extensive / sum_positive;
+				intensive_water = extensive_water / sum_positive;
+			}
+			else
+			{
+				intensive = 0;
+			}
+		}
+		add_solution(solution_ptr, extensive, intensive_water);
+	}
+	return (OK);
+}
+#ifdef SKIP_ERROR
+/* ---------------------------------------------------------------------- */
+int Phreeqc::
+add_mix(cxxMix* mix_ptr)
+/* ---------------------------------------------------------------------- */
+{
+	/*
+	 *   calls add_solution to accumulate all data in master->totals
+	 *   and other variables.
+	 */
+	LDBLE sum_fractions, intensive, extensive;
+	cxxSolution* solution_ptr;
+	int count_positive;
+	LDBLE sum_positive;
+
+	if (mix_ptr == NULL)
+		return (OK);
+	if (mix_ptr->Get_mixComps().size() == 0)
+		return (OK);
+	sum_fractions = 0.0;
+	sum_positive = 0.0;
+	count_positive = 0;
+	std::map<int, LDBLE>::const_iterator it;
+	for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
+	{
+		sum_fractions += it->second;
+		if (it->second > 0)
+		{
+			sum_positive += it->second;
+			count_positive++;
+		}
+	}
+	for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
+	{
+		solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
+		if (solution_ptr == NULL)
+		{
+			error_string = sformatf("Mix solution not found, %d.",
+				it->first);
+			error_msg(error_string, CONTINUE);
+			input_error++;
+			continue;
+		}
+		extensive = it->second;
+		intensive = extensive / sum_fractions;
+		if (count_positive < (int)mix_ptr->Get_mixComps().size())
 		{
 			if (it->second > 0)
 			{
@@ -669,6 +745,7 @@ add_mix(cxxMix *mix_ptr)
 	}
 	return (OK);
 }
+#endif
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
 add_pp_assemblage(cxxPPassemblage *pp_assemblage_ptr)

src/transport.cpp

@@ -1653,14 +1653,17 @@ init_heat_mix(int l_nmix)
 	{
 		if (implicit)
 		{
-			LDBLE viscos_f0;
+			LDBLE viscos_f;
 			l_heat_nmix = l_nmix;
 			for (i = 1; i <= count_cells + 1; i++)
 			{
 				heat_mix_array[i - 1] = heat_mix_array[i] / l_heat_nmix;	/* for implicit, m[i] has mixf with higher cell */
-				viscos_f0 = sol_D[i - 1].viscos_f0 * exp(heat_diffc / sol_D[i - 1].tk_x - heat_diffc / 298.15);
-				viscos_f0 += sol_D[i].viscos_f0 * exp(heat_diffc / sol_D[i].tk_x - heat_diffc / 298.15);
-				heat_mix_array[i - 1] *= (viscos_f0 / 2);
+				if (print_viscosity)
+				{
+					viscos_f = sol_D[i - 1].viscos_f * exp(heat_diffc / sol_D[i - 1].tk_x - heat_diffc / 298.15);
+					viscos_f += sol_D[i].viscos_f * exp(heat_diffc / sol_D[i].tk_x - heat_diffc / 298.15);
+					heat_mix_array[i - 1] *= (viscos_f / 2);
+				}
 			}
 		}
 		else
@@ -4352,9 +4355,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 					b_j *= sol_D[icell].spec[i].Dwt;
 				else
 				{
-					dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f0;
+					dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
 					dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x);
-					dum2 *= sol_D[jcell].viscos_f0;
+					dum2 *= sol_D[jcell].viscos_f;
 					b_j *= dum2;
 				}
 				if (sol_D[icell].spec[i].dw_a_v_dif)
@@ -4463,9 +4466,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 					b_i *= sol_D[jcell].spec[j].Dwt;
 				else
 				{
-					dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f0;
+					dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
 					dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x);
-					dum2 *= sol_D[icell].viscos_f0;
+					dum2 *= sol_D[icell].viscos_f;
 					b_i *= dum2;
 				}
 				if (sol_D[icell].spec[i].dw_a_v_dif)