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/Datafiles/matrices/underworld32.gz\n\
  4: and run with -f underworld32.gz\n\n";

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

  9: static PetscErrorCode replace_submats(Mat A, IS isu, IS isp)
 10: {
 12:   Mat            A11,A22,A12,A21;
 13:   Mat            nA11,nA22,nA12,nA21;
 14:   const char     *prefix;

 17:   MatCreateSubMatrix(A,isu,isu,MAT_INITIAL_MATRIX,&A11);
 18:   MatCreateSubMatrix(A,isu,isp,MAT_INITIAL_MATRIX,&A12);
 19:   MatCreateSubMatrix(A,isp,isu,MAT_INITIAL_MATRIX,&A21);
 20:   MatCreateSubMatrix(A,isp,isp,MAT_INITIAL_MATRIX,&A22);
 21:   MatDuplicate(A11,MAT_COPY_VALUES,&nA11);
 22:   MatDuplicate(A12,MAT_COPY_VALUES,&nA12);
 23:   MatDuplicate(A21,MAT_COPY_VALUES,&nA21);
 24:   MatDuplicate(A22,MAT_COPY_VALUES,&nA22);
 25:   MatGetOptionsPrefix(A11,&prefix);
 26:   MatSetOptionsPrefix(nA11,prefix);
 27:   MatGetOptionsPrefix(A22,&prefix);
 28:   MatSetOptionsPrefix(nA22,prefix);
 29:   MatNestSetSubMat(A,0,0,nA11);
 30:   MatNestSetSubMat(A,0,1,nA12);
 31:   MatNestSetSubMat(A,1,0,nA21);
 32:   MatNestSetSubMat(A,1,1,nA22);
 33:   MatDestroy(&A11);
 34:   MatDestroy(&A12);
 35:   MatDestroy(&A21);
 36:   MatDestroy(&A22);
 37:   MatDestroy(&nA11);
 38:   MatDestroy(&nA12);
 39:   MatDestroy(&nA21);
 40:   MatDestroy(&nA22);
 41:   return(0);
 42: }

 44: PetscErrorCode LSCLoadTestOperators(Mat *A11,Mat *A12,Mat *A21,Mat *A22,Vec *b1,Vec *b2)
 45: {
 46:   PetscViewer    viewer;
 48:   char           filename[PETSC_MAX_PATH_LEN];
 49:   PetscBool      flg;

 52:   MatCreate(PETSC_COMM_WORLD,A11);
 53:   MatCreate(PETSC_COMM_WORLD,A12);
 54:   MatCreate(PETSC_COMM_WORLD,A21);
 55:   MatCreate(PETSC_COMM_WORLD,A22);
 56:   MatSetOptionsPrefix(*A11,"a11_");
 57:   MatSetOptionsPrefix(*A22,"a22_");
 58:   MatSetFromOptions(*A11);
 59:   MatSetFromOptions(*A22);
 60:   VecCreate(PETSC_COMM_WORLD,b1);
 61:   VecCreate(PETSC_COMM_WORLD,b2);
 62:   /* Load matrices from a Q1-P0 discretisation of variable viscosity Stokes. The matrix blocks are packed into one file. */
 63:   PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
 64:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must provide a matrix file with -f");
 65:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 66:   MatLoad(*A11,viewer);
 67:   MatLoad(*A12,viewer);
 68:   MatLoad(*A21,viewer);
 69:   MatLoad(*A22,viewer);
 70:   VecLoad(*b1,viewer);
 71:   VecLoad(*b2,viewer);
 72:   PetscViewerDestroy(&viewer);
 73:   return(0);
 74: }

 76: PetscErrorCode LoadTestMatrices(Mat *_A,Vec *_x,Vec *_b,IS *_isu,IS *_isp)
 77: {
 78:   Vec            f,h,x,b,bX[2];
 79:   Mat            A,Auu,Aup,Apu,App,bA[2][2];
 80:   IS             is_u,is_p,bis[2];
 81:   PetscInt       lnu,lnp,nu,np,i,start_u,end_u,start_p,end_p;
 82:   VecScatter     *vscat;
 83:   PetscMPIInt    rank;

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

 90:   /* build the mat-nest */
 91:   VecGetSize(f,&nu);
 92:   VecGetSize(h,&np);

 94:   VecGetLocalSize(f,&lnu);
 95:   VecGetLocalSize(h,&lnp);

 97:   VecGetOwnershipRange(f,&start_u,&end_u);
 98:   VecGetOwnershipRange(h,&start_p,&end_p);

100:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
101:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] lnu = %D | lnp = %D \n", rank, lnu, lnp);
102:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] s_u = %D | e_u = %D \n", rank, start_u, end_u);
103:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] s_p = %D | e_p = %D \n", rank, start_p, end_p);
104:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] is_u (offset) = %D \n", rank, start_u+start_p);
105:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] is_p (offset) = %D \n", rank, start_u+start_p+lnu);
106:   PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);

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

