Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp	(revision 23908)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp	(revision 23909)
@@ -9,30 +9,73 @@
 #include "../modules/modules.h"
 #include "../analyses/analyses.h"
-
+#include <iostream>
+#include <fstream>
 
 #ifdef _HAVE_PETSC_
 
 
-void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff,Vec pf, Vec uf0,IS isv,IS isp,Parameters* parameters){/*{{{*/
+void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff, Vec pf, Vec uf0,IS isv,IS isp,Parameters* parameters){/*{{{*/
 
 	Mat                  A, B, BT;				/* Saddle point block matrices */
-	Mat						IP;						/* Preconditioner matrix */
-	Mat						IP2;
+	Mat						IP;						/* Preconditioner or mass matrix */
 	int                  nu, np;					/* No of. free nodes in velocity / pressure space */
    Vec                  p,uold,unew;			/* Solution vectors for pressure / vel. */ 
-	Vec						tmpu, tmpp, rhsu,rhsp; /* temp. vectors, arbitrary RHS in vel. / pressure space */
-	Vec						gold,gnew,wold,wnew,chi,thetaold,thetanew,eta; /* CG intermediaries */
+	Vec						tmpu,tmpu2,resu,resp,tmpp,tmpp2,rhsu,rhsp; /* temp. vectors, arbitrary RHS in vel. / pressure space */
+	Vec						gold,gnew,wold,wnew,chi; /* CG intermediaries */
 	Vec						f1,f2;					/* RHS of the global system */
-	double					rho,gamma,tmpScalar; /* Step sizes, arbitrary double */
-	KSP						kspu,kspp;				/* KSP contexts for vel. / pressure systems*/
+	double					rho,gamma,tmpScalar,tmpScalar2; /* Step sizes, arbitrary double */
+	KSP						kspu,kspip;		/* KSP contexts for vel. / pressure systems*/
 	KSPConvergedReason	reason;					/* Convergence reason for troubleshooting */
 	int						its;						/* No. of iterations for troubleshooting */
-	double					initRnorm, rnorm, TOL; /* residual norms, STOP tolerance */
+	double					initRnorm, rnorm, TOL,ELLTOL; /* residual norms, STOP tolerance */
 	PC							pcu,pcp;					/* Preconditioner contexts pertaining the KSP contexts*/
 	PetscViewer				viewer;					/* Viewer for troubleshooting */
 	IssmPDouble				t1,t2;					/* Time measurement for bottleneck analysis */
 
-	/*STOP tolerance for the rel. residual*/
-	TOL = 0.4;
+	
+	double tmp1,tmp2,tmp3;
+	int tmpi;
+	double tmp4,tmp5,tmp6,tmp7;
+
+	int noIt;
+
+
+
+	int precond = 0;
+
+	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
+	PetscTruth flag,flg;
+	#else
+	PetscBool flag,flg;
+	#endif
+
+
+	char ksp_type[50];
+	char pc_type[50];
+	int maxiter;
+
+	#if _PETSC_MINOR_<7
+	PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
+	PetscOptionsGetString(PETSC_NULL,"-pc_type",pc_type,49,&flg);
+	PetscOptionsGetReal(PETSC_NULL,"-tol",&TOL,NULL);
+	PetscOptionsGetReal(PETSC_NULL,"-elltol",&ELLTOL,NULL);
+	PetscOptionsGetInt(PETSC_NULL,"-schur_pc",&precond,NULL);
+	PetscOptionsGetInt(PETSC_NULL,"-max_iter",&maxiter,NULL);
+	#else
+	PetscOptionsGetString(NULL,PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
+	PetscOptionsGetString(NULL,PETSC_NULL,"-pc_type",pc_type,49,&flg);
+	PetscOptionsGetReal(NULL,PETSC_NULL,"-tol",&TOL,NULL);
+	PetscOptionsGetReal(NULL,PETSC_NULL,"-elltol",&ELLTOL,NULL);
+	PetscOptionsGetInt(NULL,PETSC_NULL,"-schur_pc",&precond,NULL);
+	PetscOptionsGetInt(NULL,PETSC_NULL,"-max_iter",&maxiter,NULL);
+	#endif
+
+	
+	if(precond){
+		_printf0_("Running WITH preconditioner\n");
+	}else{
+		_printf0_("Running WITHOUT preconditioner\n");
+	}
+
 
 	/*Initialize output*/
@@ -41,22 +84,15 @@
 	/* Extract block matrices from the saddle point matrix */
 	/* [ A   B ] = Kff
-    * [ B^T 0 ] 
+    * [ B^T I ]
+	 * where A is the elliptic submatrix, B^T represents the incompressibility, 
+	 * and I the Schur preconditioner (stored here, because the space was allocated either way) 
 	 *         */
-	#if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=8)
-	MatCreateSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A);
-	MatCreateSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B);
-	MatCreateSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT);
-	#else
 	MatGetSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A);
 	MatGetSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B);
 	MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT);
