code/ex29.c

   1
   2
   3 //  /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 3 ex29 -da_grid_x 600 -da_grid_y 600 -ksp_type fgmres
   4
   5
   6
   7 /*T
   8    Concepts: KSP^solving a system of linear equations
   9    Concepts: KSP^Laplacian, 2d
  10    Processors: n
  11 T*/
  12
  13 /*
  14 Added at the request of Marc Garbey.
  15
  16 Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
  17
  18    -div \rho grad u = f,  0 < x,y < 1,
  19
  20 with forcing function
  21
  22    f = e^{-x^2/\nu} e^{-y^2/\nu}
  23
  24 with Dirichlet boundary conditions
  25
  26    u = f(x,y) for x = 0, x = 1, y = 0, y = 1
  27
  28 or pure Neumman boundary conditions
  29
  30 This uses multigrid to solve the linear system
  31 */
  32
  33 static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
  34
  35 #include <petscdm.h>
  36 #include <petscdmda.h>
  37 #include <petscksp.h>
  38
  39 extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
  40 extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
  41
  42 typedef enum {DIRICHLET, NEUMANN} BCType;
  43
  44 typedef struct {
  45   PetscReal rho;
  46   PetscReal nu;
  47   BCType    bcType;
  48 } UserContext;
  49
  50 #undef __FUNCT__
  51 #define __FUNCT__ "main"
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61 int KrylovMinimize(Mat A, Vec b, Vec x) {
  62
  63
  64   //Variables
  65
  66   PetscScalar  gamma, alpha, oldgamma, beta;
  67   PetscReal norm=20, Eprecision=1e-7, cgprec=1e-40;
  68   PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
  69   PetscErrorCode ierr;
  70   PetscScalar T1, T2;
  71   KSP ksp;
  72   PetscInt total=0;
  73   PetscInt size;
  74   PetscInt Istart,Iend;
  75   PetscInt i,its;
  76   Vec       x_old, residu;
  77   Mat S, AS;
  78   PetscScalar *array;
  79   PetscInt *ind_row;
  80   Vec Alpha, p, ss, vect, r, q, Ax;
  81
  82
  83   PetscInt first=1;
  84
  85   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  86   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
  87   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  88
  89
  90
  91
  92   VecGetSize(b,&size);
  93
  94   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
  95
  96   PetscInt aa,bb;
  97   MatGetOwnershipRange(A,&aa,&bb);
  98
  99   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
 100   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
 101
 102
 103   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
 104   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
 105   ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
 106   ierr = MatSetUp(S); CHKERRQ(ierr);
 107
 108   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
 109
 110   ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
 111   ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
 112   ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
 113   ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
 114   ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
 115   ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
 116   ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
 117   ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
 118   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 119
 120   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
 121   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 122
 123
 124
 125   //indexes of row (these indexes are global)
 126   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
 127   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
 128
 129   //Initializations
 130   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
 131   ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
 132   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 133
 134
 135
 136   //GMRES WITH MINIMIZATION
 137   T1 = MPI_Wtime();
 138   while(giter<Emaxiter && norm>Eprecision ){
 139     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 140
 141
 142       //Solve
 143       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 144
 145       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 146       total += its; Iiter ++;
 147
 148
 149
 150       //Build S'
 151       ierr = VecGetArray(x, &array);
 152       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 153       VecRestoreArray(x, &array);
 154
 155
 156
 157       KSPGetResidualNorm(ksp,&norm);
 158
 159       /* //Error
 160       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 161       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 162       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 163        */
 164
 165
 166       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 167       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 168
 169
 170     }
 171
 172
 173     //minimization step
 174     if( norm>Eprecision) {
 175
 176       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 177       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 178
 179
 180       //Build AS
 181       if(first) {
 182         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 183         first=0;
 184       }
 185       else
 186         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 187
 188
 189
 190
 191       //Minimization with CGLS
 192       MatMult(AS, Alpha, r);  VecAYPX(r, -1, b); //r_0 = b-AS*x_0
 193
 194
 195       MatMultTranspose(AS, r, p); //p_0 = AS'*r_0
 196
 197
 198
 199
 200       ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
 201       VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 202       while(gamma>cgprec && iter<iterations){
 203         MatMult(AS, p, q); //q = AS*p
 204         VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
 205         alpha = gamma / alpha;
 206         VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
 207         VecAXPY(r, -alpha, q); //r -= alpha*q
 208         MatMultTranspose(AS, r, ss); // ss = AS'*r
 209         oldgamma = gamma;
 210         VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 211         beta = gamma / oldgamma;
 212         VecAYPX(p, beta, ss); //p = s+beta*p;
 213         iter ++;
 214       }
 215
 216       iter = 0; giter ++;
 217       //Minimizer
 218       MatMult(S, Alpha, x); //x = S*Alpha;
 219     }
 220
 221   }
 222   T2 = MPI_Wtime();
 223   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
 224   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
 225
 226   return 0;
 227
 228 }
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238 int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 239
 240
 241   //Variables
 242
 243   PetscScalar  alpha, beta;
 244   PetscReal norm=20, Eprecision=1e-7, tol=1e-40;
 245   PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
 246   PetscErrorCode ierr;
 247   PetscScalar T1, T2;
 248   KSP ksp;
 249   PetscInt total=0;
 250   PetscInt size;
 251   PetscInt Istart,Iend;
 252   PetscInt i,its;
 253   Vec       x_old, residu;
 254   Mat S, AS;
 255   PetscScalar *array;
 256   PetscInt *ind_row;
 257   Vec  Ax;
 258   PetscScalar norm_ksp;
 259   Vec u,v,d,uu,vt,zero_long,zero_short,x_lsqr;
 260
 261   PetscInt first=1;
 262
 263   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 264   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
 265   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 266
 267
 268
 269
 270   VecGetSize(b,&size);
 271
 272   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
 273
 274   PetscInt aa,bb;
 275   MatGetOwnershipRange(A,&aa,&bb);
 276
 277   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
 278   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
 279
 280
 281   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
 282   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
 283   ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
 284   ierr = MatSetUp(S); CHKERRQ(ierr);
 285
 286   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
 287
 288
 289
 290   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 291
 292   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
 293   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 294
 295
 296   //long vector
 297   ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
 298
 299
 300   ierr = VecDuplicate(b,&uu);CHKERRQ(ierr);
 301   ierr = VecDuplicate(b,&zero_long);CHKERRQ(ierr);
 302   ierr = VecSet(zero_long,0);CHKERRQ(ierr);
 303
 304   //small vector
 305   ierr = VecCreate(PETSC_COMM_WORLD, &v); CHKERRQ(ierr);
 306   ierr = VecSetSizes(v, PETSC_DECIDE, ColS); CHKERRQ(ierr);
 307   ierr = VecSetFromOptions(v); CHKERRQ(ierr);
 308   ierr = VecDuplicate(v,&zero_short);CHKERRQ(ierr);
 309   ierr = VecSet(zero_short,0);CHKERRQ(ierr);
 310   ierr = VecDuplicate(v,&d);CHKERRQ(ierr);
 311   ierr = VecDuplicate(v,&vt);CHKERRQ(ierr);
 312   ierr = VecDuplicate(v,&x_lsqr);CHKERRQ(ierr);
 313
 314
 315   //indexes of row (these indexes are global)
 316   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
 317   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
 318
 319   //Initializations
 320   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
 321   ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
 322   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 323
 324
 325
 326
 327   //GMRES WITH MINIMIZATION
 328   T1 = MPI_Wtime();
 329   while(giter<Emaxiter && norm>Eprecision ){
 330     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 331
 332
 333       //Solve
 334       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 335
 336       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 337       total += its; Iiter ++;
 338
 339
 340
 341       //Build S'
 342       ierr = VecGetArray(x, &array);
 343       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 344       VecRestoreArray(x, &array);
 345
 346
 347
 348       KSPGetResidualNorm(ksp,&norm);
 349       /*
 350       //Error
 351       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 352       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 353       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 354        */
 355
 356
 357       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 358       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 359
 360
 361     }
 362
 363
 364     //minimization step
 365     if( norm>Eprecision) {
 366
 367       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 368       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 369
 370
 371
 372
 373       //Build AS
 374       if(first) {
 375         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 376
 377         first=0;
 378       }
 379       else
 380         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 381
 382
 383
 384
 385
 386       //LSQR
 387       //LSQR
 388       //LSQR
 389
 390
 391
 392       PetscScalar n2b,tolb,normr,c,s,phibar,normar,norma,thet,rhot,rho,phi;
 393       PetscInt stag;
 394       tolb = tol * n2b;
 395       VecNorm(b, NORM_2, &n2b); //n2b = norm(b);
 396       ierr = VecCopy(b, u); CHKERRQ(ierr);   //u=b
 397       VecNorm(u, NORM_2, &beta); // beta=norm(u)
 398       normr=beta;
 399       if (beta != 0) {
 400         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 401       }
 402       c=1;
 403       s=0;
 404       phibar=beta;
 405       MatMultTranspose(AS, u, v); //v=A'*u
 406       ierr = VecSet(x_lsqr,0);CHKERRQ(ierr);
 407       VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 408       if (alpha != 0) {
 409         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 410       }
 411       ierr = VecSet(d,0);CHKERRQ(ierr);
 412       normar = alpha * beta;
 413       norma=0;
 414       //stag=0;
 415       for(i=0;i<iterations;i++) {
 416         MatMult(AS, v, uu); //uu=A*v
 417         VecAYPX(u, -alpha, uu); //u = uu-alpha*u;
 418         VecNorm(u, NORM_2, &beta); // beta=norm(u)
 419         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 420         norma=sqrt(norma*norma+alpha*alpha+beta*beta); // norma = norm([norma alpha beta]);
 421         thet = - s * alpha;
 422         rhot = c * alpha;
 423         rho = sqrt(rhot*rhot + beta*beta);
 424         c = rhot / rho;
 425         s = - beta / rho;
 426         phi = c * phibar;
 427         if (phi == 0) {             // stagnation of the method
 428           stag = 1;
 429         }
 430         phibar = s * phibar;
 431         VecAYPX(d,-thet,v);       //d = (v - thet * d);
 432         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
 433
 434
 435
 436         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
 437         normr = abs(s) * normr;
 438         MatMultTranspose(AS, u, vt); //vt=A'*u;
 439         VecAYPX(v,-beta,vt);  // v = vt - beta * v;
 440         VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 441         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 442         normar = alpha * abs( s * phi);
 443       }
 444
 445
 446
 447       iter = 0; giter ++;
 448       //Minimizer
 449       MatMult(S, x_lsqr, x); //x = S*x_small;
 450     }
 451
 452   }
 453   T2 = MPI_Wtime();
 454   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
 455   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
 456
 457   return 0;
 458
 459 }
 460
 461
 462
 463
 464
 465
 466
 467 int main(int argc,char **argv)
 468 {
 469   KSP            ksp;
 470   DM             da;
 471   UserContext    user;
 472   const char     *bcTypes[2] = {"dirichlet","neumann"};
 473   PetscErrorCode ierr;
 474   PetscInt       bc;
 475   Vec            b,x;
 476   Mat A;
 477   PetscScalar T1,T2,norm;
 478   Vec sol;
 479   PetscInt total;
 480
 481   PetscInitialize(&argc,&argv,(char*)0,help);
 482
 483   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 484   ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);
 485   ierr = DMDASetUniformCoordinates(da,0,1,0,1,0,0);CHKERRQ(ierr);
 486   ierr = DMDASetFieldName(da,0,"Pressure");CHKERRQ(ierr);
 487
 488   ierr        = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");
 489   user.rho    = 1.0;
 490   ierr        = PetscOptionsReal("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, NULL);CHKERRQ(ierr);
 491   user.nu     = 0.1;
 492   ierr        = PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, NULL);CHKERRQ(ierr);
 493   bc          = (PetscInt)DIRICHLET;
 494   ierr        = PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,NULL);CHKERRQ(ierr);
 495   user.bcType = (BCType)bc;
 496   ierr        = PetscOptionsEnd();
 497
 498   ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);
 499   ierr = KSPSetComputeOperators(ksp,ComputeMatrix,&user);CHKERRQ(ierr);
 500   ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
 501   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 502   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
 503
 504   ierr = KSPSetTolerances(ksp, 1e-7, 1e-7, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
 505   T1 = MPI_Wtime();
 506   ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
 507   T2 = MPI_Wtime();
 508
 509
 510
 511   ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
 512   ierr = KSPGetRhs(ksp,&b);CHKERRQ(ierr);
 513   KSPGetOperators(ksp,&A,NULL);
 514
 515   VecDuplicate(x,&sol);
 516   MatMult(A,x,sol);
 517   VecAXPY(sol,-1,b);
 518   VecNorm(sol, NORM_2, &norm);
 519   ierr = PetscPrintf(PETSC_COMM_WORLD, "Error %g\n",norm);
 520   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
 521
 522   ierr = KSPGetIterationNumber(ksp, &total); CHKERRQ(ierr);
 523   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
 524
 525  {
 526
 527     Vec x;
 528     VecDuplicate(b,&x);
 529
 530     KrylovMinimize(A, b, x);
 531
 532
 533
 534     MatMult(A,x,sol);
 535
 536
 537     VecAXPY(sol,-1,b);
 538     VecNorm(sol, NORM_2, &norm);
 539     ierr = PetscPrintf(PETSC_COMM_WORLD, "Error2 %g\n",norm);
 540
 541   }
 542
 543
 544  {
 545
 546     Vec x;
 547     VecDuplicate(b,&x);
 548
 549     KrylovMinimizeLSQR(A, b, x);
 550
 551
 552
 553     MatMult(A,x,sol);
 554
 555
 556     VecAXPY(sol,-1,b);
 557     VecNorm(sol, NORM_2, &norm);
 558     ierr = PetscPrintf(PETSC_COMM_WORLD, "Error LSQR %g\n",norm);
 559
 560   }
 561
 562
 563
 564
 565   ierr = DMDestroy(&da);CHKERRQ(ierr);
 566   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
 567   ierr = PetscFinalize();
 568
 569   return 0;
 570 }
 571
 572 #undef __FUNCT__
 573 #define __FUNCT__ "ComputeRHS"
 574 PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
 575 {
 576   UserContext    *user = (UserContext*)ctx;
 577   PetscErrorCode ierr;
 578   PetscInt       i,j,mx,my,xm,ym,xs,ys;
 579   PetscScalar    Hx,Hy;
 580   PetscScalar    **array;
 581   DM             da;
 582
 583   PetscFunctionBeginUser;
 584   ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
 585   ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
 586   Hx   = 1.0 / (PetscReal)(mx-1);
 587   Hy   = 1.0 / (PetscReal)(my-1);
 588   ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
 589   ierr = DMDAVecGetArray(da, b, &array);CHKERRQ(ierr);
 590   for (j=ys; j<ys+ym; j++) {
 591     for (i=xs; i<xs+xm; i++) {
 592       array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
 593     }
 594   }
 595   ierr = DMDAVecRestoreArray(da, b, &array);CHKERRQ(ierr);
 596   ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
 597   ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
 598
 599   /* force right hand side to be consistent for singular matrix */
 600   /* note this is really a hack, normally the model would provide you with a consistent right handside */
 601   if (user->bcType == NEUMANN) {
 602     MatNullSpace nullspace;
 603
 604     ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
 605     ierr = MatNullSpaceRemove(nullspace,b);CHKERRQ(ierr);
 606     ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
 607   }
 608   PetscFunctionReturn(0);
 609 }
 610
 611
 612 #undef __FUNCT__
 613 #define __FUNCT__ "ComputeRho"
 614 PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
 615 {
 616   PetscFunctionBeginUser;
 617   if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
 618     *rho = centerRho;
 619   } else {
 620     *rho = 1.0;
 621   }
 622   PetscFunctionReturn(0);
 623 }
 624
 625 #undef __FUNCT__
 626 #define __FUNCT__ "ComputeMatrix"
 627 PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
 628 {
 629   UserContext    *user = (UserContext*)ctx;
 630   PetscReal      centerRho;
 631   PetscErrorCode ierr;
 632   PetscInt       i,j,mx,my,xm,ym,xs,ys;
 633   PetscScalar    v[5];
 634   PetscReal      Hx,Hy,HydHx,HxdHy,rho;
 635   MatStencil     row, col[5];
 636   DM             da;
 637
 638   PetscFunctionBeginUser;
 639   ierr      = KSPGetDM(ksp,&da);CHKERRQ(ierr);
 640   centerRho = user->rho;
 641   ierr      = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
 642   Hx        = 1.0 / (PetscReal)(mx-1);
 643   Hy        = 1.0 / (PetscReal)(my-1);
 644   HxdHy     = Hx/Hy;
 645   HydHx     = Hy/Hx;
 646   ierr      = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
 647   for (j=ys; j<ys+ym; j++) {
 648     for (i=xs; i<xs+xm; i++) {
 649       row.i = i; row.j = j;
 650       ierr  = ComputeRho(i, j, mx, my, centerRho, &rho);CHKERRQ(ierr);
 651       if (i==0 || j==0 || i==mx-1 || j==my-1) {
 652         if (user->bcType == DIRICHLET) {
 653           v[0] = 2.0*rho*(HxdHy + HydHx);
 654           ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
 655         } else if (user->bcType == NEUMANN) {
 656           PetscInt numx = 0, numy = 0, num = 0;
 657           if (j!=0) {
 658             v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j-1;
 659             numy++; num++;
 660           }
 661           if (i!=0) {
 662             v[num] = -rho*HydHx;              col[num].i = i-1; col[num].j = j;
 663             numx++; num++;
 664           }
 665           if (i!=mx-1) {
 666             v[num] = -rho*HydHx;              col[num].i = i+1; col[num].j = j;
 667             numx++; num++;
 668           }
 669           if (j!=my-1) {
 670             v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j+1;
 671             numy++; num++;
 672           }
 673           v[num] = numx*rho*HydHx + numy*rho*HxdHy; col[num].i = i;   col[num].j = j;
 674           num++;
 675           ierr = MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
 676         }
 677       } else {
 678         v[0] = -rho*HxdHy;              col[0].i = i;   col[0].j = j-1;
 679         v[1] = -rho*HydHx;              col[1].i = i-1; col[1].j = j;
 680         v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
 681         v[3] = -rho*HydHx;              col[3].i = i+1; col[3].j = j;
 682         v[4] = -rho*HxdHy;              col[4].i = i;   col[4].j = j+1;
 683         ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
 684       }
 685     }
 686   }
 687   ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 688   ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 689   if (user->bcType == NEUMANN) {
 690     MatNullSpace nullspace;
 691
 692     ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
 693     ierr = MatSetNullSpace(jac,nullspace);CHKERRQ(ierr);
 694     ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
 695   }
 696   PetscFunctionReturn(0);
 697 }