Merge remote-tracking branch 'coupled/master'

phreeqcpp/GasPhase.cxx

@@ -169,9 +169,9 @@ cxxGasPhase::cxxGasPhase(std::map < int, cxxGasPhase > &entity_map,
 				this->pr_in = entity_ptr->pr_in;
 				this->temperature = entity_ptr->temperature;
 				first = false;
-				}
+			}
-				else
+			else
-				{
+			{
 				if (this->type != entity_ptr->type)
 				{
 					std::ostringstream oss;
@@ -180,9 +180,13 @@ cxxGasPhase::cxxGasPhase(std::map < int, cxxGasPhase > &entity_map,
 					return;
 				}
 
-				this->total_p += entity_ptr->total_p * it->second;
+				this->total_moles += entity_ptr->total_moles * it->second;
 				this->volume += entity_ptr->volume * it->second;
-				this->v_m += entity_ptr->v_m * it->second;
+				if (sum_fractions > 0.0)
+				{
+					this->v_m += entity_ptr->v_m * it->second / sum_fractions;
+					this->total_p += entity_ptr->total_p  * it->second / sum_fractions;
+				}
 			}
 		}
 		cxxGasPhase *gas_phase_ptr = Utilities::Rxn_find(entity_map, it->first);

phreeqcpp/Makefile.old

@@ -971,8 +971,8 @@ tester:
 #	cd ../mytest; make clean; make -k -j 1 $(SPOOL) make.out $(SPOOL2); make diff $(SPOOL) diff.out $(SPOOL2)
 	cd ../mytest $(CONCAT) make clean $(CONCAT) make -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make diff $(SPOOL) diff.out $(SPOOL2)
 	cd ../examples $(CONCAT) make -f Makefile.old clean $(CONCAT) make -f Makefile.old -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make -f Makefile.old diff $(SPOOL) diff.out $(SPOOL2)
-	svn status -q ../mytest 
+#	svn status -q ../mytest 
-	svn status -q ../examples
+#	svn status -q ../examples
 
 #ld-option
 # Usage: ldflags += $(call ld-option, -Wl$(comma)--hash-style=sysv)

phreeqcpp/Phreeqc.cpp

@@ -711,6 +711,8 @@ void Phreeqc::init(void)
 	all_cells                = 0;
 	multi_Dflag              = FALSE;
 	interlayer_Dflag         = FALSE;
+	implicit                 = FALSE;
+	max_mixf                 = 1.0;
 	default_Dw               = 0;
 	correct_Dw               = 0;
 	multi_Dpor               = 0;
@@ -1687,6 +1689,8 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc)
 	max_cells = pSrc->max_cells;
 	multi_Dflag              = pSrc->multi_Dflag;
 	interlayer_Dflag         = pSrc->interlayer_Dflag;
+	implicit                 = pSrc->implicit;
+	max_mixf                 = pSrc->max_mixf;
 	default_Dw               = pSrc->default_Dw;
 	correct_Dw               = pSrc->correct_Dw;
 	multi_Dpor               = pSrc->multi_Dpor;
@@ -2372,6 +2376,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc)
 		{
 			struct rate *rate_new = rate_copy(it->second.Get_rate());
 			it->second.Set_rate(rate_new);
+			it->second.Set_PhreeqcPtr(this);
 		}
 	}
 	

phreeqcpp/Phreeqc.h

@@ -1065,10 +1065,11 @@ public:
 	LDBLE viscosity(void);
 	LDBLE calc_vm_Cl(void);
 	int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
-	int find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant);
+	LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant);
+	void diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant);
 	int fill_spec(int cell_no);
 	void define_ct_structures(void);
-	int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec);
+	int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec, int i, int stagnant);
 	static int sort_species_name(const void *ptr1, const void *ptr2);
 	int disp_surf(LDBLE stagkin_time);
 	int diff_stag_surf(int mobile_cell);
@@ -1426,6 +1427,8 @@ protected:
 	int old_cells, max_cells, all_cells;
 	int multi_Dflag;		/* signals calc'n of multicomponent diffusion */
 	int interlayer_Dflag;	/* multicomponent diffusion and diffusion through interlayer porosity */
+	int implicit;	    /* implicit calculation of diffusion */
+	LDBLE max_mixf;     /* the maximum value of the implicit mixfactor = De * Dt / (Dx^2) */
 	LDBLE default_Dw;		/* default species diffusion coefficient in water at 25oC, m2/s */
 	int correct_Dw;         /* if true, Dw is adapted in calc_SC */
 	LDBLE multi_Dpor;		/* uniform porosity of free porewater in solid medium */

