Ignore:
Timestamp:
12/22/21 10:39:44 (3 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 26742

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        2020issm
        2121kriging
         22issm_slc
        2223issm_slr
        2324issm_ocean

  • issm/trunk/src/c/solutionsequences/solutionsequence_schurcg.cpp

    r25836 r26744  
    3333        IssmPDouble                             t1,t2;                                  /* Time measurement for bottleneck analysis */
    3434
    35        
    3635        double tmp1,tmp2,tmp3;
    3736        int tmpi;
    … …  
    3938
    4039        int noIt;
    41 
    42 
    4340
    4441        int precond = 0;
    … …  
    4946        PetscBool flag,flg;
    5047        #endif
    51 
    5248
    5349        char ksp_type[50];
    … …  
    7167        #endif
    7268
    73        
    7469        if(precond){
    7570                _printf0_("Running WITH preconditioner\n");
    … …  
    7873        }
    7974
    80 
    8175        /*Initialize output*/
    8276        Vector<IssmDouble>* out_uf=new Vector<IssmDouble>(uf0);
    83        
     77
    8478        /* Extract block matrices from the saddle point matrix */
    8579        /* [ A   B ] = Kff
    … …  
    9791        MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT);
    9892        #endif
    99        
     93
    10094        /* Extract preconditioner matrix on the pressure space*/
    10195        #if PETSC_VERSION_GT(3,8,0)
    … …  
    10599        #endif
    106100
    107 
    108101        /* Get number of velocity / pressure nodes */
    109102        MatGetSize(B,&nu,&np);
    … …  
    117110        VecGetSubVector(out_uf->pvector->vector,isv,&uold);
    118111        VecGetSubVector(out_uf->pvector->vector,isp,&p);
    119 
    120112
    121113        /* Set up intermediaries */
    … …  
    133125        VecDuplicate(p,&wnew);VecSet(wnew,0.0);
    134126        VecDuplicate(uold,&chi);VecSet(chi,0.0);
    135        
     127
    136128        /* Get global RHS (for each block sub-problem respectively)*/
    137129        VecGetSubVector(pf,isv,&f1);
    … …  
    143135        /* a(u0,v) = f1(v)-b(p0,v)  i.e.  Au0 = F1-Bp0 */
    144136        /* u0 = u_DIR on \Gamma_DIR */
    145        
     137
    146138        /* Create KSP context */
    147139        KSPCreate(IssmComm::GetComm(),&kspu);
    … …  
    174166                _error_("Suggested KSP method not implemented yet!\n");
    175167        }
    176        
     168
    177169        KSPSetInitialGuessNonzero(kspu,PETSC_TRUE);
    178        
     170
    179171        /*Strong rel. residual tolerance needed when solving for the velocity update.
    180172         * This is because ISSM uses the dimensional equations, so the initial guess chi = 0
    181173         * results in a very high initial residual (~ 1e19)*/
    182174        KSPSetTolerances(kspu,ELLTOL,PETSC_DEFAULT,PETSC_DEFAULT,maxiter); //maxiter
    183 
    184 
    185175
    186176        KSPGetPC(kspu,&pcu);
    … …  
    210200        KSPSetUp(kspu);
    211201
    212        
    213202        /* Create RHS */
    214203        /* RHS = F1-B * pold */
    … …  
    227216                VecNorm(resu,NORM_2,&tmp5);
    228217
    229 
    230         }
    231 
     218        }
    232219
    233220        /* Go solve Au0 = F1-Bp0*/
    … …  
    235222        KSPGetIterationNumber(kspu,&noIt);
    236223
    237 
    238 
    239224        if (VerboseConvergence())
    240225        {
    241        
     226
    242227        KSPGetIterationNumber(kspu,&tmpi);
    243228        VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.);
    … …  
    249234        }
    250235
    251 
    252 
    253236        /* Set up u_new */
    254237        VecDuplicate(uold,&unew);VecCopy(uold,unew);
    255238        VecAssemblyBegin(unew);VecAssemblyEnd(unew);
    256239
    257 
    258 
    259240        /* ------------------------------------------------------------- */
    260241
    261242        /*Get initial residual*/
    262243        /*(1/mu(x) * g0, q) = b(q,u0) - (f2,q)  i.e.  IP * g0 = BT * u0 - F2*/
    263        
     244
    264245        /* Create KSP context */
    265246        KSPCreate(IssmComm::GetComm(),&kspip);
    … …  
    269250        KSPSetOperators(kspip,IP,IP,DIFFERENT_NONZERO_PATTERN);
    270251        #endif
    271        
     252
    272253        /* Create RHS */
    273254        /* RHS = BT * uold - F2 */
    274255        VecScale(uold,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(uold,-1.);
    275256
    276 
    277257        /* Set KSP & PC options */
    278258        KSPSetType(kspip,KSPGMRES);
    279259        KSPSetInitialGuessNonzero(kspip,PETSC_TRUE);
    280        
    281        
     260
    282261        KSPGetPC(kspip,&pcp);
    283262        PCSetType(pcp,PCJACOBI);
    … …  
    294273        }
    295274
    296 
    297 
    298275        /* Go solve */
    299276        KSPSolve(kspip,rhsp,gold);
    300277
    301 
    302 
    303278        if (VerboseConvergence())
    304279        {
    305        
     280
    306281        KSPGetResidualNorm(kspip,&tmp1);
    307282        VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.);
    … …  
    311286        _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");   
    312287        }
    313 
    314288
    315289        /*Initial residual*/
    … …  
    318292        _printf0_("inner product norm g0: "<<initRnorm<<"\n");
    319293
    320 
    321294        /*Iterate only if inital residual large enough*/
    322295        if(initRnorm > 1e-5)
    … …  
    326299        VecDuplicate(gold,&gnew);VecCopy(gold,gnew);
    327300        VecAssemblyBegin(gnew);VecAssemblyEnd(gnew);
    328 
    329301
    330302        /* ------------------------------------------------------------ */
    … …  
    334306        VecDuplicate(gold,&wold);VecCopy(gold,wold);
    335307        VecAssemblyBegin(wold);VecAssemblyEnd(wold);
    336 
    337308
    338309        /* Count number of iterations */
    … …  
    351322                VecScale(wold,1.);MatMult(B,wold,rhsu);VecScale(wold,1.);
    352323                VecSet(chi,0.0);
    353        
    354                
     324
    355325                if(VerboseConvergence())
    356326                {
    357327                VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
    358328                VecNorm(tmpu,NORM_2,&tmp4);
    359                        
     329
    360330                KSPInitialResidual(kspu,chi,tmpu,tmpu2,resu,rhsu);
    361331                VecNorm(resu,NORM_2,&tmp5);
    … …  
    363333                }
    364334
    365                
    366335                        KSPSolve(kspu,rhsu,chi);
    367                
    368        
    369                
    370                
     336
    371337                if (VerboseConvergence())
    372338                {
    373339                VecNorm(chi,NORM_2,&tmp1);
    374340                KSPGetResidualNorm(kspu,&tmp2);
    375        
     341
    376342                VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
    377343                VecNorm(tmpu,NORM_2,&tmp3);
    378                
    379                
     344
    380345                KSPGetIterationNumber(kspu,&tmpi);
    381346                _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");
    382347                }
    383348
    384 
    385349                /* ---------------------------------------------------------- */
    386 
    387350
    388351                /*Set step size*/
    … …  
    390353                MatMult(IP,gold,tmpp);
    391354                VecDot(tmpp,wold,&rho);
    392                
     355
    393356                MatMult(BT,chi,tmpp);
    394357                VecDot(tmpp,wold,&tmpScalar);
    395358                rho = rho/tmpScalar;
    396359
    397 
    398360                /* ---------------------------------------------------------- */
    399 
    400361
    401362                /*Pressure update*/
    402363                /*p(m+1) = pm + rhom * wm*/
    403364                VecAXPY(p,-1.*rho,wold);
    404 
    405365
    406366                /*Velocity update*/
    … …  
    408368                VecWAXPY(unew,rho,chi,uold);
    409369                VecNorm(unew,NORM_2,&tmpScalar);
    410 
    411370
    412371                if (VerboseConvergence())
    … …  
    421380                _printf0_("Incompressibility norm: " << tmp7 <<"\n");
    422381                }
    423 
    424 
    425382
    426383                /* ---------------------------------------------------------- */
    … …  
    432389                VecWAXPY(rhsp,-1.*rho,tmpp,tmpp2);
    433390                KSPSolve(kspip,rhsp,gnew);
    434                
    435 
    436391
    437392                /* ---------------------------------------------------------- */
    438393
    439394                MatMult(IP,gnew,tmpp);
    440                
     395
    441396                VecDot(tmpp,gnew,&gamma);
    442397                rnorm = sqrt(gamma);
    … …  
    454409                        _printf0_("L2-Norm g: "<< rnorm << "\n");
    455410                }
    456        
     411
    457412                /*Break prematurely if solver doesn't reach desired tolerance in reasonable time frame*/
    458413                if(count > 10./TOL)   
    … …  
    461416                 break;
    462417                }
    463        
    464418
    465419                /* ---------------------------------------------------------- */
    466 
    467420
    468421                /*Directional update*/
    … …  
    475428                /* Assign new to old iterates */
    476429                VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold);
    477                
     430
    478431                count++;
    479432        }
    … …  
    487440        *puf=out_uf;
    488441
    489 
    490442        /* Cleanup */
    491443        KSPDestroy(&kspu);KSPDestroy(&kspip);
    492444
    493445        MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP);
    494        
     446
    495447        VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp);
    496448        VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi);
    … …  
    522474        double rnorm1, rnorm2;
    523475
    524 
    525476        if(VerboseModule()) _printf0_("   checking convergence\n");
    526477
    … …  
    548499        MatGetSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT);
    549500        #endif
    550        
     501
    551502                /*no. of free nodes in velocity/pressure space*/
    552503                MatGetSize(B,&n_u,&n_p);
    … …  
    560511                VecGetSubVector(uf->pvector->vector,isv,&u);
    561512                VecGetSubVector(uf->pvector->vector,isp,&p);
    562                
    563513
    564514                /*Extract values of the RHS corresponding to the first/second block*/
    … …  
    573523                VecDuplicate(p,&res2);VecSet(res2,1.0);
    574524
    575 
    576525        /*Display solver caracteristics*/
    577526        if (VerboseConvergence()){
    578                
     527
    579528                /*Calculate res1 = A*u + B*p - f1*/
    580529                VecScale(f1,-1.);MatMultAdd(A,u,f1,tmp);MatMultAdd(B,p,tmp,res1);VecScale(f1,-1.);
    581530                /*Calculate res2 = B^T * u - f2*/
    582531                VecScale(f2,-1.);MatMultAdd(BT,u,f2,res2);VecScale(f2,-1.);
    583 
    584532
    585533                /*compute norm(res1), norm(res2), norm(F) and residue*/
    … …  
    597545        VecRestoreSubVector(uf->pvector->vector,isv,&u);
    598546        VecRestoreSubVector(uf->pvector->vector,isp,&p);
    599        
     547
    600548        /*Extract values corresponding to velocity/pressure on the old solution*/
    601549        VecGetSubVector(old_uf->pvector->vector,isv,&uold);
    602550        VecGetSubVector(old_uf->pvector->vector,isp,&pold);
    603                
    604551
    605552        /*Force equilibrium (Mandatory)*/
    … …  
    609556        /*Calculate res2 = B^T * uold - f2*/
    610557        VecScale(f2,-1.);MatMultAdd(BT,uold,f2,res2);VecScale(f2,-1.);
    611        
     558
    612559        /*compute norm(res1), norm(res2), norm(F) and residue*/
    613560        VecNorm(res1,NORM_2,&rnorm1);VecNorm(res2,NORM_2,&rnorm2);
    … …  
    626573        VecRestoreSubVector(old_uf->pvector->vector,isv,&uold);
    627574        VecRestoreSubVector(old_uf->pvector->vector,isp,&pold);
    628        
    629 
    630575
    631576        //print
    … …  
    721666        int  count=0;
    722667        bool converged=false;
    723        
     668
    724669        /*Start non-linear iteration using input velocity: */
    725670        GetSolutionFromInputsx(&ug,femmodel);
    … …  
    729674        InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    730675        InputUpdateFromSolutionx(femmodel,ug);
    731        
     676
    732677        for(;;){
    733678
    … …  
    747692                PetscOptionsGetInt(NULL,PETSC_NULL,"-schur_pc",&precond,NULL);
    748693                #endif
    749        
     694
    750695                StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); 
    751        
     696
    752697                /* Construct Schur preconditioner matrix or mass matrix depending on input */
    753698                if(precond)
    … …  
    760705                        }
    761706                }else{
    762                        
     707
    763708                        for(Object* & object : femmodel->elements->objects){
    764709                                Element* element2=xDynamicCast<Element*>(object);
    … …  
    769714                }
    770715
    771                        
    772716                delete analysis;
    773        
     717
    774718                /*Obtain index sets for velocity and pressure components */
    775719                IS isv = NULL;
    … …  
    785729                Kff->Assemble();
    786730
    787        
    788731                /*Solve*/
    789732                femmodel->profiler->Start(SOLVER);
    790733                _assert_(Kff->type==PetscMatType);
    791 
    792734
    793735                SchurCGSolver(&uf,
Note: See TracChangeset for help on using the changeset viewer.