-	#endif
 	
 	/* Extract preconditioner matrix on the pressure space*/
-	#if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=8)
-	MatCreateSubMatrix(Kff,isp,isp,MAT_INITIAL_MATRIX,&IP);
-	#else
 	MatGetSubMatrix(Kff,isp,isp,MAT_INITIAL_MATRIX,&IP);
-	#endif
+
 
 	/* Get number of velocity / pressure nodes */
@@ -77,13 +113,14 @@
 	VecDuplicate(p,&f2);VecSet(f2,0.0);
 	VecDuplicate(uold,&tmpu);VecSet(tmpu,0.0);
+	VecDuplicate(uold,&tmpu2);VecSet(tmpu2,0.0);
+	VecDuplicate(uold,&resu);VecSet(resu,0.0);
 	VecDuplicate(p,&tmpp);VecSet(tmpp,0.0);
+	VecDuplicate(p,&tmpp2);VecSet(tmpp2,0.0);
 	VecDuplicate(p,&rhsp);VecSet(rhsp,0.0);
+	VecDuplicate(p,&resp);VecSet(resp,0.0);
 	VecDuplicate(uold,&rhsu);VecSet(rhsu,0.0);
 	VecDuplicate(p,&gold);VecSet(gold,0.0);
 	VecDuplicate(p,&wnew);VecSet(wnew,0.0);
 	VecDuplicate(uold,&chi);VecSet(chi,0.0);
-	VecDuplicate(p,&thetanew);VecSet(thetanew,0.0);
-	VecDuplicate(p,&thetaold);VecSet(thetaold,0.0);
-	VecDuplicate(p,&eta);VecSet(eta,0.0);
 	
 	/* Get global RHS (for each block sub-problem respectively)*/
@@ -104,10 +141,61 @@
 	KSPSetOperators(kspu,A,A,DIFFERENT_NONZERO_PATTERN);
 	#endif
-	KSPSetType(kspu,KSPCG);
+	if (strcmp(ksp_type,"gmres")==0){
+		KSPSetType(kspu,KSPGMRES);
+	}else if(strcmp(ksp_type,"pipegmres")==0){
+		KSPSetType(kspu,KSPPGMRES);
+	}else if(strcmp(ksp_type,"cg")==0){
+		KSPSetType(kspu,KSPCG);
+	}else if(strcmp(ksp_type,"pipecg")==0){
+		KSPSetType(kspu,KSPPIPECG);
+	}else if(strcmp(ksp_type,"bicg")==0){
+		KSPSetType(kspu,KSPBICG);
+	}else if(strcmp(ksp_type,"bicgstab")==0){
+		KSPSetType(kspu,KSPBCGS);
+	}else if(strcmp(ksp_type,"ibicgstab")==0){
+		KSPSetType(kspu,KSPIBCGS);
+	}else if(strcmp(ksp_type,"minres")==0){
+		KSPSetType(kspu,KSPMINRES);
+	}else if(strcmp(ksp_type,"cr")==0){
+		KSPSetType(kspu,KSPCR);
+	}else if(strcmp(ksp_type,"pipecr")==0){
+		KSPSetType(kspu,KSPPIPECR);
+	}else{
+		_error_("Suggested KSP method not implemented yet!\n");
+	}
+	
 	KSPSetInitialGuessNonzero(kspu,PETSC_TRUE);
