Actual source code: ex24.c
2: static char help[] = "Solves PDE optimization problem of ex22.c with AD for adjoint.\n\n";
4: #include petscda.h
5: #include petscpf.h
6: #include petscmg.h
7: #include petscsnes.h
8: #include petscdmmg.h
10: /*
12: Minimize F(w,u) such that G(w,u) = 0
14: L(w,u,lambda) = F(w,u) + lambda^T G(w,u)
16: w - design variables (what we change to get an optimal solution)
17: u - state variables (i.e. the PDE solution)
18: lambda - the Lagrange multipliers
20: U = (w u lambda)
22: fu, fw, flambda contain the gradient of L(w,u,lambda)
24: FU = (fw fu flambda)
26: In this example the PDE is
27: Uxx - u^2 = 2,
28: u(0) = w(0), thus this is the free parameter
29: u(1) = 0
30: the function we wish to minimize is
31: \integral u^{2}
33: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
35: Use the usual centered finite differences.
37: Note we treat the problem as non-linear though it happens to be linear
39: The lambda and u are NOT interlaced.
41: We optionally provide a preconditioner on each level from the operator
43: (1 0 0)
44: (0 J 0)
45: (0 0 J')
47:
48: */
54: typedef struct {
55: Mat J; /* Jacobian of PDE system */
56: KSP ksp; /* Solver for that Jacobian */
57: } AppCtx;
61: PetscErrorCode myPCApply(PC pc,Vec x,Vec y)
62: {
63: Vec xu,xlambda,yu,ylambda;
64: PetscScalar *xw,*yw;
66: DMMG dmmg;
67: DMComposite packer;
68: AppCtx *appctx;
71: PCShellGetContext(pc,(void**)&dmmg);
72: packer = (DMComposite)dmmg->dm;
73: appctx = (AppCtx*)dmmg->user;
74: DMCompositeGetAccess(packer,x,&xw,&xu,&xlambda);
75: DMCompositeGetAccess(packer,y,&yw,&yu,&ylambda);
76: if (yw && xw) {
77: yw[0] = xw[0];
78: }
79: KSPSolve(appctx->ksp,xu,yu);
81: KSPSolveTranspose(appctx->ksp,xlambda,ylambda);
82: /* VecCopy(xu,yu);
83: VecCopy(xlambda,ylambda); */
84: DMCompositeRestoreAccess(packer,x,&xw,&xu,&xlambda);
85: DMCompositeRestoreAccess(packer,y,&yw,&yu,&ylambda);
86: return(0);
87: }
91: PetscErrorCode myPCView(PC pc,PetscViewer v)
92: {
94: DMMG dmmg;
95: AppCtx *appctx;
98: PCShellGetContext(pc,(void**)&dmmg);
99: appctx = (AppCtx*)dmmg->user;
100: KSPView(appctx->ksp,v);
101: return(0);
102: }
106: int main(int argc,char **argv)
107: {
109: PetscInt nlevels,i,j;
110: DA da;
111: DMMG *dmmg;
112: DMComposite packer;
113: AppCtx *appctx;
114: ISColoring iscoloring;
115: PetscTruth bdp;
117: PetscInitialize(&argc,&argv,PETSC_NULL,help);
119: /* Hardwire several options; can be changed at command line */
120: PetscOptionsSetValue("-dmmg_grid_sequence",PETSC_NULL);
121: PetscOptionsSetValue("-ksp_type","fgmres");
122: PetscOptionsSetValue("-ksp_max_it","5");
123: PetscOptionsSetValue("-pc_mg_type","full");
124: PetscOptionsSetValue("-mg_coarse_ksp_type","gmres");
125: PetscOptionsSetValue("-mg_levels_ksp_type","gmres");
126: PetscOptionsSetValue("-mg_coarse_ksp_max_it","6");
127: PetscOptionsSetValue("-mg_levels_ksp_max_it","3");
128: PetscOptionsSetValue("-mat_mffd_type","wp");
129: PetscOptionsSetValue("-mat_mffd_compute_normu","no");
130: PetscOptionsSetValue("-snes_ls","basic");
131: PetscOptionsSetValue("-dmmg_jacobian_mf_fd",0);
132: /* PetscOptionsSetValue("-snes_ls","basicnonorms"); */
133: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
135: /* create DMComposite object to manage composite vector */
136: DMCompositeCreate(PETSC_COMM_WORLD,&packer);
137: DMCompositeAddArray(packer,0,1);
138: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&da);
139: DMCompositeAddDM(packer,(DM)da);
140: DMCompositeAddDM(packer,(DM)da);
141: DADestroy(da);
143: /* create nonlinear multi-level solver */
144: DMMGCreate(PETSC_COMM_WORLD,2,PETSC_NULL,&dmmg);
145: DMMGSetDM(dmmg,(DM)packer);
146: DMCompositeDestroy(packer);
148: /* Create Jacobian of PDE function for each level */
149: nlevels = DMMGGetLevels(dmmg);
150: for (i=0; i<nlevels; i++) {
151: packer = (DMComposite)dmmg[i]->dm;
152: DMCompositeGetEntries(packer,PETSC_NULL,&da,PETSC_NULL);
153: PetscNew(AppCtx,&appctx);
154: DAGetColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
155: DAGetMatrix(da,MATAIJ,&appctx->J);
156: MatSetColoring(appctx->J,iscoloring);
157: ISColoringDestroy(iscoloring);
158: DASetLocalFunction(da,(DALocalFunction1)PDEFormFunctionLocal);
159: DASetLocalAdicFunction(da,ad_PDEFormFunctionLocal);
160: dmmg[i]->user = (void*)appctx;
161: }
163: DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
164: DMMGSetFromOptions(dmmg);
166: PetscOptionsHasName(PETSC_NULL,"-bdp",&bdp);
167: if (bdp) {
168: for (i=0; i<nlevels; i++) {
169: KSP ksp;
170: PC pc,mpc;
172: appctx = (AppCtx*) dmmg[i]->user;
173: KSPCreate(PETSC_COMM_WORLD,&appctx->ksp);
174: KSPSetOptionsPrefix(appctx->ksp,"bdp_");
175: KSPSetFromOptions(appctx->ksp);
177: SNESGetKSP(dmmg[i]->snes,&ksp);
178: KSPGetPC(ksp,&pc);
179: for (j=0; j<=i; j++) {
180: PCMGGetSmoother(pc,j,&ksp);
181: KSPGetPC(ksp,&mpc);
182: PCSetType(mpc,PCSHELL);
183: PCShellSetContext(mpc,dmmg[j]);
184: PCShellSetApply(mpc,myPCApply);
185: PCShellSetView(mpc,myPCView);
186: }
187: }
188: }
190: DMMGSolve(dmmg);
192: for (i=0; i<nlevels; i++) {
193: appctx = (AppCtx*)dmmg[i]->user;
194: MatDestroy(appctx->J);
195: if (appctx->ksp) {KSPDestroy(appctx->ksp);}
196: PetscFree(appctx);
197: }
198: DMMGDestroy(dmmg);
200: PetscFinalize();
201: return 0;
202: }
203:
204: /*
205: Enforces the PDE on the grid
206: This local function acts on the ghosted version of U (accessed via DAGetLocalVector())
207: BUT the global, nonghosted version of FU
209: Process adiC(36): PDEFormFunctionLocal
210: */
213: PetscErrorCode PDEFormFunctionLocal(DALocalInfo *info,PetscScalar *u,PetscScalar *fu,PassiveScalar *w)
214: {
215: PetscInt xs = info->xs,xm = info->xm,i,mx = info->mx;
216: PetscScalar d,h;
219: d = mx-1.0;
220: h = 1.0/d;
222: for (i=xs; i<xs+xm; i++) {
223: if (i == 0) fu[i] = 2.0*d*(u[i] - w[0]) + h*u[i]*u[i];
224: else if (i == mx-1) fu[i] = 2.0*d*u[i] + h*u[i]*u[i];
225: else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];
226: }
228: PetscLogFlops(9.0*mx);
229: return 0;
230: }
232: /*
233: Evaluates FU = Gradiant(L(w,u,lambda))
235: This is the function that is usually passed to the SNESSetJacobian() or DMMGSetSNES() and
236: defines the nonlinear set of equations that are to be solved.
238: This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and
239: DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeAccess()).
241: This function uses PDEFormFunction() to enforce the PDE constraint equations and its adjoint
242: for the Lagrange multiplier equations
244: */
247: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
248: {
249: DMMG dmmg = (DMMG)dummy;
251: PetscInt xs,xm,i,N,nredundant;
252: PetscScalar *u,*w,*fw,*fu,*lambda,*flambda,d,h,h2;
253: Vec vu,vlambda,vfu,vflambda,vglambda;
254: DA da;
255: DMComposite packer = (DMComposite)dmmg->dm;
256: PetscTruth useadic = PETSC_TRUE;
257: #if defined(PETSC_HAVE_ADIC)
258: AppCtx *appctx = (AppCtx*)dmmg->user;
259: #endif
263: #if defined(PETSC_HAVE_ADIC)
264: PetscOptionsHasName(0,"-useadic",&skipadic);
265: #endif
267: DMCompositeGetEntries(packer,&nredundant,&da,PETSC_IGNORE);
268: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
269: DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
270: d = (N-1.0);
271: h = 1.0/d;
272: h2 = 2.0*h;
274: DMCompositeGetLocalVectors(packer,&w,&vu,&vlambda);
275: DMCompositeScatter(packer,U,w,vu,vlambda);
276: DMCompositeGetAccess(packer,FU,&fw,&vfu,&vflambda);
277: DMCompositeGetAccess(packer,U,0,0,&vglambda);
279: /* G() */
280: DAFormFunction1(da,vu,vfu,w);
282: #if defined(PETSC_HAVE_ADIC)
283: if (useadic) {
284: /* lambda^T G_u() */
285: DAComputeJacobian1WithAdic(da,vu,appctx->J,w);
286: if (appctx->ksp) {
287: KSPSetOperators(appctx->ksp,appctx->J,appctx->J,SAME_NONZERO_PATTERN);
288: }
289: MatMultTranspose(appctx->J,vglambda,vflambda);
290: }
291: #endif
293: DAVecGetArray(da,vu,&u);
294: DAVecGetArray(da,vfu,&fu);
295: DAVecGetArray(da,vlambda,&lambda);
296: DAVecGetArray(da,vflambda,&flambda);
298: /* L_w */
299: if (xs == 0) { /* only first processor computes this */
300: fw[0] = -2.*d*lambda[0];
301: }
303: /* lambda^T G_u() */
304: if (!useadic) {
305: for (i=xs; i<xs+xm; i++) {
306: if (i == 0) flambda[0] = 2.*d*lambda[0] - d*lambda[1] + h2*lambda[0]*u[0];
307: else if (i == 1) flambda[1] = 2.*d*lambda[1] - d*lambda[2] + h2*lambda[1]*u[1];
308: else if (i == N-1) flambda[N-1] = 2.*d*lambda[N-1] - d*lambda[N-2] + h2*lambda[N-1]*u[N-1];
309: else if (i == N-2) flambda[N-2] = 2.*d*lambda[N-2] - d*lambda[N-3] + h2*lambda[N-2]*u[N-2];
310: else flambda[i] = - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]) + h2*lambda[i]*u[i];
311: }
312: }
314: /* F_u */
315: for (i=xs; i<xs+xm; i++) {
316: if (i == 0) flambda[0] += h*u[0];
317: else if (i == 1) flambda[1] += h2*u[1];
318: else if (i == N-1) flambda[N-1] += h*u[N-1];
319: else if (i == N-2) flambda[N-2] += h2*u[N-2];
320: else flambda[i] += h2*u[i];
321: }
323: DAVecRestoreArray(da,vu,&u);
324: DAVecRestoreArray(da,vfu,&fu);
325: DAVecRestoreArray(da,vlambda,&lambda);
326: DAVecRestoreArray(da,vflambda,&flambda);
328: DMCompositeRestoreLocalVectors(packer,&w,&vu,&vlambda);
329: DMCompositeRestoreAccess(packer,FU,&fw,&vfu,&vflambda);
330: DMCompositeRestoreAccess(packer,U,0,0,&vglambda);
332: PetscLogFlops(9.0*N);
333: return(0);
334: }