![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 23197
- Timestamp:
- 08/28/18 14:10:43 (7 years ago)
- Location:
- issm/trunk-jpl/src/c/solutionsequences
- Files:
-
- 2 edited
-
solutionsequence_nonlinear.cpp (modified) (3 diffs)
-
solutionsequence_schurcg.cpp (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp
r22551 r23197 20 20 Vector<IssmDouble>* df = NULL; 21 21 Vector<IssmDouble>* ys = NULL; 22 23 Vec res; 24 double resNorm; 25 22 26 23 27 Loads* savedloads=NULL; … … 70 74 Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); 71 75 femmodel->profiler->Stop(SOLVER); 76 77 78 // VecDuplicate(uf->pvector->vector,&res); 79 // VecCopy(uf->pvector->vector,res); 80 // VecAssemblyBegin(res); 81 // VecAssemblyEnd(res); 82 // VecScale(pf->pvector->vector,-1.); 83 // MatMultAdd(Kff->pmatrix->matrix,uf->pvector->vector,pf->pvector->vector,res); 84 // VecScale(pf->pvector->vector,-1.); 85 // VecView(pf->pvector->vector,PETSC_VIEWER_STDOUT_WORLD); 86 // VecNorm(res,NORM_INFINITY,&resNorm); 87 // PetscPrintf(PETSC_COMM_WORLD,"Count = %d, Res. norm: %g\n", count,resNorm); 88 89 90 91 72 92 Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys; 73 93 … … 86 106 } 87 107 } 88 108 89 109 /*Increase count: */ 90 110 count++; -
issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp
r23162 r23197 10 10 #include "../analyses/analyses.h" 11 11 12 12 13 #ifdef _HAVE_PETSC_ 13 void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff,Mat Mff,Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/ 14 void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff,Mat Iff,Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/ 15 16 Mat A, B, BT; /* Saddle point block matrices */ 17 Mat IP; /* Preconditioner matrix */ 18 IS isv=NULL; /* Index set free velocity nodes */ 19 IS isp=NULL; /* Index set free pressure nodes */ 20 int nu, np; /* No of. free nodes in velocity / pressure space */ 21 Vec p,uold,unew; /* Solution vectors for pressure / vel. */ 22 Vec tmpu, tmpp, rhsu,rhsp; /* temp. vectors, arbitrary RHS in vel. / pressure space */ 23 Vec gold,gnew,wold,wnew,chi,thetaold,thetanew,eta; /* CG intermediaries */ 24 Vec f,BTAinvf; /* RHS of the global system */ 25 double rho,gamma,tmpScalar; /* Step sizes, arbitrary double */ 26 KSP kspu,kspp; /* KSP contexts for vel. / pressure systems*/ 27 KSPConvergedReason reason; /* Convergence reason for troubleshooting */ 28 int its; /* No. of iterations for troubleshooting */ 29 double initRnorm, rnorm, TOL; /* residual norms, STOP tolerance */ 30 PC pcu,pcp; /* Preconditioner contexts pertaining the KSP contexts*/ 31 PetscViewer viewer; /* Viewer for troubleshooting */ 32 IssmPDouble t1,t2; /* Time measurement for bottleneck analysis */ 33 34 /*STOP tolerance for the rel. residual*/ 35 TOL = 0.1; 14 36 15 37 /*Initialize output*/ 16 Vec uf = NULL; 17 18 _error_("not implemented yet"); 38 Vector<IssmDouble>* out_uf=new Vector<IssmDouble>(uf0); 39 40 /* Get velocity and pressure index sets for extraction */ 41 #if _PETSC_MAJOR_==3 42 /*Make indices out of doftypes: */ 43 if(!df)_error_("need doftypes for FS solver!\n"); 44 DofTypesToIndexSet(&isv,&isp,df,FSSolverEnum); 45 #else 46 _error_("Petsc 3.X required"); 47 #endif 48 49 50 /* Extract block matrices from the saddle point matrix */ 51 /* [ A B ] = Kff 52 * [ B^T 0 ] 53 * */ 54 MatGetSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A); 55 MatGetSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B); 56 MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT); 57 58 /* Extract preconditioner matrix on the pressure space*/ 59 MatGetSubMatrix(Iff,isp,isp,MAT_INITIAL_MATRIX,&IP); 60 61 /* Get number of velocity / pressure nodes */ 62 MatGetSize(B,&nu,&np); 63 64 /* Extract initial guesses for uold and pold */ 65 VecCreate(IssmComm::GetComm(),&p);VecSetSizes(p,PETSC_DECIDE,np);VecSetFromOptions(p); 66 VecAssemblyBegin(p);VecAssemblyEnd(p); 67 VecCreate(IssmComm::GetComm(),&uold);VecSetSizes(uold,PETSC_DECIDE,nu);VecSetFromOptions(uold); 68 VecAssemblyBegin(uold);VecAssemblyEnd(uold); 69 70 VecGetSubVector(out_uf->pvector->vector,isv,&uold); 71 VecGetSubVector(out_uf->pvector->vector,isp,&p); 72 73 74 /* 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 81 VecDuplicate(uold,&tmpu);VecSet(tmpu,0.0); 82 VecAssemblyBegin(tmpu);VecAssemblyEnd(tmpu); 83 84 VecDuplicate(p,&tmpp);VecSet(tmpp,0.0); 85 VecAssemblyBegin(tmpp);VecAssemblyEnd(tmpp); 86 87 VecDuplicate(p,&rhsp);VecSet(rhsp,0.0); 88 VecAssemblyBegin(rhsp);VecAssemblyEnd(rhsp); 89 90 VecDuplicate(uold,&rhsu);VecSet(rhsu,0.0); 91 VecAssemblyBegin(rhsu);VecAssemblyEnd(rhsu); 92 93 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 99 VecDuplicate(uold,&chi);VecSet(chi,0.0); 100 VecAssemblyBegin(chi);VecAssemblyEnd(chi); 101 102 VecDuplicate(p,&thetanew);VecSet(thetanew,0.0); 103 VecAssemblyBegin(thetanew);VecAssemblyEnd(thetanew); 104 105 VecDuplicate(p,&thetaold);VecSet(thetaold,0.0); 106 VecAssemblyBegin(thetaold);VecAssemblyEnd(thetaold); 107 108 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 114 115 /* ------------------------------------------------------------ */ 116 117 /* Generate initial value for the velocity from the pressure */ 118 /* a(u0,v) = f(v)-b(p0,v) i.e. Au0 = F-Bp0 */ 119 /* u0 = u_DIR on \Gamma_DIR */ 120 121 /* Create KSP context */ 122 KSPCreate(IssmComm::GetComm(),&kspu); 123 KSPSetOperators(kspu,A,A); 124 KSPSetType(kspu,KSPCG); 125 KSPSetInitialGuessNonzero(kspu,PETSC_TRUE); 126 KSPGetPC(kspu,&pcu); 127 PCSetType(pcu,PCJACOBI); 128 KSPSetUp(kspu); 129 130 131 /* 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*/ 136 KSPSolve(kspu,rhsu,uold); 137 138 /* Calculate B^T * A^{-1} * F for future reference */ 139 KSPSolve(kspu,f,tmpu);MatMult(BT,tmpu,BTAinvf); 140 141 142 /* Set up u_new */ 143 VecDuplicate(uold,&unew);VecCopy(uold,unew); 144 VecAssemblyBegin(unew);VecAssemblyEnd(unew); 145 146 147 148 /* ------------------------------------------------------------- */ 149 150 /*Get initial residual*/ 151 /*(1/mu(x) * g0, q) = b(q,u0) i.e. IP * g0 = BT * u0*/ 152 153 /* Create KSP context */ 154 KSPCreate(IssmComm::GetComm(),&kspp); 155 KSPSetOperators(kspp,IP,IP); 156 157 /* Create RHS */ 158 /* RHS = BT * uold */ 159 MatMult(BT,uold,rhsp); 160 161 /* Set KSP & PC options */ 162 KSPSetType(kspp,KSPCG); 163 KSPSetInitialGuessNonzero(kspp,PETSC_TRUE); 164 KSPGetPC(kspp,&pcp); 165 PCSetType(pcp,PCJACOBI); 166 /* Note: Systems in the pressure space are cheap, so we can afford a better tolerance */ 167 KSPSetTolerances(kspp,1e-10,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); 168 KSPSetUp(kspp); 169 170 /* Go solve */ 171 KSPSolve(kspp,rhsp,gold); 172 173 /*Initial residual*/ 174 VecNorm(gold,NORM_INFINITY,&initRnorm); 175 176 /* Further setup */ 177 VecDuplicate(gold,&gnew);VecCopy(gold,gnew); 178 VecAssemblyBegin(gnew);VecAssemblyEnd(gnew); 179 180 181 /* ------------------------------------------------------------ */ 182 183 /*Set initial search direction*/ 184 /*w0 = g0*/ 185 VecDuplicate(gold,&wold);VecCopy(gold,wold); 186 VecAssemblyBegin(wold);VecAssemblyEnd(wold); 187 188 /*Realizing the step size part 1: thetam */ 189 /*IP * theta = BT * uold*/ 190 MatMult(BT,uold,rhsp); 191 KSPSolve(kspp,rhsp,thetaold); 192 193 194 /* Count number of iterations */ 195 int count = 0; 196 197 /* CG iteration*/ 198 for(;;){ 199 200 /*Realizing the step size part 2: chim */ 201 /*a(chim,v) = -b(wm,v) i.e. A * chim = -B * wm */ 202 /*chim_DIR = 0*/ 203 VecScale(wold,-1.);MatMult(B,wold,rhsu);VecScale(wold,-1.); 204 KSPSolve(kspu,rhsu,chi); 205 206 /*Realizing the step size part 3: etam */ 207 MatMult(BT,chi,rhsp); 208 KSPSolve(kspp,rhsp,eta); 209 210 /* ---------------------------------------------------------- */ 211 212 213 /*Set step size*/ 214 /*rhom = [(wm)^T * IP^-1 * BT * um]/[(wm)^T * IP^-1 * BT * chim]*/ 215 VecDot(wold,thetaold,&rho); 216 VecDot(wold,eta,&tmpScalar); 217 rho = rho/tmpScalar; 218 219 220 /* ---------------------------------------------------------- */ 221 222 223 /*Pressure update*/ 224 /*p(m+1) = pm - rhom * wm*/ 225 VecAXPY(p,-1.*rho,wold); 226 227 228 /*Velocity update*/ 229 /*u(m+1) = um - rhom * chim*/ 230 VecWAXPY(unew,-1.*rho,chi,uold); 231 232 233 /* ---------------------------------------------------------- */ 234 235 /*Theta update*/ 236 /*IP * theta = BT * uold*/ 237 MatMult(BT,unew,rhsp); 238 KSPSolve(kspp,rhsp,thetanew); 239 240 241 /* ---------------------------------------------------------- */ 242 243 /*Residual update*/ 244 /*g(m+1) = gm - rhom * BT * chim*/ 245 VecWAXPY(gnew,-1.*rho,eta,gold); 246 247 /* ---------------------------------------------------------- */ 248 249 250 /*BREAK if norm(g(m+0),2) < TOL or pressure space has been full searched*/ 251 VecNorm(gnew,NORM_INFINITY,&rnorm); 252 if(rnorm < TOL*initRnorm) break; 253 if(count > np-1) break; 254 255 256 /* ---------------------------------------------------------- */ 257 258 259 /*Directional update*/ 260 /*gamma = [g(m+1)^T * theta(m+1)]/[g(m)^T * thetam]*/ 261 VecDot(gnew,thetanew,&gamma); 262 VecDot(gold,thetaold,&tmpScalar); 263 gamma = gamma/tmpScalar; 264 265 /*w(m+1) = g(m+1) + gamma * w(m)*/ 266 VecWAXPY(wnew,gamma,wold,gnew); 267 268 /* Assign new to old iterates */ 269 VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold);VecCopy(thetanew,thetaold); 270 271 count++; 272 } 273 274 /* Restore pressure and velocity sol. vectors to its global form */ 275 VecRestoreSubVector(out_uf->pvector->vector,isv,&unew); 276 VecRestoreSubVector(out_uf->pvector->vector,isp,&p); 19 277 20 278 /*return output pointer*/ 21 Vector<IssmDouble>* out_uf=new Vector<IssmDouble>(uf);22 VecFree(&uf);23 279 *puf=out_uf; 280 281 282 /* Cleanup */ 283 KSPDestroy(&kspu);KSPDestroy(&kspp); 284 285 MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP); 286 287 VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp); 288 VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi); 289 VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f);VecDestroy(&BTAinvf);VecDestroy(&eta); 290 VecDestroy(&thetanew);VecDestroy(&thetaold); 291 24 292 }/*}}}*/ 25 293 void solutionsequence_schurcg(FemModel* femmodel){/*{{{*/ … … 32 300 Vector<IssmDouble>* old_uf = NULL; 33 301 Vector<IssmDouble>* pf = NULL; 34 Vector<IssmDouble>* pf0 = NULL;35 302 Vector<IssmDouble>* df = NULL; 36 303 Vector<IssmDouble>* ys = NULL; 37 Matrix<IssmDouble>* Mff = NULL; 304 Matrix<IssmDouble>* Iff = NULL; 305 38 306 39 307 /*parameters:*/ … … 49 317 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 50 318 femmodel->UpdateConstraintsx(); 51 319 int size; 52 320 int count=0; 53 321 bool converged=false; … … 69 337 /*Get stiffness matrix and Load vector*/ 70 338 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); 71 pf0=pf->Duplicate(); pf->Copy(pf0);72 339 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); 73 340 Reduceloadx(pf, Kfs, ys); delete Kfs; … … 75 342 /*Create mass matrix*/ 76 343 int fsize; Kff->GetSize(&fsize,&fsize); 77 Mff=new Matrix<IssmDouble>(fsize,fsize,100,4);344 Iff=new Matrix<IssmDouble>(fsize,fsize,200,4); 78 345 StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); 79 346 /*Get complete stiffness matrix without penalties*/ 80 347 for(int i=0;i<femmodel->elements->Size();i++){ 81 348 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 82 ElementMatrix* Me = analysis->CreatePressureMassMatrix(element);83 if( Me) Me->AddToGlobal(Mff,NULL);84 delete Me;349 ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element); 350 if(Ie) Ie->AddToGlobal(Iff,NULL); 351 delete Ie; 85 352 } 86 Mff->Assemble();353 Iff->Assemble(); 87 354 delete analysis; 88 355 … … 90 357 femmodel->profiler->Start(SOLVER); 91 358 _assert_(Kff->type==PetscMatType); 359 92 360 SchurCGSolver(&uf, 93 361 Kff->pmatrix->matrix, 94 Mff->pmatrix->matrix,362 Iff->pmatrix->matrix, 95 363 pf->pvector->vector, 96 364 old_uf->pvector->vector, … … 98 366 femmodel->parameters); 99 367 femmodel->profiler->Stop(SOLVER); 100 delete pf0; delete Mff; 368 delete Iff; 369 101 370 102 371 /*Merge solution from f set to g set*/ … … 106 375 convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df; 107 376 count++; 377 108 378 if(count>=max_nonlinear_iterations){ 109 379 _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
Note:
See TracChangeset
for help on using the changeset viewer.