phreeqcpp/basicsubs.cpp

@@ -1470,7 +1470,8 @@ kinetics_moles(const char *kinetics_name)
 
 	error_string = sformatf( "No data for rate %s in KINETICS keyword.",
 			kinetics_name);
-	warning_msg(error_string);
+	//if (count_warnings >= 0) // appt debug cvode
+	//	warning_msg(error_string);
 	return (0);
 }
 /* ---------------------------------------------------------------------- */

phreeqcpp/class_main.cpp

@@ -277,8 +277,8 @@ write_banner(void)
 
 	/* version */
 #ifdef NPP
-	//len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8 AmpŠre");
+	//len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1 AmpŠre");
-	len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8");
+	len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1");
 #else
 	len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
 #endif
@@ -302,7 +302,7 @@ write_banner(void)
 
 	/* date */
 #ifdef NPP
-	len = sprintf(buffer, "%s", "January 17, 2019");
+	len = sprintf(buffer, "%s", "May 29, 2019");
 #else
 	len = sprintf(buffer, "%s", "@VER_DATE@");
 #endif
@@ -492,12 +492,12 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
 		screen_msg(sformatf("Database file: %s\n\n", token));
 		strcpy(db_file, token);
 #ifdef NPP
-		//output_msg(sformatf("Using PHREEQC: version 3.4.8 Ampère, compiled January 17, 2019\n"));
+		//output_msg(sformatf("Using PHREEQC: version 3.5.1 Ampère, compiled May 29, 2019\n"));
 #endif
 		output_msg(sformatf("   Input file: %s\n", in_file));
 		output_msg(sformatf("  Output file: %s\n", out_file));
 #ifdef NPP
-		output_msg(sformatf("Using PHREEQC: version 3.4.8, compiled January 17, 2019\n"));
+		output_msg(sformatf("Using PHREEQC: version 3.5.1, compiled May 29, 2019\n"));
 #endif
 		output_msg(sformatf("Database file: %s\n\n", token));
 #ifdef NPP

phreeqcpp/common/PHRQ_io.cpp

@@ -32,7 +32,7 @@ PHRQ_io(void)
 	punch_on = true;
 	error_on = true;
 	dump_on = true;
-	echo_on = true; //**appt
+	echo_on = true;
 	screen_on = true;
 	echo_destination = ECHO_OUTPUT;
 
@@ -800,7 +800,7 @@ get_line(void)
 					this->push_istream(next_stream);
 					std::ostringstream errstr;
 					errstr << "\n\tReading data from " << file_name <<" ...\n";
-					output_msg(errstr.str().c_str()); // **appt
+					output_msg(errstr.str().c_str());
 				}
 				continue;
 			}

phreeqcpp/global_structures.h

@@ -1052,12 +1052,12 @@ struct sol_D
 struct J_ij
 {
 	const char *name;
-	LDBLE tot1, tot2;        /* species change in cells i and j */
+	LDBLE tot1, tot2, charge;  /* species change in cells i and j */
 };
 struct M_S
 {
 	const char *name;
-	LDBLE tot1, tot2;
+	LDBLE tot1, tot2, charge; /* master species transport in cells i and j */
 };
 // Pitzer definitions
 typedef enum

phreeqcpp/kinetics.cpp

@@ -1312,11 +1312,11 @@ set_and_run_wrapper(int i, int use_mix, int use_kinetics, int nsaver,
 	{
 		diagonal_scale = TRUE;
 		always_full_pitzer = FALSE;
-		max_try = 15;
+		max_try = 14;
 	}
 	else
 	{
-		max_try = 15;
+		max_try = 14;
 	}
 	max_try = (max_tries < max_try) ? max_tries : max_try;
 	/*max_try = 1; */
@@ -1511,7 +1511,6 @@ restart:
 		}
 		if (j == 14)
 		{
-			//continue;
 			cxxStorageBin error_bin;
 			Use2cxxStorageBin(error_bin);
 			std::ostringstream error_input;
@@ -2479,7 +2478,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
  * Rates and moles of each reaction are calculated in calc_kinetic_reaction
  * Total number of moles in reaction is stored in kinetics[i].totals
  */
-
+	//int increase_tol = 0; // appt
 	int converge, m_iter;
 	int pr_all_save;
 	int nsaver;
@@ -2530,9 +2529,10 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
 		converge =
 			set_and_run_wrapper(i, use_mix, FALSE, nsaver, step_fraction);
 		if (converge == MASS_BALANCE)
-			error_msg
+		{
-				("Negative concentration in system. Stopping calculation.",
+			error_string = sformatf("Negative concentration in solution %d. Stopping calculation.", cell_no);
-				 STOP);
+			error_msg(error_string, STOP);
+		}
 		run_reactions_iterations += iterations;
 	}
 	else
