Tony's changes for viscosity and heat in transport. ex13_impl failed …

…in Release.

examples/ex12b

@@ -1,4 +1,4 @@
-#DATABASE ../database/phreeqc.dat
+DATABASE ../database/phreeqc.dat
 #PRINT
 #   -reset   false
 
@@ -15,14 +15,16 @@ EXCHANGE 1-31
     KX    0.048
 END
 TRANSPORT
-   -cells 30
+   -cells 10
    -lengths 0.3333333
    -shifts 1
    -flow_direction         diffusion
    -boundary_conditions    constant  closed
    -thermal_diffusion      3.0    # heat is retarded equal to Na
    -diffusion_coefficient  0.3e-9  # m^2/s
-   -time_step              1.0e+10  1 # max_mixf = 2/9 = Dt_max * De / Dx^2. Dt_max = 8.2140e+07 seconds, Number of mixes = 1e10 / 8.214e7 = 122
+# -multi_d true 0.3e-9 1 0.05 1.0 false
+# -implicit true 
+   -time_step              1.0e+9  1 # max_mixf = 2/9 = Dt_max * De / Dx^2. Dt_max = 8.2140e+07 seconds, Number of mixes = 1e10 / 8.214e7 = 122
 USER_GRAPH 1 Example 12b
    -headings Tradit:Na Cl TC Analyt
    # -headings TC Analyt
@@ -57,18 +59,35 @@ USER_GRAPH 1 Example 12b
    -end
 END
 
-Reinitialize the column...
+# Reinitialize the column...
 copy cell 31 1-30
 END
 TRANSPORT
-   -multi_d true 0.3e-9 1 0.05 0 false # will give the traditional results when tc = 25 throughout
-   -thermal_diffusion  3.0 1.33e-9  # efine the diffusion coefficient for heat to be equal to Na
-   -implicit true 3 -12             # max_mixf = 3, min_dif_LM = -12
+-shifts 1
+   -multi_d true 2.33e-9 1 0.05 0 false # will give the traditional results when tc = 25 throughout
+   -thermal_diffusion  3.0 2.33e-9  # define the diffusion coefficient for heat equal to Na
 USER_GRAPH 1 Example 12b
-   -headings MultiD&Visc&Implicit:Na Cl TC
+   -headings MultiD&Visc:Na Cl TC
    -start
   10 x = DIST
   20 PLOT_XY x, TOT("Na")*1000, symbol = Circle, line_width = 0, symbol_size = 5, color = Red
   30 PLOT_XY x, TOT("Cl")*1000, symbol = Circle, line_width = 0, symbol_size = 5, color = Green
-  40 PLOT_XY x, TC,             symbol = Circle, line_width = 0, symbol_size = 8, y-axis 2, color = Blue
+  40 PLOT_XY x, TC,             symbol = Circle, line_width = 0, symbol_size = 5, y-axis 2, color = Blue
+END
+# Reinitialize the column...
+copy cell 31 1-30
+USER_GRAPH 1; -connect_simulations false
+END
+TRANSPORT
+-shifts 1
+   -multi_d true 2.33e-9 1 0.05 0 false # will give the traditional results when tc = 25 throughout
+   -thermal_diffusion  3.0 2.33e-9  # define the diffusion coefficient for heat equal to Na
+   -implicit true 3 -12             # max_mixf = 3, min_dif_LM = -12
+USER_GRAPH 1 Example 12b
+   -headings MultiD&Visc&Implicit:Na Cl TC
+   -start
+  10 x = DIST
+  20 PLOT_XY x, TOT("Na")*1000, symbol = XCross, line_width = 0, symbol_size = 9, color = Red
+  30 PLOT_XY x, TOT("Cl")*1000, symbol = XCross, line_width = 0, symbol_size = 9, color = Green
+  40 PLOT_XY x, TC,             symbol = XCross, line_width = 0, symbol_size = 9, y-axis 2, color = Blue
 END

src/Solution.cxx

@@ -48,6 +48,7 @@ cxxSolution::cxxSolution(PHRQ_io * io)
 	this->cb = 0.0;
 	this->density = 1.0;
 	this->viscosity = 1.0;
+	this->viscos_0 = 1.0;
 	this->mass_water = 1.0;
 	this->soln_vol = 1.0;
 	this->total_alkalinity = 0.0;
@@ -82,6 +83,7 @@ cxxSolution::operator =(const cxxSolution &rhs)
 		this->total_o                    = rhs.total_o;
 		this->density                    = rhs.density;
 		this->viscosity                  = rhs.viscosity;
+		this->viscos_0                   = rhs.viscos_0;
 		this->cb                         = rhs.cb;
 		this->mass_water                 = rhs.mass_water;
 		this->soln_vol                   = rhs.soln_vol;