-	//KSPSetTolerances(kspu,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
-	//KSPMonitorSet(kspu,KSPMonitorDefault,NULL,NULL);
+	
+	/*Strong rel. residual tolerance needed when solving for the velocity update. 
+	 * This is because ISSM uses the dimensional equations, so the initial guess chi = 0
+	 * results in a very high initial residual (~ 1e19)*/
+	KSPSetTolerances(kspu,ELLTOL,PETSC_DEFAULT,PETSC_DEFAULT,maxiter); //maxiter
+
+
+
 	KSPGetPC(kspu,&pcu);
-	PCSetType(pcu,PCSOR);
+	if (strcmp(pc_type,"bjacobi")==0){
+		PCSetType(pcu,PCBJACOBI);
+	}else if(strcmp(pc_type,"asm")==0){
+		PCSetType(pcu,PCASM);
+	}else if(strcmp(pc_type,"gamg")==0){
+		PCSetType(pcu,PCGAMG);
+	}else if(strcmp(pc_type,"gasm")==0){
+		PCSetType(pcu,PCGASM);
+	}else if(strcmp(pc_type,"jacobi")==0){
+		PCSetType(pcu,PCJACOBI);
+	}else if(strcmp(pc_type,"icc")==0){
+		PCSetType(pcu,PCICC);
+	}else if(strcmp(pc_type,"ilu")==0){
+		PCSetType(pcu,PCILU);
+	}else if(strcmp(pc_type,"sor")==0){
+		PCSetType(pcu,PCSOR);
+	}else if(strcmp(pc_type,"eisenstat")==0){
+		PCSetType(pcu,PCEISENSTAT);
+	}else if(strcmp(pc_type,"none")==0){
+		PCSetType(pcu,PCNONE);
+	}else{
+	_error_("Suggested preconditioner not implemented yet!\n");
+	}
 	KSPSetUp(kspu);
 
@@ -117,7 +205,39 @@
 	VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.);
 
+	if (VerboseConvergence())
+	{
+
+		VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.);
+		VecNorm(tmpu,NORM_2,&tmp4);
+
+		VecScale(f2,-1.);MatMultAdd(BT,uold,f2,tmpp);VecScale(f2,-1.);
+		VecNorm(tmpp,NORM_2,&tmp6);
+
+		KSPInitialResidual(kspu,uold,tmpu,tmpu2,resu,rhsu);
+		VecNorm(resu,NORM_2,&tmp5);
+
+
+	}
+
+
 	/* Go solve Au0 = F1-Bp0*/
 	KSPSolve(kspu,rhsu,uold);
-	
+	KSPGetIterationNumber(kspu,&noIt);
+
+
+
+	if (VerboseConvergence())
+	{
+	
+	KSPGetIterationNumber(kspu,&tmpi);
+	VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.);
+	VecNorm(tmpu,NORM_2,&tmp2);
+	KSPGetResidualNorm(kspu,&tmp1);
+
+	_printf0_("||Au_00 - rhsu||_euc: " << tmp4 <<"\n||P(-1)(Au_00 - rhsu)||_euc: " << tmp5 << "\n ||Au_n0 - rhsu||_euc: " << tmp2<< "\n||P(-1)(Au_n0 - rhsu)||_euc: " << tmp1 << "\nIteration number: "<<tmpi<<"\n"); 	
+	_printf0_("||BTu_00 - f2||_euc: " << tmp6 <<"\n");
+	}
+
+
 
 	/* Set up u_new */
@@ -133,30 +253,64 @@
 	
 	/* Create KSP context */
