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*/