Actual source code: ex22.c
2: static char help[] = "Solves PDE optimization problem.\n\n";
4: #include <petscdmda.h>
5: #include <petscdmcomposite.h>
6: #include <petscpf.h>
7: #include <petscsnes.h>
8: #include <petscdmmg.h>
10: /*
12: w - design variables (what we change to get an optimal solution)
13: u - state variables (i.e. the PDE solution)
14: lambda - the Lagrange multipliers
16: U = (w [u_0 lambda_0 u_1 lambda_1 .....])
18: fu, fw, flambda contain the gradient of L(w,u,lambda)
20: FU = (fw [fu_0 flambda_0 .....])
22: In this example the PDE is
23: Uxx = 2,
24: u(0) = w(0), thus this is the free parameter
25: u(1) = 0
26: the function we wish to minimize is
27: \integral u^{2}
29: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
31: Use the usual centered finite differences.
33: Note we treat the problem as non-linear though it happens to be linear
35: See ex21.c for the same code, but that does NOT interlaces the u and the lambda
37: The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
38: */
40: typedef struct {
41: PetscViewer u_lambda_viewer;
42: PetscViewer fu_lambda_viewer;
43: } UserCtx;
48: /*
49: Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the
50: smoother on all levels. This is because (1) in the matrix free case no matrix entries are
51: available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal
52: entry for the control variable is zero which means default SOR will not work.
54: */
55: char common_options[] = "-dmmg_grid_sequence \
56: -dmmg_nlevels 5 \
57: -mg_levels_pc_type none \
58: -mg_coarse_pc_type none \
59: -pc_mg_type full \
60: -mg_coarse_ksp_type gmres \
61: -mg_levels_ksp_type gmres \
62: -mg_coarse_ksp_max_it 6 \
63: -mg_levels_ksp_max_it 3";
65: char matrix_free_options[] = "-mat_mffd_compute_normu no \
66: -mat_mffd_type wp \
67: -dmmg_jacobian_mf_fd";
69: /*
70: Currently only global coloring is supported with DMComposite
71: */
72: char matrix_based_options[] = "-dmmg_iscoloring_type global \
73: -mg_coarse_ksp_type gmres \
74: -ksp_type fgmres ";
76: /*
77: The -use_matrix_based version does not work! This is because the DMComposite code cannot determine the nonzero
78: pattern of the Jacobian since the coupling between the boundary condition (array variable) and DMDA variables is problem
79: dependent. To get the explicit Jacobian correct you would need to use the DMCompositeSetCoupling() to indicate the extra nonzero
80: pattern and run with -dmmg_coloring_from_mat.
81: */
86: int main(int argc,char **argv)
87: {
89: UserCtx user;
90: DM da;
91: DMMG *dmmg;
92: DM packer;
93: PetscBool use_matrix_based = PETSC_FALSE,use_monitor = PETSC_FALSE;
94: PetscInt i;
96: PetscInitialize(&argc,&argv,PETSC_NULL,help);
97: PetscOptionsSetFromOptions();
99: /* Hardwire several options; can be changed at command line */
100: PetscOptionsInsertString(common_options);
101: PetscOptionsGetBool(PETSC_NULL,"-use_matrix_based",&use_matrix_based,PETSC_IGNORE);
102: if (use_matrix_based) {
103: PetscOptionsInsertString(matrix_based_options);
104: } else {
105: PetscOptionsInsertString(matrix_free_options);
106: }
107: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
108: PetscOptionsGetBool(PETSC_NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);
110: /* Create a global vector that includes a single redundant array and two da arrays */
111: DMCompositeCreate(PETSC_COMM_WORLD,&packer);
112: DMCompositeAddArray(packer,0,1);
113: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-5,2,1,PETSC_NULL,&da);
114: DMCompositeAddDM(packer,(DM)da);
117: /* create nonlinear multi-level solver */
118: DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
119: DMMGSetDM(dmmg,(DM)packer);
120: DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
121: DMMGSetFromOptions(dmmg);
123: if (use_monitor) {
124: /* create graphics windows */
125: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);
126: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);
127: for (i=0; i<DMMGGetLevels(dmmg); i++) {
128: SNESMonitorSet(dmmg[i]->snes,Monitor,dmmg[i],0);
129: }
130: }
132: DMMGSolve(dmmg);
133: DMMGDestroy(dmmg);
135: DMDestroy(&da);
136: DMDestroy(&packer);
137: if (use_monitor) {
138: PetscViewerDestroy(&user.u_lambda_viewer);
139: PetscViewerDestroy(&user.fu_lambda_viewer);
140: }
142: PetscFinalize();
143: return 0;
144: }
146: typedef struct {
147: PetscScalar u;
148: PetscScalar lambda;
149: } ULambda;
153: /*
154: Evaluates FU = Gradiant(L(w,u,lambda))
156: This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and
157: DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()).
159: */
160: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
161: {
162: DMMG dmmg = (DMMG)dummy;
164: PetscInt xs,xm,i,N,nredundant;
165: ULambda *u_lambda,*fu_lambda;
166: PetscScalar d,h,*w,*fw;
167: Vec vu_lambda,vfu_lambda;
168: DM da;
169: DM packer = (DM)dmmg->dm;
172: DMCompositeGetEntries(packer,&nredundant,&da);
173: DMCompositeGetLocalVectors(packer,&w,&vu_lambda);
174: DMCompositeScatter(packer,U,w,vu_lambda);
175: DMCompositeGetAccess(packer,FU,&fw,&vfu_lambda);
177: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
178: DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
179: DMDAVecGetArray(da,vu_lambda,&u_lambda);
180: DMDAVecGetArray(da,vfu_lambda,&fu_lambda);
181: d = N-1.0;
182: h = 1.0/d;
184: /* derivative of L() w.r.t. w */
185: if (xs == 0) { /* only first processor computes this */
186: fw[0] = -2.0*d*u_lambda[0].lambda;
187: }
189: /* derivative of L() w.r.t. u */
190: for (i=xs; i<xs+xm; i++) {
191: if (i == 0) fu_lambda[0].lambda = h*u_lambda[0].u + 2.*d*u_lambda[0].lambda - d*u_lambda[1].lambda;
192: 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;
193: 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;
194: 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;
195: 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);
196: }
198: /* derivative of L() w.r.t. lambda */
199: for (i=xs; i<xs+xm; i++) {
200: if (i == 0) fu_lambda[0].u = 2.0*d*(u_lambda[0].u - w[0]);
201: else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;
202: 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);
203: }
205: DMDAVecRestoreArray(da,vu_lambda,&u_lambda);
206: DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda);
207: DMCompositeRestoreLocalVectors(packer,&w,&vu_lambda);
208: DMCompositeRestoreAccess(packer,FU,&fw,&vfu_lambda);
209: PetscLogFlops(13.0*N);
210: return(0);
211: }
215: /*
216: Computes the exact solution
217: */
218: PetscErrorCode u_solution(void *dummy,PetscInt n,const PetscScalar *x,PetscScalar *u)
219: {
220: PetscInt i;
222: for (i=0; i<n; i++) {
223: u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
224: }
225: return(0);
226: }
230: PetscErrorCode ExactSolution(DM packer,Vec U)
231: {
232: PF pf;
233: Vec x,u_global;
234: PetscScalar *w;
235: DM da;
237: PetscInt m;
240: DMCompositeGetEntries(packer,&m,&da);
242: PFCreate(PETSC_COMM_WORLD,1,2,&pf);
243: PFSetType(pf,PFQUICK,(void*)u_solution);
244: DMDAGetCoordinates(da,&x);
245: if (!x) {
246: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
247: DMDAGetCoordinates(da,&x);
248: }
249: DMCompositeGetAccess(packer,U,&w,&u_global,0);
250: if (w) w[0] = .25;
251: PFApplyVec(pf,x,u_global);
252: PFDestroy(&pf);
253: DMCompositeRestoreAccess(packer,U,&w,&u_global,0);
254: return(0);
255: }
259: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
260: {
261: DMMG dmmg = (DMMG)dummy;
262: UserCtx *user = (UserCtx*)dmmg->user;
264: PetscInt m,N;
265: PetscScalar *w,*dw;
266: Vec u_lambda,U,F,Uexact;
267: DM packer = (DM)dmmg->dm;
268: PetscReal norm;
269: DM da;
272: SNESGetSolution(snes,&U);
273: DMCompositeGetAccess(packer,U,&w,&u_lambda);
274: VecView(u_lambda,user->u_lambda_viewer);
275: DMCompositeRestoreAccess(packer,U,&w,&u_lambda);
277: SNESGetFunction(snes,&F,0,0);
278: DMCompositeGetAccess(packer,F,&w,&u_lambda);
279: /* VecView(u_lambda,user->fu_lambda_viewer); */
280: DMCompositeRestoreAccess(packer,U,&w,&u_lambda);
282: DMCompositeGetEntries(packer,&m,&da);
283: DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
284: VecDuplicate(U,&Uexact);
285: ExactSolution(packer,Uexact);
286: VecAXPY(Uexact,-1.0,U);
287: DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda);
288: VecStrideNorm(u_lambda,0,NORM_2,&norm);
289: norm = norm/sqrt(N-1.);
290: if (dw) PetscPrintf(dmmg->comm,"Norm of error %G Error at x = 0 %G\n",norm,PetscRealPart(dw[0]));
291: VecView(u_lambda,user->fu_lambda_viewer);
292: DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda);
293: VecDestroy(&Uexact);
294: return(0);
295: }