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