@@ -2604,9 +2604,10 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
 			else
 				converge = set_and_run_wrapper(i, use_mix, FALSE, i, step_fraction);
 			if (converge == MASS_BALANCE)
-				error_msg
+			{
-					("Negative concentration in system. Stopping calculation.",
+				error_string = sformatf("Negative concentration in solution %d. Stopping calculation.", cell_no);
-					 STOP);
+				error_msg(error_string, STOP);
+			}
 			saver();
 			pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, i);
 			ss_assemblage_ptr = Utilities::Rxn_find(Rxn_ss_assemblage_map, i);
@@ -2676,6 +2677,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
 			   cvode_mem = CVodeMalloc(n_reactions, f, 0.0, y, ADAMS, FUNCTIONAL, SV, &reltol, abstol, NULL, NULL, FALSE, iopt, ropt, machEnv);
 			   iopt[MXSTEP] is maximum number of steps that CVODE tries.
 			 */
+			//iopt[SLDET] = TRUE; // appt
 			iopt[MXSTEP] = kinetics_ptr->Get_cvode_steps();
 			iopt[MAXORD] = kinetics_ptr->Get_cvode_order();
 			kinetics_cvode_mem =
@@ -2705,30 +2707,49 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
 			 */
 			m_iter = 0;
 			sum_t = 0;
-		  RESTART:
+		RESTART:
 			while (flag != SUCCESS)
 			{
 				sum_t += cvode_last_good_time;
-				if (state != TRANSPORT)
 				{
-					error_string = sformatf(
+					error_string = sformatf("CV_ODE: Time: %8.2e s. Delta t: %8.2e s. Calls: %d.", (double)(sum_t), (double) cvode_last_good_time, m_iter);
-						"CVode incomplete at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d, continuing...\n",
+					status(0, error_string, true);
-						(int)iopt[NST], cell_no, (double)sum_t, m_iter + 1);
-					warning_msg(error_string);
 				}
+				//if (state != TRANSPORT)
+				//{
+				//	error_string = sformatf(
+				//		"CVode incomplete at cvode_steps %d. Cell: %d. Time: %8.2e s. Cvode calls: %d, continuing...\n",
+				//		(int)iopt[NST], cell_no, (double)sum_t, m_iter + 1);
+				//	warning_msg(error_string);
+				//}
 #ifdef DEBUG_KINETICS
 				if (m_iter > 5)
 					dump_kinetics_stderr(cell_no);
 #endif
 
+				//if (m_iter > 0.5 * kinetics_ptr->Get_bad_step_max() &&
+				//	(cvode_last_good_time < 1e-6 || cvode_last_good_time < 1e-6 * tout)) // appt
+				//{
+				//	if (increase_tol < 3)
+				//	{
+				//		increase_tol += 1;
+				//		for (size_t j = 0; j < kinetics_ptr->Get_kinetics_comps().size(); j++)
+				//		{
+				//			cxxKineticsComp * kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]);
+				//			LDBLE tr = kinetics_comp_ptr->Get_tol() * 10.0;
+				//			kinetics_comp_ptr->Set_tol(tr);
+				//			tr += 0;
+				//		}
+				//	}
+				//}
 				cvode_last_good_time = 0;
 				if (++m_iter >= kinetics_ptr->Get_bad_step_max())
 				{
 					m_temp = (LDBLE *) free_check_null(m_temp);
 					m_original = (LDBLE *) free_check_null(m_original);
 					error_string = sformatf(
-						"Restart of integration at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d.",
+						"CVode is at maximum calls: %d. Cell: %d. Time: %8.2e s\nERROR: Please increase the maximum calls with -bad_step_max.",
-						(int)iopt[NST], cell_no, (double)sum_t, m_iter - 1);
+						m_iter, cell_no, (double)sum_t);
 					error_msg(error_string, STOP);
 				}
 				tout1 = tout - sum_t;
@@ -2811,6 +2832,11 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
 			free_cvode();
 			use.Set_mix_in(use_save.Get_mix_in());
 			use.Set_mix_ptr(use_save.Get_mix_ptr());
