![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 23909
- Timestamp:
- 05/07/19 00:50:30 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp
r23587 r23909 9 9 #include "../modules/modules.h" 10 10 #include "../analyses/analyses.h" 11 11 #include <iostream> 12 #include <fstream> 12 13 13 14 #ifdef _HAVE_PETSC_ 14 15 15 16 16 void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff, Vec pf, Vec uf0,IS isv,IS isp,Parameters* parameters){/*{{{*/17 void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff, Vec pf, Vec uf0,IS isv,IS isp,Parameters* parameters){/*{{{*/ 17 18 18 19 Mat A, B, BT; /* Saddle point block matrices */ 19 Mat IP; /* Preconditioner matrix */ 20 Mat IP2; 20 Mat IP; /* Preconditioner or mass matrix */ 21 21 int nu, np; /* No of. free nodes in velocity / pressure space */ 22 22 Vec p,uold,unew; /* Solution vectors for pressure / vel. */ 23 Vec tmpu, tmpp,rhsu,rhsp; /* temp. vectors, arbitrary RHS in vel. / pressure space */24 Vec gold,gnew,wold,wnew,chi ,thetaold,thetanew,eta; /* CG intermediaries */23 Vec tmpu,tmpu2,resu,resp,tmpp,tmpp2,rhsu,rhsp; /* temp. vectors, arbitrary RHS in vel. / pressure space */ 24 Vec gold,gnew,wold,wnew,chi; /* CG intermediaries */ 25 25 Vec f1,f2; /* RHS of the global system */ 26 double rho,gamma,tmpScalar ; /* Step sizes, arbitrary double */27 KSP kspu,ksp p;/* KSP contexts for vel. / pressure systems*/26 double rho,gamma,tmpScalar,tmpScalar2; /* Step sizes, arbitrary double */ 27 KSP kspu,kspip; /* KSP contexts for vel. / pressure systems*/ 28 28 KSPConvergedReason reason; /* Convergence reason for troubleshooting */ 29 29 int its; /* No. of iterations for troubleshooting */ 30 double initRnorm, rnorm, TOL ; /* residual norms, STOP tolerance */30 double initRnorm, rnorm, TOL,ELLTOL; /* residual norms, STOP tolerance */ 31 31 PC pcu,pcp; /* Preconditioner contexts pertaining the KSP contexts*/ 32 32 PetscViewer viewer; /* Viewer for troubleshooting */ 33 33 IssmPDouble t1,t2; /* Time measurement for bottleneck analysis */ 34 34 35 /*STOP tolerance for the rel. residual*/ 36 TOL = 0.4; 35 36 double tmp1,tmp2,tmp3; 37 int tmpi; 38 double tmp4,tmp5,tmp6,tmp7; 39 40 int noIt; 41 42 43 44 int precond = 0; 45 46 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) 47 PetscTruth flag,flg; 48 #else 49 PetscBool flag,flg; 50 #endif 51 52 53 char ksp_type[50]; 54 char pc_type[50]; 55 int maxiter; 56 57 #if _PETSC_MINOR_<7 58 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg); 59 PetscOptionsGetString(PETSC_NULL,"-pc_type",pc_type,49,&flg); 60 PetscOptionsGetReal(PETSC_NULL,"-tol",&TOL,NULL); 61 PetscOptionsGetReal(PETSC_NULL,"-elltol",&ELLTOL,NULL); 62 PetscOptionsGetInt(PETSC_NULL,"-schur_pc",&precond,NULL); 63 PetscOptionsGetInt(PETSC_NULL,"-max_iter",&maxiter,NULL); 64 #else 65 PetscOptionsGetString(NULL,PETSC_NULL,"-ksp_type",ksp_type,49,&flg); 66 PetscOptionsGetString(NULL,PETSC_NULL,"-pc_type",pc_type,49,&flg); 67 PetscOptionsGetReal(NULL,PETSC_NULL,"-tol",&TOL,NULL); 68 PetscOptionsGetReal(NULL,PETSC_NULL,"-elltol",&ELLTOL,NULL); 69 PetscOptionsGetInt(NULL,PETSC_NULL,"-schur_pc",&precond,NULL); 70 PetscOptionsGetInt(NULL,PETSC_NULL,"-max_iter",&maxiter,NULL); 71 #endif 72 73 74 if(precond){ 75 _printf0_("Running WITH preconditioner\n"); 76 }else{ 77 _printf0_("Running WITHOUT preconditioner\n"); 78 } 79 37 80 38 81 /*Initialize output*/ … … 41 84 /* Extract block matrices from the saddle point matrix */ 42 85 /* [ A B ] = Kff 43 * [ B^T 0 ] 86 * [ B^T I ] 87 * where A is the elliptic submatrix, B^T represents the incompressibility, 88 * and I the Schur preconditioner (stored here, because the space was allocated either way) 44 89 * */ 45 #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=8)46 MatCreateSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A);47 MatCreateSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B);48 MatCreateSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT);49 #else50 90 MatGetSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A); 51 91 MatGetSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B); 52 92 MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT); 53 #endif54 93 55 94 /* Extract preconditioner matrix on the pressure space*/ 56 #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=8)57 MatCreateSubMatrix(Kff,isp,isp,MAT_INITIAL_MATRIX,&IP);58 #else59 95 MatGetSubMatrix(Kff,isp,isp,MAT_INITIAL_MATRIX,&IP); 60 #endif 96 61 97 62 98 /* Get number of velocity / pressure nodes */ … … 77 113 VecDuplicate(p,&f2);VecSet(f2,0.0); 78 114 VecDuplicate(uold,&tmpu);VecSet(tmpu,0.0); 115 VecDuplicate(uold,&tmpu2);VecSet(tmpu2,0.0); 116 VecDuplicate(uold,&resu);VecSet(resu,0.0); 79 117 VecDuplicate(p,&tmpp);VecSet(tmpp,0.0); 118 VecDuplicate(p,&tmpp2);VecSet(tmpp2,0.0); 80 119 VecDuplicate(p,&rhsp);VecSet(rhsp,0.0); 120 VecDuplicate(p,&resp);VecSet(resp,0.0); 81 121 VecDuplicate(uold,&rhsu);VecSet(rhsu,0.0); 82 122 VecDuplicate(p,&gold);VecSet(gold,0.0); 83 123 VecDuplicate(p,&wnew);VecSet(wnew,0.0); 84 124 VecDuplicate(uold,&chi);VecSet(chi,0.0); 85 VecDuplicate(p,&thetanew);VecSet(thetanew,0.0);86 VecDuplicate(p,&thetaold);VecSet(thetaold,0.0);87 VecDuplicate(p,&eta);VecSet(eta,0.0);88 125 89 126 /* Get global RHS (for each block sub-problem respectively)*/ … … 104 141 KSPSetOperators(kspu,A,A,DIFFERENT_NONZERO_PATTERN); 105 142 #endif 106 KSPSetType(kspu,KSPCG); 143 if (strcmp(ksp_type,"gmres")==0){ 144 KSPSetType(kspu,KSPGMRES); 145 }else if(strcmp(ksp_type,"pipegmres")==0){ 146 KSPSetType(kspu,KSPPGMRES); 147 }else if(strcmp(ksp_type,"cg")==0){ 148 KSPSetType(kspu,KSPCG); 149 }else if(strcmp(ksp_type,"pipecg")==0){ 150 KSPSetType(kspu,KSPPIPECG); 151 }else if(strcmp(ksp_type,"bicg")==0){ 152 KSPSetType(kspu,KSPBICG); 153 }else if(strcmp(ksp_type,"bicgstab")==0){ 154 KSPSetType(kspu,KSPBCGS); 155 }else if(strcmp(ksp_type,"ibicgstab")==0){ 156 KSPSetType(kspu,KSPIBCGS); 157 }else if(strcmp(ksp_type,"minres")==0){ 158 KSPSetType(kspu,KSPMINRES); 159 }else if(strcmp(ksp_type,"cr")==0){ 160 KSPSetType(kspu,KSPCR); 161 }else if(strcmp(ksp_type,"pipecr")==0){ 162 KSPSetType(kspu,KSPPIPECR); 163 }else{ 164 _error_("Suggested KSP method not implemented yet!\n"); 165 } 166 107 167 KSPSetInitialGuessNonzero(kspu,PETSC_TRUE); 108 //KSPSetTolerances(kspu,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); 109 //KSPMonitorSet(kspu,KSPMonitorDefault,NULL,NULL); 168 169 /*Strong rel. residual tolerance needed when solving for the velocity update. 170 * This is because ISSM uses the dimensional equations, so the initial guess chi = 0 171 * results in a very high initial residual (~ 1e19)*/ 172 KSPSetTolerances(kspu,ELLTOL,PETSC_DEFAULT,PETSC_DEFAULT,maxiter); //maxiter 173 174 175 110 176 KSPGetPC(kspu,&pcu); 111 PCSetType(pcu,PCSOR); 177 if (strcmp(pc_type,"bjacobi")==0){ 178 PCSetType(pcu,PCBJACOBI); 179 }else if(strcmp(pc_type,"asm")==0){ 180 PCSetType(pcu,PCASM); 181 }else if(strcmp(pc_type,"gamg")==0){ 182 PCSetType(pcu,PCGAMG); 183 }else if(strcmp(pc_type,"gasm")==0){ 184 PCSetType(pcu,PCGASM); 185 }else if(strcmp(pc_type,"jacobi")==0){ 186 PCSetType(pcu,PCJACOBI); 187 }else if(strcmp(pc_type,"icc")==0){ 188 PCSetType(pcu,PCICC); 189 }else if(strcmp(pc_type,"ilu")==0){ 190 PCSetType(pcu,PCILU); 191 }else if(strcmp(pc_type,"sor")==0){ 192 PCSetType(pcu,PCSOR); 193 }else if(strcmp(pc_type,"eisenstat")==0){ 194 PCSetType(pcu,PCEISENSTAT); 195 }else if(strcmp(pc_type,"none")==0){ 196 PCSetType(pcu,PCNONE); 197 }else{ 198 _error_("Suggested preconditioner not implemented yet!\n"); 199 } 112 200 KSPSetUp(kspu); 113 201 … … 117 205 VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.); 118 206 207 if (VerboseConvergence()) 208 { 209 210 VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.); 211 VecNorm(tmpu,NORM_2,&tmp4); 212 213 VecScale(f2,-1.);MatMultAdd(BT,uold,f2,tmpp);VecScale(f2,-1.); 214 VecNorm(tmpp,NORM_2,&tmp6); 215 216 KSPInitialResidual(kspu,uold,tmpu,tmpu2,resu,rhsu); 217 VecNorm(resu,NORM_2,&tmp5); 218 219 220 } 221 222 119 223 /* Go solve Au0 = F1-Bp0*/ 120 224 KSPSolve(kspu,rhsu,uold); 121 225 KSPGetIterationNumber(kspu,&noIt); 226 227 228 229 if (VerboseConvergence()) 230 { 231 232 KSPGetIterationNumber(kspu,&tmpi); 233 VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.); 234 VecNorm(tmpu,NORM_2,&tmp2); 235 KSPGetResidualNorm(kspu,&tmp1); 236 237 _printf0_("||Au_00 - rhsu||_euc: " << tmp4 <<"\n||P(-1)(Au_00 - rhsu)||_euc: " << tmp5 << "\n ||Au_n0 - rhsu||_euc: " << tmp2<< "\n||P(-1)(Au_n0 - rhsu)||_euc: " << tmp1 << "\nIteration number: "<<tmpi<<"\n"); 238 _printf0_("||BTu_00 - f2||_euc: " << tmp6 <<"\n"); 239 } 240 241 122 242 123 243 /* Set up u_new */ … … 133 253 134 254 /* Create KSP context */ 135 KSPCreate(IssmComm::GetComm(),&ksp p);255 KSPCreate(IssmComm::GetComm(),&kspip); 136 256 #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=5) 137 KSPSetOperators(ksp p,IP,IP);257 KSPSetOperators(kspip,IP,IP); 138 258 #else 139 KSPSetOperators(ksp p,IP,IP,DIFFERENT_NONZERO_PATTERN);259 KSPSetOperators(kspip,IP,IP,DIFFERENT_NONZERO_PATTERN); 140 260 #endif 141 261 142 262 /* Create RHS */ 143 263 /* RHS = BT * uold - F2 */ 144 VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.); 264 VecScale(uold,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(uold,-1.); 265 145 266 146 267 /* Set KSP & PC options */ 147 KSPSetType(kspp,KSPCG); 148 KSPSetInitialGuessNonzero(kspp,PETSC_TRUE); 149 KSPGetPC(kspp,&pcp); 268 KSPSetType(kspip,KSPGMRES); 269 KSPSetInitialGuessNonzero(kspip,PETSC_TRUE); 270 271 272 KSPGetPC(kspip,&pcp); 150 273 PCSetType(pcp,PCJACOBI); 151 /* Note: Systems in the pressure space are cheap, so we can afford a better tolerance */ 152 KSPSetTolerances(kspp,1e-10,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); 153 KSPSetUp(kspp); 154 274 /* Note: Systems in the pressure space are cheap, so we can afford a good tolerance */ 275 KSPSetTolerances(kspip,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); 276 KSPSetUp(kspip); 277 278 if (VerboseConvergence()) 279 { 280 KSPInitialResidual(kspip,gold,tmpp,tmpp2,resp,rhsp); 281 VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.); 282 VecNorm(resp,NORM_2,&tmp5); 283 VecNorm(tmpp,NORM_2,&tmp4); 284 } 285 286 287 155 288 /* Go solve */ 156 KSPSolve(kspp,rhsp,gold); 157 289 KSPSolve(kspip,rhsp,gold); 290 291 292 293 if (VerboseConvergence()) 294 { 295 296 KSPGetResidualNorm(kspip,&tmp1); 297 VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.); 298 VecNorm(tmpp,NORM_2,&tmp2); 299 300 KSPGetIterationNumber(kspip,&tmpi); 301 _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"); 302 } 303 304 158 305 /*Initial residual*/ 159 VecNorm(gold,NORM_INFINITY,&initRnorm); 160 306 MatMult(IP,gold,tmpp);VecDot(gold,tmpp,&initRnorm); 307 initRnorm = sqrt(initRnorm); 308 _printf0_("inner product norm g0: "<<initRnorm<<"\n"); 309 310 311 /*Iterate only if inital residual large enough*/ 312 if(initRnorm > 1e-5) 313 { 314 161 315 /* Further setup */ 162 316 VecDuplicate(gold,&gnew);VecCopy(gold,gnew); … … 171 325 VecAssemblyBegin(wold);VecAssemblyEnd(wold); 172 326 173 /*Realizing the step size part 1: thetam */174 /*IP * theta = BT * uold - F2*/175 VecScale(f2,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(f2,-1.);176 KSPSolve(kspp,rhsp,thetaold);177 178 327 179 328 /* Count number of iterations */ … … 182 331 /* CG iteration*/ 183 332 for(;;){ 333 if(VerboseConvergence()) 334 { 335 _printf0_("\n###### Step " << count << " ######\n"); 336 } 184 337 185 338 /*Realizing the step size part 2: chim */ 186 /*a(chim,v) = -b(wm,v) i.e. A * chim = -B * wm */339 /*a(chim,v) = b(wm,v) i.e. A * chim = B * wm */ 187 340 /*chim_DIR = 0*/ 188 VecScale(wold,-1.);MatMult(B,wold,rhsu);VecScale(wold,-1.); 189 KSPSolve(kspu,rhsu,chi); 190 191 /*Realizing the step size part 3: etam */ 192 MatMult(BT,chi,rhsp); 193 KSPSolve(kspp,rhsp,eta); 194 341 VecScale(wold,1.);MatMult(B,wold,rhsu);VecScale(wold,1.); 342 VecSet(chi,0.0); 343 344 345 if(VerboseConvergence()) 346 { 347 VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.); 348 VecNorm(tmpu,NORM_2,&tmp4); 349 350 KSPInitialResidual(kspu,chi,tmpu,tmpu2,resu,rhsu); 351 VecNorm(resu,NORM_2,&tmp5); 352 353 } 354 355 356 KSPSolve(kspu,rhsu,chi); 357 358 359 360 361 if (VerboseConvergence()) 362 { 363 VecNorm(chi,NORM_2,&tmp1); 364 KSPGetResidualNorm(kspu,&tmp2); 365 366 VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.); 367 VecNorm(tmpu,NORM_2,&tmp3); 368 369 370 KSPGetIterationNumber(kspu,&tmpi); 371 _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"); 372 } 373 374 195 375 /* ---------------------------------------------------------- */ 196 376 … … 198 378 /*Set step size*/ 199 379 /*rhom = [(wm)^T * IP^-1 * (BT * um - F2)]/[(wm)^T * IP^-1 * BT * chim]*/ 200 VecDot(wold,thetaold,&rho); 201 VecDot(wold,eta,&tmpScalar); 380 MatMult(IP,gold,tmpp); 381 VecDot(tmpp,wold,&rho); 382 383 MatMult(BT,chi,tmpp); 384 VecDot(tmpp,wold,&tmpScalar); 202 385 rho = rho/tmpScalar; 203 386 … … 207 390 208 391 /*Pressure update*/ 209 /*p(m+1) = pm -rhom * wm*/392 /*p(m+1) = pm + rhom * wm*/ 210 393 VecAXPY(p,-1.*rho,wold); 211 394 … … 213 396 /*Velocity update*/ 214 397 /*u(m+1) = um - rhom * chim*/ 215 VecWAXPY(unew,-1.*rho,chi,uold); 216 217 218 /* ---------------------------------------------------------- */ 219 220 /*Theta update*/ 221 /*IP * theta = BT * uold - F2*/ 222 VecScale(f2,-1.);MatMultAdd(BT,unew,f2,rhsp);VecScale(f2,-1.); 223 KSPSolve(kspp,rhsp,thetanew); 398 VecWAXPY(unew,rho,chi,uold); 399 VecNorm(unew,NORM_2,&tmpScalar); 400 401 402 if (VerboseConvergence()) 403 { 404 VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.); 405 VecScale(rhsu,-1.);MatMultAdd(A,unew,rhsu,tmpu);VecScale(rhsu,-1.); 406 VecNorm(tmpu,NORM_2,&tmp6); 407 408 VecScale(f2,-1);MatMultAdd(BT,unew,f2,tmpp);VecScale(f2,-1); 409 VecNorm(tmpp,NORM_2,&tmp7); 410 _printf0_("Momentum balance norm: " << tmp6 <<"\n"); 411 _printf0_("Incompressibility norm: " << tmp7 <<"\n"); 412 } 413 224 414 225 415 … … 228 418 /*Residual update*/ 229 419 /*g(m+1) = gm - rhom * BT * chim*/ 230 VecWAXPY(gnew,-1.*rho,eta,gold); 420 MatMult(BT,chi,tmpp); 421 MatMult(IP,gold,tmpp2); 422 VecWAXPY(rhsp,-1.*rho,tmpp,tmpp2); 423 KSPSolve(kspip,rhsp,gnew); 424 425 231 426 232 427 /* ---------------------------------------------------------- */ 233 428 234 235 /*BREAK if norm(g(m+0),2) < TOL or pressure space has been full searched*/ 236 VecNorm(gnew,NORM_INFINITY,&rnorm); 429 MatMult(IP,gnew,tmpp); 430 431 VecDot(tmpp,gnew,&gamma); 432 rnorm = sqrt(gamma); 433 434 /*BREAK if norm(g(m+0),2) < TOL or pressure space has been fully searched*/ 237 435 if(rnorm < TOL*initRnorm) 436 { 437 break; 438 }else if(rnorm > 1000*initRnorm) 439 { 440 _printf0_("L2-Norm g: "<< rnorm << "\n"); 441 _printf0_("Solver diverged and ends prematurely.\n"); 238 442 break; 239 else if(rnorm > 100*initRnorm) 240 _error_("Solver diverged. This shouldn't happen\n"); 241 //else 242 // PetscPrintf(PETSC_COMM_WORLD,"rel. residual at step %d: %g, at TOL = %g\n",count,rnorm/initRnorm,TOL); 243 244 245 246 247 if(count > np-1) break; 443 }else{ 444 _printf0_("L2-Norm g: "<< rnorm << "\n"); 445 } 446 447 /*Break prematurely if solver doesn't reach desired tolerance in reasonable time frame*/ 448 if(count > 10./TOL) 449 { 450 _printf0_("Ended after " << ceil(5./TOL) << " CG iterations\n"); 451 break; 452 } 248 453 249 454 … … 252 457 253 458 /*Directional update*/ 254 /*gamma = [g(m+1)^T * theta(m+1)]/[g(m)^T * thetam]*/ 255 VecDot(gnew,thetanew,&gamma); 256 VecDot(gold,thetaold,&tmpScalar); 459 MatMult(IP,gold,tmpp); 460 VecDot(tmpp,gold,&tmpScalar); 257 461 gamma = gamma/tmpScalar; 258 462 259 /*w(m+1) = g(m+1) + gamma * w(m)*/260 463 VecWAXPY(wnew,gamma,wold,gnew); 261 464 262 465 /* Assign new to old iterates */ 263 VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold); VecCopy(thetanew,thetaold);466 VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold); 264 467 265 468 count++; 266 469 } 267 470 } 268 471 269 472 /* Restore pressure and velocity sol. vectors to its global form */ … … 276 479 277 480 /* Cleanup */ 278 KSPDestroy(&kspu);KSPDestroy(&ksp p);481 KSPDestroy(&kspu);KSPDestroy(&kspip); 279 482 280 483 MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP); … … 282 485 VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp); 283 486 VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi); 284 VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f1);VecDestroy(&f2); VecDestroy(&eta);285 VecDestroy(&t hetanew);VecDestroy(&thetaold);487 VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f1);VecDestroy(&f2); 488 VecDestroy(&tmpu2);VecDestroy(&tmpp2);VecDestroy(&resu);VecDestroy(&resp); 286 489 287 490 }/*}}}*/ … … 323 526 * Kff = | | 324 527 * [B^T IP] 325 528 /* To calculate the residual, only the necessary blocks need to be extracted */ 326 529 327 530 /*Extract A, B, B^T */ 328 #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=8)329 MatCreateSubMatrix(Kff->pmatrix->matrix,isv,isv,MAT_INITIAL_MATRIX,&A);330 MatCreateSubMatrix(Kff->pmatrix->matrix,isv,isp,MAT_INITIAL_MATRIX,&B);331 MatCreateSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT);332 #else333 531 MatGetSubMatrix(Kff->pmatrix->matrix,isv,isv,MAT_INITIAL_MATRIX,&A); 334 532 MatGetSubMatrix(Kff->pmatrix->matrix,isv,isp,MAT_INITIAL_MATRIX,&B); 335 533 MatGetSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT); 336 #endif337 534 338 535 /*no. of free nodes in velocity/pressure space*/ … … 493 690 494 691 692 495 693 /*parameters:*/ 496 694 int max_nonlinear_iterations; 695 int configuration_type; 696 int precond; 497 697 IssmDouble eps_res,eps_rel,eps_abs; 498 698 … … 502 702 femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum); 503 703 femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum); 704 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 504 705 femmodel->UpdateConstraintsx(); 505 706 int size; … … 507 708 bool converged=false; 508 709 710 509 711 /*Start non-linear iteration using input velocity: */ 510 712 GetSolutionFromInputsx(&ug,femmodel); … … 515 717 InputUpdateFromSolutionx(femmodel,ug); 516 718 719 517 720 for(;;){ 518 721 … … 521 724 delete ug; 522 725 726 727 523 728 /*Get stiffness matrix and Load vector*/ 524 729 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); … … 526 731 Reduceloadx(pf, Kfs, ys); delete Kfs; 527 732 528 /*Create mass matrix*/ 529 StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); 530 /*Get complete stiffness matrix without penalties*/ 531 for(int i=0;i<femmodel->elements->Size();i++){ 532 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 533 ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element); 534 if(Ie) Ie->AddToGlobal(Kff,NULL); 535 delete Ie; 536 } 537 Kff->Assemble(); 733 /*Create pressure matrix of choice*/ 734 PetscOptionsGetInt(PETSC_NULL,"-schur_pc",&precond,NULL); 735 736 StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); 737 738 /* Construct Schur preconditioner matrix or mass matrix depending on input */ 739 if(precond) 740 { 741 for(int i=0;i<femmodel->elements->Size();i++){ 742 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 743 ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element); 744 if(Ie) Ie->AddToGlobal(Kff,NULL); 745 delete Ie; 746 } 747 }else{ 748 749 for(int i=0;i<femmodel->elements->Size();i++){ 750 Element* element2=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 751 ElementMatrix* Ie2 = analysis->CreatePressureMassMatrix(element2); 752 if(Ie2) Ie2->AddToGlobal(Kff,NULL); 753 delete Ie2; 754 } 755 } 756 757 538 758 delete analysis; 539 759 540 760 /*Obtain index sets for velocity and pressure components */ 541 761 IS isv = NULL; 542 762 IS isp = NULL; 543 763 #if _PETSC_MAJOR_==3 544 /*Make indices out of doftypes: */ 545 if(!(df->pvector->vector))_error_("need doftypes for FS solver!\n"); 764 765 /*Make indices out of doftypes: */ 766 if(!(df->pvector->vector))_error_("need doftypes for FS solver!\n"); 546 767 DofTypesToIndexSet(&isv,&isp,df->pvector->vector,FSSolverEnum); 547 768 #else 548 769 _error_("Petsc 3.X required"); 549 770 #endif 550 551 771 Kff->Assemble(); 772 773 552 774 /*Solve*/ 553 775 femmodel->profiler->Start(SOLVER); 554 776 _assert_(Kff->type==PetscMatType); 555 777 778 556 779 SchurCGSolver(&uf, 557 780 Kff->pmatrix->matrix, … … 562 785 femmodel->parameters); 563 786 femmodel->profiler->Stop(SOLVER); 564 565 787 566 788 /*Merge solution from f set to g set*/ … … 575 797 converged=true; 576 798 } 799 577 800 InputUpdateFromConstantx(femmodel,converged,ConvergedEnum); 578 801 InputUpdateFromSolutionx(femmodel,ug);
Note:
See TracChangeset
for help on using the changeset viewer.
