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-6, cgprec=1e-40;
  68   PetscInt giter=0, ColS=8, 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, 30); CHKERRQ(ierr);
 132   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 133
 134
 135
 136   //GMRES WITH MINIMIZATION
 137   T1 = MPI_Wtime();
 138   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
 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=8, col=0, Emaxiter=50000000, iter=0, iterations=15, 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, 30); CHKERRQ(ierr);
 323   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 324
 325
 326
 327
 328   //GMRES WITH MINIMIZATION
 329   T1 = MPI_Wtime();
 330   ierr = KSPSetUp(ksp); CHKERRQ(ierr);
 331   while(giter<Emaxiter && norm>Eprecision ){
 332     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 333
 334
 335       //Solve
 336       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 337
 338       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 339       total += its; Iiter ++;
 340
 341
 342
 343       //Build S'
 344       ierr = VecGetArray(x, &array);
 345       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 346       VecRestoreArray(x, &array);
 347
 348
 349
 350       KSPGetResidualNorm(ksp,&norm);
 351       /*
 352       //Error
 353       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 354       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 355       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 356        */
 357
 358
 359       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 360       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 361
 362
 363     }
 364
 365
 366     //minimization step
 367     if( norm>Eprecision) {
 368
 369       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 370       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 371
 372
 373
 374
 375       //Build AS
 376       if(first) {
 377         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 378
 379         first=0;
 380       }
 381       else
 382         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 383
 384
 385
 386
 387
 388       //LSQR
 389       //LSQR
 390       //LSQR
 391
 392
 393
 394       PetscScalar n2b,tolb,normr,c,s,phibar,normar,norma,thet,rhot,rho,phi;
 395       PetscInt stag;
 396       tolb = tol * n2b;
 397       VecNorm(b, NORM_2, &n2b); //n2b = norm(b);
 398       ierr = VecCopy(b, u); CHKERRQ(ierr);   //u=b
 399       VecNorm(u, NORM_2, &beta); // beta=norm(u)
 400       normr=beta;
 401       if (beta != 0) {
 402         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 403       }
 404       c=1;
 405       s=0;
 406       phibar=beta;
 407       MatMultTranspose(AS, u, v); //v=A'*u
 408       ierr = VecSet(x_lsqr,0);CHKERRQ(ierr);
 409       VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 410       if (alpha != 0) {
 411         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 412       }
 413       ierr = VecSet(d,0);CHKERRQ(ierr);
 414       normar = alpha * beta;
 415       norma=0;
 416       //stag=0;
 417       for(i=0;i<iterations;i++) {
 418         MatMult(AS, v, uu); //uu=A*v
 419         VecAYPX(u, -alpha, uu); //u = uu-alpha*u;
 420         VecNorm(u, NORM_2, &beta); // beta=norm(u)
 421         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 422         norma=sqrt(norma*norma+alpha*alpha+beta*beta); // norma = norm([norma alpha beta]);
 423         thet = - s * alpha;
 424         rhot = c * alpha;
 425         rho = sqrt(rhot*rhot + beta*beta);
 426         c = rhot / rho;
 427         s = - beta / rho;
 428         phi = c * phibar;
 429         if (phi == 0) {             // stagnation of the method
 430           stag = 1;
 431         }
 432         phibar = s * phibar;
 433         VecAYPX(d,-thet,v);       //d = (v - thet * d);
 434         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
 435
 436
 437
 438         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
 439         normr = abs(s) * normr;
 440         MatMultTranspose(AS, u, vt); //vt=A'*u;
 441         VecAYPX(v,-beta,vt);  // v = vt - beta * v;
 442         VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 443         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 444         normar = alpha * abs( s * phi);
 445       }
 446
 447
 448
 449       iter = 0; giter ++;
 450       //Minimizer
 451       MatMult(S, x_lsqr, x); //x = S*x_small;
 452     }
 453
 454   }
 455   T2 = MPI_Wtime();
 456   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
 457   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
 458
 459   return 0;
 460
 461 }
 462
 463
 464
 465
 466
 467
 468
 469 int main(int argc,char **argv)
 470 {
 471   KSP            ksp;
 472   DM             da;
 473   UserContext    user;
 474   const char     *bcTypes[2] = {"dirichlet","neumann"};
 475   PetscErrorCode ierr;
 476   PetscInt       bc;
 477   Vec            b,x;
 478   Mat A;
 479   PetscScalar T1,T2,norm;
 480   Vec sol;
 481   PetscInt total;
 482
 483   PetscInitialize(&argc,&argv,(char*)0,help);
 484
 485
 486   PetscMPIInt size;
 487   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 488   PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d\n",size);
 489
 490   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 491   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);
 492   ierr = DMDASetUniformCoordinates(da,0,1,0,1,0,0);CHKERRQ(ierr);
 493   ierr = DMDASetFieldName(da,0,"Pressure");CHKERRQ(ierr);
 494
 495   ierr        = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");
 496   user.rho    = 1.0;
 497   ierr        = PetscOptionsReal("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, NULL);CHKERRQ(ierr);
 498   user.nu     = 0.1;
 499   ierr        = PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, NULL);CHKERRQ(ierr);
 500   bc          = (PetscInt)DIRICHLET;
 501   ierr        = PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,NULL);CHKERRQ(ierr);
 502   user.bcType = (BCType)bc;
 503   ierr        = PetscOptionsEnd();
 504
 505   ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);
 506   ierr = KSPSetComputeOperators(ksp,ComputeMatrix,&user);CHKERRQ(ierr);
 507   ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
 508   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 509
 510
 511
 512   PC             pc;
 513   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 514   KSPGetPC(ksp, &pc);
 515   PCType         type;
 516   PCGetType(pc, &type);
 517
 518   PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
 519
 520
 521   ierr = KSPSetTolerances(ksp, 1e-6, 1e-6, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
 522
 523   T1 = MPI_Wtime();
 524     ierr = KSPSetUp(ksp);CHKERRQ(ierr);
 525   ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
 526   T2 = MPI_Wtime();
 527
 528
 529
 530   ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
 531   ierr = KSPGetRhs(ksp,&b);CHKERRQ(ierr);
 532   KSPGetOperators(ksp,&A,NULL);
 533
 534   VecDuplicate(x,&sol);
 535   MatMult(A,x,sol);
 536   VecAXPY(sol,-1,b);
 537   VecNorm(sol, NORM_2, &norm);
 538   ierr = PetscPrintf(PETSC_COMM_WORLD, "Error %g\n",norm);
 539   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
 540
 541   ierr = KSPGetIterationNumber(ksp, &total); CHKERRQ(ierr);
 542   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
 543
 544  {
 545
 546     Vec x;
 547     VecDuplicate(b,&x);
 548
 549     KrylovMinimize(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, "Error2 %g\n",norm);
 559
 560   }
 561
 562
 563  {
 564
 565     Vec x;
 566     VecDuplicate(b,&x);
 567
 568     KrylovMinimizeLSQR(A, b, x);
 569
 570
 571
 572     MatMult(A,x,sol);
 573
 574
 575     VecAXPY(sol,-1,b);
 576     VecNorm(sol, NORM_2, &norm);
 577     ierr = PetscPrintf(PETSC_COMM_WORLD, "Error LSQR %g\n",norm);
 578
 579   }
 580
 581
 582
 583
 584   ierr = DMDestroy(&da);CHKERRQ(ierr);
 585   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
 586   ierr = PetscFinalize();
 587
 588   return 0;
 589 }
 590
 591 #undef __FUNCT__
 592 #define __FUNCT__ "ComputeRHS"
 593 PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
 594 {
 595   UserContext    *user = (UserContext*)ctx;
 596   PetscErrorCode ierr;
 597   PetscInt       i,j,mx,my,xm,ym,xs,ys;
 598   PetscScalar    Hx,Hy;
 599   PetscScalar    **array;
 600   DM             da;
 601
 602   PetscFunctionBeginUser;
 603   ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
 604   ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
 605   Hx   = 1.0 / (PetscReal)(mx-1);
 606   Hy   = 1.0 / (PetscReal)(my-1);
 607   ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
 608   ierr = DMDAVecGetArray(da, b, &array);CHKERRQ(ierr);
 609   for (j=ys; j<ys+ym; j++) {
 610     for (i=xs; i<xs+xm; i++) {
 611       array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
 612     }
 613   }
 614   ierr = DMDAVecRestoreArray(da, b, &array);CHKERRQ(ierr);
 615   ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
 616   ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
 617
 618   /* force right hand side to be consistent for singular matrix */
 619   /* note this is really a hack, normally the model would provide you with a consistent right handside */
 620   if (user->bcType == NEUMANN) {
 621     MatNullSpace nullspace;
 622
 623     ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
 624     ierr = MatNullSpaceRemove(nullspace,b);CHKERRQ(ierr);
 625     ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
 626   }
 627   PetscFunctionReturn(0);
 628 }
 629
 630
 631 #undef __FUNCT__
 632 #define __FUNCT__ "ComputeRho"
 633 PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
 634 {
 635   PetscFunctionBeginUser;
 636   if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
 637     *rho = centerRho;
 638   } else {
 639     *rho = 1.0;
 640   }
 641   PetscFunctionReturn(0);
 642 }
 643
 644 #undef __FUNCT__
 645 #define __FUNCT__ "ComputeMatrix"
 646 PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
 647 {
 648   UserContext    *user = (UserContext*)ctx;
 649   PetscReal      centerRho;
 650   PetscErrorCode ierr;
 651   PetscInt       i,j,mx,my,xm,ym,xs,ys;
 652   PetscScalar    v[5];
 653   PetscReal      Hx,Hy,HydHx,HxdHy,rho;
 654   MatStencil     row, col[5];
 655   DM             da;
 656
 657   PetscFunctionBeginUser;
 658   ierr      = KSPGetDM(ksp,&da);CHKERRQ(ierr);
 659   centerRho = user->rho;
 660   ierr      = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
 661   Hx        = 1.0 / (PetscReal)(mx-1);
 662   Hy        = 1.0 / (PetscReal)(my-1);
 663   HxdHy     = Hx/Hy;
 664   HydHx     = Hy/Hx;
 665   ierr      = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
 666   for (j=ys; j<ys+ym; j++) {
 667     for (i=xs; i<xs+xm; i++) {
 668       row.i = i; row.j = j;
 669       ierr  = ComputeRho(i, j, mx, my, centerRho, &rho);CHKERRQ(ierr);
 670       if (i==0 || j==0 || i==mx-1 || j==my-1) {
 671         if (user->bcType == DIRICHLET) {
 672           v[0] = 2.0*rho*(HxdHy + HydHx);
 673           ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
 674         } else if (user->bcType == NEUMANN) {
 675           PetscInt numx = 0, numy = 0, num = 0;
 676           if (j!=0) {
 677             v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j-1;
 678             numy++; num++;
 679           }
 680           if (i!=0) {
 681             v[num] = -rho*HydHx;              col[num].i = i-1; col[num].j = j;
 682             numx++; num++;
 683           }
 684           if (i!=mx-1) {
 685             v[num] = -rho*HydHx;              col[num].i = i+1; col[num].j = j;
 686             numx++; num++;
 687           }
 688           if (j!=my-1) {
 689             v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j+1;
 690             numy++; num++;
 691           }
 692           v[num] = numx*rho*HydHx + numy*rho*HxdHy; col[num].i = i;   col[num].j = j;
 693           num++;
 694           ierr = MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
 695         }
 696       } else {
 697         v[0] = -rho*HxdHy;              col[0].i = i;   col[0].j = j-1;
 698         v[1] = -rho*HydHx;              col[1].i = i-1; col[1].j = j;
 699         v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
 700         v[3] = -rho*HydHx;              col[3].i = i+1; col[3].j = j;
 701         v[4] = -rho*HxdHy;              col[4].i = i;   col[4].j = j+1;
 702         ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
 703       }
 704     }
 705   }
 706   ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 707   ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 708   if (user->bcType == NEUMANN) {
 709     MatNullSpace nullspace;
 710
 711     ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
 712     ierr = MatSetNullSpace(jac,nullspace);CHKERRQ(ierr);
 713     ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
 714   }
 715   PetscFunctionReturn(0);
 716 }