+
+			error_string = sformatf("CV_ODE: Final Delta t: %8.2e s. Calls: %d.             ", (double)cvode_last_good_time, m_iter);
+			status(0, error_string, true);
+
+			//status(0, NULL);
 		}
 
 		rate_sim_time = rate_sim_time_start + kin_time;

phreeqcpp/phqalloc.cpp

@@ -142,6 +142,7 @@ PHRQ_calloc(size_t num, size_t size
 	assert((s_pTail == NULL) || (s_pTail->pNext == NULL));
 
 	p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num);
+	//p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt
 
 	if (p == NULL)
 		return NULL;

phreeqcpp/pitzer.cpp

@@ -2555,7 +2555,6 @@ gammas_pz(bool exch_a_f)
 			case 9:			/* activity water */
 				break;
 			case 4:			/* Exchange */
-
 				/*
 				 *   Find CEC
 				 *   z contains valence of cation for exchange species, alk contains cec
@@ -2580,7 +2579,6 @@ gammas_pz(bool exch_a_f)
 				/*
 				 *   All other species
 				 */
-
 				/* modific 29 july 2005... */
 				if (s_x[i]->equiv != 0 && s_x[i]->alk > 0)
 				{
@@ -2595,7 +2593,6 @@ gammas_pz(bool exch_a_f)
 							continue;
 						coef = s_x[i]->rxn_x->token[j].coef;
 						s_x[i]->lg += coef * s_x[i]->rxn_x->token[j].s->lg;
-						/*s_x[i]->dg += coef * s_x[i]->rxn_x->token[j].s->dg;*//* remove? */
 					}
 				}
 				if (s_x[i]->a_f && s_x[i]->primary == NULL && s_x[i]->moles)

phreeqcpp/readtr.cpp

@@ -94,9 +94,10 @@ read_transport(void)
 		"porosities",			/* 42 */
 		"porosity",				/* 43 */
 		"fix_current",			/* 44 */
-		"current"			    /* 45 */
+		"current",			    /* 45 */
+		"implicit"			    /* 46 */
 	};
-	int count_opt_list = 46;
+	int count_opt_list = 47;
 
 	strcpy(file_name, "phreeqc.dmp");
 	/*
@@ -682,6 +683,31 @@ read_transport(void)
 			}
 			opt_save = OPTION_DEFAULT;
 			break;
+		case 46:				/* implicit diffusion */
+			copy_token(token, &next_char, &l);
+			str_tolower(token);
+			if (strstr(token, "f") == token)
+				implicit = FALSE;
+			else if (strstr(token, "t") == token)
+				implicit = TRUE;
+			else
+			{
+				input_error++;
+				error_msg
+				("Expected flag for implicit diffusion calc`s: 'true' or 'false'.",
+					CONTINUE);
+			}
+			if (copy_token(token, &next_char, &l) == DIGIT)
+			{
+				sscanf(token, SCANFORMAT, &max_mixf);
+			}
+			else
+			{
+				//warning_msg("Expected the maximal value for the mixfactor (= D * Dt / Dx^2) in implicit calc`s of diffusion.");
+				max_mixf = 1.0;
+			}
+			opt_save = OPTION_DEFAULT;
+			break;
 		}
 		if (return_value == EOF || return_value == KEYWORD)
 			break;

phreeqcpp/tidy.cpp

@@ -1435,7 +1435,17 @@ tidy_inverse(void)
 			inverse[i].elts[j].master->in = TRUE;
 		}
 		s_eminus->primary->in = TRUE;	/* Include electrons */
+		if (master_alk_ptr)
+		{
 		master_alk_ptr->in = TRUE;	/* Include alkalinity */
+		}
+		else
+		{
+			input_error++;
+			error_string = sformatf(
+				"Alkalinity must be defined in SOLUTION_MASTER_SPECIES to be able to use INVERSE_MODELING.");
+			error_msg(error_string, CONTINUE);
+		}
 /*
  *   Unmark primary and mark secondary master species for redox elements
  */

phreeqcpp/utilities.cpp

@@ -6,7 +6,6 @@
 #include "Solution.h"
 #include <time.h>
 
-
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
 add_elt_list(struct elt_list *elt_list_ptr, LDBLE coef)
@@ -1516,7 +1515,7 @@ status(int count, const char *str, bool rk_string)
 			}
 			else
 			{
-				screen_string = sformatf("%-15s%-27s%1s%37s", sim_str, state_str, spin_str, stdstr.c_str());
+				screen_string = sformatf("%-15s%-27s%1s%45s", sim_str, state_str, spin_str, stdstr.c_str());
 				status_string = screen_string;
 			}
 		}