-	KSPCreate(IssmComm::GetComm(),&kspp);
+	KSPCreate(IssmComm::GetComm(),&kspip);
 	#if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=5)
-	KSPSetOperators(kspp,IP,IP);
+	KSPSetOperators(kspip,IP,IP);
 	#else
-	KSPSetOperators(kspp,IP,IP,DIFFERENT_NONZERO_PATTERN);
+	KSPSetOperators(kspip,IP,IP,DIFFERENT_NONZERO_PATTERN);
 	#endif
 	
 	/* Create RHS */
 	/* RHS = BT * uold - F2 */
-	VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.);
+	VecScale(uold,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(uold,-1.);
+
 
 	/* Set KSP & PC options */
-	KSPSetType(kspp,KSPCG);
-	KSPSetInitialGuessNonzero(kspp,PETSC_TRUE);
-	KSPGetPC(kspp,&pcp);
+	KSPSetType(kspip,KSPGMRES);
+	KSPSetInitialGuessNonzero(kspip,PETSC_TRUE);
+	
+	
+	KSPGetPC(kspip,&pcp);
 	PCSetType(pcp,PCJACOBI);
-	/* Note: Systems in the pressure space are cheap, so we can afford a better tolerance */
-	KSPSetTolerances(kspp,1e-10,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
-	KSPSetUp(kspp);
-	
+	/* Note: Systems in the pressure space are cheap, so we can afford a good tolerance */
+	KSPSetTolerances(kspip,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
+	KSPSetUp(kspip);
+
+	if (VerboseConvergence())
+	{
+		KSPInitialResidual(kspip,gold,tmpp,tmpp2,resp,rhsp);
+		VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.);
+		VecNorm(resp,NORM_2,&tmp5);
+		VecNorm(tmpp,NORM_2,&tmp4);
+	}
+
+
+
 	/* Go solve */
-	KSPSolve(kspp,rhsp,gold);
-	
+	KSPSolve(kspip,rhsp,gold);
+
+
+
+	if (VerboseConvergence())
+	{
+	
+	KSPGetResidualNorm(kspip,&tmp1);
+	VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.);
+	VecNorm(tmpp,NORM_2,&tmp2);
+
+	KSPGetIterationNumber(kspip,&tmpi);
+	_printf0_("||IP*g00 - rhsp||_euc: " << tmp4 <<"\n||Jac(-1)*(IP*g00-rhsp)||_euc: " << tmp5 << "\n||IP*gn0-rhsp||_euc: " << tmp2<< "\n||Jac(-1)*(IP*gn0-rhsp)||_euc: " << tmp1 << "\nIteration number: "<<tmpi<<"\n"); 	
+	}
+
+
 	/*Initial residual*/
-	VecNorm(gold,NORM_INFINITY,&initRnorm);
-	
+	MatMult(IP,gold,tmpp);VecDot(gold,tmpp,&initRnorm);
+	initRnorm = sqrt(initRnorm);
+	_printf0_("inner product norm g0: "<<initRnorm<<"\n");
+
+
+	/*Iterate only if inital residual large enough*/
+	if(initRnorm > 1e-5)
+	{
+
 	/* Further setup */
 	VecDuplicate(gold,&gnew);VecCopy(gold,gnew);
@@ -171,9 +325,4 @@
 	VecAssemblyBegin(wold);VecAssemblyEnd(wold);
 
-	/*Realizing the step size part 1: thetam */
-	/*IP * theta = BT * uold - F2*/
-	VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.);
-	KSPSolve(kspp,rhsp,thetaold);
-
 
 	/* Count number of iterations */
@@ -182,15 +331,46 @@
 	/* CG iteration*/
 	for(;;){
+		if(VerboseConvergence())
+		{
+		 _printf0_("\n###### Step " << count << " ######\n");
+		}
 
 		/*Realizing the step size part 2: chim */
-		/*a(chim,v) = -b(wm,v)  i.e.  A * chim = -B * wm */
+		/*a(chim,v) = b(wm,v)  i.e.  A * chim = B * wm */
 		/*chim_DIR = 0*/
-		VecScale(wold,-1.);MatMult(B,wold,rhsu);VecScale(wold,-1.);
-		KSPSolve(kspu,rhsu,chi);
-
-		/*Realizing the step size part 3: etam */
-		MatMult(BT,chi,rhsp);
-		KSPSolve(kspp,rhsp,eta);
-	
+		VecScale(wold,1.);MatMult(B,wold,rhsu);VecScale(wold,1.);
+		VecSet(chi,0.0);
+	
+		
+		if(VerboseConvergence())
+		{
+		VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
+		VecNorm(tmpu,NORM_2,&tmp4);
+			
+		KSPInitialResidual(kspu,chi,tmpu,tmpu2,resu,rhsu);
+		VecNorm(resu,NORM_2,&tmp5);
+
+		}
+
+		
+			KSPSolve(kspu,rhsu,chi);
+		
+	
+		
+		
+		if (VerboseConvergence())
+		{
+		VecNorm(chi,NORM_2,&tmp1);
+		KSPGetResidualNorm(kspu,&tmp2);
+	
+		VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
+		VecNorm(tmpu,NORM_2,&tmp3);
+		
+		
+		KSPGetIterationNumber(kspu,&tmpi);
+		_printf0_("||chi_nk||_euc: "<< tmp1 << "\n||A*chi_0k - rhsu||_euc: "<<tmp4<< "\n||P(-1)*(A*chi_0k-rhsu)||_euc: " << tmp5 << "\n||A*chi_nk - rhsu||_euc: "<< tmp3 <<"\n||P(-1)*(A*chi_nk - rhsu)||_euc: " << tmp2 <<"\nIteration Number: " << tmpi <<"\n");
+		}
+
+
 		/* ---------------------------------------------------------- */
 
@@ -198,6 +378,9 @@
 		/*Set step size*/
 		/*rhom = [(wm)^T * IP^-1 * (BT * um - F2)]/[(wm)^T * IP^-1 * BT * chim]*/
-		VecDot(wold,thetaold,&rho);
-		VecDot(wold,eta,&tmpScalar);
+		MatMult(IP,gold,tmpp);
+		VecDot(tmpp,wold,&rho);
+		
+		MatMult(BT,chi,tmpp);
+		VecDot(tmpp,wold,&tmpScalar);
 		rho = rho/tmpScalar;
 
@@ -207,5 +390,5 @@
 
 		/*Pressure update*/
-		/*p(m+1) = pm - rhom * wm*/
+		/*p(m+1) = pm + rhom * wm*/
 		VecAXPY(p,-1.*rho,wold);
 
@@ -213,13 +396,20 @@
 		/*Velocity update*/
 		/*u(m+1) = um - rhom * chim*/
-		VecWAXPY(unew,-1.*rho,chi,uold);
-
-
-		/* ---------------------------------------------------------- */
-
-		/*Theta update*/
-		/*IP * theta = BT * uold - F2*/
-		VecScale(f2,-1.);MatMultAdd(BT,unew,f2,rhsp);VecScale(f2,-1.);
-		KSPSolve(kspp,rhsp,thetanew);
+		VecWAXPY(unew,rho,chi,uold);
+		VecNorm(unew,NORM_2,&tmpScalar);
+
+
+		if (VerboseConvergence())
+		{
+		VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.);
+		VecScale(rhsu,-1.);MatMultAdd(A,unew,rhsu,tmpu);VecScale(rhsu,-1.);
+		VecNorm(tmpu,NORM_2,&tmp6);
+
+		VecScale(f2,-1);MatMultAdd(BT,unew,f2,tmpp);VecScale(f2,-1);
+		VecNorm(tmpp,NORM_2,&tmp7);
+		_printf0_("Momentum balance norm: " << tmp6 <<"\n");
+		_printf0_("Incompressibility norm: " << tmp7 <<"\n");
+		}
+
 
 
@@ -228,22 +418,37 @@
 		/*Residual update*/
 		/*g(m+1) = gm - rhom * BT * chim*/
-		VecWAXPY(gnew,-1.*rho,eta,gold);
+		MatMult(BT,chi,tmpp);
+		MatMult(IP,gold,tmpp2);
+		VecWAXPY(rhsp,-1.*rho,tmpp,tmpp2);
+		KSPSolve(kspip,rhsp,gnew);
+		
+
 
 		/* ---------------------------------------------------------- */
 
-
-		/*BREAK if norm(g(m+0),2) < TOL or pressure space has been full searched*/
-		VecNorm(gnew,NORM_INFINITY,&rnorm);
+		MatMult(IP,gnew,tmpp);
+		
+		VecDot(tmpp,gnew,&gamma);
+		rnorm = sqrt(gamma);
+
+		/*BREAK if norm(g(m+0),2) < TOL or pressure space has been fully searched*/
 		if(rnorm < TOL*initRnorm) 
+		{
+			break;
+		}else if(rnorm > 1000*initRnorm)
+		{
+		 _printf0_("L2-Norm g: "<< rnorm << "\n");
+		 _printf0_("Solver diverged and ends prematurely.\n");
 		 break;
-		else if(rnorm > 100*initRnorm)
-		 _error_("Solver diverged. This shouldn't happen\n");
-		//else
-		// PetscPrintf(PETSC_COMM_WORLD,"rel. residual at step %d: %g, at TOL = %g\n",count,rnorm/initRnorm,TOL);
-	
-
-
-
-		if(count > np-1) break;
+		}else{
+			_printf0_("L2-Norm g: "<< rnorm << "\n");
+		}
+	
+		/*Break prematurely if solver doesn't reach desired tolerance in reasonable time frame*/
+		if(count > 10./TOL)   
+		{ 
+		 _printf0_("Ended after " << ceil(5./TOL) << " CG iterations\n");
+		 break;
+		}
 	
 
@@ -252,18 +457,16 @@
 
 		/*Directional update*/
-		/*gamma = [g(m+1)^T * theta(m+1)]/[g(m)^T * thetam]*/
-		VecDot(gnew,thetanew,&gamma);
-		VecDot(gold,thetaold,&tmpScalar);
+		MatMult(IP,gold,tmpp);
+		VecDot(tmpp,gold,&tmpScalar);
 		gamma = gamma/tmpScalar;
 
-		/*w(m+1) = g(m+1) + gamma * w(m)*/
 		VecWAXPY(wnew,gamma,wold,gnew);
 
 		/* Assign new to old iterates */
-		VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold);VecCopy(thetanew,thetaold);
+		VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold);
 		
 		count++;
 	}
