Ice Sheet System Model  4.18
Code documentation
solutionsequence_schurcg.cpp
Go to the documentation of this file.
1 
5 #include "./solutionsequences.h"
6 #include "../toolkits/toolkits.h"
7 #include "../classes/classes.h"
8 #include "../shared/shared.h"
9 #include "../modules/modules.h"
10 #include "../analyses/analyses.h"
11 #include <iostream>
12 #include <fstream>
13 
14 #ifdef _HAVE_PETSC_
15 
16 
17 void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff, Vec pf, Vec uf0,IS isv,IS isp,Parameters* parameters){/*{{{*/
18 
19  Mat A, B, BT; /* Saddle point block matrices */
20  Mat IP; /* Preconditioner or mass matrix */
21  int nu, np; /* No of. free nodes in velocity / pressure space */
22  Vec p,uold,unew; /* Solution vectors for pressure / vel. */
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  Vec f1,f2; /* RHS of the global system */
26  double rho,gamma,tmpScalar,tmpScalar2; /* Step sizes, arbitrary double */
27  KSP kspu,kspip; /* KSP contexts for vel. / pressure systems*/
28  KSPConvergedReason reason; /* Convergence reason for troubleshooting */
29  int its; /* No. of iterations for troubleshooting */
30  double initRnorm, rnorm, TOL,ELLTOL; /* residual norms, STOP tolerance */
31  PC pcu,pcp; /* Preconditioner contexts pertaining the KSP contexts*/
32  PetscViewer viewer; /* Viewer for troubleshooting */
33  IssmPDouble t1,t2; /* Time measurement for bottleneck analysis */
34 
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 
80 
81  /*Initialize output*/
82  Vector<IssmDouble>* out_uf=new Vector<IssmDouble>(uf0);
83 
84  /* Extract block matrices from the saddle point matrix */
85  /* [ A B ] = Kff
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)
89  * */
90  #if _PETSC_MINOR_>8
91  MatCreateSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A);
92  MatCreateSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B);
93  MatCreateSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT);
94  #else
95  MatGetSubMatrix(Kff,isv,isv,MAT_INITIAL_MATRIX,&A);
96  MatGetSubMatrix(Kff,isv,isp,MAT_INITIAL_MATRIX,&B);
97  MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT);
98  #endif
99 
100  /* Extract preconditioner matrix on the pressure space*/
101  #if _PETSC_MINOR_>8
102  MatCreateSubMatrix(Kff,isp,isp,MAT_INITIAL_MATRIX,&IP);
103  #else
104  MatGetSubMatrix(Kff,isp,isp,MAT_INITIAL_MATRIX,&IP);
105  #endif
106 
107 
108  /* Get number of velocity / pressure nodes */
109  MatGetSize(B,&nu,&np);
110 
111  /* Extract initial guesses for uold and pold */
112  VecCreate(IssmComm::GetComm(),&p);VecSetSizes(p,PETSC_DECIDE,np);VecSetFromOptions(p);
113  VecAssemblyBegin(p);VecAssemblyEnd(p);
114  VecCreate(IssmComm::GetComm(),&uold);VecSetSizes(uold,PETSC_DECIDE,nu);VecSetFromOptions(uold);
115  VecAssemblyBegin(uold);VecAssemblyEnd(uold);
116 
117  VecGetSubVector(out_uf->pvector->vector,isv,&uold);
118  VecGetSubVector(out_uf->pvector->vector,isp,&p);
119 
120 
121  /* Set up intermediaries */
122  VecDuplicate(uold,&f1);VecSet(f1,0.0);
123  VecDuplicate(p,&f2);VecSet(f2,0.0);
124  VecDuplicate(uold,&tmpu);VecSet(tmpu,0.0);
125  VecDuplicate(uold,&tmpu2);VecSet(tmpu2,0.0);
126  VecDuplicate(uold,&resu);VecSet(resu,0.0);
127  VecDuplicate(p,&tmpp);VecSet(tmpp,0.0);
128  VecDuplicate(p,&tmpp2);VecSet(tmpp2,0.0);
129  VecDuplicate(p,&rhsp);VecSet(rhsp,0.0);
130  VecDuplicate(p,&resp);VecSet(resp,0.0);
131  VecDuplicate(uold,&rhsu);VecSet(rhsu,0.0);
132  VecDuplicate(p,&gold);VecSet(gold,0.0);
133  VecDuplicate(p,&wnew);VecSet(wnew,0.0);
134  VecDuplicate(uold,&chi);VecSet(chi,0.0);
135 
136  /* Get global RHS (for each block sub-problem respectively)*/
137  VecGetSubVector(pf,isv,&f1);
138  VecGetSubVector(pf,isp,&f2);
139 
140  /* ------------------------------------------------------------ */
141 
142  /* Generate initial value for the velocity from the pressure */
143  /* a(u0,v) = f1(v)-b(p0,v) i.e. Au0 = F1-Bp0 */
144  /* u0 = u_DIR on \Gamma_DIR */
145 
146  /* Create KSP context */
147  KSPCreate(IssmComm::GetComm(),&kspu);
148  #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=5)
149  KSPSetOperators(kspu,A,A);
150  #else
151  KSPSetOperators(kspu,A,A,DIFFERENT_NONZERO_PATTERN);
152  #endif
153  if (strcmp(ksp_type,"gmres")==0){
154  KSPSetType(kspu,KSPGMRES);
155  }else if(strcmp(ksp_type,"pipegmres")==0){
156  KSPSetType(kspu,KSPPGMRES);
157  }else if(strcmp(ksp_type,"cg")==0){
158  KSPSetType(kspu,KSPCG);
159  }else if(strcmp(ksp_type,"pipecg")==0){
160  KSPSetType(kspu,KSPPIPECG);
161  }else if(strcmp(ksp_type,"bicg")==0){
162  KSPSetType(kspu,KSPBICG);
163  }else if(strcmp(ksp_type,"bicgstab")==0){
164  KSPSetType(kspu,KSPBCGS);
165  }else if(strcmp(ksp_type,"ibicgstab")==0){
166  KSPSetType(kspu,KSPIBCGS);
167  }else if(strcmp(ksp_type,"minres")==0){
168  KSPSetType(kspu,KSPMINRES);
169  }else if(strcmp(ksp_type,"cr")==0){
170  KSPSetType(kspu,KSPCR);
171  }else if(strcmp(ksp_type,"pipecr")==0){
172  KSPSetType(kspu,KSPPIPECR);
173  }else{
174  _error_("Suggested KSP method not implemented yet!\n");
175  }
176 
177  KSPSetInitialGuessNonzero(kspu,PETSC_TRUE);
178 
179  /*Strong rel. residual tolerance needed when solving for the velocity update.
180  * This is because ISSM uses the dimensional equations, so the initial guess chi = 0
181  * results in a very high initial residual (~ 1e19)*/
182  KSPSetTolerances(kspu,ELLTOL,PETSC_DEFAULT,PETSC_DEFAULT,maxiter); //maxiter
183 
184 
185 
186  KSPGetPC(kspu,&pcu);
187  if (strcmp(pc_type,"bjacobi")==0){
188  PCSetType(pcu,PCBJACOBI);
189  }else if(strcmp(pc_type,"asm")==0){
190  PCSetType(pcu,PCASM);
191  }else if(strcmp(pc_type,"gamg")==0){
192  PCSetType(pcu,PCGAMG);
193  }else if(strcmp(pc_type,"gasm")==0){
194  PCSetType(pcu,PCGASM);
195  }else if(strcmp(pc_type,"jacobi")==0){
196  PCSetType(pcu,PCJACOBI);
197  }else if(strcmp(pc_type,"icc")==0){
198  PCSetType(pcu,PCICC);
199  }else if(strcmp(pc_type,"ilu")==0){
200  PCSetType(pcu,PCILU);
201  }else if(strcmp(pc_type,"sor")==0){
202  PCSetType(pcu,PCSOR);
203  }else if(strcmp(pc_type,"eisenstat")==0){
204  PCSetType(pcu,PCEISENSTAT);
205  }else if(strcmp(pc_type,"none")==0){
206  PCSetType(pcu,PCNONE);
207  }else{
208  _error_("Suggested preconditioner not implemented yet!\n");
209  }
210  KSPSetUp(kspu);
211 
212 
213  /* Create RHS */
214  /* RHS = F1-B * pold */
215  VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.);
216 
217  if (VerboseConvergence())
218  {
219 
220  VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.);
221  VecNorm(tmpu,NORM_2,&tmp4);
222 
223  VecScale(f2,-1.);MatMultAdd(BT,uold,f2,tmpp);VecScale(f2,-1.);
224  VecNorm(tmpp,NORM_2,&tmp6);
225 
226  KSPInitialResidual(kspu,uold,tmpu,tmpu2,resu,rhsu);
227  VecNorm(resu,NORM_2,&tmp5);
228 
229 
230  }
231 
232 
233  /* Go solve Au0 = F1-Bp0*/
234  KSPSolve(kspu,rhsu,uold);
235  KSPGetIterationNumber(kspu,&noIt);
236 
237 
238 
239  if (VerboseConvergence())
240  {
241 
242  KSPGetIterationNumber(kspu,&tmpi);
243  VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.);
244  VecNorm(tmpu,NORM_2,&tmp2);
245  KSPGetResidualNorm(kspu,&tmp1);
246 
247  _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");
248  _printf0_("||BTu_00 - f2||_euc: " << tmp6 <<"\n");
249  }
250 
251 
252 
253  /* Set up u_new */
254  VecDuplicate(uold,&unew);VecCopy(uold,unew);
255  VecAssemblyBegin(unew);VecAssemblyEnd(unew);
256 
257 
258 
259  /* ------------------------------------------------------------- */
260 
261  /*Get initial residual*/
262  /*(1/mu(x) * g0, q) = b(q,u0) - (f2,q) i.e. IP * g0 = BT * u0 - F2*/
263 
264  /* Create KSP context */
265  KSPCreate(IssmComm::GetComm(),&kspip);
266  #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=5)
267  KSPSetOperators(kspip,IP,IP);
268  #else
269  KSPSetOperators(kspip,IP,IP,DIFFERENT_NONZERO_PATTERN);
270  #endif
271 
272  /* Create RHS */
273  /* RHS = BT * uold - F2 */
274  VecScale(uold,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(uold,-1.);
275 
276 
277  /* Set KSP & PC options */
278  KSPSetType(kspip,KSPGMRES);
279  KSPSetInitialGuessNonzero(kspip,PETSC_TRUE);
280 
281 
282  KSPGetPC(kspip,&pcp);
283  PCSetType(pcp,PCJACOBI);
284  /* Note: Systems in the pressure space are cheap, so we can afford a good tolerance */
285  KSPSetTolerances(kspip,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
286  KSPSetUp(kspip);
287 
288  if (VerboseConvergence())
289  {
290  KSPInitialResidual(kspip,gold,tmpp,tmpp2,resp,rhsp);
291  VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.);
292  VecNorm(resp,NORM_2,&tmp5);
293  VecNorm(tmpp,NORM_2,&tmp4);
294  }
295 
296 
297 
298  /* Go solve */
299  KSPSolve(kspip,rhsp,gold);
300 
301 
302 
303  if (VerboseConvergence())
304  {
305 
306  KSPGetResidualNorm(kspip,&tmp1);
307  VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.);
308  VecNorm(tmpp,NORM_2,&tmp2);
309 
310  KSPGetIterationNumber(kspip,&tmpi);
311  _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");
312  }
313 
314 
315  /*Initial residual*/
316  MatMult(IP,gold,tmpp);VecDot(gold,tmpp,&initRnorm);
317  initRnorm = sqrt(initRnorm);
318  _printf0_("inner product norm g0: "<<initRnorm<<"\n");
319 
320 
321  /*Iterate only if inital residual large enough*/
322  if(initRnorm > 1e-5)
323  {
324 
325  /* Further setup */
326  VecDuplicate(gold,&gnew);VecCopy(gold,gnew);
327  VecAssemblyBegin(gnew);VecAssemblyEnd(gnew);
328 
329 
330  /* ------------------------------------------------------------ */
331 
332  /*Set initial search direction*/
333  /*w0 = g0*/
334  VecDuplicate(gold,&wold);VecCopy(gold,wold);
335  VecAssemblyBegin(wold);VecAssemblyEnd(wold);
336 
337 
338  /* Count number of iterations */
339  int count = 0;
340 
341  /* CG iteration*/
342  for(;;){
343  if(VerboseConvergence())
344  {
345  _printf0_("\n###### Step " << count << " ######\n");
346  }
347 
348  /*Realizing the step size part 2: chim */
349  /*a(chim,v) = b(wm,v) i.e. A * chim = B * wm */
350  /*chim_DIR = 0*/
351  VecScale(wold,1.);MatMult(B,wold,rhsu);VecScale(wold,1.);
352  VecSet(chi,0.0);
353 
354 
355  if(VerboseConvergence())
356  {
357  VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
358  VecNorm(tmpu,NORM_2,&tmp4);
359 
360  KSPInitialResidual(kspu,chi,tmpu,tmpu2,resu,rhsu);
361  VecNorm(resu,NORM_2,&tmp5);
362 
363  }
364 
365 
366  KSPSolve(kspu,rhsu,chi);
367 
368 
369 
370 
371  if (VerboseConvergence())
372  {
373  VecNorm(chi,NORM_2,&tmp1);
374  KSPGetResidualNorm(kspu,&tmp2);
375 
376  VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
377  VecNorm(tmpu,NORM_2,&tmp3);
378 
379 
380  KSPGetIterationNumber(kspu,&tmpi);
381  _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");
382  }
383 
384 
385  /* ---------------------------------------------------------- */
386 
387 
388  /*Set step size*/
389  /*rhom = [(wm)^T * IP^-1 * (BT * um - F2)]/[(wm)^T * IP^-1 * BT * chim]*/
390  MatMult(IP,gold,tmpp);
391  VecDot(tmpp,wold,&rho);
392 
393  MatMult(BT,chi,tmpp);
394  VecDot(tmpp,wold,&tmpScalar);
395  rho = rho/tmpScalar;
396 
397 
398  /* ---------------------------------------------------------- */
399 
400 
401  /*Pressure update*/
402  /*p(m+1) = pm + rhom * wm*/
403  VecAXPY(p,-1.*rho,wold);
404 
405 
406  /*Velocity update*/
407  /*u(m+1) = um - rhom * chim*/
408  VecWAXPY(unew,rho,chi,uold);
409  VecNorm(unew,NORM_2,&tmpScalar);
410 
411 
412  if (VerboseConvergence())
413  {
414  VecScale(p,-1.);MatMultAdd(B,p,f1,rhsu);VecScale(p,-1.);
415  VecScale(rhsu,-1.);MatMultAdd(A,unew,rhsu,tmpu);VecScale(rhsu,-1.);
416  VecNorm(tmpu,NORM_2,&tmp6);
417 
418  VecScale(f2,-1);MatMultAdd(BT,unew,f2,tmpp);VecScale(f2,-1);
419  VecNorm(tmpp,NORM_2,&tmp7);
420  _printf0_("Momentum balance norm: " << tmp6 <<"\n");
421  _printf0_("Incompressibility norm: " << tmp7 <<"\n");
422  }
423 
424 
425 
426  /* ---------------------------------------------------------- */
427 
428  /*Residual update*/
429  /*g(m+1) = gm - rhom * BT * chim*/
430  MatMult(BT,chi,tmpp);
431  MatMult(IP,gold,tmpp2);
432  VecWAXPY(rhsp,-1.*rho,tmpp,tmpp2);
433  KSPSolve(kspip,rhsp,gnew);
434 
435 
436 
437  /* ---------------------------------------------------------- */
438 
439  MatMult(IP,gnew,tmpp);
440 
441  VecDot(tmpp,gnew,&gamma);
442  rnorm = sqrt(gamma);
443 
444  /*BREAK if norm(g(m+0),2) < TOL or pressure space has been fully searched*/
445  if(rnorm < TOL*initRnorm)
446  {
447  break;
448  }else if(rnorm > 1000*initRnorm)
449  {
450  _printf0_("L2-Norm g: "<< rnorm << "\n");
451  _printf0_("Solver diverged and ends prematurely.\n");
452  break;
453  }else{
454  _printf0_("L2-Norm g: "<< rnorm << "\n");
455  }
456 
457  /*Break prematurely if solver doesn't reach desired tolerance in reasonable time frame*/
458  if(count > 10./TOL)
459  {
460  _printf0_("Ended after " << ceil(5./TOL) << " CG iterations\n");
461  break;
462  }
463 
464 
465  /* ---------------------------------------------------------- */
466 
467 
468  /*Directional update*/
469  MatMult(IP,gold,tmpp);
470  VecDot(tmpp,gold,&tmpScalar);
471  gamma = gamma/tmpScalar;
472 
473  VecWAXPY(wnew,gamma,wold,gnew);
474 
475  /* Assign new to old iterates */
476  VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold);
477 
478  count++;
479  }
480  }
481 
482  /* Restore pressure and velocity sol. vectors to its global form */
483  VecRestoreSubVector(out_uf->pvector->vector,isv,&unew);
484  VecRestoreSubVector(out_uf->pvector->vector,isp,&p);
485 
486  /*return output pointer*/
487  *puf=out_uf;
488 
489 
490  /* Cleanup */
491  KSPDestroy(&kspu);KSPDestroy(&kspip);
492 
493  MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP);
494 
495  VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp);
496  VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi);
497  VecDestroy(&tmpp);VecDestroy(&tmpu);VecDestroy(&f1);VecDestroy(&f2);
498  VecDestroy(&tmpu2);VecDestroy(&tmpp2);VecDestroy(&resu);VecDestroy(&resp);
499 
500 }/*}}}*/
501 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){/*{{{*/
502 
503  /*output*/
504  bool converged=false;
505 
506  /*intermediary*/
507  //Vector<IssmDouble>* KU=NULL;
508  //Vector<IssmDouble>* KUF=NULL;
509  //Vector<IssmDouble>* KUold=NULL;
510  //Vector<IssmDouble>* KUoldF=NULL;
511  Vector<IssmDouble>* duf=NULL;
512  IssmDouble ndu,nduinf,nu;
513  IssmDouble nKUF;
514  IssmDouble nKUoldF;
515  IssmDouble nF;
516  IssmDouble solver_residue,res;
517  int analysis_type;
518 
519  Mat A, B, BT;
520  Vec u,p,uold,pold,f1,f2,tmp,res1,res2;
521  int n_u,n_p;
522  double rnorm1, rnorm2;
523 
524 
525  if(VerboseModule()) _printf0_(" checking convergence\n");
526 
527  /*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from
528  * the get go: */
529  if(uf->IsEmpty()){
530  *pconverged=true;
531  return;
532  }
533 
534  /* Note: SchurCG also constructs the Schur preconditioner and stores it in the free block of Kff
535  * [A B]
536  * Kff = | |
537  * [B^T IP]
538  * To calculate the residual, only the necessary blocks need to be extracted */
539 
540  /*Extract A, B, B^T */
541  #if _PETSC_MINOR_>8
542  MatCreateSubMatrix(Kff->pmatrix->matrix,isv,isv,MAT_INITIAL_MATRIX,&A);
543  MatCreateSubMatrix(Kff->pmatrix->matrix,isv,isp,MAT_INITIAL_MATRIX,&B);
544  MatCreateSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT);
545  #else
546  MatGetSubMatrix(Kff->pmatrix->matrix,isv,isv,MAT_INITIAL_MATRIX,&A);
547  MatGetSubMatrix(Kff->pmatrix->matrix,isv,isp,MAT_INITIAL_MATRIX,&B);
548  MatGetSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT);
549  #endif
550 
551  /*no. of free nodes in velocity/pressure space*/
552  MatGetSize(B,&n_u,&n_p);
553 
554  /*Extract values corresponding to the free velocity/pressure nodes*/
555  VecCreate(IssmComm::GetComm(),&p);VecSetSizes(p,PETSC_DECIDE,n_p);VecSetFromOptions(p);
556  VecAssemblyBegin(p);VecAssemblyEnd(p);
557  VecCreate(IssmComm::GetComm(),&u);VecSetSizes(u,PETSC_DECIDE,n_u);VecSetFromOptions(u);
558  VecAssemblyBegin(u);VecAssemblyEnd(u);
559 
560  VecGetSubVector(uf->pvector->vector,isv,&u);
561  VecGetSubVector(uf->pvector->vector,isp,&p);
562 
563 
564  /*Extract values of the RHS corresponding to the first/second block*/
565  VecDuplicate(u,&f1);VecSet(f1,1.0);
566  VecDuplicate(p,&f2);VecSet(f2,1.0);
567  VecGetSubVector(pf->pvector->vector,isv,&f1);
568  VecGetSubVector(pf->pvector->vector,isp,&f2);
569 
570  /*Allocate intermediaries*/
571  VecDuplicate(u,&res1);VecSet(res1,1.0);
572  VecDuplicate(u,&tmp);VecSet(tmp,1.0);
573  VecDuplicate(p,&res2);VecSet(res2,1.0);
574 
575 
576  /*Display solver caracteristics*/
577  if (VerboseConvergence()){
578 
579  /*Calculate res1 = A*u + B*p - f1*/
580  VecScale(f1,-1.);MatMultAdd(A,u,f1,tmp);MatMultAdd(B,p,tmp,res1);VecScale(f1,-1.);
581  /*Calculate res2 = B^T * u - f2*/
582  VecScale(f2,-1.);MatMultAdd(BT,u,f2,res2);VecScale(f2,-1.);
583 
584 
585  /*compute norm(res1), norm(res2), norm(F) and residue*/
586  VecNorm(res1,NORM_2,&rnorm1);VecNorm(res2,NORM_2,&rnorm2);
587  nKUF=sqrt(rnorm1*rnorm1 + rnorm2*rnorm2);
588  nF=pf->Norm(NORM_TWO);
589  solver_residue=nKUF/nF;
590  _printf0_("\n" << " solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n");
591  if(xIsNan<IssmDouble>(solver_residue)){
592  //Kff->Echo();
593  }
594 
595  }
596  /*clean up*/
597  VecRestoreSubVector(uf->pvector->vector,isv,&u);
598  VecRestoreSubVector(uf->pvector->vector,isp,&p);
599 
600  /*Extract values corresponding to velocity/pressure on the old solution*/
601  VecGetSubVector(old_uf->pvector->vector,isv,&uold);
602  VecGetSubVector(old_uf->pvector->vector,isp,&pold);
603 
604 
605  /*Force equilibrium (Mandatory)*/
606 
607  /*Calculate res1 = A*uold + B*pold - f1*/
608  VecScale(f1,-1.);MatMultAdd(A,uold,f1,tmp);MatMultAdd(B,pold,tmp,res1);VecScale(f1,-1.);
609  /*Calculate res2 = B^T * uold - f2*/
610  VecScale(f2,-1.);MatMultAdd(BT,uold,f2,res2);VecScale(f2,-1.);
611 
612  /*compute norm(res1), norm(res2), norm(F) and residue*/
613  VecNorm(res1,NORM_2,&rnorm1);VecNorm(res2,NORM_2,&rnorm2);
614  nKUoldF=sqrt(rnorm1*rnorm1 + rnorm2*rnorm2);
615  nF=pf->Norm(NORM_TWO);
616  res=nKUoldF/nF;
617  if (xIsNan<IssmDouble>(res)){
618  _printf0_("norm nf = " << nF << "f and norm kuold = " << nKUoldF << "f\n");
619  _error_("mechanical equilibrium convergence criterion is NaN!");
620  }
621 
622  MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);
623  VecRestoreSubVector(pf->pvector->vector,isv,&f1);
624  VecRestoreSubVector(pf->pvector->vector,isp,&f2);
625  VecDestroy(&res1);VecDestroy(&res2);VecDestroy(&tmp);
626  VecRestoreSubVector(old_uf->pvector->vector,isv,&uold);
627  VecRestoreSubVector(old_uf->pvector->vector,isp,&pold);
628 
629 
630 
631  //print
632  if(res<eps_res){
633  if(VerboseConvergence()) _printf0_(setw(50)<<left<<" mechanical equilibrium convergence criterion"<<res*100<< " < "<<eps_res*100<<" %\n");
634  converged=true;
635  }
636  else{
637  if(VerboseConvergence()) _printf0_(setw(50)<<left<<" mechanical equilibrium convergence criterion"<<res*100<<" > "<<eps_res*100<<" %\n");
638  converged=false;
639  }
640 
641  /*Relative criterion (optional)*/
642  if (!xIsNan<IssmDouble>(eps_rel) || (VerboseConvergence())){
643 
644  //compute norm(du)/norm(u)
645  duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
646  ndu=duf->Norm(NORM_TWO); nu=old_uf->Norm(NORM_TWO);
647 
648  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
649 
650  //clean up
651  delete duf;
652 
653  //print
654  if (!xIsNan<IssmDouble>(eps_rel)){
655  if((ndu/nu)<eps_rel){
656  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
657  }
658  else{
659  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
660  converged=false;
661  }
662  }
663  else _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " %\n");
664 
665  }
666 
667  /*Absolute criterion (Optional) = max(du)*/
668  if (!xIsNan<IssmDouble>(eps_abs) || (VerboseConvergence())){
669 
670  //compute max(du)
671  duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
672  ndu=duf->Norm(NORM_TWO); nduinf=duf->Norm(NORM_INF);
673  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
674 
675  //clean up
676  delete duf;
677 
678  //print
679  if (!xIsNan<IssmDouble>(eps_abs)){
680  if ((nduinf)<eps_abs){
681  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: max(du)" << nduinf << " < " << eps_abs << "\n");
682  }
683  else{
684  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: max(du)" << nduinf << " > " << eps_abs << "\n");
685  converged=false;
686  }
687  }
688  else _printf0_(setw(50) << left << " Convergence criterion: max(du)" << nduinf << "\n");
689 
690  }
691 
692  /*assign output*/
693  *pconverged=converged;
694 }/*}}}*/
695 void solutionsequence_schurcg(FemModel* femmodel){/*{{{*/
696 
697  /*intermediary: */
698  Matrix<IssmDouble>* Kff = NULL;
699  Matrix<IssmDouble>* Kfs = NULL;
700  Vector<IssmDouble>* ug = NULL;
701  Vector<IssmDouble>* uf = NULL;
702  Vector<IssmDouble>* old_uf = NULL;
703  Vector<IssmDouble>* pf = NULL;
704  Vector<IssmDouble>* df = NULL;
705  Vector<IssmDouble>* ys = NULL;
706 
707  /*parameters:*/
708  int max_nonlinear_iterations;
709  int configuration_type;
710  int precond;
711  IssmDouble eps_res,eps_rel,eps_abs;
712 
713  /*Recover parameters: */
714  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
715  femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
716  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
717  femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
718  femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
719  femmodel->UpdateConstraintsx();
720  int size;
721  int count=0;
722  bool converged=false;
723 
724  /*Start non-linear iteration using input velocity: */
725  GetSolutionFromInputsx(&ug,femmodel);
726  Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
727 
728  /*Update once again the solution to make sure that vx and vxold are similar*/
729  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
730  InputUpdateFromSolutionx(femmodel,ug);
731 
732  for(;;){
733 
734  /*save pointer to old velocity*/
735  delete old_uf; old_uf=uf;
736  delete ug;
737 
738  /*Get stiffness matrix and Load vector*/
739  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
740  CreateNodalConstraintsx(&ys,femmodel->nodes);
741  Reduceloadx(pf, Kfs, ys); delete Kfs;
742 
743  /*Create pressure matrix of choice*/
744  #if _PETSC_MINOR_<7
745  PetscOptionsGetInt(PETSC_NULL,"-schur_pc",&precond,NULL);
746  #else
747  PetscOptionsGetInt(NULL,PETSC_NULL,"-schur_pc",&precond,NULL);
748  #endif
749 
750  StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
751 
752  /* Construct Schur preconditioner matrix or mass matrix depending on input */
753  if(precond)
754  {
755  for(int i=0;i<femmodel->elements->Size();i++){
756  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
757  ElementMatrix* Ie = analysis->CreateSchurPrecondMatrix(element);
758  if(Ie) Ie->AddToGlobal(Kff,NULL);
759  delete Ie;
760  }
761  }else{
762 
763  for(int i=0;i<femmodel->elements->Size();i++){
764  Element* element2=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
765  ElementMatrix* Ie2 = analysis->CreatePressureMassMatrix(element2);
766  if(Ie2) Ie2->AddToGlobal(Kff,NULL);
767  delete Ie2;
768  }
769  }
770 
771 
772  delete analysis;
773 
774  /*Obtain index sets for velocity and pressure components */
775  IS isv = NULL;
776  IS isp = NULL;
777  #if _PETSC_MAJOR_==3
778 
779  /*Make indices out of doftypes: */
780  if(!(df->pvector->vector))_error_("need doftypes for FS solver!\n");
781  DofTypesToIndexSet(&isv,&isp,df->pvector->vector,FSSolverEnum);
782  #else
783  _error_("Petsc 3.X required");
784  #endif
785  Kff->Assemble();
786 
787 
788  /*Solve*/
789  femmodel->profiler->Start(SOLVER);
790  _assert_(Kff->type==PetscMatType);
791 
792 
793  SchurCGSolver(&uf,
794  Kff->pmatrix->matrix,
795  pf->pvector->vector,
796  old_uf->pvector->vector,
797  isv,
798  isp,
799  femmodel->parameters);
800  femmodel->profiler->Stop(SOLVER);
801 
802  /*Merge solution from f set to g set*/
803  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
804 
805  /*Check for convergence and update inputs accordingly*/
806  convergence_schurcg(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs,isv,isp); delete Kff; delete pf; delete df;
807  count++;
808 
809  if(count>=max_nonlinear_iterations){
810  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
811  converged=true;
812  }
813 
814  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
815  InputUpdateFromSolutionx(femmodel,ug);
816 
817  /*Increase count: */
818  if(converged==true){
819  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,StressbalanceConvergenceNumStepsEnum,count));
820  break;
821  }
822  }
823 
824  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count << "\n");
825 
826  /*clean-up*/
827  delete uf;
828  delete ug;
829  delete old_uf;
830 
831 }/*}}}*/
832 
833 #else
834 void solutionsequence_schurcg(FemModel* femmodel){_error_("PETSc needs to be installed");}
835 #endif
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Matrix< IssmDouble >
SOLVER
#define SOLVER
Definition: Profiler.h:16
StressbalanceAnalysis
Definition: StressbalanceAnalysis.h:11
Mergesolutionfromftogx
void Mergesolutionfromftogx(Vector< IssmDouble > **pug, Vector< IssmDouble > *uf, Vector< IssmDouble > *ys, Nodes *nodes, Parameters *parameters, bool flag_ys0)
Definition: Mergesolutionfromftogx.cpp:8
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
VerboseConvergence
bool VerboseConvergence(void)
Definition: Verbosity.cpp:26
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
ConvergedEnum
@ ConvergedEnum
Definition: EnumDefinitions.h:514
solutionsequence_schurcg
void solutionsequence_schurcg(FemModel *femmodel)
Definition: solutionsequence_schurcg.cpp:834
StressbalanceAnalysis::CreatePressureMassMatrix
ElementMatrix * CreatePressureMassMatrix(Element *element)
Definition: StressbalanceAnalysis.cpp:3340
StressbalanceAnalysis::CreateSchurPrecondMatrix
ElementMatrix * CreateSchurPrecondMatrix(Element *element)
Definition: StressbalanceAnalysis.cpp:3389
FemModel::results
Results * results
Definition: FemModel.h:48
StressbalanceConvergenceNumStepsEnum
@ StressbalanceConvergenceNumStepsEnum
Definition: EnumDefinitions.h:1286
Matrix::Assemble
void Assemble(void)
Definition: Matrix.h:185
NORM_INF
@ NORM_INF
Definition: toolkitsenums.h:15
Vector::Duplicate
Vector< doubletype > * Duplicate(void)
Definition: Vector.h:230
FSSolverEnum
@ FSSolverEnum
Definition: EnumDefinitions.h:1061
VerboseModule
bool VerboseModule(void)
Definition: Verbosity.cpp:23
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
StressbalanceMaxiterEnum
@ StressbalanceMaxiterEnum
Definition: EnumDefinitions.h:407
Vector::Norm
doubletype Norm(NormMode norm_type)
Definition: Vector.h:342
Element
Definition: Element.h:41
Matrix::type
int type
Definition: Matrix.h:31
FemModel::UpdateConstraintsx
void UpdateConstraintsx(void)
Definition: FemModel.cpp:3027
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
solutionsequences.h
StressbalanceAbstolEnum
@ StressbalanceAbstolEnum
Definition: EnumDefinitions.h:404
StressbalanceReltolEnum
@ StressbalanceReltolEnum
Definition: EnumDefinitions.h:410
ConfigurationTypeEnum
@ ConfigurationTypeEnum
Definition: EnumDefinitions.h:101
Results::AddResult
int AddResult(ExternalResult *result)
Definition: Results.cpp:33
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
Profiler::Stop
void Stop(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:179
DofTypesToIndexSet
void DofTypesToIndexSet(IS *pisv, IS *pisp, Vec df, int typeenum)
Definition: PetscSolver.cpp:168
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Reducevectorgtofx
void Reducevectorgtofx(Vector< IssmDouble > **puf, Vector< IssmDouble > *ug, Nodes *nodes, Parameters *parameters)
Definition: Reducevectorgtofx.cpp:8
FemModel
Definition: FemModel.h:31
Reduceloadx
void Reduceloadx(Vector< IssmDouble > *pf, Matrix< IssmDouble > *Kfs, Vector< IssmDouble > *y_s, bool flag_ys0)
Definition: Reduceloadx.cpp:14
Vector::IsEmpty
bool IsEmpty(void)
Definition: Vector.h:196
GenericExternalResult
Definition: GenericExternalResult.h:21
Vector::AYPX
void AYPX(Vector *X, doubletype a)
Definition: Vector.h:267
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Profiler::Start
void Start(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:139
InputUpdateFromSolutionx
void InputUpdateFromSolutionx(FemModel *femmodel, Vector< IssmDouble > *solution)
Definition: InputUpdateFromSolutionx.cpp:9
StressbalanceRestolEnum
@ StressbalanceRestolEnum
Definition: EnumDefinitions.h:412
GetSolutionFromInputsx
void GetSolutionFromInputsx(Vector< IssmDouble > **psolution, FemModel *femmodel)
Definition: GetSolutionFromInputsx.cpp:9
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ElementMatrix::AddToGlobal
void AddToGlobal(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: ElementMatrix.cpp:271
Vector::Copy
void Copy(Vector *to)
Definition: Vector.h:319
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
CreateNodalConstraintsx
void CreateNodalConstraintsx(Vector< IssmDouble > **pys, Nodes *nodes)
Definition: CreateNodalConstraintsx.cpp:10
SystemMatricesx
void SystemMatricesx(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **ppf, Vector< IssmDouble > **pdf, IssmDouble *pkmax, FemModel *femmodel, bool isAllocated)
Definition: SystemMatricesx.cpp:10
ElementMatrix
Definition: ElementMatrix.h:19
FemModel::profiler
Profiler * profiler
Definition: FemModel.h:42
Vector< IssmDouble >
PetscMatType
@ PetscMatType
Definition: Matrix.h:22
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:22 for Ice Sheet System Model by   doxygen 1.8.19