Actual source code: ex22.c

  2: static char help[] = "Solves PDE optimization problem.\n\n";

 4:  #include petscda.h
 5:  #include petscpf.h
 6:  #include petscsnes.h
 7:  #include petscdmmg.h

  9: /*

 11:        w - design variables (what we change to get an optimal solution)
 12:        u - state variables (i.e. the PDE solution)
 13:        lambda - the Lagrange multipliers

 15:             U = (w [u_0 lambda_0 u_1 lambda_1 .....])

 17:        fu, fw, flambda contain the gradient of L(w,u,lambda)

 19:             FU = (fw [fu_0 flambda_0 .....])

 21:        In this example the PDE is 
 22:                              Uxx = 2, 
 23:                             u(0) = w(0), thus this is the free parameter
 24:                             u(1) = 0
 25:        the function we wish to minimize is 
 26:                             \integral u^{2}

 28:        The exact solution for u is given by u(x) = x*x - 1.25*x + .25

 30:        Use the usual centered finite differences.

 32:        Note we treat the problem as non-linear though it happens to be linear

 34:        See ex21.c for the same code, but that does NOT interlaces the u and the lambda

 36:        The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
 37: */

 39: typedef struct {
 40:   PetscViewer  u_lambda_viewer;
 41:   PetscViewer  fu_lambda_viewer;
 42: } UserCtx;


 47: /*
 48:     Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the 
 49:   smoother on all levels. This is because (1) in the matrix free case no matrix entries are 
 50:   available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal
 51:   entry for the control variable is zero which means default SOR will not work.

 53: */
 54: char  common_options[]      = "-dmmg_grid_sequence \
 55:                                -dmmg_nlevels 5 \
 56:                                -mg_levels_pc_type none \
 57:                                -mg_coarse_pc_type none \
 58:                                -pc_mg_type full \
 59:                                -mg_coarse_ksp_type gmres \
 60:                                -mg_levels_ksp_type gmres \
 61:                                -mg_coarse_ksp_max_it 6 \
 62:                                -mg_levels_ksp_max_it 3";

 64: char  matrix_free_options[] = "-mat_mffd_compute_normu no \
 65:                                -mat_mffd_type wp \
 66:                                -dmmg_jacobian_mf_fd";

 68: /*
 69:     Currently only global coloring is supported with DMComposite
 70: */
 71: char  matrix_based_options[] = "-dmmg_iscoloring_type global";




 78: int main(int argc,char **argv)
 79: {
 81:   UserCtx        user;
 82:   DA             da;
 83:   DMMG           *dmmg;
 84:   DMComposite    packer;
 85:   PetscTruth     use_matrix_based = PETSC_FALSE,use_monitor = PETSC_FALSE;
 86:   PetscInt       i;

 88:   PetscInitialize(&argc,&argv,PETSC_NULL,help);
 89:   PetscOptionsSetFromOptions();

 91:   /* Hardwire several options; can be changed at command line */
 92:   PetscOptionsInsertString(common_options);
 93:   PetscOptionsGetTruth(PETSC_NULL,"-use_matrix_based",&use_matrix_based,PETSC_IGNORE);
 94:   if (use_matrix_based) {
 95:     PetscOptionsInsertString(matrix_based_options);
 96:   } else {
 97:     PetscOptionsInsertString(matrix_free_options);
 98:   }
 99:   PetscOptionsInsert(&argc,&argv,PETSC_NULL);
100:   PetscOptionsGetTruth(PETSC_NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);

102:   /* Create a global vector that includes a single redundant array and two da arrays */
103:   DMCompositeCreate(PETSC_COMM_WORLD,&packer);
104:   DMCompositeAddArray(packer,0,1);
105:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,2,1,PETSC_NULL,&da);
106:   DMCompositeAddDM(packer,(DM)da);


109:   /* create nonlinear multi-level solver */
110:   DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
111:   DMMGSetDM(dmmg,(DM)packer);
112:   DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
113:   DMMGSetFromOptions(dmmg);

115:   if (use_monitor) {
116:     /* create graphics windows */
117:     PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);
118:     PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);
119:     for (i=0; i<DMMGGetLevels(dmmg); i++) {
120:       SNESMonitorSet(dmmg[i]->snes,Monitor,dmmg[i],0);
121:     }
122:   }

124:   DMMGSolve(dmmg);
125:   DMMGDestroy(dmmg);

127:   DADestroy(da);
128:   DMCompositeDestroy(packer);
129:   if (use_monitor) {
130:     PetscViewerDestroy(user.u_lambda_viewer);
131:     PetscViewerDestroy(user.fu_lambda_viewer);
132:   }

134:   PetscFinalize();
135:   return 0;
136: }

