Actual source code: ex22.c

petsc-master 2017-02-22
Report Typos and Errors

  2: static const char help[] = "Solves PDE optimization problem using full-space method, interlaces state and adjoint variables.\n\n";

  4:  #include <petscdm.h>
  5:  #include <petscdmda.h>
  6:  #include <petscdmredundant.h>
  7:  #include <petscdmcomposite.h>
  8:  #include <petscpf.h>
  9:  #include <petscsnes.h>
 10:  #include <petsc/private/dmimpl.h>

 12: /*

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

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

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

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

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

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

 33:        Use the usual centered finite differences.

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

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

 39:        The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
 40: */

 42: typedef struct {
 43:   PetscViewer u_lambda_viewer;
 44:   PetscViewer fu_lambda_viewer;
 45: } UserCtx;

 47: extern PetscErrorCode ComputeFunction(SNES,Vec,Vec,void*);
 48: extern PetscErrorCode ComputeJacobian_MF(SNES,Vec,Mat,Mat,void*);
 49: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);

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

 57: */
 58: char common_options[] = "-ksp_type fgmres\
 59:                          -snes_grid_sequence 2 \
 60:                          -pc_type mg\
 61:                          -mg_levels_pc_type none \
 62:                          -mg_coarse_pc_type none \
 63:                          -pc_mg_type full \
 64:                          -mg_coarse_ksp_type gmres \
 65:                          -mg_levels_ksp_type gmres \
 66:                          -mg_coarse_ksp_max_it 6 \
 67:                          -mg_levels_ksp_max_it 3";

 69: char matrix_free_options[] = "-mat_mffd_compute_normu no \
 70:                               -mat_mffd_type wp";

 72: extern PetscErrorCode DMCreateMatrix_MF(DM,Mat*);

 74: int main(int argc,char **argv)
 75: {
 77:   UserCtx        user;
 78:   DM             red,da;
 79:   SNES           snes;
 80:   DM             packer;
 81:   PetscBool      use_monitor = PETSC_FALSE;

 83:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
 84:   PetscOptionsSetFromOptions(NULL);

 86:   /* Hardwire several options; can be changed at command line */
 87:   PetscOptionsInsertString(NULL,common_options);
 88:   PetscOptionsInsertString(NULL,matrix_free_options);
 89:   PetscOptionsInsert(NULL,&argc,&argv,NULL);
 90:   PetscOptionsGetBool(NULL,NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);

 92:   /* Create a global vector that includes a single redundant array and two da arrays */
 93:   DMCompositeCreate(PETSC_COMM_WORLD,&packer);
 94:   DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red);
 95:   DMSetOptionsPrefix(red,"red_");
 96:   DMCompositeAddDM(packer,red);
 97:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,2,1,NULL,&da);
 98:   DMSetOptionsPrefix(red,"da_");
 99:   DMSetFromOptions(da);
100:   DMSetUp(da);
101:   DMCompositeAddDM(packer,(DM)da);
102:   DMSetApplicationContext(packer,&user);

104:   packer->ops->creatematrix = DMCreateMatrix_MF;

106:   /* create nonlinear multi-level solver */
107:   SNESCreate(PETSC_COMM_WORLD,&snes);
108:   SNESSetDM(snes,packer);
109:   SNESSetFunction(snes,NULL,ComputeFunction,NULL);
110:   SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL);

112:   SNESSetFromOptions(snes);

114:   if (use_monitor) {
115:     /* create graphics windows */
116:     PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);
117:     PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);
118:     SNESMonitorSet(snes,Monitor,0,0);
119:   }

121:   SNESSolve(snes,NULL,NULL);
122:   SNESDestroy(&snes);

124:   DMDestroy(&red);
125:   DMDestroy(&da);
126:   DMDestroy(&packer);
127:   if (use_monitor) {
128:     PetscViewerDestroy(&user.u_lambda_viewer);
129:     PetscViewerDestroy(&user.fu_lambda_viewer);
130:   }
131:   PetscFinalize();
132:   return ierr;
133: }

135: typedef struct {
136:   PetscScalar u;
137:   PetscScalar lambda;
138: } ULambda;

