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: }