111:   bis[0]   = is_u; bis[1]   = is_p;
112:   bA[0][0] = Auu;  bA[0][1] = Aup;
113:   bA[1][0] = Apu;  bA[1][1] = App;
114:   MatCreateNest(PETSC_COMM_WORLD,2,bis,2,bis,&bA[0][0],&A);
115:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
116:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

118:   /* Pull f,h into b */
119:   MatCreateVecs(A,&b,&x);
120:   bX[0] = f;  bX[1] = h;
121:   PetscMalloc1(2,&vscat);
122:   for (i=0; i<2; i++) {
123:     VecScatterCreate(b,bis[i],bX[i],NULL,&vscat[i]);
124:     VecScatterBegin(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);
125:     VecScatterEnd(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);
126:     VecScatterDestroy(&vscat[i]);
127:   }

129:   PetscFree(vscat);
130:   MatDestroy(&Auu);
131:   MatDestroy(&Aup);
132:   MatDestroy(&Apu);
133:   MatDestroy(&App);
134:   VecDestroy(&f);
135:   VecDestroy(&h);

137:   *_isu = is_u;
138:   *_isp = is_p;
139:   *_A   = A;
140:   *_x   = x;
141:   *_b   = b;
142:   return(0);
143: }

145: PetscErrorCode port_lsd_bfbt(void)
146: {
147:   Mat            A,P;
148:   Vec            x,b;
149:   KSP            ksp_A;
150:   PC             pc_A;
151:   IS             isu,isp;
152:   PetscBool      test_fs = PETSC_FALSE;

156:   LoadTestMatrices(&A,&x,&b,&isu,&isp);
157:   PetscOptionsGetBool(NULL,NULL,"-test_fs",&test_fs,NULL);
158:   if (!test_fs) {
159:     PetscObjectReference((PetscObject)A);
160:     P = A;
161:   } else {
162:     MatDuplicate(A,MAT_COPY_VALUES,&P);
163:   }
164:   KSPCreate(PETSC_COMM_WORLD,&ksp_A);
165:   KSPSetOptionsPrefix(ksp_A,"fc_");
166:   KSPSetOperators(ksp_A,A,P);

168:   KSPSetFromOptions(ksp_A);
169:   KSPGetPC(ksp_A,&pc_A);
170:   PetscObjectTypeCompare((PetscObject)pc_A,PCLU,&test_fs);
171:   if (test_fs) {
172:     MatDestroy(&P);
173:     MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&P);
174:     KSPSetOperators(ksp_A,A,P);
175:     KSPSetFromOptions(ksp_A);
176:     KSPSolve(ksp_A,b,x);
177:   } else {
178:     PCFieldSplitSetBlockSize(pc_A,2);
179:     PCFieldSplitSetIS(pc_A,"velocity",isu);
180:     PCFieldSplitSetIS(pc_A,"pressure",isp);
181:     KSPSolve(ksp_A,b,x);

183:     /* Pull u,p out of x */
184:     {
185:       PetscInt    loc;
186:       PetscReal   max,norm;
187:       PetscScalar sum;
188:       Vec         uvec,pvec;
189:       VecScatter  uscat,pscat;
190:       Mat         A11,A22;

192:       /* grab matrices and create the compatable u,p vectors */
193:       MatCreateSubMatrix(A,isu,isu,MAT_INITIAL_MATRIX,&A11);
194:       MatCreateSubMatrix(A,isp,isp,MAT_INITIAL_MATRIX,&A22);

196:       MatCreateVecs(A11,&uvec,NULL);
197:       MatCreateVecs(A22,&pvec,NULL);

199:       /* perform the scatter from x -> (u,p) */
200:       VecScatterCreate(x,isu,uvec,NULL,&uscat);
201:       VecScatterBegin(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);
202:       VecScatterEnd(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);

204:       VecScatterCreate(x,isp,pvec,NULL,&pscat);
205:       VecScatterBegin(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);
206:       VecScatterEnd(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);

208:       PetscPrintf(PETSC_COMM_WORLD,"-- vector vector values --\n");
209:       VecMin(uvec,&loc,&max);
210:       PetscPrintf(PETSC_COMM_WORLD,"  Min(u)  = %1.6f [loc=%D]\n",(double)max,loc);
211:       VecMax(uvec,&loc,&max);
212:       PetscPrintf(PETSC_COMM_WORLD,"  Max(u)  = %1.6f [loc=%D]\n",(double)max,loc);
213:       VecNorm(uvec,NORM_2,&norm);
214:       PetscPrintf(PETSC_COMM_WORLD,"  Norm(u) = %1.6f \n",(double)norm);
215:       VecSum(uvec,&sum);
216:       PetscPrintf(PETSC_COMM_WORLD,"  Sum(u)  = %1.6f \n",(double)PetscRealPart(sum));

218:       PetscPrintf(PETSC_COMM_WORLD,"-- pressure vector values --\n");
219:       VecMin(pvec,&loc,&max);
220:       PetscPrintf(PETSC_COMM_WORLD,"  Min(p)  = %1.6f [loc=%D]\n",(double)max,loc);
221:       VecMax(pvec,&loc,&max);
222:       PetscPrintf(PETSC_COMM_WORLD,"  Max(p)  = %1.6f [loc=%D]\n",(double)max,loc);
223:       VecNorm(pvec,NORM_2,&norm);
224:       PetscPrintf(PETSC_COMM_WORLD,"  Norm(p) = %1.6f \n",(double)norm);
225:       VecSum(pvec,&sum);
226:       PetscPrintf(PETSC_COMM_WORLD,"  Sum(p)  = %1.6f \n",(double)PetscRealPart(sum));

228:       PetscPrintf(PETSC_COMM_WORLD,"-- Full vector values --\n");
229:       VecMin(x,&loc,&max);
230:       PetscPrintf(PETSC_COMM_WORLD,"  Min(u,p)  = %1.6f [loc=%D]\n",(double)max,loc);
231:       VecMax(x,&loc,&max);
232:       PetscPrintf(PETSC_COMM_WORLD,"  Max(u,p)  = %1.6f [loc=%D]\n",(double)max,loc);
233:       VecNorm(x,NORM_2,&norm);
234:       PetscPrintf(PETSC_COMM_WORLD,"  Norm(u,p) = %1.6f \n",(double)norm);
235:       VecSum(x,&sum);
236:       PetscPrintf(PETSC_COMM_WORLD,"  Sum(u,p)  = %1.6f \n",(double)PetscRealPart(sum));

238:       VecScatterDestroy(&uscat);
239:       VecScatterDestroy(&pscat);
240:       VecDestroy(&uvec);
241:       VecDestroy(&pvec);
242:       MatDestroy(&A11);
243:       MatDestroy(&A22);
244:     }

246:     /* test second solve by changing the mat associated to the MATNEST blocks */
247:     {
248:       replace_submats(A,isu,isp);
249:       replace_submats(P,isu,isp);
250:       KSPSolve(ksp_A,b,x);
251:     }
252:   }