-
+	}
 
 	/* Restore pressure and velocity sol. vectors to its global form */
@@ -276,5 +479,5 @@
 
 	/* Cleanup */
-	KSPDestroy(&kspu);KSPDestroy(&kspp);
+	KSPDestroy(&kspu);KSPDestroy(&kspip);
 
 	MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP);
@@ -282,6 +485,6 @@
 	VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp);
 	VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi);
-	VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f1);VecDestroy(&f2);VecDestroy(&eta);
-	VecDestroy(&thetanew);VecDestroy(&thetaold);
+	VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f1);VecDestroy(&f2);
+	VecDestroy(&tmpu2);VecDestroy(&tmpp2);VecDestroy(&resu);VecDestroy(&resp);
 
 }/*}}}*/
@@ -323,16 +526,10 @@
 	* Kff =  |      |
 	*			[B^T IP]
-   * To calculate the residual, only the necessary blocks need to be extracted */
+  /* To calculate the residual, only the necessary blocks need to be extracted */
 
 		/*Extract A, B, B^T */
-		#if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=8)
-		MatCreateSubMatrix(Kff->pmatrix->matrix,isv,isv,MAT_INITIAL_MATRIX,&A);
-		MatCreateSubMatrix(Kff->pmatrix->matrix,isv,isp,MAT_INITIAL_MATRIX,&B);
-		MatCreateSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT);
-		#else
 		MatGetSubMatrix(Kff->pmatrix->matrix,isv,isv,MAT_INITIAL_MATRIX,&A);
 		MatGetSubMatrix(Kff->pmatrix->matrix,isv,isp,MAT_INITIAL_MATRIX,&B);
 		MatGetSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT);
