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(DMMG dmmg,Vec x,Vec y)
62: {
63:   Vec            xu,xlambda,yu,ylambda;
64:   PetscScalar    *xw,*yw;
66:   VecPack        packer = (VecPack)dmmg->dm;
67:   AppCtx         *appctx = (AppCtx*)dmmg->user;

70:   VecPackGetAccess(packer,x,&xw,&xu,&xlambda);
71:   VecPackGetAccess(packer,y,&yw,&yu,&ylambda);
72:   if (yw && xw) {
73:     yw[0] = xw[0];
74:   }
75:   KSPSolve(appctx->ksp,xu,yu);

77:   KSPSolveTranspose(appctx->ksp,xlambda,ylambda);
78:   /*  VecCopy(xu,yu);
79:       VecCopy(xlambda,ylambda); */
80:   VecPackRestoreAccess(packer,x,&xw,&xu,&xlambda);
81:   VecPackRestoreAccess(packer,y,&yw,&yu,&ylambda);
82:   return(0);
83: }

87: PetscErrorCode myPCView(DMMG dmmg,PetscViewer v)
88: {
90:   AppCtx         *appctx = (AppCtx*)dmmg->user;

93:   KSPView(appctx->ksp,v);
94:   return(0);
95: }

99: int main(int argc,char **argv)
100: {
102:   PetscInt       nlevels,i,j;
103:   DA             da;
104:   DMMG           *dmmg;
105:   VecPack        packer;
106:   AppCtx         *appctx;
107:   ISColoring     iscoloring;
108:   PetscTruth     bdp;

110:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

112:   /* Hardwire several options; can be changed at command line */
113:   PetscOptionsSetValue("-dmmg_grid_sequence",PETSC_NULL);
114:   PetscOptionsSetValue("-ksp_type","fgmres");
115:   PetscOptionsSetValue("-ksp_max_it","5");
116:   PetscOptionsSetValue("-pc_mg_type","full");
117:   PetscOptionsSetValue("-mg_coarse_ksp_type","gmres");
118:   PetscOptionsSetValue("-mg_levels_ksp_type","gmres");
119:   PetscOptionsSetValue("-mg_coarse_ksp_max_it","6");
120:   PetscOptionsSetValue("-mg_levels_ksp_max_it","3");
121:   PetscOptionsSetValue("-snes_mf_type","wp");
122:   PetscOptionsSetValue("-snes_mf_compute_normu","no");
123:   PetscOptionsSetValue("-snes_ls","basic");
124:   PetscOptionsSetValue("-dmmg_jacobian_mf_fd",0);
125:   /* PetscOptionsSetValue("-snes_ls","basicnonorms"); */
126:   PetscOptionsInsert(&argc,&argv,PETSC_NULL);

128:   /* create VecPack object to manage composite vector */
129:   VecPackCreate(PETSC_COMM_WORLD,&packer);
131:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&da);

136:   /* create nonlinear multi-level solver */
137:   DMMGCreate(PETSC_COMM_WORLD,2,PETSC_NULL,&dmmg);
138:   DMMGSetDM(dmmg,(DM)packer);
139:   VecPackDestroy(packer);

141:   /* Create Jacobian of PDE function for each level */
142:   nlevels = DMMGGetLevels(dmmg);
143:   for (i=0; i<nlevels; i++) {
144:     packer = (VecPack)dmmg[i]->dm;
145:     VecPackGetEntries(packer,PETSC_NULL,&da,PETSC_NULL);
146:     PetscNew(AppCtx,&appctx);
147:     DAGetColoring(da,IS_COLORING_GHOSTED,&iscoloring);
148:     DAGetMatrix(da,MATAIJ,&appctx->J);
149:     MatSetColoring(appctx->J,iscoloring);
150:     ISColoringDestroy(iscoloring);
151:     DASetLocalFunction(da,(DALocalFunction1)PDEFormFunctionLocal);
153:     dmmg[i]->user = (void*)appctx;
154:   }

156:   DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);

158:   PetscOptionsHasName(PETSC_NULL,"-bdp",&bdp);
159:   if (bdp) {
160:     for (i=0; i<nlevels; i++) {
161:       KSP  ksp;
162:       PC   pc,mpc;

164:       appctx = (AppCtx*) dmmg[i]->user;
165:       KSPCreate(PETSC_COMM_WORLD,&appctx->ksp);
166:       KSPSetOptionsPrefix(appctx->ksp,"bdp_");
167:       KSPSetFromOptions(appctx->ksp);

169:       SNESGetKSP(dmmg[i]->snes,&ksp);
170:       KSPGetPC(ksp,&pc);
171:       for (j=0; j<=i; j++) {
172:         PCMGGetSmoother(pc,j,&ksp);
173:         KSPGetPC(ksp,&mpc);
174:         PCSetType(mpc,PCSHELL);
175:         PCShellSetContext(mpc,dmmg[j]);
176:         PCShellSetApply(mpc,(PetscErrorCode (*)(void*,Vec,Vec))myPCApply);
177:         PCShellSetView(mpc,(PetscErrorCode (*)(void*,PetscViewer))myPCView);
178:       }
179:     }
180:   }

182:   DMMGSolve(dmmg);

184:   for (i=0; i<nlevels; i++) {
185:     appctx = (AppCtx*)dmmg[i]->user;
186:     MatDestroy(appctx->J);
187:     if (appctx->ksp) {KSPDestroy(appctx->ksp);}
188:     PetscFree(appctx);
189:   }
190:   DMMGDestroy(dmmg);

