Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18350)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18351)
@@ -159,199 +159,224 @@
 
 }/*}}}*/
-#endif
-void solutionsequence_fct(FemModel* femmodel){
-
-	/*intermediary: */
-	Vector<IssmDouble>*  Ml = NULL;
-	Matrix<IssmDouble>*  K  = NULL;
-	Matrix<IssmDouble>*  Mc = NULL;
-	Vector<IssmDouble>*  ug = NULL;
-	Vector<IssmDouble>*  uf = NULL;
-
-	IssmDouble theta,deltat,dmax;
-	int        dof,ncols,ncols2,rstart,rend;
-	int        configuration_type;
-	double     d;
-	int*       cols  = NULL;
-	int*       cols2 = NULL;
-	double*    vals  = NULL;
-	double*    vals2 = NULL;
-
-	/*Create analysis*/
-	MasstransportAnalysis* analysis = new MasstransportAnalysis();
-
-	/*Recover parameters: */
-	femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
-	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-	femmodel->UpdateConstraintsx();
-	theta = 0.5;
-
-	/*Create lumped mass matrix*/
-	analysis->LumpedMassMatrix(&Ml,femmodel);
-	analysis->MassMatrix(&Mc,femmodel);
-	analysis->FctKMatrix(&K,NULL,femmodel);
-
-	/*Convert matrices to PETSc matrices*/
-	Mat D_petsc  = NULL;
-	Mat LHS      = NULL;
-	Vec RHS      = NULL;
-	Vec u        = NULL;
-	Mat K_petsc  = K->pmatrix->matrix;
-	Vec Ml_petsc = Ml->pvector->vector;
-	Mat Mc_petsc = Mc->pmatrix->matrix;
-
-	/*Create D Matrix*/
-	#ifdef _HAVE_PETSC_
-	CreateDMatrix(&D_petsc,K_petsc);
-
-	/*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
-	CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
-
-	/*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
-	GetSolutionFromInputsx(&ug,femmodel);
-	Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
-	delete ug;
-	CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
-	delete uf;
-
-	/*Go solve!*/
-	SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
-	MatFree(&LHS);
-	VecFree(&RHS);
-
-	/*Richardson to calculate udot*/
-	Vec udot = NULL;
+void RichardsonUdot(Vec* pudot,Vec u,Vec Ml,Mat K,Mat Mc){/*{{{*/
+	/*Use Richardson's formulato get udot using 5 steps and an initial guess of 0
+	 *
+	 * udot_new = udot_old + Ml^-1 (K^(n+1) u - Mc udot_old)
+	 *
+	 * */
+
+	/*Intermediaries*/
+	Vec udot  = NULL;
+	Vec temp1 = NULL;
+	Vec temp2 = NULL;
+
+	/*Initialize vectors*/
 	VecDuplicate(u,&udot);
+	VecDuplicate(u,&temp1);
+	VecDuplicate(u,&temp2);
+
+	/*Initial guess*/
 	VecZeroEntries(udot);
-	Vec temp1 = NULL; VecDuplicate(u,&temp1);
-	Vec temp2 = NULL; VecDuplicate(u,&temp2);
+	
+	/*Richardson's iterations*/
 	for(int i=0;i<5;i++){
-		/*udot_new = udot_old + Ml^-1 (K^(n+1) u -- Mc udot_old)*/
-		MatMult(Mc_petsc,udot,temp1);
-		MatMult(K_petsc, u,   temp2);
-		VecAXPBY(temp2,-1.,1.,temp1); // temp2 = (K^(n+1) u -- Mc udot_old)
-		VecPointwiseDivide(temp1,temp2,Ml_petsc); //temp1 = Ml^-1 temp2
+		MatMult(Mc,udot,temp1);
+		MatMult(K, u,   temp2);
+		VecAXPBY(temp2,-1.,1.,temp1);       // temp2 = (K^(n+1) u -- Mc udot_old)
+		VecPointwiseDivide(temp1,temp2,Ml); // temp1 = Ml^-1 temp2
 		VecAXPBY(udot,1.,1.,temp1);
 	}
+
+	/*Clean up and assign output pointer*/
 	VecFree(&temp1);
 	VecFree(&temp2);
-
-	/*Serialize u and udot*/
-	IssmDouble* udot_serial = NULL;
-	IssmDouble* u_serial    = NULL;
-	IssmDouble* ml_serial    = NULL;
-	VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
-	VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
-	VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
-
-	/*Anti diffusive fluxes*/
-	Vec Ri_plus  = NULL;
-	Vec Ri_minus = NULL;
-	double uiLmax = 3.;
-	double uiLmin = 2.;
-	VecDuplicate(u,&Ri_plus);
-	VecDuplicate(u,&Ri_minus);
-	MatGetOwnershipRange(K_petsc,&rstart,&rend);
-	for(int row=rstart; row<rend; row++){
-		double Pi_plus  = 0.;
-		double Pi_minus = 0.;
-		MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-		_assert_(ncols==ncols2);
-		for(int j=0; j<ncols; j++) {
-			_assert_(cols[j]==cols2[j]);
-			d = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
-			if(row!=cols[j]){
-				if(d>0.){
-					Pi_plus  += d;
+	*pudot=udot;
+
+}/*}}}*/
+#endif
+	void solutionsequence_fct(FemModel* femmodel){
+
+		/*intermediary: */
+		Vector<IssmDouble>*  Ml = NULL;
+		Matrix<IssmDouble>*  K  = NULL;
+		Matrix<IssmDouble>*  Mc = NULL;
+		Vector<IssmDouble>*  ug = NULL;
+		Vector<IssmDouble>*  uf = NULL;
+
+		IssmDouble theta,deltat,dmax;
+		int        dof,ncols,ncols2,rstart,rend;
+		int        configuration_type;
+		double     d;
+		int*       cols  = NULL;
+		int*       cols2 = NULL;
+		double*    vals  = NULL;
+		double*    vals2 = NULL;
+
+		/*Create analysis*/
+		MasstransportAnalysis* analysis = new MasstransportAnalysis();
+
+		/*Recover parameters: */
+		femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
+		femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+		femmodel->UpdateConstraintsx();
+		theta = 0.5;
+
+		/*Create lumped mass matrix*/
+		analysis->LumpedMassMatrix(&Ml,femmodel);
+		analysis->MassMatrix(&Mc,femmodel);
+		analysis->FctKMatrix(&K,NULL,femmodel);
+
+		/*Convert matrices to PETSc matrices*/
+		Mat D_petsc  = NULL;
+		Mat LHS      = NULL;
+		Vec RHS      = NULL;
+		Vec u        = NULL;
+		Vec udot     = NULL;
+		Mat K_petsc  = K->pmatrix->matrix;
+		Vec Ml_petsc = Ml->pvector->vector;
+		Mat Mc_petsc = Mc->pmatrix->matrix;
+
+		/*Create D Matrix*/
+		#ifdef _HAVE_PETSC_
+		CreateDMatrix(&D_petsc,K_petsc);
+
+		/*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
+		CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
+
+		/*Get previous solution u^n*/
+		GetSolutionFromInputsx(&ug,femmodel);
+		Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
+		delete ug;
+
+		/*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
+		CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
+		delete uf;
+
+		/*Go solve!*/
+		SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
+		MatFree(&LHS);
+		VecFree(&RHS);
+
+		/*Richardson to calculate udot*/
+		RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
+
+		/*Serialize u and udot*/
+		IssmDouble* udot_serial = NULL;
+		IssmDouble* u_serial    = NULL;
+		IssmDouble* ml_serial    = NULL;
+		VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
+		VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
+		VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
+
+		/*Anti diffusive fluxes*/
+		Vec Ri_plus  = NULL;
+		Vec Ri_minus = NULL;
+		double uiLmax = 3.;
+		double uiLmin = 2.;
+		VecDuplicate(u,&Ri_plus);
+		VecDuplicate(u,&Ri_minus);
+		MatGetOwnershipRange(K_petsc,&rstart,&rend);
+		for(int row=rstart; row<rend; row++){
+			double Pi_plus  = 0.;
+			double Pi_minus = 0.;
+			MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
+			MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+			_assert_(ncols==ncols2);
+			for(int j=0; j<ncols; j++) {
+				_assert_(cols[j]==cols2[j]);
+				d = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
+				if(row!=cols[j]){
+					if(d>0.){
+						Pi_plus  += d;
+					}
+					else{
+						Pi_minus += d;
+					}
+				}
+			}
+
+			/*Compute Qis and Ris*/
+			double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u_serial[row]);
+			double Qi_minus = ml_serial[row]/deltat*(uiLmin - u_serial[row]);
+			d = 1.;
+			if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
+			VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
+			d = 1.;
+			if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
+			VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
+
+			MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
+			MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		}
+		VecAssemblyBegin(Ri_plus);
+		VecAssemblyEnd(  Ri_plus);
+		VecAssemblyBegin(Ri_minus);
+		VecAssemblyEnd(  Ri_minus);
+
+		/*Serialize Ris*/
+		IssmDouble* Ri_plus_serial  = NULL;
+		IssmDouble* Ri_minus_serial = NULL;
+		VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
+		VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
+		VecFree(&Ri_plus);
+		VecFree(&Ri_minus);
+
+		/*Build fbar*/
+		Vec Fbar = NULL;
+		VecDuplicate(u,&Fbar);
+		for(int row=rstart; row<rend; row++){
+			MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
+			MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+			_assert_(ncols==ncols2);
+			d = 0.;
+			for(int j=0; j<ncols; j++) {
+				_assert_(cols[j]==cols2[j]);
+				if(row==cols[j]) continue;
+				double f_ij = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
+				if(f_ij>0){
+					d += min(Ri_plus_serial[row],Ri_minus_serial[cols[j]]) * f_ij;
 				}
 				else{
-					Pi_minus += d;
+					d += min(Ri_minus_serial[row],Ri_plus_serial[cols[j]]) * f_ij;
 				}
 			}
-		}
-
-		/*Compute Qis and Ris*/
-		double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u_serial[row]);
-		double Qi_minus = ml_serial[row]/deltat*(uiLmin - u_serial[row]);
-		d = 1.;
-		if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
-		VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
-		d = 1.;
-		if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
-		VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
-
-		MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-	}
-	VecAssemblyBegin(Ri_plus);
-	VecAssemblyEnd(  Ri_plus);
-	VecAssemblyBegin(Ri_minus);
-	VecAssemblyEnd(  Ri_minus);
-
-	/*Serialize Ris*/
-	IssmDouble* Ri_plus_serial  = NULL;
-	IssmDouble* Ri_minus_serial = NULL;
-	VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
-	VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
-	VecFree(&Ri_plus);
-	VecFree(&Ri_minus);
-
-	/*Build fbar*/
-	Vec Fbar = NULL;
-	VecDuplicate(u,&Fbar);
-	for(int row=rstart; row<rend; row++){
-		MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-		_assert_(ncols==ncols2);
-		d = 0.;
-		for(int j=0; j<ncols; j++) {
-			_assert_(cols[j]==cols2[j]);
-			if(row==cols[j]) continue;
-			double f_ij = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
-			if(f_ij>0){
-				d += min(Ri_plus_serial[row],Ri_minus_serial[cols[j]]) * f_ij;
-			}
-			else{
-				d += min(Ri_minus_serial[row],Ri_plus_serial[cols[j]]) * f_ij;
-			}
-		}
-		VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
-		MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-	}
-	VecAssemblyBegin(Fbar);
-	VecAssemblyEnd(  Fbar);
-
-	MatFree(&D_petsc);
-	delete Mc;
-	xDelete<IssmDouble>(udot_serial);
-	xDelete<IssmDouble>(u_serial);
-	xDelete<IssmDouble>(ml_serial);
-	xDelete<IssmDouble>(Ri_plus_serial);
-	xDelete<IssmDouble>(Ri_minus_serial);
-
-	/*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
-	VecDuplicate(u,&temp1);
-	VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2
-	VecAXPBY(udot,1.,1.,temp1);
-	VecAXPY(u,deltat,temp1);
-	VecFree(&Fbar);
-	VecFree(&udot);
-	VecFree(&temp1);
-
-	uf =new Vector<IssmDouble>(u);
-	VecFree(&u);
-
-	InputUpdateFromSolutionx(femmodel,uf); 
-	delete uf;
-
-	#else
-	_error_("PETSc needs to be installed");
-	#endif
-
-	delete Ml;
-	delete K;
-	delete analysis;
-
-}
+			VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
+			MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
+			MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		}
+		VecAssemblyBegin(Fbar);
+		VecAssemblyEnd(  Fbar);
+
+		MatFree(&D_petsc);
+		delete Mc;
+		xDelete<IssmDouble>(udot_serial);
+		xDelete<IssmDouble>(u_serial);
+		xDelete<IssmDouble>(ml_serial);
+		xDelete<IssmDouble>(Ri_plus_serial);
+		xDelete<IssmDouble>(Ri_minus_serial);
+
+		/*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
+		Vec temp1 = NULL; VecDuplicate(u,&temp1);
+		VecDuplicate(u,&temp1);
+		VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2
+		VecAXPBY(udot,1.,1.,temp1);
+		VecAXPY(u,deltat,temp1);
+		VecFree(&Fbar);
+		VecFree(&udot);
+		VecFree(&temp1);
+
+		uf =new Vector<IssmDouble>(u);
+		VecFree(&u);
+
+		InputUpdateFromSolutionx(femmodel,uf); 
+		delete uf;
+
+		#else
+		_error_("PETSc needs to be installed");
+		#endif
+
+		delete Ml;
+		delete K;
+		delete analysis;
+
+	}