Changeset 23205

Timestamp:

08/31/18 12:39:46 (7 years ago)

Author:

wester

Message:

CHG: Reduced SchurCG storage req. by storing preconditioner in Kff. Added adjusted routine to check for convergence.

Location:

issm/trunk-jpl/src/c/solutionsequences

Files:

• convergence.cpp (modified) (3 diffs)
• solutionsequence_schurcg.cpp (modified) (10 diffs)

issm/trunk-jpl/src/c/solutionsequences/convergence.cpp

@@ -7 +7 @@
 #include "../shared/shared.h"
 
-void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){
+void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){/*{{{*/
 
         /*output*/
@@ -34 +34 @@
         }
 
+
         /*Display solver caracteristics*/
         if (VerboseConvergence()){
+
+        
 
                 /*compute KUF = KU - F = K*U - F*/
@@ -136 +139 @@
         /*assign output*/
         *pconverged=converged;
-}
+}/*}}}*/
+
+

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp

@@ -12 +12 @@
 
 #ifdef _HAVE_PETSC_
-void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff,Mat Iff,Vec pf, Vec uf0,Vec df,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 */
-        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 +39 @@
         Vector<IssmDouble>* out_uf=new Vector<IssmDouble>(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 +48 @@
         
         /* 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 +221 @@
                 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 +267 @@
         VecDestroy(&thetanew);VecDestroy(&thetaold);
 
+}/*}}}*/
+void convergence_schurcg(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs,IS isv,IS isp){/*{{{*/
+
+        /*output*/
+        bool converged=false;
+
+        /*intermediary*/
+        //Vector<IssmDouble>* KU=NULL;
+        //Vector<IssmDouble>* KUF=NULL;
+        //Vector<IssmDouble>* KUold=NULL;
+        //Vector<IssmDouble>* KUoldF=NULL;
+        Vector<IssmDouble>* 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<IssmDouble>(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<IssmDouble>(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_res){
+                if(VerboseConvergence()) _printf0_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<< " < "<<eps_res*100<<" %\n");
+                converged=true;
+        }
+        else{ 
+                if(VerboseConvergence()) _printf0_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<<" > "<<eps_res*100<<" %\n");
+                converged=false;
+        }
+
+        /*Relative criterion (optional)*/
+        if (!xIsNan<IssmDouble>(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<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+
+                //clean up
+                delete duf;
+
+                //print
+                if (!xIsNan<IssmDouble>(eps_rel)){
+                        if((ndu/nu)<eps_rel){
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
+                        }
+                        else{ 
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << 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<IssmDouble>(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<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+
+                //clean up
+                delete duf;
+
+                //print
+                if (!xIsNan<IssmDouble>(eps_abs)){
+                        if ((nduinf)<eps_abs){
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf << " < " << eps_abs << "\n");
+                        }
+                        else{
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << 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 +467 @@
         Vector<IssmDouble>* df  = NULL;
         Vector<IssmDouble>* ys  = NULL;
-        Matrix<IssmDouble>* Iff = NULL;
 
 
@@ -324 +505 @@
 
                 /*Create mass matrix*/
-                int fsize; Kff->GetSize(&fsize,&fsize);
-                Iff=new Matrix<IssmDouble>(fsize,fsize,300,4);
                 StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
                 /*Get complete stiffness matrix without penalties*/
@@ -331 +510 @@
                         Element* element=xDynamicCast<Element*>(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 +534 @@
                 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 +571 @@
 
 }/*}}}*/
+
 #else
 void solutionsequence_schurcg(FemModel* femmodel){_error_("PETSc needs to be installed");}