Actual source code: ex3.c

petsc-master 2020-06-03
Report Typos and Errors
  1: static char help[] = "Reduced formulation of the mother problem of PDE-constrained optimisation.\n";

We solve the mother problem

min 1/2 || y - y_d ||^2_2 + \alpha/2 || u ||^2_2

subject to

- \laplace y = u on \Omega
y = 0 on \Gamma

where u is in L^2 and y is in H^1_0.

We formulate the reduced problem solely in terms of the control
by using the state equation to express y in terms of u, and then
apply LMVM/BLMVM to the resulting reduced problem.

Mesh independence is achieved by configuring the Riesz map for the control
space.

Example meshes where the Riesz map is crucial can be downloaded from the
http://people.maths.ox.ac.uk/~farrellp/meshes.tar.gz

Contributed by: Patrick Farrell <patrick.farrell@maths.ox.ac.uk>

Run with e.g.:
./ex3 -laplace_ksp_type cg -laplace_pc_type hypre -mat_lmvm_ksp_type cg -mat_lmvm_pc_type gamg -laplace_ksp_monitor_true_residual -tao_monitor -petscspace_degree 1 -tao_converged_reason -tao_gatol 1.0e-9 -dm_view hdf5:solution.h5 -sol_view hdf5:solution.h5::append -use_riesz 1 -f $DATAFILESPATH/meshes/mesh-1.h5

and visualise in paraview with ../../../../petsc_gen_xdmf.py solution.h5.

Toggle the Riesz map (-use_riesz 0) to see the difference setting the Riesz maps makes.

TODO: broken for parallel runs
 37:  #include <petsc.h>
 38:  #include <petscfe.h>
 39:  #include <petscviewerhdf5.h>

 41: typedef struct {
 42:   DM  dm;
 43:   Mat mass;
 44:   Vec data;
 45:   Vec state;
 46:   Vec tmp1;
 47:   Vec tmp2;
 48:   Vec adjoint;
 49:   Mat laplace;
 50:   KSP ksp_laplace;
 51:   PetscInt  num_bc_dofs;
 52:   PetscInt* bc_indices;
 53:   PetscScalar* bc_values;
 54:   PetscBool use_riesz;
 55: } AppCtx;

 57: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 58: {
 59:   DM             distributedMesh = NULL;
 61:   const PetscInt dim = 2;
 62:   char filename[2048];
 63:   PetscBool flg;

 66:   PetscOptionsBegin(comm, "", "Poisson mother problem options", "DMPLEX");
 67:   filename[0] = '\0';
 68:   user->use_riesz = PETSC_TRUE;

 70:   PetscOptionsBool("-use_riesz", "Use the Riesz map to achieve mesh independence", "ex2.c", user->use_riesz, &user->use_riesz, NULL);
 71:   PetscOptionsString("-f", "filename to read", "ex2.c", filename, filename, sizeof(filename), &flg);
 72:   PetscOptionsEnd();

 74:   if (!flg) {
 75:     DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
 76:   } else {
 77: #if defined(PETSC_HAVE_HDF5)
 78:     const PetscInt vertices_per_cell = 3;
 79:     PetscViewer    viewer;
 80:     Vec            coordinates;
 81:     Vec            topology;
 82:     PetscInt       numCells;
 83:     PetscInt       numVertices;
 84:     PetscScalar*   coords;
 85:     PetscScalar*   topo_f;
 86:     PetscInt*      cells;
 87:     PetscInt       j;
 88:     DMLabel        label;

 90:     /* Read in FEniCS HDF5 output */
 91:     PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer);

 93:     /* create Vecs to read in the data from H5 */
 94:     VecCreate(comm, &coordinates);
 95:     PetscObjectSetName((PetscObject)coordinates, "coordinates");
 96:     VecSetBlockSize(coordinates, dim);
 97:     VecCreate(comm, &topology);
 98:     PetscObjectSetName((PetscObject)topology, "topology");
 99:     VecSetBlockSize(topology, vertices_per_cell);

101:     /* navigate to the right group */
102:     PetscViewerHDF5PushGroup(viewer, "/Mesh/mesh");

104:     /* Read the Vecs */
105:     VecLoad(coordinates, viewer);
106:     VecLoad(topology, viewer);

108:     /* do some ugly calculations */
109:     VecGetSize(topology, &numCells);
110:     numCells = numCells / vertices_per_cell;
111:     VecGetSize(coordinates, &numVertices);
112:     numVertices = numVertices / dim;

114:     VecGetArray(coordinates, &coords);
115:     VecGetArray(topology, &topo_f);
116:     /* and now we have to convert the double representation to integers to pass over, argh */
117:     PetscMalloc1(numCells*vertices_per_cell, &cells);
118:     for (j = 0; j < numCells*vertices_per_cell; j++) {
119:       cells[j] = (PetscInt) topo_f[j];
120:     }

122:     /* Now create the DM */
123:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, vertices_per_cell, PETSC_TRUE, cells, dim, coords, dm);
124:     /* Check for flipped first cell */
125:     {
126:       PetscReal v0[3], J[9], invJ[9], detJ;

128:       DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ);
129:       if (detJ < 0) {
130:         DMPlexReverseCell(*dm, 0);
131:         DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ);
132:         if (detJ < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Something is wrong");
133:       }
134:     }
135:     DMPlexOrient(*dm);
136:     DMCreateLabel(*dm, "marker");
137:     DMGetLabel(*dm, "marker", &label);
138:     DMPlexMarkBoundaryFaces(*dm, 1, label);
139:     DMPlexLabelComplete(*dm, label);