140: /*
141:       Evaluates FU = Gradiant(L(w,u,lambda))

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

146: */
147: PetscErrorCode ComputeFunction(SNES snes,Vec U,Vec FU,void *ctx)
148: {
150:   PetscInt       xs,xm,i,N;
151:   ULambda        *u_lambda,*fu_lambda;
152:   PetscScalar    d,h,*w,*fw;
153:   Vec            vw,vfw,vu_lambda,vfu_lambda;
154:   DM             packer,red,da;

157:   VecGetDM(U, &packer);
158:   DMCompositeGetEntries(packer,&red,&da);
159:   DMCompositeGetLocalVectors(packer,&vw,&vu_lambda);
160:   DMCompositeScatter(packer,U,vw,vu_lambda);
161:   DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda);

163:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
164:   DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
165:   VecGetArray(vw,&w);
166:   VecGetArray(vfw,&fw);
167:   DMDAVecGetArray(da,vu_lambda,&u_lambda);
168:   DMDAVecGetArray(da,vfu_lambda,&fu_lambda);
169:   d    = N-1.0;
170:   h    = 1.0/d;

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

177:   /* derivative of L() w.r.t. u */
178:   for (i=xs; i<xs+xm; i++) {
179:     if      (i == 0)   fu_lambda[0].lambda   =    h*u_lambda[0].u   + 2.*d*u_lambda[0].lambda   - d*u_lambda[1].lambda;
180:     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;
181:     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;
182:     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;
183:     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);
184:   }

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

193:   VecRestoreArray(vw,&w);
194:   VecRestoreArray(vfw,&fw);
195:   DMDAVecRestoreArray(da,vu_lambda,&u_lambda);
196:   DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda);
197:   DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda);
198:   DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda);
199:   PetscLogFlops(13.0*N);
200:   return(0);
201: }

203: /*
204:     Computes the exact solution
205: */
206: PetscErrorCode u_solution(void *dummy,PetscInt n,const PetscScalar *x,PetscScalar *u)
207: {
208:   PetscInt i;

211:   for (i=0; i<n; i++) u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
212:   return(0);
213: }

215: PetscErrorCode ExactSolution(DM packer,Vec U)
216: {
217:   PF             pf;
218:   Vec            x,u_global;
219:   PetscScalar    *w;
220:   DM             da;
222:   PetscInt       m;

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

227:   PFCreate(PETSC_COMM_WORLD,1,2,&pf);
228:   /* The cast through PETSC_UINTPTR_T is so that compilers will warn about casting to void * from void(*)(void) */
229:   PFSetType(pf,PFQUICK,(void*)(PETSC_UINTPTR_T)u_solution);
230:   DMGetCoordinates(da,&x);
231:   if (!x) {
232:     DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
233:     DMGetCoordinates(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:   DMCompositeRestoreAccess(packer,U,&w,&u_global,0);
240:   return(0);
241: }

243: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
244: {
245:   UserCtx        *user;
247:   PetscInt       m,N;
248:   PetscScalar    *w,*dw;
249:   Vec            u_lambda,U,F,Uexact;
250:   DM             packer;
251:   PetscReal      norm;
252:   DM             da;

255:   SNESGetDM(snes,&packer);
256:   DMGetApplicationContext(packer,&user);
257:   SNESGetSolution(snes,&U);
258:   DMCompositeGetAccess(packer,U,&w,&u_lambda);
259:   VecView(u_lambda,user->u_lambda_viewer);
260:   DMCompositeRestoreAccess(packer,U,&w,&u_lambda);

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

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

282: PetscErrorCode DMCreateMatrix_MF(DM packer,Mat *A)
283: {
285:   Vec            t;
286:   PetscInt       m;

289:   DMGetGlobalVector(packer,&t);
290:   VecGetLocalSize(t,&m);
291:   DMRestoreGlobalVector(packer,&t);
292:   MatCreateMFFD(PETSC_COMM_WORLD,m,m,PETSC_DETERMINE,PETSC_DETERMINE,A);
293:   MatSetUp(*A);
294:   return(0);
295: }

297: PetscErrorCode ComputeJacobian_MF(SNES snes,Vec x,Mat A,Mat B,void *ctx)
298: {

302:   MatMFFDSetFunction(A,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);
303:   MatMFFDSetBase(A,x,NULL);
304:   return(0);
305: }