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