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