Actual source code: ex11.c

petsc-3.10.5 2019-03-28
Report Typos and Errors
  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: PetscErrorCode LSCLoadTestOperators(Mat *A11,Mat *A12,Mat *A21,Mat *A22,Vec *b1,Vec *b2)
 10: {
 11:   PetscViewer    viewer;
 13:   char           filename[PETSC_MAX_PATH_LEN];
 14:   PetscBool      flg;

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

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

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

 55:   /* build the mat-nest */
 56:   VecGetSize(f,&nu);
 57:   VecGetSize(h,&np);

 59:   VecGetLocalSize(f,&lnu);
 60:   VecGetLocalSize(h,&lnp);

 62:   VecGetOwnershipRange(f,&start_u,&end_u);
 63:   VecGetOwnershipRange(h,&start_p,&end_p);

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

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

 76:   bis[0]   = is_u; bis[1]   = is_p;
 77:   bA[0][0] = Auu;  bA[0][1] = Aup;
 78:   bA[1][0] = Apu;  bA[1][1] = App;
 79:   MatCreateNest(PETSC_COMM_WORLD,2,bis,2,bis,&bA[0][0],&A);
 80:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 81:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 83:   /* Pull f,h into b */
 84:   MatCreateVecs(A,&b,&x);
 85:   bX[0] = f;  bX[1] = h;
 86:   PetscMalloc1(2,&vscat);
 87:   for (i=0; i<2; i++) {
 88:     VecScatterCreate(b,bis[i],bX[i],NULL,&vscat[i]);
 89:     VecScatterBegin(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);
 90:   }
 91:   for (i=0; i<2; i++) {
 92:     VecScatterEnd(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);
 93:   }

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

107:   *_isu = is_u;
108:   *_isp = is_p;
109:   *_A   = A;
110:   *_x   = x;
111:   *_b   = b;
112:   return(0);
113: }


116: PetscErrorCode port_lsd_bfbt(void)
117: {
118:   Mat            A;
119:   Vec            x,b;
120:   KSP            ksp_A;
121:   PC             pc_A;
122:   IS             isu,isp;

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

128:   KSPCreate(PETSC_COMM_WORLD,&ksp_A);
129:   KSPSetOptionsPrefix(ksp_A,"fc_");
130:   KSPSetOperators(ksp_A,A,A);

132:   KSPGetPC(ksp_A,&pc_A);
133:   PCSetType(pc_A,PCFIELDSPLIT);
134:   PCFieldSplitSetBlockSize(pc_A,2);
135:   PCFieldSplitSetIS(pc_A,"velocity",isu);
136:   PCFieldSplitSetIS(pc_A,"pressure",isp);

138:   KSPSetFromOptions(ksp_A);
139:   KSPSolve(ksp_A,b,x);

141:   /* Pull u,p out of x */
142:   {
143:     PetscInt    loc;
144:     PetscReal   max,norm;
145:     PetscScalar sum;
146:     Vec         uvec,pvec;
147:     VecScatter  uscat,pscat;
148:     Mat         A11,A22;

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

154:     MatCreateVecs(A11,&uvec,NULL);
155:     MatCreateVecs(A22,&pvec,NULL);

157:     /* perform the scatter from x -> (u,p) */
158:     VecScatterCreate(x,isu,uvec,NULL,&uscat);
159:     VecScatterBegin(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);
160:     VecScatterEnd(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);

162:     VecScatterCreate(x,isp,pvec,NULL,&pscat);
163:     VecScatterBegin(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);
164:     VecScatterEnd(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);

166:     PetscPrintf(PETSC_COMM_WORLD,"-- vector vector values --\n");
167:     VecMin(uvec,&loc,&max);
168:     PetscPrintf(PETSC_COMM_WORLD,"  Min(u)  = %1.6f [loc=%D]\n",(double)max,loc);
169:     VecMax(uvec,&loc,&max);
170:     PetscPrintf(PETSC_COMM_WORLD,"  Max(u)  = %1.6f [loc=%D]\n",(double)max,loc);
171:     VecNorm(uvec,NORM_2,&norm);
172:     PetscPrintf(PETSC_COMM_WORLD,"  Norm(u) = %1.6f \n",(double)norm);
173:     VecSum(uvec,&sum);
174:     PetscPrintf(PETSC_COMM_WORLD,"  Sum(u)  = %1.6f \n",(double)PetscRealPart(sum));

176:     PetscPrintf(PETSC_COMM_WORLD,"-- pressure vector values --\n");
177:     VecMin(pvec,&loc,&max);
178:     PetscPrintf(PETSC_COMM_WORLD,"  Min(p)  = %1.6f [loc=%D]\n",(double)max,loc);
179:     VecMax(pvec,&loc,&max);
180:     PetscPrintf(PETSC_COMM_WORLD,"  Max(p)  = %1.6f [loc=%D]\n",(double)max,loc);
181:     VecNorm(pvec,NORM_2,&norm);
182:     PetscPrintf(PETSC_COMM_WORLD,"  Norm(p) = %1.6f \n",(double)norm);
183:     VecSum(pvec,&sum);
184:     PetscPrintf(PETSC_COMM_WORLD,"  Sum(p)  = %1.6f \n",(double)PetscRealPart(sum));

186:     PetscPrintf(PETSC_COMM_WORLD,"-- Full vector values --\n");
187:     VecMin(x,&loc,&max);
188:     PetscPrintf(PETSC_COMM_WORLD,"  Min(u,p)  = %1.6f [loc=%D]\n",(double)max,loc);
189:     VecMax(x,&loc,&max);
190:     PetscPrintf(PETSC_COMM_WORLD,"  Max(u,p)  = %1.6f [loc=%D]\n",(double)max,loc);
191:     VecNorm(x,NORM_2,&norm);
192:     PetscPrintf(PETSC_COMM_WORLD,"  Norm(u,p) = %1.6f \n",(double)norm);
193:     VecSum(x,&sum);
194:     PetscPrintf(PETSC_COMM_WORLD,"  Sum(u,p)  = %1.6f \n",(double)PetscRealPart(sum));

196:     VecScatterDestroy(&uscat);
197:     VecScatterDestroy(&pscat);
198:     VecDestroy(&uvec);
199:     VecDestroy(&pvec);
200:     MatDestroy(&A11);
201:     MatDestroy(&A22);
202:   }

204:   KSPDestroy(&ksp_A);
205:   MatDestroy(&A);
206:   VecDestroy(&x);
207:   VecDestroy(&b);
208:   ISDestroy(&isu);
209:   ISDestroy(&isp);
210:   return(0);
211: }


214: int main(int argc,char **argv)
215: {

218:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
219:   port_lsd_bfbt();
220:   PetscFinalize();
221:   return ierr;
222: }

224: /*TEST

226:     test:
227:       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
228:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)

230:     test:
231:       suffix: 2
232:       nsize: 4
233:       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
234:       requires: datafilespath double  !complex !define(PETSC_USE_64BIT_INDICES)

236: TEST*/