Actual source code: ex11.c
petsc-3.3-p6 2013-02-11
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: }