Actual source code: ex24.c

  1: /*$Id: ex24.c,v 1.25 2001/08/07 03:04:16 balay Exp $*/

  3: static char help[] = "Solves PDE optimization problem of ex22.c with AD for adjoint.nn";

 5:  #include petscda.h
 6:  #include petscpf.h
 7:  #include petscmg.h
 8:  #include petscsnes.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: */


 51: extern int FormFunction(SNES,Vec,Vec,void*);
 52: extern int PDEFormFunctionLocal(DALocalInfo*,PetscScalar*,PetscScalar*,PassiveScalar*);

 54: typedef struct {
 55:   Mat        J;           /* Jacobian of PDE system */
 56:   SLES       sles;        /* Solver for that Jacobian */
 57: } AppCtx;

 59: #undef __FUNCT__
 61: int myPCApply(DMMG dmmg,Vec x,Vec y)
 62: {
 63:   Vec          xu,xlambda,yu,ylambda;
 64:   PetscScalar  *xw,*yw;
 65:   int          ierr;
 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:   SLESSolve(appctx->sles,xu,yu,PETSC_IGNORE);
 76:   SLESSolveTranspose(appctx->sles,xlambda,ylambda,PETSC_IGNORE);
 77:   /*  VecCopy(xu,yu);
 78:       VecCopy(xlambda,ylambda); */
 79:   VecPackRestoreAccess(packer,x,&xw,&xu,&xlambda);
 80:   VecPackRestoreAccess(packer,y,&yw,&yu,&ylambda);
 81:   return(0);
 82: }

 84: #undef __FUNCT__
 86: int myPCView(DMMG dmmg,PetscViewer v)
 87: {
 88:   int     ierr;
 89:   AppCtx  *appctx = (AppCtx*)dmmg->user;

 92:   SLESView(appctx->sles,v);
 93:   return(0);
 94: }

 96: #undef __FUNCT__
 98: int main(int argc,char **argv)
 99: {
100:   int        ierr,nlevels,i,j;
101:   DA         da;
102:   DMMG       *dmmg;
103:   VecPack    packer;
104:   AppCtx     *appctx;
105:   ISColoring iscoloring;
106:   PetscTruth bdp;

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

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

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

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

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

155:   DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);

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

163:       appctx = (AppCtx*) dmmg[i]->user;
164:       ierr   = SLESCreate(PETSC_COMM_WORLD,&appctx->sles);
165:       ierr   = SLESSetOptionsPrefix(appctx->sles,"bdp_");
166:       ierr   = SLESSetFromOptions(appctx->sles);

168:       SNESGetSLES(dmmg[i]->snes,&sles);
169:       SLESGetPC(sles,&pc);
170:       for (j=0; j<=i; j++) {
171:         MGGetSmoother(pc,j,&sles);
172:         SLESGetPC(sles,&mpc);
173:         PCSetType(mpc,PCSHELL);
174:         PCShellSetApply(mpc,(int (*)(void*,Vec,Vec))myPCApply,dmmg[j]);
175:         PCShellSetView(mpc,(int (*)(void*,PetscViewer))myPCView);
176:       }
177:     }
178:   }

180:   DMMGSolve(dmmg);

182:   /* VecView(DMMGGetx(dmmg),PETSC_VIEWER_SOCKET_WORLD); */
183:   for (i=0; i<nlevels; i++) {
184:     appctx = (AppCtx*)dmmg[i]->user;
185:     ierr   = MatDestroy(appctx->J);
186:     if (appctx->sles) {SLESDestroy(appctx->sles);}
187:     ierr   = PetscFree(appctx);
188:   }
189:   DMMGDestroy(dmmg);

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

200:      Process adiC: PDEFormFunctionLocal
201: */
202: #undef __FUNCT__
204: int PDEFormFunctionLocal(DALocalInfo *info,PetscScalar *u,PetscScalar *fu,PassiveScalar *w)
205: {
206:   int          xs = info->xs,xm = info->xm,i,mx = info->mx;
207:   PetscScalar  d,h;

209:   d    = mx-1.0;
210:   h    = 1.0/d;

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

218:   PetscLogFlops(9*mx);
219:   return 0;
220: }