141:     PetscViewerDestroy(&viewer);
142:     VecRestoreArray(coordinates, &coords);
143:     VecRestoreArray(topology, &topo_f);
144:     PetscFree(cells);
145:     VecDestroy(&coordinates);
146:     VecDestroy(&topology);
147: #else
148:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Reconfigure PETSc with --download-hdf5");
149: #endif
150:   }

152:   PetscObjectSetName((PetscObject) *dm, "Mesh");
153:   DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
154:   if (distributedMesh) {
155:     DMDestroy(dm);
156:     *dm  = distributedMesh;
157:   }

159:   DMSetFromOptions(*dm);
160:   DMViewFromOptions(*dm, NULL, "-dm_view");
161:   return(0);
162: }

164: void mass_kernel(PetscInt dim, PetscInt Nf, PetscInt NfAux,
165:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
166:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
167:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
168: {
169:   g0[0] = 1.0;
170: }

172: void laplace_kernel(PetscInt dim, PetscInt Nf, PetscInt NfAux,
173:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
174:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
175:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
176: {
177:   PetscInt d;
178:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
179: }

181: /* data we seek to match */
182: PetscErrorCode data_kernel(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *y, void *ctx)
183: {
184:   *y = 1.0/(2*PETSC_PI*PETSC_PI) * PetscSinReal(PETSC_PI*x[0]) * PetscSinReal(PETSC_PI*x[1]);
185:   /* the associated control is sin(pi*x[0])*sin(pi*x[1]) */
186:   return 0;
187: }
188: PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
189: {
190:   *u = 0.0;
191:   return 0;
192: }

194: PetscErrorCode CreateCtx(DM dm, AppCtx* user)
195: {

198:   DM             dm_mass;
199:   DM             dm_laplace;
200:   PetscDS        prob_mass;
201:   PetscDS        prob_laplace;
202:   PetscFE        fe;
203:   DMLabel        label;
204:   PetscSection   section;
205:   PetscInt       n, k, p, d;
206:   PetscInt       dof, off;
207:   IS             is;
208:   const PetscInt* points;
209:   const PetscInt dim = 2;
210:   const PetscInt id  = 1;
211:   PetscErrorCode (**wtf)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);


215:   /* make the data we seek to match */
216:   PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, PETSC_TRUE, NULL, 4, &fe);

218:   DMSetField(dm, 0, NULL, (PetscObject) fe);
219:   DMCreateDS(dm);
220:   DMCreateGlobalVector(dm, &user->data);

222:   /* ugh, this is hideous */
223:   /* y_d = interpolate(Expression("sin(x[0]) + .."), V) */
224:   PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &wtf);
225:   wtf[0] = data_kernel;
226:   DMProjectFunction(dm, 0.0, wtf, NULL, INSERT_VALUES, user->data);
227:   PetscFree(wtf);

229:   /* assemble(inner(u, v)*dx), almost */
230:   DMClone(dm, &dm_mass);
231:   DMCopyDisc(dm, dm_mass);
232:   DMSetNumFields(dm_mass, 1);
233:   DMPlexCopyCoordinates(dm, dm_mass); /* why do I have to do this separately? */
234:   DMGetDS(dm_mass, &prob_mass);
235:   PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL, NULL, NULL);
236:   PetscDSSetDiscretization(prob_mass, 0, (PetscObject) fe);
237:   DMCreateMatrix(dm_mass, &user->mass);
238:   DMPlexSNESComputeJacobianFEM(dm_mass, user->data, user->mass, user->mass, NULL);
239:   MatSetOption(user->mass, MAT_SYMMETRIC, PETSC_TRUE);
240:   DMDestroy(&dm_mass);

