Actual source code: ex22.c
2: static char help[] = "Solves PDE optimization problem.\n\n";
4: #include petscda.h
5: #include petscpf.h
6: #include petscsnes.h
7: #include petscdmmg.h
9: /*
11: w - design variables (what we change to get an optimal solution)
12: u - state variables (i.e. the PDE solution)
13: lambda - the Lagrange multipliers
15: U = (w [u_0 lambda_0 u_1 lambda_1 .....])
17: fu, fw, flambda contain the gradient of L(w,u,lambda)
19: FU = (fw [fu_0 flambda_0 .....])
21: In this example the PDE is
22: Uxx = 2,
23: u(0) = w(0), thus this is the free parameter
24: u(1) = 0
25: the function we wish to minimize is
26: \integral u^{2}
28: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
30: Use the usual centered finite differences.
32: Note we treat the problem as non-linear though it happens to be linear
34: See ex21.c for the same code, but that does NOT interlaces the u and the lambda
36: The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
37: */
39: typedef struct {
40: PetscViewer u_lambda_viewer;
41: PetscViewer fu_lambda_viewer;
42: } UserCtx;
47: /*
48: Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the
49: smoother on all levels. This is because (1) in the matrix free case no matrix entries are
50: available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal
51: entry for the control variable is zero which means default SOR will not work.
53: */
54: char common_options[] = "-dmmg_grid_sequence \
55: -dmmg_nlevels 5 \
56: -mg_levels_pc_type none \
57: -mg_coarse_pc_type none \
58: -pc_mg_type full \
59: -mg_coarse_ksp_type gmres \
60: -mg_levels_ksp_type gmres \
61: -mg_coarse_ksp_max_it 6 \
62: -mg_levels_ksp_max_it 3";
64: char matrix_free_options[] = "-mat_mffd_compute_normu no \
65: -mat_mffd_type wp \
66: -dmmg_jacobian_mf_fd";
68: /*
69: Currently only global coloring is supported with DMComposite
70: */
71: char matrix_based_options[] = "-dmmg_iscoloring_type global";
73: /*
74: The -use_matrix_based version does not work! This is because the DMComposite code cannot determine the nonzero
75: pattern of the Jacobian since the coupling between the boundary condition (array variable) and DA variables is problem
76: dependent. To get the explicit Jacobian correct you would need to use the DMCompositeSetCoupling() to indicate the extra nonzero
77: pattern and run with -dmmg_coloring_from_mat.
78: */
83: int main(int argc,char **argv)
84: {
86: UserCtx user;
87: DA da;
88: DMMG *dmmg;
89: DMComposite packer;
90: PetscTruth use_matrix_based = PETSC_FALSE,use_monitor = PETSC_FALSE;
91: PetscInt i;
93: PetscInitialize(&argc,&argv,PETSC_NULL,help);
94: PetscOptionsSetFromOptions();
96: /* Hardwire several options; can be changed at command line */
97: PetscOptionsInsertString(common_options);
98: PetscOptionsGetTruth(PETSC_NULL,"-use_matrix_based",&use_matrix_based,PETSC_IGNORE);
99: if (use_matrix_based) {
100: PetscOptionsInsertString(matrix_based_options);
101: } else {
102: PetscOptionsInsertString(matrix_free_options);
103: }
104: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
105: PetscOptionsGetTruth(PETSC_NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);
107: /* Create a global vector that includes a single redundant array and two da arrays */
108: DMCompositeCreate(PETSC_COMM_WORLD,&packer);
109: DMCompositeAddArray(packer,0,1);
110: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,2,1,PETSC_NULL,&da);
111: DMCompositeAddDM(packer,(DM)da);
114: /* create nonlinear multi-level solver */
115: DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
116: DMMGSetDM(dmmg,(DM)packer);
117: DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
118: DMMGSetFromOptions(dmmg);
120: if (use_monitor) {
121: /* create graphics windows */
122: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);
123: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);
124: for (i=0; i<DMMGGetLevels(dmmg); i++) {
125: SNESMonitorSet(dmmg[i]->snes,Monitor,dmmg[i],0);
126: }
127: }
129: DMMGSolve(dmmg);
130: DMMGDestroy(dmmg);
132: DADestroy(da);
133: DMCompositeDestroy(packer);
134: if (use_monitor) {
135: PetscViewerDestroy(user.u_lambda_viewer);
136: PetscViewerDestroy(user.fu_lambda_viewer);
137: }
139: PetscFinalize();
140: return 0;
141: }
143: typedef struct {
144: PetscScalar u;
145: PetscScalar lambda;
146: } ULambda;
147:
148: /*
149: Evaluates FU = Gradiant(L(w,u,lambda))
151: This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and
152: DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()).
154: */
155: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
156: {
157: DMMG dmmg = (DMMG)dummy;
159: PetscInt xs,xm,i,N,nredundant;
160: ULambda *u_lambda,*fu_lambda;
161: PetscScalar d,h,*w,*fw;
162: Vec vu_lambda,vfu_lambda;
163: DA da;
164: DMComposite packer = (DMComposite)dmmg->dm;
167: DMCompositeGetEntries(packer,&nredundant,&da);
168: DMCompositeGetLocalVectors(packer,&w,&vu_lambda);
169: DMCompositeScatter(packer,U,w,vu_lambda);
170: DMCompositeGetAccess(packer,FU,&fw,&vfu_lambda);
172: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
173: DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
174: DAVecGetArray(da,vu_lambda,&u_lambda);
175: DAVecGetArray(da,vfu_lambda,&fu_lambda);
176: d = N-1.0;
177: h = 1.0/d;
179: /* derivative of L() w.r.t. w */
180: if (xs == 0) { /* only first processor computes this */
181: fw[0] = -2.0*d*u_lambda[0].lambda;
182: }
184: /* derivative of L() w.r.t. u */
185: for (i=xs; i<xs+xm; i++) {
186: if (i == 0) fu_lambda[0].lambda = h*u_lambda[0].u + 2.*d*u_lambda[0].lambda - d*u_lambda[1].lambda;
187: else if (i == 1) fu_lambda[1].lambda = 2.*h*u_lambda[1].u + 2.*d*u_lambda[1].lambda - d*u_lambda[2].lambda;
188: else if (i == N-1) fu_lambda[N-1].lambda = h*u_lambda[N-1].u + 2.*d*u_lambda[N-1].lambda - d*u_lambda[N-2].lambda;
189: else if (i == N-2) fu_lambda[N-2].lambda = 2.*h*u_lambda[N-2].u + 2.*d*u_lambda[N-2].lambda - d*u_lambda[N-3].lambda;
190: else fu_lambda[i].lambda = 2.*h*u_lambda[i].u - d*(u_lambda[i+1].lambda - 2.0*u_lambda[i].lambda + u_lambda[i-1].lambda);
191: }
193: /* derivative of L() w.r.t. lambda */
194: for (i=xs; i<xs+xm; i++) {
195: if (i == 0) fu_lambda[0].u = 2.0*d*(u_lambda[0].u - w[0]);
196: else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;
197: else fu_lambda[i].u = -(d*(u_lambda[i+1].u - 2.0*u_lambda[i].u + u_lambda[i-1].u) - 2.0*h);
198: }
200: DAVecRestoreArray(da,vu_lambda,&u_lambda);
201: DAVecRestoreArray(da,vfu_lambda,&fu_lambda);
202: DMCompositeRestoreLocalVectors(packer,&w,&vu_lambda);
203: DMCompositeRestoreAccess(packer,FU,&fw,&vfu_lambda);
204: PetscLogFlops(13.0*N);
205: return(0);
206: }
208: /*
209: Computes the exact solution
210: */
211: PetscErrorCode u_solution(void *dummy,PetscInt n,PetscScalar *x,PetscScalar *u)
212: {
213: PetscInt i;
215: for (i=0; i<n; i++) {
216: u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
217: }
218: return(0);
219: }
221: PetscErrorCode ExactSolution(DMComposite packer,Vec U)
222: {
223: PF pf;
224: Vec x,u_global;
225: PetscScalar *w;
226: DA da;
228: PetscInt m;
231: DMCompositeGetEntries(packer,&m,&da);
233: PFCreate(PETSC_COMM_WORLD,1,2,&pf);
234: PFSetType(pf,PFQUICK,(void*)u_solution);
235: DAGetCoordinates(da,&x);
236: if (!x) {
237: DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
238: DAGetCoordinates(da,&x);
239: }
240: DMCompositeGetAccess(packer,U,&w,&u_global,0);
241: if (w) w[0] = .25;
242: PFApplyVec(pf,x,u_global);
243: PFDestroy(pf);
244: VecDestroy(x);
245: DMCompositeRestoreAccess(packer,U,&w,&u_global,0);
246: return(0);
247: }
250: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
251: {
252: DMMG dmmg = (DMMG)dummy;
253: UserCtx *user = (UserCtx*)dmmg->user;
255: PetscInt m,N;
256: PetscScalar *w,*dw;
257: Vec u_lambda,U,F,Uexact;
258: DMComposite packer = (DMComposite)dmmg->dm;
259: PetscReal norm;
260: DA da;
263: SNESGetSolution(snes,&U);
264: DMCompositeGetAccess(packer,U,&w,&u_lambda);
265: VecView(u_lambda,user->u_lambda_viewer);
266: DMCompositeRestoreAccess(packer,U,&w,&u_lambda);
268: SNESGetFunction(snes,&F,0,0);
269: DMCompositeGetAccess(packer,F,&w,&u_lambda);
270: /* VecView(u_lambda,user->fu_lambda_viewer); */
271: DMCompositeRestoreAccess(packer,U,&w,&u_lambda);
273: DMCompositeGetEntries(packer,&m,&da);
274: DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
275: VecDuplicate(U,&Uexact);
276: ExactSolution(packer,Uexact);
277: VecAXPY(Uexact,-1.0,U);
278: DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda);
279: VecStrideNorm(u_lambda,0,NORM_2,&norm);
280: norm = norm/sqrt(N-1.);
281: if (dw) PetscPrintf(dmmg->comm,"Norm of error %G Error at x = 0 %G\n",norm,PetscRealPart(dw[0]));
282: VecView(u_lambda,user->fu_lambda_viewer);
283: DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda);
284: VecDestroy(Uexact);
285: return(0);
286: }