Actual source code: ex11.c
petsc-3.4.2 2013-07-02

      1: static const char help[] = "Solves a Q1-P0 Stokes problem from Underworld.\n\
      2: \n\
      3: You can obtain a sample matrix from http://ftp.mcs.anl.gov/pub/petsc/matrices/underworld32.gz\n\
      4: and run with -f underworld32.gz\n\n";

      6: #include <petscksp.h>
      7: #include <petscdmda.h>

     11: PetscErrorCode LSCLoadTestOperators(Mat *A11,Mat *A12,Mat *A21,Mat *A22,Vec *b1,Vec *b2)
     12: {
     13:   PetscViewer    viewer;
     15:   char           filename[PETSC_MAX_PATH_LEN];
     16:   PetscBool      flg;

     19:   MatCreate(PETSC_COMM_WORLD,A11);
     20:   MatCreate(PETSC_COMM_WORLD,A12);
     21:   MatCreate(PETSC_COMM_WORLD,A21);
     22:   MatCreate(PETSC_COMM_WORLD,A22);
     23:   MatSetOptionsPrefix(*A11,"a11_");
     24:   MatSetOptionsPrefix(*A22,"a22_");
     25:   MatSetFromOptions(*A11);
     26:   MatSetFromOptions(*A22);
     27:   VecCreate(PETSC_COMM_WORLD,b1);
     28:   VecCreate(PETSC_COMM_WORLD,b2);
     29:   /* Load matrices from a Q1-P0 discretisation of variable viscosity Stokes. The matrix blocks are packed into one file. */
     30:   PetscOptionsGetString(NULL,"-f",filename,sizeof(filename),&flg);
     31:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must provide a matrix file with -f");
     32:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
     33:   MatLoad(*A11,viewer);
     34:   MatLoad(*A12,viewer);
     35:   MatLoad(*A21,viewer);
     36:   MatLoad(*A22,viewer);
     37:   VecLoad(*b1,viewer);
     38:   VecLoad(*b2,viewer);
     39:   PetscViewerDestroy(&viewer);
     40:   return(0);
     41: }

     45: PetscErrorCode LoadTestMatrices(Mat *_A,Vec *_x,Vec *_b,IS *_isu,IS *_isp)
     46: {
     47:   Vec            f,h,x,b,bX[2];
     48:   Mat            A,Auu,Aup,Apu,App,bA[2][2];
     49:   IS             is_u,is_p,bis[2];
     50:   PetscInt       lnu,lnp,nu,np,i,start_u,end_u,start_p,end_p;
     51:   VecScatter     *vscat;
     52:   PetscMPIInt    rank;

     56:   /* fetch test matrices and vectors */
     57:   LSCLoadTestOperators(&Auu,&Aup,&Apu,&App,&f,&h);

     59:   /* build the mat-nest */
     60:   VecGetSize(f,&nu);
     61:   VecGetSize(h,&np);

     63:   VecGetLocalSize(f,&lnu);
     64:   VecGetLocalSize(h,&lnp);

     66:   VecGetOwnershipRange(f,&start_u,&end_u);
     67:   VecGetOwnershipRange(h,&start_p,&end_p);

     69:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     70:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] lnu = %D | lnp = %D \n", rank, lnu, lnp);
     71:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] s_u = %D | e_u = %D \n", rank, start_u, end_u);
     72:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] s_p = %D | e_p = %D \n", rank, start_p, end_p);
     73:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] is_u (offset) = %D \n", rank, start_u+start_p);
     74:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] is_p (offset) = %D \n", rank, start_u+start_p+lnu);
     75:   PetscSynchronizedFlush(PETSC_COMM_WORLD);

     77:   ISCreateStride(PETSC_COMM_WORLD,lnu,start_u+start_p,1,&is_u);
     78:   ISCreateStride(PETSC_COMM_WORLD,lnp,start_u+start_p+lnu,1,&is_p);

     80:   bis[0]   = is_u; bis[1]   = is_p;
     81:   bA[0][0] = Auu;  bA[0][1] = Aup;
     82:   bA[1][0] = Apu;  bA[1][1] = App;
     83:   MatCreateNest(PETSC_COMM_WORLD,2,bis,2,bis,&bA[0][0],&A);
     84:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     85:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

     87:   /* Pull f,h into b */
     88:   MatGetVecs(A,&b,&x);
     89:   bX[0] = f;  bX[1] = h;
     90:   PetscMalloc(sizeof(VecScatter)*2,&vscat);
     91:   for (i=0; i<2; i++) {
     92:     VecScatterCreate(b,bis[i],bX[i],NULL,&vscat[i]);
     93:     VecScatterBegin(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);
     94:   }
     95:   for (i=0; i<2; i++) {
     96:     VecScatterEnd(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);
     97:   }

     99:   /* tidy up */
    100:   for (i=0; i<2; i++) {
    101:     VecScatterDestroy(&vscat[i]);
    102:   }
    103:   PetscFree(vscat);
    104:   MatDestroy(&Auu);
    105:   MatDestroy(&Aup);
    106:   MatDestroy(&Apu);
    107:   MatDestroy(&App);
    108:   VecDestroy(&f);
    109:   VecDestroy(&h);

    111:   *_isu = is_u;
    112:   *_isp = is_p;
    113:   *_A   = A;
    114:   *_x   = x;
    115:   *_b   = b;
    116:   return(0);
    117: }

    122: PetscErrorCode port_lsd_bfbt(void)
    123: {
    124:   Mat            A;
    125:   Vec            x,b;
    126:   KSP            ksp_A;
    127:   PC             pc_A;
    128:   IS             isu,isp;

    132:   LoadTestMatrices(&A,&x,&b,&isu,&isp);

    134:   KSPCreate(PETSC_COMM_WORLD,&ksp_A);
    135:   KSPSetOptionsPrefix(ksp_A,"fc_");
    136:   KSPSetOperators(ksp_A,A,A,SAME_NONZERO_PATTERN);

    138:   KSPGetPC(ksp_A,&pc_A);
    139:   PCSetType(pc_A,PCFIELDSPLIT);
    140:   PCFieldSplitSetBlockSize(pc_A,2);
    141:   PCFieldSplitSetIS(pc_A,"velocity",isu);
    142:   PCFieldSplitSetIS(pc_A,"pressure",isp);

    144:   KSPSetFromOptions(ksp_A);
    145:   KSPSolve(ksp_A,b,x);

    147:   /* Pull u,p out of x */
    148:   {
    149:     PetscInt    loc;
    150:     PetscReal   max,norm;
    151:     PetscScalar sum;
    152:     Vec         uvec,pvec;
    153:     VecScatter  uscat,pscat;
    154:     Mat         A11,A22;

    156:     /* grab matrices and create the compatable u,p vectors */
    157:     MatGetSubMatrix(A,isu,isu,MAT_INITIAL_MATRIX,&A11);
    158:     MatGetSubMatrix(A,isp,isp,MAT_INITIAL_MATRIX,&A22);

    160:     MatGetVecs(A11,&uvec,NULL);
    161:     MatGetVecs(A22,&pvec,NULL);

    163:     /* perform the scatter from x -> (u,p) */
    164:     VecScatterCreate(x,isu,uvec,NULL,&uscat);
    165:     VecScatterBegin(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);
    166:     VecScatterEnd(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);

    168:     VecScatterCreate(x,isp,pvec,NULL,&pscat);
    169:     VecScatterBegin(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);
    170:     VecScatterEnd(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);

    172:     PetscPrintf(PETSC_COMM_WORLD,"-- vector vector values --\n");
    173:     VecMin(uvec,&loc,&max);
    174:     PetscPrintf(PETSC_COMM_WORLD,"  Min(u)  = %1.6F [loc=%d]\n",max,loc);
    175:     VecMax(uvec,&loc,&max);
    176:     PetscPrintf(PETSC_COMM_WORLD,"  Max(u)  = %1.6F [loc=%d]\n",max,loc);
    177:     VecNorm(uvec,NORM_2,&norm);
    178:     PetscPrintf(PETSC_COMM_WORLD,"  Norm(u) = %1.6F \n",norm);
    179:     VecSum(uvec,&sum);
    180:     PetscPrintf(PETSC_COMM_WORLD,"  Sum(u)  = %1.6F \n",PetscRealPart(sum));

    182:     PetscPrintf(PETSC_COMM_WORLD,"-- pressure vector values --\n");
    183:     VecMin(pvec,&loc,&max);
    184:     PetscPrintf(PETSC_COMM_WORLD,"  Min(p)  = %1.6F [loc=%d]\n",max,loc);
    185:     VecMax(pvec,&loc,&max);
    186:     PetscPrintf(PETSC_COMM_WORLD,"  Max(p)  = %1.6F [loc=%d]\n",max,loc);
    187:     VecNorm(pvec,NORM_2,&norm);
    188:     PetscPrintf(PETSC_COMM_WORLD,"  Norm(p) = %1.6F \n",norm);
    189:     VecSum(pvec,&sum);
    190:     PetscPrintf(PETSC_COMM_WORLD,"  Sum(p)  = %1.6F \n",PetscRealPart(sum));

    192:     PetscPrintf(PETSC_COMM_WORLD,"-- Full vector values --\n");
    193:     VecMin(x,&loc,&max);
    194:     PetscPrintf(PETSC_COMM_WORLD,"  Min(u,p)  = %1.6F [loc=%d]\n",max,loc);
    195:     VecMax(x,&loc,&max);
    196:     PetscPrintf(PETSC_COMM_WORLD,"  Max(u,p)  = %1.6F [loc=%d]\n",max,loc);
    197:     VecNorm(x,NORM_2,&norm);
    198:     PetscPrintf(PETSC_COMM_WORLD,"  Norm(u,p) = %1.6F \n",norm);
    199:     VecSum(x,&sum);
    200:     PetscPrintf(PETSC_COMM_WORLD,"  Sum(u,p)  = %1.6F \n",PetscRealPart(sum));

    202:     VecScatterDestroy(&uscat);
    203:     VecScatterDestroy(&pscat);
    204:     VecDestroy(&uvec);
    205:     VecDestroy(&pvec);
    206:     MatDestroy(&A11);
    207:     MatDestroy(&A22);
    208:   }

    210:   KSPDestroy(&ksp_A);
    211:   MatDestroy(&A);
    212:   VecDestroy(&x);
    213:   VecDestroy(&b);
    214:   ISDestroy(&isu);
    215:   ISDestroy(&isp);
    216:   return(0);
    217: }

    222: int main(int argc,char **argv)
    223: {

    226:   PetscInitialize(&argc,&argv,0,help);
    227:   port_lsd_bfbt();
    228:   PetscFinalize();
    229:   return 0;
    230: }