Actual source code: ex11.c

  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(PETSC_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;


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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

161:     MatGetVecs(A11,&uvec,PETSC_NULL);
162:     MatGetVecs(A22,&pvec,PETSC_NULL);

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

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

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

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

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

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

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


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

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