242:   /* inner(grad(u), grad(v))*dx with homogeneous Dirichlet boundary conditions */
243:   DMClone(dm, &dm_laplace);
244:   DMCopyDisc(dm, dm_laplace);
245:   DMSetNumFields(dm_laplace, 1);
246:   DMPlexCopyCoordinates(dm, dm_laplace);
247:   DMGetDS(dm_laplace, &prob_laplace);
248:   PetscDSSetJacobian(prob_laplace, 0, 0, NULL, NULL, NULL, laplace_kernel);
249:   PetscDSSetDiscretization(prob_laplace, 0, (PetscObject) fe);
250:   DMCreateMatrix(dm_laplace, &user->laplace);
251:   DMPlexSNESComputeJacobianFEM(dm_laplace, user->data, user->laplace, user->laplace, NULL);

253:   /* Code from Matt to get the indices associated with the boundary dofs */
254:   DMAddBoundary(dm_laplace, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) zero, 1, &id, NULL);
255:   DMGetLocalSection(dm_laplace, &section);
256:   DMGetLabel(dm_laplace, "marker", &label);
257:   DMLabelGetStratumSize(label, 1, &n);
258:   DMLabelGetStratumIS(label, 1, &is);
259:   ISGetIndices(is, &points);
260:   user->num_bc_dofs = 0;
261:   for (p = 0; p < n; ++p) {
262:     PetscSectionGetDof(section, points[p], &dof);
263:     user->num_bc_dofs += dof;
264:   }
265:   PetscMalloc1(user->num_bc_dofs, &user->bc_indices);
266:   for (p = 0, k = 0; p < n; ++p) {
267:     PetscSectionGetDof(section, points[p], &dof);
268:     PetscSectionGetOffset(section, points[p], &off);
269:     for (d = 0; d < dof; ++d) user->bc_indices[k++] = off+d;
270:   }
271:   ISRestoreIndices(is, &points);
272:   ISDestroy(&is);
273:   DMDestroy(&dm_laplace);

275:   /* This is how I handle boundary conditions. I can't figure out how to get
276:      plex to play with the way I want to impose the BCs. This loses symmetry,
277:      but not in a disastrous way. If someone can improve it, please do! */
278:   MatZeroRows(user->laplace, user->num_bc_dofs, user->bc_indices, 1.0, NULL, NULL);
279:   PetscCalloc1(user->num_bc_dofs, &user->bc_values);

281:   /* also create the KSP for solving the Laplace system */
282:   KSPCreate(PETSC_COMM_WORLD, &user->ksp_laplace);
283:   KSPSetOperators(user->ksp_laplace, user->laplace, user->laplace);
284:   KSPSetOptionsPrefix(user->ksp_laplace, "laplace_");
285:   KSPSetFromOptions(user->ksp_laplace);

287:   /* A bit of setting up the user context */
288:   user->dm = dm;
289:   VecDuplicate(user->data, &user->state);
290:   VecDuplicate(user->data, &user->adjoint);
291:   VecDuplicate(user->data, &user->tmp1);
292:   VecDuplicate(user->data, &user->tmp2);

294:   PetscFEDestroy(&fe);

296:   return(0);
297: }

299: PetscErrorCode DestroyCtx(AppCtx* user)
300: {


305:   MatDestroy(&user->mass);
306:   MatDestroy(&user->laplace);
307:   KSPDestroy(&user->ksp_laplace);
308:   VecDestroy(&user->data);
309:   VecDestroy(&user->state);
310:   VecDestroy(&user->adjoint);
311:   VecDestroy(&user->tmp1);
312:   VecDestroy(&user->tmp2);
313:   PetscFree(user->bc_indices);
314:   PetscFree(user->bc_values);

316:   return(0);
317: }

