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 -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-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=20, 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
 432         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
 433         normr = abs(s) * normr;
 434         MatMultTranspose(AS, u, vt); //vt=A'*u;
 435         VecAYPX(v,-beta,vt);  // v = vt - beta * v;
 436         VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 437         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 438         normar = alpha * abs( s * phi);
 439       }
 440
 441
 442
 443       iter = 0; giter ++;
 444       //Minimizer
 445       MatMult(S, x_lsqr, x); //x = S*x_small;
 446     }
 447
 448   }
 449   T2 = MPI_Wtime();
 450   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
 451   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
 452
 453   return 0;
 454
 455 }
 456
 457
 458
 459
 460
 461
 462
 463 int main(int argc,char **argv)
 464 {
 465   KSP            ksp;
 466   DM             da;
 467   UserContext    user;
 468   const char     *bcTypes[2] = {"dirichlet","neumann"};
 469   PetscErrorCode ierr;
 470   PetscInt       bc;
 471   Vec            b,x;
 472   Mat A;
 473   PetscScalar T1,T2,norm;
 474   Vec sol;
 475   PetscInt total;
 476
 477   PetscInitialize(&argc,&argv,(char*)0,help);
 478
 479   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 480   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);
 481   ierr = DMDASetUniformCoordinates(da,0,1,0,1,0,0);CHKERRQ(ierr);
 482   ierr = DMDASetFieldName(da,0,"Pressure");CHKERRQ(ierr);
 483
 484   ierr        = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");
 485   user.rho    = 1.0;
 486   ierr        = PetscOptionsReal("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, NULL);CHKERRQ(ierr);
 487   user.nu     = 0.1;
 488   ierr        = PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, NULL);CHKERRQ(ierr);
 489   bc          = (PetscInt)DIRICHLET;
 490   ierr        = PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,NULL);CHKERRQ(ierr);
 491   user.bcType = (BCType)bc;
 492   ierr        = PetscOptionsEnd();
 493
 494   ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);
 495   ierr = KSPSetComputeOperators(ksp,ComputeMatrix,&user);CHKERRQ(ierr);
 496   ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
 497   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 498   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
 499
 500   ierr = KSPSetTolerances(ksp, 1e-5, 1e-5, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
 501   T1 = MPI_Wtime();
 502   ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
 503   T2 = MPI_Wtime();
 504
 505
 506
 507   ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
 508   ierr = KSPGetRhs(ksp,&b);CHKERRQ(ierr);
 509   KSPGetOperators(ksp,&A,NULL);
 510
 511   VecDuplicate(x,&sol);
 512   MatMult(A,x,sol);
 513   VecAXPY(sol,-1,b);
 514   VecNorm(sol, NORM_2, &norm);
 515   ierr = PetscPrintf(PETSC_COMM_WORLD, "Error %g\n",norm);
 516   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
 517
 518   ierr = KSPGetIterationNumber(ksp, &total); CHKERRQ(ierr);
 519   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
 520
 521  {
 522
 523     Vec x;
 524     VecDuplicate(b,&x);
 525
 526     KrylovMinimize(A, b, x);
 527
 528
 529
 530     MatMult(A,x,sol);
 531
 532
 533     VecAXPY(sol,-1,b);
 534     VecNorm(sol, NORM_2, &norm);
 535     ierr = PetscPrintf(PETSC_COMM_WORLD, "Error2 %g\n",norm);
 536
 537   }
 538
 539
 540  {
 541
 542     Vec x;
 543     VecDuplicate(b,&x);
 544
 545     KrylovMinimizeLSQR(A, b, x);
 546
 547
 548
 549     MatMult(A,x,sol);
 550
 551
 552     VecAXPY(sol,-1,b);
 553     VecNorm(sol, NORM_2, &norm);
 554     ierr = PetscPrintf(PETSC_COMM_WORLD, "Error LSQR %g\n",norm);
 555
 556   }
 557
 558
 559
 560
 561   ierr = DMDestroy(&da);CHKERRQ(ierr);
 562   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
 563   ierr = PetscFinalize();
 564
 565   return 0;
 566 }
 567
 568 #undef __FUNCT__
 569 #define __FUNCT__ "ComputeRHS"
 570 PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
 571 {
 572   UserContext    *user = (UserContext*)ctx;
 573   PetscErrorCode ierr;
 574   PetscInt       i,j,mx,my,xm,ym,xs,ys;
 575   PetscScalar    Hx,Hy;
 576   PetscScalar    **array;
 577   DM             da;
 578
 579   PetscFunctionBeginUser;
 580   ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
 581   ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
 582   Hx   = 1.0 / (PetscReal)(mx-1);
 583   Hy   = 1.0 / (PetscReal)(my-1);
 584   ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
 585   ierr = DMDAVecGetArray(da, b, &array);CHKERRQ(ierr);
 586   for (j=ys; j<ys+ym; j++) {
 587     for (i=xs; i<xs+xm; i++) {
 588       array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
 589     }
 590   }
 591   ierr = DMDAVecRestoreArray(da, b, &array);CHKERRQ(ierr);
 592   ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
 593   ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
 594
 595   /* force right hand side to be consistent for singular matrix */
 596   /* note this is really a hack, normally the model would provide you with a consistent right handside */
 597   if (user->bcType == NEUMANN) {
 598     MatNullSpace nullspace;
 599
 600     ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
 601     ierr = MatNullSpaceRemove(nullspace,b);CHKERRQ(ierr);
 602     ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
 603   }
 604   PetscFunctionReturn(0);
 605 }
 606
 607
 608 #undef __FUNCT__
 609 #define __FUNCT__ "ComputeRho"
 610 PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
 611 {
 612   PetscFunctionBeginUser;
 613   if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
 614     *rho = centerRho;
 615   } else {
 616     *rho = 1.0;
 617   }
 618   PetscFunctionReturn(0);
 619 }
 620
 621 #undef __FUNCT__
 622 #define __FUNCT__ "ComputeMatrix"
 623 PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
 624 {
 625   UserContext    *user = (UserContext*)ctx;
 626   PetscReal      centerRho;
 627   PetscErrorCode ierr;
 628   PetscInt       i,j,mx,my,xm,ym,xs,ys;
 629   PetscScalar    v[5];
 630   PetscReal      Hx,Hy,HydHx,HxdHy,rho;
 631   MatStencil     row, col[5];
 632   DM             da;
 633
 634   PetscFunctionBeginUser;
 635   ierr      = KSPGetDM(ksp,&da);CHKERRQ(ierr);
 636   centerRho = user->rho;
 637   ierr      = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
 638   Hx        = 1.0 / (PetscReal)(mx-1);
 639   Hy        = 1.0 / (PetscReal)(my-1);
 640   HxdHy     = Hx/Hy;
 641   HydHx     = Hy/Hx;
 642   ierr      = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
 643   for (j=ys; j<ys+ym; j++) {
 644     for (i=xs; i<xs+xm; i++) {
 645       row.i = i; row.j = j;
 646       ierr  = ComputeRho(i, j, mx, my, centerRho, &rho);CHKERRQ(ierr);
 647       if (i==0 || j==0 || i==mx-1 || j==my-1) {
 648         if (user->bcType == DIRICHLET) {
 649           v[0] = 2.0*rho*(HxdHy + HydHx);
 650           ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
 651         } else if (user->bcType == NEUMANN) {
 652           PetscInt numx = 0, numy = 0, num = 0;
 653           if (j!=0) {
 654             v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j-1;
 655             numy++; num++;
 656           }
 657           if (i!=0) {
 658             v[num] = -rho*HydHx;              col[num].i = i-1; col[num].j = j;
 659             numx++; num++;
 660           }
 661           if (i!=mx-1) {
 662             v[num] = -rho*HydHx;              col[num].i = i+1; col[num].j = j;
 663             numx++; num++;
 664           }
 665           if (j!=my-1) {
 666             v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j+1;
 667             numy++; num++;
 668           }
 669           v[num] = numx*rho*HydHx + numy*rho*HxdHy; col[num].i = i;   col[num].j = j;
 670           num++;
 671           ierr = MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
 672         }
 673       } else {
 674         v[0] = -rho*HxdHy;              col[0].i = i;   col[0].j = j-1;
 675         v[1] = -rho*HydHx;              col[1].i = i-1; col[1].j = j;
 676         v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
 677         v[3] = -rho*HydHx;              col[3].i = i+1; col[3].j = j;
 678         v[4] = -rho*HxdHy;              col[4].i = i;   col[4].j = j+1;
 679         ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
 680       }
 681     }
 682   }
 683   ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 684   ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 685   if (user->bcType == NEUMANN) {
 686     MatNullSpace nullspace;
 687
 688     ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
 689     ierr = MatSetNullSpace(jac,nullspace);CHKERRQ(ierr);
 690     ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
 691   }
 692   PetscFunctionReturn(0);
 693 }