Actual source code: ex24.c
2: static char help[] = "Solves PDE optimization problem of ex22.c with AD for adjoint.\n\n";
4: #include <petscdmda.h>
5: #include <petscdmredundant.h>
6: #include <petscdmcomposite.h>
7: #include <petscpf.h>
8: #include <petscpcmg.h>
9: #include <petscsnes.h>
10: #include <petscdmmg.h>
12: /*
14: Minimize F(w,u) such that G(w,u) = 0
16: L(w,u,lambda) = F(w,u) + lambda^T G(w,u)
18: w - design variables (what we change to get an optimal solution)
19: u - state variables (i.e. the PDE solution)
20: lambda - the Lagrange multipliers
22: U = (w u lambda)
24: fu, fw, flambda contain the gradient of L(w,u,lambda)
26: FU = (fw fu flambda)
28: In this example the PDE is
29: Uxx - u^2 = 2,
30: u(0) = w(0), thus this is the free parameter
31: u(1) = 0
32: the function we wish to minimize is
33: \integral u^{2}
35: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
37: Use the usual centered finite differences.
39: Note we treat the problem as non-linear though it happens to be linear
41: The lambda and u are NOT interlaced.
43: We optionally provide a preconditioner on each level from the operator
45: (1 0 0)
46: (0 J 0)
47: (0 0 J')
49:
50: */
53: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
54: extern PetscErrorCode PDEFormFunctionLocal(DMDALocalInfo*,PetscScalar*,PetscScalar*,PassiveScalar*);
56: typedef struct {
57: Mat J; /* Jacobian of PDE system */
58: KSP ksp; /* Solver for that Jacobian */
59: } AppCtx;
63: PetscErrorCode myPCApply(PC pc,Vec x,Vec y)
64: {
65: Vec xu,xlambda,yu,ylambda;
66: PetscScalar *xw,*yw;
68: DMMG dmmg;
69: DM packer;
70: AppCtx *appctx;
73: PCShellGetContext(pc,(void**)&dmmg);
74: packer = dmmg->dm;
75: appctx = (AppCtx*)dmmg->user;
76: DMCompositeGetAccess(packer,x,&xw,&xu,&xlambda);
77: DMCompositeGetAccess(packer,y,&yw,&yu,&ylambda);
78: if (yw && xw) {
79: yw[0] = xw[0];
80: }
81: KSPSolve(appctx->ksp,xu,yu);
83: KSPSolveTranspose(appctx->ksp,xlambda,ylambda);
84: /* VecCopy(xu,yu);
85: VecCopy(xlambda,ylambda); */
86: DMCompositeRestoreAccess(packer,x,&xw,&xu,&xlambda);
87: DMCompositeRestoreAccess(packer,y,&yw,&yu,&ylambda);
88: return(0);
89: }
93: PetscErrorCode myPCView(PC pc,PetscViewer v)
94: {
96: DMMG dmmg;
97: AppCtx *appctx;
100: PCShellGetContext(pc,(void**)&dmmg);
101: appctx = (AppCtx*)dmmg->user;
102: KSPView(appctx->ksp,v);
103: return(0);
104: }
108: int main(int argc,char **argv)
109: {
111: PetscInt nlevels,i,j;
112: DM red,da;
113: DMMG *dmmg;
114: DM packer;
115: AppCtx *appctx;
116: ISColoring iscoloring;
117: PetscBool bdp;
119: PetscInitialize(&argc,&argv,PETSC_NULL,help);
121: /* Hardwire several options; can be changed at command line */
122: PetscOptionsSetValue("-dmmg_grid_sequence",PETSC_NULL);
123: PetscOptionsSetValue("-ksp_type","fgmres");
124: PetscOptionsSetValue("-ksp_max_it","5");
125: PetscOptionsSetValue("-pc_mg_type","full");
126: PetscOptionsSetValue("-mg_coarse_ksp_type","gmres");
127: PetscOptionsSetValue("-mg_levels_ksp_type","gmres");
128: PetscOptionsSetValue("-mg_coarse_ksp_max_it","6");
129: PetscOptionsSetValue("-mg_levels_ksp_max_it","3");
130: PetscOptionsSetValue("-mat_mffd_type","wp");
131: PetscOptionsSetValue("-mat_mffd_compute_normu","no");
132: PetscOptionsSetValue("-snes_ls","basic");
133: PetscOptionsSetValue("-dmmg_jacobian_mf_fd",0);
134: /* PetscOptionsSetValue("-snes_ls","basicnonorms"); */
135: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
137: /* create DMComposite object to manage composite vector */
138: DMCompositeCreate(PETSC_COMM_WORLD,&packer);
139: DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red);
140: DMCompositeAddDM(packer,red);
141: DMDestroy(&red);
142: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-5,1,1,PETSC_NULL,&da);
143: DMCompositeAddDM(packer,(DM)da);
144: DMCompositeAddDM(packer,(DM)da);
145: DMDestroy(&da);
147: /* create nonlinear multi-level solver */
148: DMMGCreate(PETSC_COMM_WORLD,2,PETSC_NULL,&dmmg);
149: DMMGSetDM(dmmg,(DM)packer);
150: DMDestroy(&packer);
152: /* Create Jacobian of PDE function for each level */
153: nlevels = DMMGGetLevels(dmmg);
154: for (i=0; i<nlevels; i++) {
155: packer = dmmg[i]->dm;
156: DMCompositeGetEntries(packer,PETSC_NULL,&da,PETSC_NULL);
157: PetscNew(AppCtx,&appctx);
158: DMCreateColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
159: DMCreateMatrix(da,MATAIJ,&appctx->J);
160: MatSetColoring(appctx->J,iscoloring);
161: ISColoringDestroy(&iscoloring);
162: DMDASetLocalFunction(da,(DMDALocalFunction1)PDEFormFunctionLocal);
163: DMDASetLocalAdicFunction(da,ad_PDEFormFunctionLocal);
164: dmmg[i]->user = (void*)appctx;
165: }
167: DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
168: DMMGSetFromOptions(dmmg);
170: PetscOptionsHasName(PETSC_NULL,"-bdp",&bdp);
171: if (bdp) {
172: for (i=0; i<nlevels; i++) {
173: KSP ksp;
174: PC pc,mpc;
176: appctx = (AppCtx*) dmmg[i]->user;
177: KSPCreate(PETSC_COMM_WORLD,&appctx->ksp);
178: KSPSetOptionsPrefix(appctx->ksp,"bdp_");
179: KSPSetFromOptions(appctx->ksp);
181: SNESGetKSP(dmmg[i]->snes,&ksp);
182: KSPGetPC(ksp,&pc);
183: for (j=0; j<=i; j++) {
184: PCMGGetSmoother(pc,j,&ksp);
185: KSPGetPC(ksp,&mpc);
186: PCSetType(mpc,PCSHELL);
187: PCShellSetContext(mpc,dmmg[j]);
188: PCShellSetApply(mpc,myPCApply);
189: PCShellSetView(mpc,myPCView);
190: }
191: }
192: }
194: DMMGSolve(dmmg);
196: for (i=0; i<nlevels; i++) {
197: appctx = (AppCtx*)dmmg[i]->user;
198: MatDestroy(&appctx->J);
199: if (appctx->ksp) {KSPDestroy(&appctx->ksp);}
200: PetscFree(appctx);
201: }
202: DMMGDestroy(dmmg);
204: PetscFinalize();
205: return 0;
206: }
207:
208: /*
209: Enforces the PDE on the grid
210: This local function acts on the ghosted version of U (accessed via DMGetLocalVector())
211: BUT the global, nonghosted version of FU
213: Process adiC(36): PDEFormFunctionLocal
214: */
217: PetscErrorCode PDEFormFunctionLocal(DMDALocalInfo *info,PetscScalar *u,PetscScalar *fu,PassiveScalar *w)
218: {
219: PetscInt xs = info->xs,xm = info->xm,i,mx = info->mx;
220: PetscScalar d,h;
223: d = mx-1.0;
224: h = 1.0/d;
226: for (i=xs; i<xs+xm; i++) {
227: if (i == 0) fu[i] = 2.0*d*(u[i] - w[0]) + h*u[i]*u[i];
228: else if (i == mx-1) fu[i] = 2.0*d*u[i] + h*u[i]*u[i];
229: else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];
230: }
232: PetscLogFlops(9.0*mx);
233: return 0;
234: }
236: /*
237: Evaluates FU = Gradiant(L(w,u,lambda))
239: This is the function that is usually passed to the SNESSetJacobian() or DMMGSetSNES() and
240: defines the nonlinear set of equations that are to be solved.
242: This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and
243: DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeAccess()).
245: This function uses PDEFormFunction() to enforce the PDE constraint equations and its adjoint
246: for the Lagrange multiplier equations
248: */
251: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
252: {
253: DMMG dmmg = (DMMG)dummy;
255: PetscInt xs,xm,i,N;
256: PetscScalar *u,*w,*fw,*fu,*lambda,*flambda,d,h,h2;
257: Vec vw,vu,vlambda,vfw,vfu,vflambda,vglambda;
258: DM red,da;
259: DM packer = (DM)dmmg->dm;
260: PetscBool useadic = PETSC_TRUE;
261: #if defined(PETSC_HAVE_ADIC)
262: AppCtx *appctx = (AppCtx*)dmmg->user;
263: #endif
267: #if defined(PETSC_HAVE_ADIC)
268: PetscOptionsHasName(0,"-useadic",&skipadic);
269: #endif
271: DMCompositeGetEntries(packer,&red,&da,PETSC_IGNORE);
272: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
273: DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
274: d = (N-1.0);
275: h = 1.0/d;
276: h2 = 2.0*h;
278: DMCompositeGetLocalVectors(packer,&vw,&vu,&vlambda);
279: DMCompositeScatter(packer,U,vw,vu,vlambda);
280: DMCompositeGetAccess(packer,FU,&vfw,&vfu,&vflambda);
281: DMCompositeGetAccess(packer,U,PETSC_NULL,PETSC_NULL,&vglambda);
283: /* G() */
284: DMDAComputeFunction1(da,vu,vfu,vw);
286: #if defined(PETSC_HAVE_ADIC)
287: if (useadic) {
288: /* lambda^T G_u() */
289: DMDAComputeJacobian1WithAdic(da,vu,appctx->J,vw);
290: if (appctx->ksp) {
291: KSPSetOperators(appctx->ksp,appctx->J,appctx->J,SAME_NONZERO_PATTERN);
292: }
293: MatMultTranspose(appctx->J,vglambda,vflambda);
294: }
295: #endif
297: VecGetArray(vw,&w);
298: VecGetArray(vfw,&fw);
299: DMDAVecGetArray(da,vu,&u);
300: DMDAVecGetArray(da,vfu,&fu);
301: DMDAVecGetArray(da,vlambda,&lambda);
302: DMDAVecGetArray(da,vflambda,&flambda);
304: /* L_w */
305: if (xs == 0) { /* only first processor computes this */
306: fw[0] = -2.*d*lambda[0];
307: }
309: /* lambda^T G_u() */
310: if (!useadic) {
311: for (i=xs; i<xs+xm; i++) {
312: if (i == 0) flambda[0] = 2.*d*lambda[0] - d*lambda[1] + h2*lambda[0]*u[0];
313: else if (i == 1) flambda[1] = 2.*d*lambda[1] - d*lambda[2] + h2*lambda[1]*u[1];
314: else if (i == N-1) flambda[N-1] = 2.*d*lambda[N-1] - d*lambda[N-2] + h2*lambda[N-1]*u[N-1];
315: else if (i == N-2) flambda[N-2] = 2.*d*lambda[N-2] - d*lambda[N-3] + h2*lambda[N-2]*u[N-2];
316: else flambda[i] = - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]) + h2*lambda[i]*u[i];
317: }
318: }
320: /* F_u */
321: for (i=xs; i<xs+xm; i++) {
322: if (i == 0) flambda[0] += h*u[0];
323: else if (i == 1) flambda[1] += h2*u[1];
324: else if (i == N-1) flambda[N-1] += h*u[N-1];
325: else if (i == N-2) flambda[N-2] += h2*u[N-2];
326: else flambda[i] += h2*u[i];
327: }
329: VecRestoreArray(vw,&w);
330: VecRestoreArray(vfw,&fw);
331: DMDAVecRestoreArray(da,vu,&u);
332: DMDAVecRestoreArray(da,vfu,&fu);
333: DMDAVecRestoreArray(da,vlambda,&lambda);
334: DMDAVecRestoreArray(da,vflambda,&flambda);
336: DMCompositeRestoreLocalVectors(packer,&vw,&vu,&vlambda);
337: DMCompositeRestoreAccess(packer,FU,&vfw,&vfu,&vflambda);
338: DMCompositeRestoreAccess(packer,U,PETSC_NULL,PETSC_NULL,&vglambda);
340: PetscLogFlops(9.0*N);
341: return(0);
342: }