319: PetscErrorCode ReducedFunctionGradient(Tao tao, Vec u, PetscReal* func, Vec g, void* userv)
320: {
322:   AppCtx* user = (AppCtx*) userv;
323:   const PetscReal alpha = 1.0e-6; /* regularisation parameter */
324:   PetscReal inner;


328:   MatMult(user->mass, u, user->tmp1);
329:   VecDot(u, user->tmp1, &inner);               /* regularisation contribution to */
330:   *func = alpha * 0.5 * inner;                                      /* the functional                 */

332:   VecSet(g, 0.0);
333:   VecAXPY(g, alpha, user->tmp1);               /* regularisation contribution to the gradient */

335:   /* Now compute the forward state. */
336:   VecSetValues(user->tmp1, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES);
337:   VecAssemblyBegin(user->tmp1);
338:   VecAssemblyEnd(user->tmp1);
339:   KSPSolve(user->ksp_laplace, user->tmp1, user->state); /* forward solve */

341:   /* Now compute the adjoint state also. */
342:   VecCopy(user->state, user->tmp1);
343:   VecAXPY(user->tmp1, -1.0, user->data);
344:   MatMult(user->mass, user->tmp1, user->tmp2);
345:   VecDot(user->tmp1, user->tmp2, &inner);      /* misfit contribution to */
346:   *func += 0.5 * inner;                                             /* the functional         */

348:   VecSetValues(user->tmp2, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES);
349:   VecAssemblyBegin(user->tmp2);
350:   VecAssemblyEnd(user->tmp2);
351:   KSPSolve(user->ksp_laplace, user->tmp2, user->adjoint); /* adjoint solve */

353:   /* And bring it home with the gradient. */
354:   MatMult(user->mass, user->adjoint, user->tmp1);
355:   VecAXPY(g, 1.0, user->tmp1);                 /* adjoint contribution to the gradient */

357:   return(0);
358: }

360: int main(int argc, char **argv)
361: {
362:   DM             dm;
363:   Tao            tao;
364:   Vec            u, lb, ub;
365:   AppCtx         user;

368:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
369:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
370:   CreateCtx(dm, &user);

372:   DMCreateGlobalVector(dm, &u);
373:   VecSet(u, 0.0);
374:   VecDuplicate(u, &lb);
375:   VecDuplicate(u, &ub);
376:   VecSet(lb, 0.0); /* satisfied at the minimum anyway */
377:   VecSet(ub, 0.8); /* a nontrivial upper bound */

379:   TaoCreate(PETSC_COMM_WORLD, &tao);
380:   TaoSetInitialVector(tao, u);
381:   TaoSetObjectiveAndGradientRoutine(tao, ReducedFunctionGradient, &user);
382:   TaoSetVariableBounds(tao, lb, ub);
383:   TaoSetType(tao, TAOBLMVM);
384:   TaoSetFromOptions(tao);

386:   if (user.use_riesz) {
387:     TaoLMVMSetH0(tao, user.mass);       /* crucial for mesh independence */
388:     TaoSetGradientNorm(tao, user.mass);
389:   }

391:   TaoSolve(tao);

393:   DMViewFromOptions(dm, NULL, "-dm_view");
394:   VecViewFromOptions(u, NULL, "-sol_view");

396:   TaoDestroy(&tao);
397:   DMDestroy(&dm);
398:   VecDestroy(&u);
399:   VecDestroy(&lb);
400:   VecDestroy(&ub);
401:   DestroyCtx(&user);
402:   PetscFinalize();
403:   return ierr;
404: }

406: /*TEST

408:     build:
409:       requires: !complex !single

411:     test:
412:       requires: hdf5 double datafilespath !define(PETSC_USE_64BIT_INDICES) hypre
413:       args: -laplace_ksp_type cg -laplace_pc_type hypre -laplace_ksp_monitor_true_residual -laplace_ksp_max_it 4 -mat_lmvm_ksp_type cg -mat_lmvm_pc_type gamg -tao_monitor -petscspace_degree 1 -tao_converged_reason -tao_gatol 1.0e-9 -dm_view hdf5:solution.h5 -sol_view hdf5:solution.h5::append -use_riesz 1 -f $DATAFILESPATH/meshes/mesh-1.h5
414:       filter: sed -e "s/-nan/nan/g"

416:     test:
417:       suffix: guess_pod
418:       requires: double triangle
419:       args: -laplace_ksp_type cg -laplace_pc_type gamg  -laplace_pc_gamg_esteig_ksp_type cg -laplace_ksp_converged_reason -mat_lmvm_ksp_type cg -mat_lmvm_pc_type gamg -tao_monitor -petscspace_degree 1 -tao_converged_reason -dm_refine 3 -laplace_ksp_guess_type pod -tao_gatol 1e-6 -mat_lmvm_pc_gamg_esteig_ksp_type cg
420:       filter: sed -e "s/-nan/nan/g"

422: TEST*/