222: /*
223:       Evaluates FU = Gradiant(L(w,u,lambda))

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

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

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

234: */
235: #undef __FUNCT__
237: int FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
238: {
239:   DMMG         dmmg = (DMMG)dummy;
240:   int          ierr,xs,xm,i,N,nredundant;
241:   PetscScalar  *u,*w,*fw,*fu,*lambda,*flambda,d,h,h2;
242:   Vec          vu,vlambda,vfu,vflambda,vglambda;
243:   DA           da;
244:   VecPack      packer = (VecPack)dmmg->dm;
245:   AppCtx       *appctx = (AppCtx*)dmmg->user;
246:   PetscTruth   skipadic;

249:   PetscOptionsHasName(0,"-skipadic",&skipadic);

251:   VecPackGetEntries(packer,&nredundant,&da,PETSC_IGNORE);
252:   DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
253:   DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
254:   d    = (N-1.0);
255:   h    = 1.0/d;
256:   h2   = 2.0*h;

258:   VecPackGetLocalVectors(packer,&w,&vu,&vlambda);
259:   VecPackScatter(packer,U,w,vu,vlambda);
260:   VecPackGetAccess(packer,FU,&fw,&vfu,&vflambda);
261:   VecPackGetAccess(packer,U,0,0,&vglambda);

263:   /* G() */
264:   DAFormFunction1(da,vu,vfu,w);
265:   if (!skipadic) {
266:     /* lambda^T G_u() */
267:     DAComputeJacobian1WithAdic(da,vu,appctx->J,w);
268:     if (appctx->sles) {
269:       SLESSetOperators(appctx->sles,appctx->J,appctx->J,SAME_NONZERO_PATTERN);
270:     }
271:     MatMultTranspose(appctx->J,vglambda,vflambda);
272:   }

274:   DAVecGetArray(da,vu,(void**)&u);
275:   DAVecGetArray(da,vfu,(void**)&fu);
276:   DAVecGetArray(da,vlambda,(void**)&lambda);
277:   DAVecGetArray(da,vflambda,(void**)&flambda);

279:   /* L_w */
280:   if (xs == 0) { /* only first processor computes this */
281:     fw[0] = -2.*d*lambda[0];
282:   }

284:   /* lambda^T G_u() */
285:   if (skipadic) {
286:     for (i=xs; i<xs+xm; i++) {
287:       if      (i == 0)   flambda[0]   = 2.*d*lambda[0]   - d*lambda[1] + h2*lambda[0]*u[0];
288:       else if (i == 1)   flambda[1]   = 2.*d*lambda[1]   - d*lambda[2] + h2*lambda[1]*u[1];
289:       else if (i == N-1) flambda[N-1] = 2.*d*lambda[N-1] - d*lambda[N-2] + h2*lambda[N-1]*u[N-1];
290:       else if (i == N-2) flambda[N-2] = 2.*d*lambda[N-2] - d*lambda[N-3] + h2*lambda[N-2]*u[N-2];
291:       else               flambda[i]   = - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]) + h2*lambda[i]*u[i];
292:     }
293:   }

295:   /* F_u */
296:   for (i=xs; i<xs+xm; i++) {
297:     if      (i == 0)   flambda[0]   +=    h*u[0];
298:     else if (i == 1)   flambda[1]   +=    h2*u[1];
299:     else if (i == N-1) flambda[N-1] +=    h*u[N-1];
300:     else if (i == N-2) flambda[N-2] +=    h2*u[N-2];
301:     else               flambda[i]   +=    h2*u[i];
302:   }

304:   DAVecRestoreArray(da,vu,(void**)&u);
305:   DAVecRestoreArray(da,vfu,(void**)&fu);
306:   DAVecRestoreArray(da,vlambda,(void**)&lambda);
307:   DAVecRestoreArray(da,vflambda,(void**)&flambda);

309:   VecPackRestoreLocalVectors(packer,&w,&vu,&vlambda);
310:   VecPackRestoreAccess(packer,FU,&fw,&vfu,&vflambda);
311:   VecPackRestoreAccess(packer,U,0,0,&vglambda);

313:   PetscLogFlops(9*N);
314:   return(0);
315: }