Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 23196) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 23197) @@ -20,4 +20,8 @@ Vector* df = NULL; Vector* ys = NULL; + + Vec res; + double resNorm; + Loads* savedloads=NULL; @@ -70,4 +74,20 @@ Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); femmodel->profiler->Stop(SOLVER); + + +// VecDuplicate(uf->pvector->vector,&res); +// VecCopy(uf->pvector->vector,res); +// VecAssemblyBegin(res); +// VecAssemblyEnd(res); +// VecScale(pf->pvector->vector,-1.); +// MatMultAdd(Kff->pmatrix->matrix,uf->pvector->vector,pf->pvector->vector,res); +// VecScale(pf->pvector->vector,-1.); +// VecView(pf->pvector->vector,PETSC_VIEWER_STDOUT_WORLD); +// VecNorm(res,NORM_INFINITY,&resNorm); +// PetscPrintf(PETSC_COMM_WORLD,"Count = %d, Res. norm: %g\n", count,resNorm); + + + + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys; @@ -86,5 +106,5 @@ } } - + /*Increase count: */ count++; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 23196) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 23197) @@ -10,16 +10,284 @@ #include "../analyses/analyses.h" + #ifdef _HAVE_PETSC_ -void SchurCGSolver(Vector** puf,Mat Kff,Mat Mff,Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/ +void SchurCGSolver(Vector** puf,Mat Kff,Mat Iff,Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/ + + Mat A, B, BT; /* Saddle point block matrices */ + Mat IP; /* Preconditioner matrix */ + IS isv=NULL; /* Index set free velocity nodes */ + IS isp=NULL; /* Index set free pressure nodes */ + 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 f,BTAinvf; /* RHS of the global system */ + double rho,gamma,tmpScalar; /* Step sizes, arbitrary double */ + KSP kspu,kspp; /* 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 */ + 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.1; /*Initialize output*/ - Vec uf = NULL; - - _error_("not implemented yet"); + Vector* out_uf=new Vector(uf0); + + /* Get velocity and pressure index sets for extraction */ + #if _PETSC_MAJOR_==3 + /*Make indices out of doftypes: */ + if(!df)_error_("need doftypes for FS solver!\n"); + DofTypesToIndexSet(&isv,&isp,df,FSSolverEnum); + #else + _error_("Petsc 3.X required"); + #endif + + + /* Extract block matrices from the saddle point matrix */ + /* [ A B ] = Kff + * [ B^T 0 ] + * */ + MatGetSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A); + MatGetSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B); + MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT); + + /* Extract preconditioner matrix on the pressure space*/ + MatGetSubMatrix(Iff,isp,isp,MAT_INITIAL_MATRIX,&IP); + + /* Get number of velocity / pressure nodes */ + MatGetSize(B,&nu,&np); + + /* Extract initial guesses for uold and pold */ + VecCreate(IssmComm::GetComm(),&p);VecSetSizes(p,PETSC_DECIDE,np);VecSetFromOptions(p); + VecAssemblyBegin(p);VecAssemblyEnd(p); + VecCreate(IssmComm::GetComm(),&uold);VecSetSizes(uold,PETSC_DECIDE,nu);VecSetFromOptions(uold); + VecAssemblyBegin(uold);VecAssemblyEnd(uold); + + VecGetSubVector(out_uf->pvector->vector,isv,&uold); + VecGetSubVector(out_uf->pvector->vector,isp,&p); + + + /* 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,&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(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); + + + /* ------------------------------------------------------------ */ + + /* Generate initial value for the velocity from the pressure */ + /* a(u0,v) = f(v)-b(p0,v) i.e. Au0 = F-Bp0 */ + /* u0 = u_DIR on \Gamma_DIR */ + + /* Create KSP context */ + KSPCreate(IssmComm::GetComm(),&kspu); + KSPSetOperators(kspu,A,A); + KSPSetType(kspu,KSPCG); + KSPSetInitialGuessNonzero(kspu,PETSC_TRUE); + KSPGetPC(kspu,&pcu); + PCSetType(pcu,PCJACOBI); + KSPSetUp(kspu); + + + /* Create RHS */ + /* RHS = F-B * pold */ + VecScale(p,-1.);MatMultAdd(B,p,f,rhsu);VecScale(p,-1.); + + /* Go solve Au0 = F-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 */ + VecDuplicate(uold,&unew);VecCopy(uold,unew); + VecAssemblyBegin(unew);VecAssemblyEnd(unew); + + + + /* ------------------------------------------------------------- */ + + /*Get initial residual*/ + /*(1/mu(x) * g0, q) = b(q,u0) i.e. IP * g0 = BT * u0*/ + + /* Create KSP context */ + KSPCreate(IssmComm::GetComm(),&kspp); + KSPSetOperators(kspp,IP,IP); + + /* Create RHS */ + /* RHS = BT * uold */ + MatMult(BT,uold,rhsp); + + /* Set KSP & PC options */ + KSPSetType(kspp,KSPCG); + KSPSetInitialGuessNonzero(kspp,PETSC_TRUE); + KSPGetPC(kspp,&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); + + /* Go solve */ + KSPSolve(kspp,rhsp,gold); + + /*Initial residual*/ + VecNorm(gold,NORM_INFINITY,&initRnorm); + + /* Further setup */ + VecDuplicate(gold,&gnew);VecCopy(gold,gnew); + VecAssemblyBegin(gnew);VecAssemblyEnd(gnew); + + + /* ------------------------------------------------------------ */ + + /*Set initial search direction*/ + /*w0 = g0*/ + VecDuplicate(gold,&wold);VecCopy(gold,wold); + VecAssemblyBegin(wold);VecAssemblyEnd(wold); + + /*Realizing the step size part 1: thetam */ + /*IP * theta = BT * uold*/ + MatMult(BT,uold,rhsp); + KSPSolve(kspp,rhsp,thetaold); + + + /* Count number of iterations */ + int count = 0; + + /* CG iteration*/ + for(;;){ + + /*Realizing the step size part 2: chim */ + /*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); + + /* ---------------------------------------------------------- */ + + + /*Set step size*/ + /*rhom = [(wm)^T * IP^-1 * BT * um]/[(wm)^T * IP^-1 * BT * chim]*/ + VecDot(wold,thetaold,&rho); + VecDot(wold,eta,&tmpScalar); + rho = rho/tmpScalar; + + + /* ---------------------------------------------------------- */ + + + /*Pressure update*/ + /*p(m+1) = pm - rhom * wm*/ + VecAXPY(p,-1.*rho,wold); + + + /*Velocity update*/ + /*u(m+1) = um - rhom * chim*/ + VecWAXPY(unew,-1.*rho,chi,uold); + + + /* ---------------------------------------------------------- */ + + /*Theta update*/ + /*IP * theta = BT * uold*/ + MatMult(BT,unew,rhsp); + KSPSolve(kspp,rhsp,thetanew); + + + /* ---------------------------------------------------------- */ + + /*Residual update*/ + /*g(m+1) = gm - rhom * BT * chim*/ + VecWAXPY(gnew,-1.*rho,eta,gold); + + /* ---------------------------------------------------------- */ + + + /*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(count > np-1) break; + + + /* ---------------------------------------------------------- */ + + + /*Directional update*/ + /*gamma = [g(m+1)^T * theta(m+1)]/[g(m)^T * thetam]*/ + VecDot(gnew,thetanew,&gamma); + VecDot(gold,thetaold,&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); + + count++; + } + + /* Restore pressure and velocity sol. vectors to its global form */ + VecRestoreSubVector(out_uf->pvector->vector,isv,&unew); + VecRestoreSubVector(out_uf->pvector->vector,isp,&p); /*return output pointer*/ - Vector* out_uf=new Vector(uf); - VecFree(&uf); *puf=out_uf; + + + /* Cleanup */ + KSPDestroy(&kspu);KSPDestroy(&kspp); + + MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP); + + 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(&thetanew);VecDestroy(&thetaold); + }/*}}}*/ void solutionsequence_schurcg(FemModel* femmodel){/*{{{*/ @@ -32,8 +300,8 @@ Vector* old_uf = NULL; Vector* pf = NULL; - Vector* pf0 = NULL; Vector* df = NULL; Vector* ys = NULL; - Matrix* Mff = NULL; + Matrix* Iff = NULL; + /*parameters:*/ @@ -49,5 +317,5 @@ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); femmodel->UpdateConstraintsx(); - + int size; int count=0; bool converged=false; @@ -69,5 +337,4 @@ /*Get stiffness matrix and Load vector*/ SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); - pf0=pf->Duplicate(); pf->Copy(pf0); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf, Kfs, ys); delete Kfs; @@ -75,14 +342,14 @@ /*Create mass matrix*/ int fsize; Kff->GetSize(&fsize,&fsize); - Mff=new Matrix(fsize,fsize,100,4); + Iff=new Matrix(fsize,fsize,200,4); StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); /*Get complete stiffness matrix without penalties*/ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - ElementMatrix* Me = analysis->CreatePressureMassMatrix(element); - if(Me) Me->AddToGlobal(Mff,NULL); - delete Me; + ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element); + if(Ie) Ie->AddToGlobal(Iff,NULL); + delete Ie; } - Mff->Assemble(); + Iff->Assemble(); delete analysis; @@ -90,7 +357,8 @@ femmodel->profiler->Start(SOLVER); _assert_(Kff->type==PetscMatType); + SchurCGSolver(&uf, Kff->pmatrix->matrix, - Mff->pmatrix->matrix, + Iff->pmatrix->matrix, pf->pvector->vector, old_uf->pvector->vector, @@ -98,5 +366,6 @@ femmodel->parameters); femmodel->profiler->Stop(SOLVER); - delete pf0; delete Mff; + delete Iff; + /*Merge solution from f set to g set*/ @@ -106,4 +375,5 @@ convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df; count++; + if(count>=max_nonlinear_iterations){ _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");