src/Solution.h

@@ -51,6 +51,8 @@ class cxxSolution:public cxxNumKeyword
 	void Set_density(LDBLE l_density)                   {this->density = l_density;}
 	LDBLE Get_viscosity() const                         { return this->viscosity; }
 	void Set_viscosity(LDBLE l_viscos)                  { this->viscosity = l_viscos; }
+	LDBLE Get_viscos_0() const                          { return this->viscos_0; }
+	void Set_viscos_0(LDBLE l_viscos_0)                 { this->viscos_0 = l_viscos_0; }
 	LDBLE Get_mass_water() const                        {return this->mass_water;}
 	void Set_mass_water(LDBLE l_mass_water)             {this->mass_water = l_mass_water;}
 	LDBLE Get_total_alkalinity() const                  {return this->total_alkalinity;}
@@ -134,6 +136,7 @@ class cxxSolution:public cxxNumKeyword
 	LDBLE mass_water;
 	LDBLE density;
 	LDBLE viscosity;
+	LDBLE viscos_0;
 	LDBLE soln_vol;
 	LDBLE total_alkalinity;
 	cxxNameDouble totals;

src/global_structures.h

@@ -109,7 +109,7 @@
 #define TRANSPORT		 8
 #define PHAST		     9
 
-/* constants in mass balance */
+/* constraints in mass balance */
 #define EITHER 0
 #define DISSOLVE 1
 #define PRECIPITATE -1
@@ -1486,7 +1486,9 @@ class spec
 		c = 0;
 		// charge number
 		z = 0;
-		// temperature corrected free water diffusion coefficient, m2/s
+		// free water diffusion coefficient, m2/s
+		Dw = 0;
+		// temperature and viscosity corrected free water diffusion coefficient, m2/s
 		Dwt = 0;
 		// temperature factor for Dw
 		dw_t = 0;
@@ -1503,6 +1505,7 @@ class spec
 	LDBLE lg;
 	LDBLE c;
 	LDBLE z;
+	LDBLE Dw;
 	LDBLE Dwt;
 	LDBLE dw_t;
 	LDBLE dw_a_v_dif;
@@ -1521,17 +1524,17 @@ class sol_D
 		count_exch_spec = 0;
 		// total moles of X-, max X- in transport step in sol_D[1], tk
 		exch_total = 0, x_max = 0, tk_x = 0;
-		// (tk_x * viscos_0_25) / (298 * viscos_0) 
-		viscos_f0 = 0;
-		// (viscos_0) / (298 * viscos) 
-		viscos_f = 0;
+		// viscos_0 at I = 0
+		viscos_0 = 0;
+		// viscosity of solution
+		viscos = 0;
 		spec = NULL;
 		spec_size = 0;
 	}
 	int count_spec;
 	int count_exch_spec;
 	LDBLE exch_total, x_max, tk_x;
-	LDBLE viscos_f0, viscos_f;
+	LDBLE viscos_0, viscos;
 	class spec* spec;
 	int spec_size;
 };

src/kinetics.cpp

@@ -1581,7 +1581,8 @@ set_and_run(int i, int use_mix, int use_kinetics, int nsaver,
 	sum_species();
 	viscos = viscosity(NULL);
 	use.Get_solution_ptr()->Set_viscosity(viscos);
-	if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL)
+	use.Get_solution_ptr()->Set_viscos_0(viscos_0);
+	if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL && use.Get_surface_ptr()->Get_calc_viscosity())
 		use.Get_surface_ptr()->Set_DDL_viscosity(viscosity(use.Get_surface_ptr()));
 	return (converge);
 }

src/mainsubs.cpp

@@ -438,6 +438,7 @@ initial_solutions(int print)
 			sum_species();
 			viscos = viscosity(NULL);
 			use.Get_solution_ptr()->Set_viscosity(viscos);
+			use.Get_solution_ptr()->Set_viscos_0(viscos_0);
 			if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL)
 				use.Get_surface_ptr()->Set_DDL_viscosity(viscosity(use.Get_surface_ptr()));
 			add_isotopes(solution_ref);
@@ -1267,6 +1268,7 @@ xsolution_save(int n_user)
 	// the subroutine is called at the start of a new simulation, and the following 2 go wrong since s_x is not updated 
 	temp_solution.Set_density(density_x);
 	temp_solution.Set_viscosity(viscos);
+	temp_solution.Set_viscos_0(viscos_0);
 	temp_solution.Set_total_h(total_h_x);
 	temp_solution.Set_total_o(total_o_x);
 	temp_solution.Set_cb(cb_x);	/* cb_x does not include surface charge after sum_species */