Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 23203) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 23204) @@ -22,5 +22,5 @@ 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 f,BTAinvf; /* RHS of the global system */ + 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*/ @@ -33,5 +33,5 @@ /*STOP tolerance for the rel. residual*/ - TOL = 0.1; + TOL = 0.4; /*Initialize output*/ @@ -73,48 +73,25 @@ /* Set up intermediaries */ - VecDuplicate(uold,&f);VecSet(f,0.0); - VecAssemblyBegin(f);VecAssemblyEnd(f); - - VecDuplicate(p,&BTAinvf);VecSet(BTAinvf,0.0); - VecAssemblyBegin(BTAinvf);VecAssemblyEnd(BTAinvf); - + VecDuplicate(uold,&f1);VecSet(f1,0.0); + VecDuplicate(p,&f2);VecSet(f2,0.0); VecDuplicate(uold,&tmpu);VecSet(tmpu,0.0); - VecAssemblyBegin(tmpu);VecAssemblyEnd(tmpu); - VecDuplicate(p,&tmpp);VecSet(tmpp,0.0); - VecAssemblyBegin(tmpp);VecAssemblyEnd(tmpp); - VecDuplicate(p,&rhsp);VecSet(rhsp,0.0); - VecAssemblyBegin(rhsp);VecAssemblyEnd(rhsp); - VecDuplicate(uold,&rhsu);VecSet(rhsu,0.0); - VecAssemblyBegin(rhsu);VecAssemblyEnd(rhsu); - VecDuplicate(p,&gold);VecSet(gold,0.0); - VecAssemblyBegin(gold);VecAssemblyEnd(gold); - - VecDuplicate(gold,&wnew);VecSet(wnew,0.0); - VecAssemblyBegin(wnew);VecAssemblyEnd(wnew); - + VecDuplicate(p,&wnew);VecSet(wnew,0.0); VecDuplicate(uold,&chi);VecSet(chi,0.0); - VecAssemblyBegin(chi);VecAssemblyEnd(chi); - VecDuplicate(p,&thetanew);VecSet(thetanew,0.0); - VecAssemblyBegin(thetanew);VecAssemblyEnd(thetanew); - VecDuplicate(p,&thetaold);VecSet(thetaold,0.0); - VecAssemblyBegin(thetaold);VecAssemblyEnd(thetaold); - VecDuplicate(p,&eta);VecSet(eta,0.0); - VecAssemblyBegin(eta);VecAssemblyEnd(eta); - - /* Get global RHS (restricted to the velocity space */ - VecGetSubVector(pf,isv,&f); - + + /* Get global RHS (for each block sub-problem respectively)*/ + VecGetSubVector(pf,isv,&f1); + VecGetSubVector(pf,isp,&f2); /* ------------------------------------------------------------ */ /* Generate initial value for the velocity from the pressure */ - /* a(u0,v) = f(v)-b(p0,v) i.e. Au0 = F-Bp0 */ + /* a(u0,v) = f1(v)-b(p0,v) i.e. Au0 = F1-Bp0 */ /* u0 = u_DIR on \Gamma_DIR */ @@ -124,19 +101,18 @@ KSPSetType(kspu,KSPCG); KSPSetInitialGuessNonzero(kspu,PETSC_TRUE); + //KSPSetTolerances(kspu,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); + //KSPMonitorSet(kspu,KSPMonitorDefault,NULL,NULL); KSPGetPC(kspu,&pcu); - PCSetType(pcu,PCJACOBI); + PCSetType(pcu,PCSOR); KSPSetUp(kspu); /* Create RHS */ - /* RHS = F-B * pold */ - VecScale(p,-1.);MatMultAdd(B,p,f,rhsu);VecScale(p,-1.); - - /* Go solve Au0 = F-Bp0*/ + /* RHS = F1-B * pold */ + VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.); + + /* Go solve Au0 = F1-Bp0*/ KSPSolve(kspu,rhsu,uold); - /* Calculate B^T * A^{-1} * F for future reference */ - KSPSolve(kspu,f,tmpu);MatMult(BT,tmpu,BTAinvf); - /* Set up u_new */ @@ -149,5 +125,5 @@ /*Get initial residual*/ - /*(1/mu(x) * g0, q) = b(q,u0) i.e. IP * g0 = BT * u0*/ + /*(1/mu(x) * g0, q) = b(q,u0) - (f2,q) i.e. IP * g0 = BT * u0 - F2*/ /* Create KSP context */ @@ -156,6 +132,6 @@ /* Create RHS */ - /* RHS = BT * uold */ - MatMult(BT,uold,rhsp); + /* RHS = BT * uold - F2 */ + VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.); /* Set KSP & PC options */ @@ -187,6 +163,6 @@ /*Realizing the step size part 1: thetam */ - /*IP * theta = BT * uold*/ - MatMult(BT,uold,rhsp); + /*IP * theta = BT * uold - F2*/ + VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.); KSPSolve(kspp,rhsp,thetaold); @@ -212,5 +188,5 @@ /*Set step size*/ - /*rhom = [(wm)^T * IP^-1 * BT * um]/[(wm)^T * IP^-1 * BT * chim]*/ + /*rhom = [(wm)^T * IP^-1 * (BT * um - F2)]/[(wm)^T * IP^-1 * BT * chim]*/ VecDot(wold,thetaold,&rho); VecDot(wold,eta,&tmpScalar); @@ -234,6 +210,6 @@ /*Theta update*/ - /*IP * theta = BT * uold*/ - MatMult(BT,unew,rhsp); + /*IP * theta = BT * uold - F2*/ + VecScale(f2,-1.);MatMultAdd(BT,unew,f2,rhsp);VecScale(f2,-1.); KSPSolve(kspp,rhsp,thetanew); @@ -250,7 +226,13 @@ /*BREAK if norm(g(m+0),2) < TOL or pressure space has been full searched*/ VecNorm(gnew,NORM_INFINITY,&rnorm); - if(rnorm < TOL*initRnorm) break; + if(rnorm < TOL*initRnorm) + 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; - + /* ---------------------------------------------------------- */ @@ -272,4 +254,5 @@ } + /* Restore pressure and velocity sol. vectors to its global form */ VecRestoreSubVector(out_uf->pvector->vector,isv,&unew); @@ -287,5 +270,5 @@ 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(&f);VecDestroy(&BTAinvf);VecDestroy(&eta); + VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f1);VecDestroy(&f2);VecDestroy(&eta); VecDestroy(&thetanew);VecDestroy(&thetaold); @@ -342,5 +325,5 @@ /*Create mass matrix*/ int fsize; Kff->GetSize(&fsize,&fsize); - Iff=new Matrix(fsize,fsize,200,4); + Iff=new Matrix(fsize,fsize,300,4); StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); /*Get complete stiffness matrix without penalties*/ @@ -367,6 +350,5 @@ femmodel->profiler->Stop(SOLVER); delete Iff; - - + /*Merge solution from f set to g set*/ Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;