254:   KSPDestroy(&ksp_A);
255:   MatDestroy(&P);
256:   MatDestroy(&A);
257:   VecDestroy(&x);
258:   VecDestroy(&b);
259:   ISDestroy(&isu);
260:   ISDestroy(&isp);
261:   return(0);
262: }

264: int main(int argc,char **argv)
265: {

268:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
269:   port_lsd_bfbt();
270:   PetscFinalize();
271:   return ierr;
272: }

274: /*TEST

276:     test:
277:       args: -f ${DATAFILESPATH}/matrices/underworld32.gz -fc_ksp_view -fc_ksp_monitor_short -fc_ksp_type fgmres -fc_ksp_max_it 4000 -fc_ksp_diagonal_scale -fc_pc_type fieldsplit -fc_pc_fieldsplit_type SCHUR -fc_pc_fieldsplit_schur_fact_type UPPER -fc_fieldsplit_velocity_ksp_type cg -fc_fieldsplit_velocity_pc_type cholesky -fc_fieldsplit_velocity_pc_factor_mat_ordering_type nd -fc_fieldsplit_pressure_ksp_max_it 100 -fc_fieldsplit_pressure_ksp_constant_null_space -fc_fieldsplit_pressure_ksp_monitor_short -fc_fieldsplit_pressure_pc_type lsc -fc_fieldsplit_pressure_lsc_ksp_type cg -fc_fieldsplit_pressure_lsc_ksp_max_it 100 -fc_fieldsplit_pressure_lsc_ksp_constant_null_space -fc_fieldsplit_pressure_lsc_ksp_converged_reason -fc_fieldsplit_pressure_lsc_pc_type icc -test_fs {{0 1}separate output} -fc_pc_fieldsplit_off_diag_use_amat {{0 1}separate output} -fc_pc_fieldsplit_diag_use_amat {{0 1}separate output}
278:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)

280:     test:
281:       suffix: 2
282:       nsize: 4
283:       args: -f ${DATAFILESPATH}/matrices/underworld32.gz -fc_ksp_view -fc_ksp_monitor_short -fc_ksp_type fgmres -fc_ksp_max_it 4000 -fc_ksp_diagonal_scale -fc_pc_type fieldsplit -fc_pc_fieldsplit_type SCHUR -fc_pc_fieldsplit_schur_fact_type UPPER -fc_fieldsplit_velocity_ksp_type cg -fc_fieldsplit_velocity_ksp_rtol 1.0e-6 -fc_fieldsplit_velocity_pc_type bjacobi -fc_fieldsplit_velocity_sub_pc_type cholesky -fc_fieldsplit_velocity_sub_pc_factor_mat_ordering_type nd -fc_fieldsplit_pressure_ksp_type fgmres -fc_fieldsplit_pressure_ksp_constant_null_space -fc_fieldsplit_pressure_ksp_monitor_short -fc_fieldsplit_pressure_pc_type lsc -fc_fieldsplit_pressure_lsc_ksp_type cg -fc_fieldsplit_pressure_lsc_ksp_rtol 1.0e-2 -fc_fieldsplit_pressure_lsc_ksp_constant_null_space -fc_fieldsplit_pressure_lsc_ksp_converged_reason -fc_fieldsplit_pressure_lsc_pc_type bjacobi -fc_fieldsplit_pressure_lsc_sub_pc_type icc -test_fs {{0 1}separate output} -fc_pc_fieldsplit_off_diag_use_amat {{0 1}separate output} -fc_pc_fieldsplit_diag_use_amat {{0 1}separate output}
284:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)

286:     test:
287:       suffix: 3
288:       nsize: 2
289:       args: -f ${DATAFILESPATH}/matrices/underworld32.gz -fc_ksp_view_pre -fc_pc_type lu
290:       requires: datafilespath mumps double !complex !defined(PETSC_USE_64BIT_INDICES)

292: TEST*/