-		#endif
 	
 		/*no. of free nodes in velocity/pressure space*/
@@ -493,6 +690,9 @@
 
 
+
 	/*parameters:*/
 	int max_nonlinear_iterations;
+	int configuration_type;
+	int precond;
 	IssmDouble eps_res,eps_rel,eps_abs;
 
@@ -502,4 +702,5 @@
 	femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
 	femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
 	femmodel->UpdateConstraintsx();
 	int size;
@@ -507,4 +708,5 @@
 	bool converged=false;
 
+	
 	/*Start non-linear iteration using input velocity: */
 	GetSolutionFromInputsx(&ug,femmodel);
@@ -515,4 +717,5 @@
 	InputUpdateFromSolutionx(femmodel,ug);
 
+	
 	for(;;){
 
@@ -521,4 +724,6 @@
 		delete ug;
 
+
+		
 		/*Get stiffness matrix and Load vector*/
 		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
@@ -526,32 +731,50 @@
 		Reduceloadx(pf, Kfs, ys); delete Kfs;
 
-		/*Create mass matrix*/
-		StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
-		/*Get complete stiffness matrix without penalties*/
-		for(int i=0;i<femmodel->elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element);
-			if(Ie) Ie->AddToGlobal(Kff,NULL);
-			delete Ie;
-		}
-		Kff->Assemble();
+		/*Create pressure matrix of choice*/
+		PetscOptionsGetInt(PETSC_NULL,"-schur_pc",&precond,NULL);
+	
+		StressbalanceAnalysis* analysis = new StressbalanceAnalysis();	
+	
+		/* Construct Schur preconditioner matrix or mass matrix depending on input */
+		if(precond)
+		{
+			for(int i=0;i<femmodel->elements->Size();i++){
+				Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+				ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element);
+				if(Ie) Ie->AddToGlobal(Kff,NULL);
+				delete Ie;
+			}
+		}else{
+			
+			for(int i=0;i<femmodel->elements->Size();i++){
+				Element* element2=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+				ElementMatrix* Ie2 = analysis->CreatePressureMassMatrix(element2);
+				if(Ie2) Ie2->AddToGlobal(Kff,NULL);
+				delete Ie2;
+			}
+		}
+
+	   		
 		delete analysis;
-
+	
 		/*Obtain index sets for velocity and pressure components */
 		IS isv = NULL;
 		IS isp = NULL;
 		#if _PETSC_MAJOR_==3
-			/*Make indices out of doftypes: */
-			if(!(df->pvector->vector))_error_("need doftypes for FS solver!\n");
+
+		/*Make indices out of doftypes: */
+		if(!(df->pvector->vector))_error_("need doftypes for FS solver!\n");
 			DofTypesToIndexSet(&isv,&isp,df->pvector->vector,FSSolverEnum);
 		#else
 			_error_("Petsc 3.X required");
 		#endif
-
-
+		Kff->Assemble();
+
+	
 		/*Solve*/
 		femmodel->profiler->Start(SOLVER);
 		_assert_(Kff->type==PetscMatType); 
-		
+
+
 		SchurCGSolver(&uf,
 					Kff->pmatrix->matrix,
@@ -562,5 +785,4 @@
 					femmodel->parameters);
 		femmodel->profiler->Stop(SOLVER);
-	
 
 		/*Merge solution from f set to g set*/
@@ -575,4 +797,5 @@
 			converged=true;
 		}
+
 		InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
 		InputUpdateFromSolutionx(femmodel,ug);