Changeset 23204


Ignore:
Timestamp:
08/30/18 11:14:27 (7 years ago)
Author:
wester
Message:

CHG: Fixed bug in the SchurCG solver for nonzero Dirichlet B.C.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r23197 r23204  
    2222        Vec                                             tmpu, tmpp, rhsu,rhsp; /* temp. vectors, arbitrary RHS in vel. / pressure space */
    2323        Vec                                             gold,gnew,wold,wnew,chi,thetaold,thetanew,eta; /* CG intermediaries */
    24         Vec                                             f,BTAinvf;                              /* RHS of the global system */
     24        Vec                                             f1,f2;                                  /* RHS of the global system */
    2525        double                                  rho,gamma,tmpScalar; /* Step sizes, arbitrary double */
    2626        KSP                                             kspu,kspp;                              /* KSP contexts for vel. / pressure systems*/
    … …  
    3333
    3434        /*STOP tolerance for the rel. residual*/
    35         TOL = 0.1;
     35        TOL = 0.4;
    3636
    3737        /*Initialize output*/
    … …  
    7373
    7474        /* Set up intermediaries */
    75         VecDuplicate(uold,&f);VecSet(f,0.0);
    76         VecAssemblyBegin(f);VecAssemblyEnd(f);
    77 
    78         VecDuplicate(p,&BTAinvf);VecSet(BTAinvf,0.0);
    79         VecAssemblyBegin(BTAinvf);VecAssemblyEnd(BTAinvf);
    80 
     75        VecDuplicate(uold,&f1);VecSet(f1,0.0);
     76        VecDuplicate(p,&f2);VecSet(f2,0.0);
    8177        VecDuplicate(uold,&tmpu);VecSet(tmpu,0.0);
    82         VecAssemblyBegin(tmpu);VecAssemblyEnd(tmpu);
    83 
    8478        VecDuplicate(p,&tmpp);VecSet(tmpp,0.0);
    85         VecAssemblyBegin(tmpp);VecAssemblyEnd(tmpp);
    86 
    8779        VecDuplicate(p,&rhsp);VecSet(rhsp,0.0);
    88         VecAssemblyBegin(rhsp);VecAssemblyEnd(rhsp);
    89 
    9080        VecDuplicate(uold,&rhsu);VecSet(rhsu,0.0);
    91         VecAssemblyBegin(rhsu);VecAssemblyEnd(rhsu);
    92 
    9381        VecDuplicate(p,&gold);VecSet(gold,0.0);
    94         VecAssemblyBegin(gold);VecAssemblyEnd(gold);
    95 
    96         VecDuplicate(gold,&wnew);VecSet(wnew,0.0);
    97         VecAssemblyBegin(wnew);VecAssemblyEnd(wnew);
    98 
     82        VecDuplicate(p,&wnew);VecSet(wnew,0.0);
    9983        VecDuplicate(uold,&chi);VecSet(chi,0.0);
    100         VecAssemblyBegin(chi);VecAssemblyEnd(chi);
    101        
    10284        VecDuplicate(p,&thetanew);VecSet(thetanew,0.0);
    103         VecAssemblyBegin(thetanew);VecAssemblyEnd(thetanew);
    104 
    10585        VecDuplicate(p,&thetaold);VecSet(thetaold,0.0);
    106         VecAssemblyBegin(thetaold);VecAssemblyEnd(thetaold);
    107 
    10886        VecDuplicate(p,&eta);VecSet(eta,0.0);
    109         VecAssemblyBegin(eta);VecAssemblyEnd(eta);
    110        
    111         /* Get global RHS (restricted to the velocity space */
    112         VecGetSubVector(pf,isv,&f);
    113 
     87       
     88        /* Get global RHS (for each block sub-problem respectively)*/
     89        VecGetSubVector(pf,isv,&f1);
     90        VecGetSubVector(pf,isp,&f2);
    11491
    11592   /* ------------------------------------------------------------ */
    11693
    11794        /* Generate initial value for the velocity from the pressure */
    118         /* a(u0,v) = f(v)-b(p0,v)  i.e.  Au0 = F-Bp0 */
     95        /* a(u0,v) = f1(v)-b(p0,v)  i.e.  Au0 = F1-Bp0 */
    11996        /* u0 = u_DIR on \Gamma_DIR */
    12097       
    … …  
    124101        KSPSetType(kspu,KSPCG);
    125102        KSPSetInitialGuessNonzero(kspu,PETSC_TRUE);
     103        //KSPSetTolerances(kspu,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     104        //KSPMonitorSet(kspu,KSPMonitorDefault,NULL,NULL);
    126105        KSPGetPC(kspu,&pcu);
    127         PCSetType(pcu,PCJACOBI);
     106        PCSetType(pcu,PCSOR);
    128107        KSPSetUp(kspu);
    129108
    130109       
    131110        /* Create RHS */
    132         /* RHS = F-B * pold */
    133         VecScale(p,-1.);MatMultAdd(B,p,f,rhsu);VecScale(p,-1.);
    134 
    135         /* Go solve Au0 = F-Bp0*/
     111        /* RHS = F1-B * pold */
     112        VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.);
     113
     114        /* Go solve Au0 = F1-Bp0*/
    136115        KSPSolve(kspu,rhsu,uold);
    137116       
    138         /* Calculate B^T * A^{-1} * F for future reference */
    139         KSPSolve(kspu,f,tmpu);MatMult(BT,tmpu,BTAinvf);
    140 
    141117
    142118        /* Set up u_new */
    … …  
    149125
    150126        /*Get initial residual*/
    151         /*(1/mu(x) * g0, q) = b(q,u0)  i.e.  IP * g0 = BT * u0*/
     127        /*(1/mu(x) * g0, q) = b(q,u0) - (f2,q)  i.e.  IP * g0 = BT * u0 - F2*/
    152128       
    153129        /* Create KSP context */
    … …  
    156132       
    157133        /* Create RHS */
    158         /* RHS = BT * uold */
    159         MatMult(BT,uold,rhsp);
     134        /* RHS = BT * uold - F2 */
     135        VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.);
    160136
    161137        /* Set KSP & PC options */
    … …  
    187163
    188164        /*Realizing the step size part 1: thetam */
    189         /*IP * theta = BT * uold*/
    190         MatMult(BT,uold,rhsp);
     165        /*IP * theta = BT * uold - F2*/
     166        VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.);
    191167        KSPSolve(kspp,rhsp,thetaold);
    192168
    … …  
    212188
    213189                /*Set step size*/
    214                 /*rhom = [(wm)^T * IP^-1 * BT * um]/[(wm)^T * IP^-1 * BT * chim]*/
     190                /*rhom = [(wm)^T * IP^-1 * (BT * um - F2)]/[(wm)^T * IP^-1 * BT * chim]*/
    215191                VecDot(wold,thetaold,&rho);
    216192                VecDot(wold,eta,&tmpScalar);
    … …  
    234210
    235211                /*Theta update*/
    236                 /*IP * theta = BT * uold*/
    237                 MatMult(BT,unew,rhsp);
     212                /*IP * theta = BT * uold - F2*/
     213                VecScale(f2,-1.);MatMultAdd(BT,unew,f2,rhsp);VecScale(f2,-1.);
    238214                KSPSolve(kspp,rhsp,thetanew);
    239215
    … …  
    250226                /*BREAK if norm(g(m+0),2) < TOL or pressure space has been full searched*/
    251227                VecNorm(gnew,NORM_INFINITY,&rnorm);
    252                 if(rnorm < TOL*initRnorm) break;
     228                if(rnorm < TOL*initRnorm)
     229                 break;
     230                else if(rnorm > 100*initRnorm)
     231                 _error_("Solver diverged. This shouldn't happen\n");
     232                else
     233                 PetscPrintf(PETSC_COMM_WORLD,"rel. residual at step %d: %g, at TOL = %g\n",count,rnorm/initRnorm,TOL);
     234               
    253235                if(count > np-1) break;
    254 
     236       
    255237
    256238                /* ---------------------------------------------------------- */
    … …  
    272254        }
    273255
     256
    274257        /* Restore pressure and velocity sol. vectors to its global form */
    275258        VecRestoreSubVector(out_uf->pvector->vector,isv,&unew);
    … …  
    287270        VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp);
    288271        VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi);
    289         VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f);VecDestroy(&BTAinvf);VecDestroy(&eta);
     272        VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f1);VecDestroy(&f2);VecDestroy(&eta);
    290273        VecDestroy(&thetanew);VecDestroy(&thetaold);
    291274
    … …  
    342325                /*Create mass matrix*/
    343326                int fsize; Kff->GetSize(&fsize,&fsize);
    344                 Iff=new Matrix<IssmDouble>(fsize,fsize,200,4);
     327                Iff=new Matrix<IssmDouble>(fsize,fsize,300,4);
    345328                StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
    346329                /*Get complete stiffness matrix without penalties*/
    … …  
    367350                femmodel->profiler->Stop(SOLVER);
    368351                delete Iff;
    369 
    370 
     352               
    371353                /*Merge solution from f set to g set*/
    372354                Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
Note: See TracChangeset for help on using the changeset viewer.