Actual source code: ex11.c
petsc-3.12.0 2019-09-29
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: VecScatterEnd(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);
91: VecScatterDestroy(&vscat[i]);
92: }
94: PetscFree(vscat);
95: MatDestroy(&Auu);
96: MatDestroy(&Aup);
97: MatDestroy(&Apu);
98: MatDestroy(&App);
99: VecDestroy(&f);
100: VecDestroy(&h);
102: *_isu = is_u;
103: *_isp = is_p;
104: *_A = A;
105: *_x = x;
106: *_b = b;
107: return(0);
108: }
111: PetscErrorCode port_lsd_bfbt(void)
112: {
113: Mat A;
114: Vec x,b;
115: KSP ksp_A;
116: PC pc_A;
117: IS isu,isp;
121: LoadTestMatrices(&A,&x,&b,&isu,&isp);
123: KSPCreate(PETSC_COMM_WORLD,&ksp_A);
124: KSPSetOptionsPrefix(ksp_A,"fc_");
125: KSPSetOperators(ksp_A,A,A);
127: KSPGetPC(ksp_A,&pc_A);
128: PCSetType(pc_A,PCFIELDSPLIT);
129: PCFieldSplitSetBlockSize(pc_A,2);
130: PCFieldSplitSetIS(pc_A,"velocity",isu);
131: PCFieldSplitSetIS(pc_A,"pressure",isp);
133: KSPSetFromOptions(ksp_A);
134: KSPSolve(ksp_A,b,x);
136: /* Pull u,p out of x */
137: {
138: PetscInt loc;
139: PetscReal max,norm;
140: PetscScalar sum;
141: Vec uvec,pvec;
142: VecScatter uscat,pscat;
143: Mat A11,A22;
145: /* grab matrices and create the compatable u,p vectors */
146: MatCreateSubMatrix(A,isu,isu,MAT_INITIAL_MATRIX,&A11);
147: MatCreateSubMatrix(A,isp,isp,MAT_INITIAL_MATRIX,&A22);
149: MatCreateVecs(A11,&uvec,NULL);
150: MatCreateVecs(A22,&pvec,NULL);
152: /* perform the scatter from x -> (u,p) */
153: VecScatterCreate(x,isu,uvec,NULL,&uscat);
154: VecScatterBegin(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);
155: VecScatterEnd(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);
157: VecScatterCreate(x,isp,pvec,NULL,&pscat);
158: VecScatterBegin(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);
159: VecScatterEnd(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);
161: PetscPrintf(PETSC_COMM_WORLD,"-- vector vector values --\n");
162: VecMin(uvec,&loc,&max);
163: PetscPrintf(PETSC_COMM_WORLD," Min(u) = %1.6f [loc=%D]\n",(double)max,loc);
164: VecMax(uvec,&loc,&max);
165: PetscPrintf(PETSC_COMM_WORLD," Max(u) = %1.6f [loc=%D]\n",(double)max,loc);
166: VecNorm(uvec,NORM_2,&norm);
167: PetscPrintf(PETSC_COMM_WORLD," Norm(u) = %1.6f \n",(double)norm);
168: VecSum(uvec,&sum);
169: PetscPrintf(PETSC_COMM_WORLD," Sum(u) = %1.6f \n",(double)PetscRealPart(sum));
171: PetscPrintf(PETSC_COMM_WORLD,"-- pressure vector values --\n");
172: VecMin(pvec,&loc,&max);
173: PetscPrintf(PETSC_COMM_WORLD," Min(p) = %1.6f [loc=%D]\n",(double)max,loc);
174: VecMax(pvec,&loc,&max);
175: PetscPrintf(PETSC_COMM_WORLD," Max(p) = %1.6f [loc=%D]\n",(double)max,loc);
176: VecNorm(pvec,NORM_2,&norm);
177: PetscPrintf(PETSC_COMM_WORLD," Norm(p) = %1.6f \n",(double)norm);
178: VecSum(pvec,&sum);
179: PetscPrintf(PETSC_COMM_WORLD," Sum(p) = %1.6f \n",(double)PetscRealPart(sum));
181: PetscPrintf(PETSC_COMM_WORLD,"-- Full vector values --\n");
182: VecMin(x,&loc,&max);
183: PetscPrintf(PETSC_COMM_WORLD," Min(u,p) = %1.6f [loc=%D]\n",(double)max,loc);
184: VecMax(x,&loc,&max);
185: PetscPrintf(PETSC_COMM_WORLD," Max(u,p) = %1.6f [loc=%D]\n",(double)max,loc);
186: VecNorm(x,NORM_2,&norm);
187: PetscPrintf(PETSC_COMM_WORLD," Norm(u,p) = %1.6f \n",(double)norm);
188: VecSum(x,&sum);
189: PetscPrintf(PETSC_COMM_WORLD," Sum(u,p) = %1.6f \n",(double)PetscRealPart(sum));
191: VecScatterDestroy(&uscat);
192: VecScatterDestroy(&pscat);
193: VecDestroy(&uvec);
194: VecDestroy(&pvec);
195: MatDestroy(&A11);
196: MatDestroy(&A22);
197: }
199: KSPDestroy(&ksp_A);
200: MatDestroy(&A);
201: VecDestroy(&x);
202: VecDestroy(&b);
203: ISDestroy(&isu);
204: ISDestroy(&isp);
205: return(0);
206: }
209: int main(int argc,char **argv)
210: {
213: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
214: port_lsd_bfbt();
215: PetscFinalize();
216: return ierr;
217: }
219: /*TEST
221: test:
222: 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
223: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
225: test:
226: suffix: 2
227: nsize: 4
228: 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
229: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
231: TEST*/