Index: /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp (revision 23204) +++ /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp (revision 23205) @@ -7,5 +7,5 @@ #include "../shared/shared.h" -void convergence(bool* pconverged, Matrix* Kff,Vector* pf,Vector* uf,Vector* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){ +void convergence(bool* pconverged, Matrix* Kff,Vector* pf,Vector* uf,Vector* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){/*{{{*/ /*output*/ @@ -34,6 +34,9 @@ } + /*Display solver caracteristics*/ if (VerboseConvergence()){ + + /*compute KUF = KU - F = K*U - F*/ @@ -136,3 +139,5 @@ /*assign output*/ *pconverged=converged; -} +}/*}}}*/ + + Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 23204) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 23205) @@ -12,10 +12,11 @@ #ifdef _HAVE_PETSC_ -void SchurCGSolver(Vector** puf,Mat Kff,Mat Iff,Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/ + + +void SchurCGSolver(Vector** 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 */ - IS isv=NULL; /* Index set free velocity nodes */ - IS isp=NULL; /* Index set free pressure nodes */ + Mat IP2; int nu, np; /* No of. free nodes in velocity / pressure space */ Vec p,uold,unew; /* Solution vectors for pressure / vel. */ @@ -38,14 +39,4 @@ 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 @@ -57,6 +48,6 @@ /* Extract preconditioner matrix on the pressure space*/ - MatGetSubMatrix(Iff,isp,isp,MAT_INITIAL_MATRIX,&IP); - + MatGetSubMatrix(Kff,isp,isp,MAT_INITIAL_MATRIX,&IP); + /* Get number of velocity / pressure nodes */ MatGetSize(B,&nu,&np); @@ -230,7 +221,10 @@ 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); - + //else + // PetscPrintf(PETSC_COMM_WORLD,"rel. residual at step %d: %g, at TOL = %g\n",count,rnorm/initRnorm,TOL); + + + + if(count > np-1) break; @@ -273,4 +267,192 @@ VecDestroy(&thetanew);VecDestroy(&thetaold); +}/*}}}*/ +void convergence_schurcg(bool* pconverged, Matrix* Kff,Vector* pf,Vector* uf,Vector* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs,IS isv,IS isp){/*{{{*/ + + /*output*/ + bool converged=false; + + /*intermediary*/ + //Vector* KU=NULL; + //Vector* KUF=NULL; + //Vector* KUold=NULL; + //Vector* KUoldF=NULL; + Vector* duf=NULL; + IssmDouble ndu,nduinf,nu; + IssmDouble nKUF; + IssmDouble nKUoldF; + IssmDouble nF; + IssmDouble solver_residue,res; + int analysis_type; + + Mat A, B, BT; + Vec u,p,uold,pold,f1,f2,tmp,res1,res2; + int n_u,n_p; + double rnorm1, rnorm2; + + + if(VerboseModule()) _printf0_(" checking convergence\n"); + + /*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from + * the get go: */ + if(uf->IsEmpty()){ + *pconverged=true; + return; + } + + /* Note: SchurCG also constructs the Schur preconditioner and stores it in the free block of Kff */ + /* [A B] + * Kff = | | + * [B^T IP] + /* To calculate the residual, only the necessary blocks need to be extracted */ + + /*Extract A, B, B^T */ + 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); + + /*no. of free nodes in velocity/pressure space*/ + MatGetSize(B,&n_u,&n_p); + + /*Extract values corresponding to the free velocity/pressure nodes*/ + VecCreate(IssmComm::GetComm(),&p);VecSetSizes(p,PETSC_DECIDE,n_p);VecSetFromOptions(p); + VecAssemblyBegin(p);VecAssemblyEnd(p); + VecCreate(IssmComm::GetComm(),&u);VecSetSizes(u,PETSC_DECIDE,n_u);VecSetFromOptions(u); + VecAssemblyBegin(u);VecAssemblyEnd(u); + + VecGetSubVector(uf->pvector->vector,isv,&u); + VecGetSubVector(uf->pvector->vector,isp,&p); + + + /*Extract values of the RHS corresponding to the first/second block*/ + VecDuplicate(u,&f1);VecSet(f1,1.0); + VecDuplicate(p,&f2);VecSet(f2,1.0); + VecGetSubVector(pf->pvector->vector,isv,&f1); + VecGetSubVector(pf->pvector->vector,isp,&f2); + + /*Allocate intermediaries*/ + VecDuplicate(u,&res1);VecSet(res1,1.0); + VecDuplicate(u,&tmp);VecSet(tmp,1.0); + VecDuplicate(p,&res2);VecSet(res2,1.0); + + + /*Display solver caracteristics*/ + if (VerboseConvergence()){ + + /*Calculate res1 = A*u + B*p - f1*/ + VecScale(f1,-1.);MatMultAdd(A,u,f1,tmp);MatMultAdd(B,p,tmp,res1);VecScale(f1,-1.); + /*Calculate res2 = B^T * u - f2*/ + VecScale(f2,-1.);MatMultAdd(BT,u,f2,res2);VecScale(f2,-1.); + + + /*compute norm(res1), norm(res2), norm(F) and residue*/ + VecNorm(res1,NORM_2,&rnorm1);VecNorm(res2,NORM_2,&rnorm2); + nKUF=sqrt(rnorm1*rnorm1 + rnorm2*rnorm2); + nF=pf->Norm(NORM_TWO); + solver_residue=nKUF/nF; + _printf0_("\n" << " solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n"); + if(xIsNan(solver_residue)){ + //Kff->Echo(); + } + + } + /*clean up*/ + VecRestoreSubVector(uf->pvector->vector,isv,&u); + VecRestoreSubVector(uf->pvector->vector,isp,&p); + + /*Extract values corresponding to velocity/pressure on the old solution*/ + VecGetSubVector(old_uf->pvector->vector,isv,&uold); + VecGetSubVector(old_uf->pvector->vector,isp,&pold); + + + /*Force equilibrium (Mandatory)*/ + + /*Calculate res1 = A*uold + B*pold - f1*/ + VecScale(f1,-1.);MatMultAdd(A,uold,f1,tmp);MatMultAdd(B,pold,tmp,res1);VecScale(f1,-1.); + /*Calculate res2 = B^T * uold - f2*/ + VecScale(f2,-1.);MatMultAdd(BT,uold,f2,res2);VecScale(f2,-1.); + + /*compute norm(res1), norm(res2), norm(F) and residue*/ + VecNorm(res1,NORM_2,&rnorm1);VecNorm(res2,NORM_2,&rnorm2); + nKUoldF=sqrt(rnorm1*rnorm1 + rnorm2*rnorm2); + nF=pf->Norm(NORM_TWO); + res=nKUoldF/nF; + if (xIsNan(res)){ + _printf0_("norm nf = " << nF << "f and norm kuold = " << nKUoldF << "f\n"); + _error_("mechanical equilibrium convergence criterion is NaN!"); + } + + MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT); + VecRestoreSubVector(pf->pvector->vector,isv,&f1); + VecRestoreSubVector(pf->pvector->vector,isp,&f2); + VecDestroy(&res1);VecDestroy(&res2);VecDestroy(&tmp); + VecRestoreSubVector(old_uf->pvector->vector,isv,&uold); + VecRestoreSubVector(old_uf->pvector->vector,isp,&pold); + + + + //print + if(res "<(eps_rel) || (VerboseConvergence())){ + + //compute norm(du)/norm(u) + duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0); + ndu=duf->Norm(NORM_TWO); nu=old_uf->Norm(NORM_TWO); + + if (xIsNan(ndu) || xIsNan(nu)) _error_("convergence criterion is NaN!"); + + //clean up + delete duf; + + //print + if (!xIsNan(eps_rel)){ + if((ndu/nu) " << eps_rel*100 << " %\n"); + converged=false; + } + } + else _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " %\n"); + + } + + /*Absolute criterion (Optional) = max(du)*/ + if (!xIsNan(eps_abs) || (VerboseConvergence())){ + + //compute max(du) + duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0); + ndu=duf->Norm(NORM_TWO); nduinf=duf->Norm(NORM_INF); + if (xIsNan(ndu) || xIsNan(nu)) _error_("convergence criterion is NaN!"); + + //clean up + delete duf; + + //print + if (!xIsNan(eps_abs)){ + if ((nduinf) " << eps_abs << "\n"); + converged=false; + } + } + else _printf0_(setw(50) << left << " Convergence criterion: max(du)" << nduinf << "\n"); + + } + + /*assign output*/ + *pconverged=converged; }/*}}}*/ void solutionsequence_schurcg(FemModel* femmodel){/*{{{*/ @@ -285,5 +467,4 @@ Vector* df = NULL; Vector* ys = NULL; - Matrix* Iff = NULL; @@ -324,6 +505,4 @@ /*Create mass matrix*/ - int fsize; Kff->GetSize(&fsize,&fsize); - Iff=new Matrix(fsize,fsize,300,4); StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); /*Get complete stiffness matrix without penalties*/ @@ -331,9 +510,21 @@ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element); - if(Ie) Ie->AddToGlobal(Iff,NULL); + if(Ie) Ie->AddToGlobal(Kff,NULL); delete Ie; } - Iff->Assemble(); + Kff->Assemble(); 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"); + DofTypesToIndexSet(&isv,&isp,df->pvector->vector,FSSolverEnum); + #else + _error_("Petsc 3.X required"); + #endif + /*Solve*/ @@ -343,17 +534,17 @@ SchurCGSolver(&uf, Kff->pmatrix->matrix, - Iff->pmatrix->matrix, pf->pvector->vector, old_uf->pvector->vector, - df->pvector->vector, + isv, + isp, femmodel->parameters); femmodel->profiler->Stop(SOLVER); - delete Iff; - + + /*Merge solution from f set to g set*/ Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys; /*Check for convergence and update inputs accordingly*/ - convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df; + convergence_schurcg(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs,isv,isp); delete Kff; delete pf; delete df; count++; @@ -380,4 +571,5 @@ }/*}}}*/ + #else void solutionsequence_schurcg(FemModel* femmodel){_error_("PETSc needs to be installed");}