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";

 73: /*
 74:      The -use_matrix_based version does not work! This is because the DMComposite code cannot determine the nonzero
 75:   pattern of the Jacobian since the coupling between the boundary condition (array variable) and DA variables is problem 
 76:   dependent. To get the explicit Jacobian correct you would need to use the DMCompositeSetCoupling() to indicate the extra nonzero 
 77:   pattern and run with -dmmg_coloring_from_mat.
 78: */


 83: int main(int argc,char **argv)
 84: {
 86:   UserCtx        user;
 87:   DA             da;
 88:   DMMG           *dmmg;
 89:   DMComposite    packer;
 90:   PetscTruth     use_matrix_based = PETSC_FALSE,use_monitor = PETSC_FALSE;
 91:   PetscInt       i;

 93:   PetscInitialize(&argc,&argv,PETSC_NULL,help);
 94:   PetscOptionsSetFromOptions();

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

107:   /* Create a global vector that includes a single redundant array and two da arrays */
108:   DMCompositeCreate(PETSC_COMM_WORLD,&packer);
109:   DMCompositeAddArray(packer,0,1);
110:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,2,1,PETSC_NULL,&da);
111:   DMCompositeAddDM(packer,(DM)da);


114:   /* create nonlinear multi-level solver */
115:   DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
116:   DMMGSetDM(dmmg,(DM)packer);
117:   DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
118:   DMMGSetFromOptions(dmmg);

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

129:   DMMGSolve(dmmg);
130:   DMMGDestroy(dmmg);

132:   DADestroy(da);
133:   DMCompositeDestroy(packer);
134:   if (use_monitor) {
135:     PetscViewerDestroy(user.u_lambda_viewer);
136:     PetscViewerDestroy(user.fu_lambda_viewer);
137:   }

139:   PetscFinalize();
140:   return 0;
141: }

143: typedef struct {
144:   PetscScalar u;
145:   PetscScalar lambda;
146: } ULambda;
147: 
148: /*
149:       Evaluates FU = Gradiant(L(w,u,lambda))

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

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

167:   DMCompositeGetEntries(packer,&nredundant,&da);
168:   DMCompositeGetLocalVectors(packer,&w,&vu_lambda);
169:   DMCompositeScatter(packer,U,w,vu_lambda);
170:   DMCompositeGetAccess(packer,FU,&fw,&vfu_lambda);

172:   DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
173:   DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
174:   DAVecGetArray(da,vu_lambda,&u_lambda);
175:   DAVecGetArray(da,vfu_lambda,&fu_lambda);
176:   d    = N-1.0;
177:   h    = 1.0/d;

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

184:   /* derivative of L() w.r.t. u */
185:   for (i=xs; i<xs+xm; i++) {
186:     if      (i == 0)   fu_lambda[0].lambda   =    h*u_lambda[0].u   + 2.*d*u_lambda[0].lambda   - d*u_lambda[1].lambda;
187:     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;
188:     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;
189:     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;
190:     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);
191:   }

193:   /* derivative of L() w.r.t. lambda */
194:   for (i=xs; i<xs+xm; i++) {
195:     if      (i == 0)   fu_lambda[0].u   = 2.0*d*(u_lambda[0].u - w[0]);
196:     else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;
197:     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);
198:   }

200:   DAVecRestoreArray(da,vu_lambda,&u_lambda);
201:   DAVecRestoreArray(da,vfu_lambda,&fu_lambda);
202:   DMCompositeRestoreLocalVectors(packer,&w,&vu_lambda);
203:   DMCompositeRestoreAccess(packer,FU,&fw,&vfu_lambda);
204:   PetscLogFlops(13.0*N);
205:   return(0);
206: }

208: /* 
209:     Computes the exact solution
210: */
211: PetscErrorCode u_solution(void *dummy,PetscInt n,PetscScalar *x,PetscScalar *u)
212: {
213:   PetscInt i;
215:   for (i=0; i<n; i++) {
216:     u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
217:   }
218:   return(0);
219: }

221: PetscErrorCode ExactSolution(DMComposite packer,Vec U)
222: {
223:   PF             pf;
224:   Vec            x,u_global;
225:   PetscScalar    *w;
226:   DA             da;
228:   PetscInt       m;

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

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


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

263:   SNESGetSolution(snes,&U);
264:   DMCompositeGetAccess(packer,U,&w,&u_lambda);
265:   VecView(u_lambda,user->u_lambda_viewer);
266:   DMCompositeRestoreAccess(packer,U,&w,&u_lambda);

268:   SNESGetFunction(snes,&F,0,0);
269:   DMCompositeGetAccess(packer,F,&w,&u_lambda);
270:   /* VecView(u_lambda,user->fu_lambda_viewer); */
271:   DMCompositeRestoreAccess(packer,U,&w,&u_lambda);

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