138: typedef struct {
139:   PetscScalar u;
140:   PetscScalar lambda;
141: } ULambda;
142: 
143: /*
144:       Evaluates FU = Gradiant(L(w,u,lambda))

146:      This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and
147:    DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()).

149: */
150: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
151: {
152:   DMMG           dmmg = (DMMG)dummy;
154:   PetscInt       xs,xm,i,N,nredundant;
155:   ULambda        *u_lambda,*fu_lambda;
156:   PetscScalar    d,h,*w,*fw;
157:   Vec            vu_lambda,vfu_lambda;
158:   DA             da;
159:   DMComposite        packer = (DMComposite)dmmg->dm;

162:   DMCompositeGetEntries(packer,&nredundant,&da);
163:   DMCompositeGetLocalVectors(packer,&w,&vu_lambda);
164:   DMCompositeScatter(packer,U,w,vu_lambda);
165:   DMCompositeGetAccess(packer,FU,&fw,&vfu_lambda);

167:   DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
168:   DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
169:   DAVecGetArray(da,vu_lambda,&u_lambda);
170:   DAVecGetArray(da,vfu_lambda,&fu_lambda);
171:   d    = N-1.0;
172:   h    = 1.0/d;

174:   /* derivative of L() w.r.t. w */
175:   if (xs == 0) { /* only first processor computes this */
176:     fw[0] = -2.0*d*u_lambda[0].lambda;
177:   }

179:   /* derivative of L() w.r.t. u */
180:   for (i=xs; i<xs+xm; i++) {
181:     if      (i == 0)   fu_lambda[0].lambda   =    h*u_lambda[0].u   + 2.*d*u_lambda[0].lambda   - d*u_lambda[1].lambda;
182:     else if (i == 1)   fu_lambda[1].lambda   = 2.*h*u_lambda[1].u   + 2.*d*u_lambda[1].lambda   - d*u_lambda[2].lambda;
183:     else if (i == N-1) fu_lambda[N-1].lambda =    h*u_lambda[N-1].u + 2.*d*u_lambda[N-1].lambda - d*u_lambda[N-2].lambda;
184:     else if (i == N-2) fu_lambda[N-2].lambda = 2.*h*u_lambda[N-2].u + 2.*d*u_lambda[N-2].lambda - d*u_lambda[N-3].lambda;
185:     else               fu_lambda[i].lambda   = 2.*h*u_lambda[i].u   - d*(u_lambda[i+1].lambda - 2.0*u_lambda[i].lambda + u_lambda[i-1].lambda);
186:   }

188:   /* derivative of L() w.r.t. lambda */
189:   for (i=xs; i<xs+xm; i++) {
190:     if      (i == 0)   fu_lambda[0].u   = 2.0*d*(u_lambda[0].u - w[0]);
191:     else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;
192:     else               fu_lambda[i].u   = -(d*(u_lambda[i+1].u - 2.0*u_lambda[i].u + u_lambda[i-1].u) - 2.0*h);
193:   }

195:   DAVecRestoreArray(da,vu_lambda,&u_lambda);
196:   DAVecRestoreArray(da,vfu_lambda,&fu_lambda);
197:   DMCompositeRestoreLocalVectors(packer,&w,&vu_lambda);
198:   DMCompositeRestoreAccess(packer,FU,&fw,&vfu_lambda);
199:   PetscLogFlops(13*N);
200:   return(0);
201: }

203: /* 
204:     Computes the exact solution
205: */
206: PetscErrorCode u_solution(void *dummy,PetscInt n,PetscScalar *x,PetscScalar *u)
207: {
208:   PetscInt i;
210:   for (i=0; i<n; i++) {
211:     u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
212:   }
213:   return(0);
214: }

216: PetscErrorCode ExactSolution(DMComposite packer,Vec U)
217: {
218:   PF             pf;
219:   Vec            x,u_global;
220:   PetscScalar    *w;
221:   DA             da;
223:   PetscInt       m;

226:   DMCompositeGetEntries(packer,&m,&da);

228:   PFCreate(PETSC_COMM_WORLD,1,2,&pf);
229:   PFSetType(pf,PFQUICK,(void*)u_solution);
230:   DAGetCoordinates(da,&x);
231:   if (!x) {
232:     DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
233:     DAGetCoordinates(da,&x);
234:   }
235:   DMCompositeGetAccess(packer,U,&w,&u_global,0);
236:   if (w) w[0] = .25;
237:   PFApplyVec(pf,x,u_global);
238:   PFDestroy(pf);
239:   VecDestroy(x);
240:   DMCompositeRestoreAccess(packer,U,&w,&u_global,0);
241:   return(0);
242: }


245: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
246: {
247:   DMMG           dmmg = (DMMG)dummy;
248:   UserCtx        *user = (UserCtx*)dmmg->user;
250:   PetscInt       m,N;
251:   PetscScalar    *w,*dw;
252:   Vec            u_lambda,U,F,Uexact;
253:   DMComposite        packer = (DMComposite)dmmg->dm;
254:   PetscReal      norm;
255:   DA             da;

258:   SNESGetSolution(snes,&U);
259:   DMCompositeGetAccess(packer,U,&w,&u_lambda);
260:   VecView(u_lambda,user->u_lambda_viewer);
261:   DMCompositeRestoreAccess(packer,U,&w,&u_lambda);

263:   SNESGetFunction(snes,&F,0,0);
264:   DMCompositeGetAccess(packer,F,&w,&u_lambda);
265:   /* VecView(u_lambda,user->fu_lambda_viewer); */
266:   DMCompositeRestoreAccess(packer,U,&w,&u_lambda);

268:   DMCompositeGetEntries(packer,&m,&da);
269:   DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
270:   VecDuplicate(U,&Uexact);
271:   ExactSolution(packer,Uexact);
272:   VecAXPY(Uexact,-1.0,U);
273:   DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda);
274:   VecStrideNorm(u_lambda,0,NORM_2,&norm);
275:   norm = norm/sqrt(N-1.);
276:   if (dw) PetscPrintf(dmmg->comm,"Norm of error %G Error at x = 0 %G\n",norm,PetscRealPart(dw[0]));
277:   VecView(u_lambda,user->fu_lambda_viewer);
278:   DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda);
279:   VecDestroy(Uexact);
280:   return(0);
281: }