192:   PetscFinalize();
193:   return 0;
194: }
195:
196: /*
197:      Enforces the PDE on the grid
198:      This local function acts on the ghosted version of U (accessed via DAGetLocalVector())
199:      BUT the global, nonghosted version of FU

202: */
205: PetscErrorCode PDEFormFunctionLocal(DALocalInfo *info,PetscScalar *u,PetscScalar *fu,PassiveScalar *w)
206: {
207:   PetscInt       xs = info->xs,xm = info->xm,i,mx = info->mx;
208:   PetscScalar    d,h;

211:   d    = mx-1.0;
212:   h    = 1.0/d;

214:   for (i=xs; i<xs+xm; i++) {
215:     if      (i == 0)    fu[i]   = 2.0*d*(u[i] - w[0]) + h*u[i]*u[i];
216:     else if (i == mx-1) fu[i]   = 2.0*d*u[i] + h*u[i]*u[i];
217:     else                fu[i]   = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];
218:   }

220:   PetscLogFlops(9*mx);
221:   return 0;
222: }

224: /*

227:       This is the function that is usually passed to the SNESSetJacobian() or DMMGSetSNES() and
228:     defines the nonlinear set of equations that are to be solved.

230:      This local function acts on the ghosted version of U (accessed via VecPackGetLocalVectors() and
231:    VecPackScatter()) BUT the global, nonghosted version of FU (via VecPackAccess()).

233:      This function uses PDEFormFunction() to enforce the PDE constraint equations and its adjoint
234:    for the Lagrange multiplier equations

236: */
239: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
240: {
241:   DMMG           dmmg = (DMMG)dummy;
243:   PetscInt       xs,xm,i,N,nredundant;
244:   PetscScalar    *u,*w,*fw,*fu,*lambda,*flambda,d,h,h2;
245:   Vec            vu,vlambda,vfu,vflambda,vglambda;
246:   DA             da;
247:   VecPack        packer = (VecPack)dmmg->dm;
250:   AppCtx         *appctx = (AppCtx*)dmmg->user;
251: #endif

257: #endif

259:   VecPackGetEntries(packer,&nredundant,&da,PETSC_IGNORE);
260:   DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
261:   DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
262:   d    = (N-1.0);
263:   h    = 1.0/d;
264:   h2   = 2.0*h;

266:   VecPackGetLocalVectors(packer,&w,&vu,&vlambda);
267:   VecPackScatter(packer,U,w,vu,vlambda);
268:   VecPackGetAccess(packer,FU,&fw,&vfu,&vflambda);
269:   VecPackGetAccess(packer,U,0,0,&vglambda);

271:   /* G() */
272:   DAFormFunction1(da,vu,vfu,w);

276:     /* lambda^T G_u() */
278:     if (appctx->ksp) {
279:       KSPSetOperators(appctx->ksp,appctx->J,appctx->J,SAME_NONZERO_PATTERN);
280:     }
281:     MatMultTranspose(appctx->J,vglambda,vflambda);
282:   }
283: #endif

285:   DAVecGetArray(da,vu,&u);
286:   DAVecGetArray(da,vfu,&fu);
287:   DAVecGetArray(da,vlambda,&lambda);
288:   DAVecGetArray(da,vflambda,&flambda);

290:   /* L_w */
291:   if (xs == 0) { /* only first processor computes this */
292:     fw[0] = -2.*d*lambda[0];
293:   }

295:   /* lambda^T G_u() */
297:     for (i=xs; i<xs+xm; i++) {
298:       if      (i == 0)   flambda[0]   = 2.*d*lambda[0]   - d*lambda[1] + h2*lambda[0]*u[0];
299:       else if (i == 1)   flambda[1]   = 2.*d*lambda[1]   - d*lambda[2] + h2*lambda[1]*u[1];
300:       else if (i == N-1) flambda[N-1] = 2.*d*lambda[N-1] - d*lambda[N-2] + h2*lambda[N-1]*u[N-1];
301:       else if (i == N-2) flambda[N-2] = 2.*d*lambda[N-2] - d*lambda[N-3] + h2*lambda[N-2]*u[N-2];
302:       else               flambda[i]   = - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]) + h2*lambda[i]*u[i];
303:     }
304:   }

306:   /* F_u */
307:   for (i=xs; i<xs+xm; i++) {
308:     if      (i == 0)   flambda[0]   +=    h*u[0];
309:     else if (i == 1)   flambda[1]   +=    h2*u[1];
310:     else if (i == N-1) flambda[N-1] +=    h*u[N-1];
311:     else if (i == N-2) flambda[N-2] +=    h2*u[N-2];
312:     else               flambda[i]   +=    h2*u[i];
313:   }

315:   DAVecRestoreArray(da,vu,&u);
316:   DAVecRestoreArray(da,vfu,&fu);
317:   DAVecRestoreArray(da,vlambda,&lambda);
318:   DAVecRestoreArray(da,vflambda,&flambda);

320:   VecPackRestoreLocalVectors(packer,&w,&vu,&vlambda);
321:   VecPackRestoreAccess(packer,FU,&fw,&vfu,&vflambda);
322:   VecPackRestoreAccess(packer,U,0,0,&vglambda);

324:   PetscLogFlops(9*N);
325:   return(0);
326: }

```