1: static char help[] = "Reduced formulation of the mother problem of PDE-constrained optimisation.\n";

 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:   PetscBool      flg;
 60:   char           filename[2048];

 64:   filename[0] = '\0';
 65:   user->use_riesz = PETSC_TRUE;

 67:   PetscOptionsBegin(comm, "", "Poisson mother problem options", "DMPLEX");
 68:   PetscOptionsBool("-use_riesz", "Use the Riesz map to achieve mesh independence", "ex3.c", user->use_riesz, &user->use_riesz, NULL);
 69:   PetscOptionsString("-f", "filename to read", "ex3.c", filename, filename, sizeof(filename), &flg);
 70:   PetscOptionsEnd();

 72:   if (!flg) {
 73:     DMCreate(comm, dm);
 74:     DMSetType(*dm, DMPLEX);
 75:   } else {
 76:     /* TODO Eliminate this in favor of DMLoad() in new code */
 77: #if defined(PETSC_HAVE_HDF5)
 78:     const PetscInt vertices_per_cell = 3;
 79:     PetscViewer    viewer;
 80:     Vec            coordinates;
 81:     Vec            topology;
 82:     PetscInt       dim = 2, 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:     DMPlexCreateFromCellListPetsc(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:         DMPlexOrientPoint(*dm, 0, -1);
131:         DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ);
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:   }
151:   DMSetFromOptions(*dm);
152:   DMViewFromOptions(*dm, NULL, "-dm_view");
153:   return 0;
154: }

156: void mass_kernel(PetscInt dim, PetscInt Nf, PetscInt NfAux,
157:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
158:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
159:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
160: {
161:   g0[0] = 1.0;
162: }

164: void laplace_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 g3[])
168: {
169:   PetscInt d;
170:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
171: }

173: /* data we seek to match */
174: PetscErrorCode data_kernel(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *y, void *ctx)
175: {
176:   *y = 1.0/(2*PETSC_PI*PETSC_PI) * PetscSinReal(PETSC_PI*x[0]) * PetscSinReal(PETSC_PI*x[1]);
177:   /* the associated control is sin(pi*x[0])*sin(pi*x[1]) */
178:   return 0;
179: }
180: PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
181: {
182:   *u = 0.0;
183:   return 0;
184: }

186: PetscErrorCode CreateCtx(DM dm, AppCtx* user)
187: {

189:   DM             dm_mass;
190:   DM             dm_laplace;
191:   PetscDS        prob_mass;
192:   PetscDS        prob_laplace;
193:   PetscFE        fe;
194:   DMLabel        label;
195:   PetscSection   section;
196:   PetscInt       n, k, p, d;
197:   PetscInt       dof, off;
198:   IS             is;
199:   const PetscInt* points;
200:   const PetscInt dim = 2;
201:   const PetscInt id  = 1;
202:   PetscErrorCode (**wtf)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);


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

209:   DMSetField(dm, 0, NULL, (PetscObject) fe);
210:   DMCreateDS(dm);
211:   DMCreateGlobalVector(dm, &user->data);

213:   /* ugh, this is hideous */
214:   /* y_d = interpolate(Expression("sin(x[0]) + .."), V) */
215:   PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &wtf);
216:   wtf[0] = data_kernel;
217:   DMProjectFunction(dm, 0.0, wtf, NULL, INSERT_VALUES, user->data);
218:   PetscFree(wtf);

220:   /* assemble(inner(u, v)*dx), almost */
221:   DMClone(dm, &dm_mass);
222:   DMCopyDisc(dm, dm_mass);
223:   DMSetNumFields(dm_mass, 1);
224:   DMPlexCopyCoordinates(dm, dm_mass); /* why do I have to do this separately? */
225:   DMGetDS(dm_mass, &prob_mass);
226:   PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL, NULL, NULL);
227:   PetscDSSetDiscretization(prob_mass, 0, (PetscObject) fe);
228:   DMCreateMatrix(dm_mass, &user->mass);
229:   DMPlexSNESComputeJacobianFEM(dm_mass, user->data, user->mass, user->mass, NULL);
230:   MatSetOption(user->mass, MAT_SYMMETRIC, PETSC_TRUE);
231:   DMDestroy(&dm_mass);

233:   /* inner(grad(u), grad(v))*dx with homogeneous Dirichlet boundary conditions */
234:   DMClone(dm, &dm_laplace);
235:   DMCopyDisc(dm, dm_laplace);
236:   DMSetNumFields(dm_laplace, 1);
237:   DMPlexCopyCoordinates(dm, dm_laplace);
238:   DMGetDS(dm_laplace, &prob_laplace);
239:   PetscDSSetJacobian(prob_laplace, 0, 0, NULL, NULL, NULL, laplace_kernel);
240:   PetscDSSetDiscretization(prob_laplace, 0, (PetscObject) fe);
241:   DMCreateMatrix(dm_laplace, &user->laplace);
242:   DMPlexSNESComputeJacobianFEM(dm_laplace, user->data, user->laplace, user->laplace, NULL);

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

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

272:   /* also create the KSP for solving the Laplace system */
273:   KSPCreate(PETSC_COMM_WORLD, &user->ksp_laplace);
274:   KSPSetOperators(user->ksp_laplace, user->laplace, user->laplace);
275:   KSPSetOptionsPrefix(user->ksp_laplace, "laplace_");
276:   KSPSetFromOptions(user->ksp_laplace);

278:   /* A bit of setting up the user context */
279:   user->dm = dm;
280:   VecDuplicate(user->data, &user->state);
281:   VecDuplicate(user->data, &user->adjoint);
282:   VecDuplicate(user->data, &user->tmp1);
283:   VecDuplicate(user->data, &user->tmp2);

285:   PetscFEDestroy(&fe);

287:   return 0;
288: }

290: PetscErrorCode DestroyCtx(AppCtx* user)
291: {

294:   MatDestroy(&user->mass);
295:   MatDestroy(&user->laplace);
296:   KSPDestroy(&user->ksp_laplace);
297:   VecDestroy(&user->data);
298:   VecDestroy(&user->state);
299:   VecDestroy(&user->adjoint);
300:   VecDestroy(&user->tmp1);
301:   VecDestroy(&user->tmp2);
302:   PetscFree(user->bc_indices);
303:   PetscFree(user->bc_values);

305:   return 0;
306: }

308: PetscErrorCode ReducedFunctionGradient(Tao tao, Vec u, PetscReal* func, Vec g, void* userv)
309: {
310:   AppCtx* user = (AppCtx*) userv;
311:   const PetscReal alpha = 1.0e-6; /* regularisation parameter */
312:   PetscReal inner;


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

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

323:   /* Now compute the forward state. */
324:   VecSetValues(user->tmp1, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES);
325:   VecAssemblyBegin(user->tmp1);
326:   VecAssemblyEnd(user->tmp1);
327:   KSPSolve(user->ksp_laplace, user->tmp1, user->state); /* forward solve */

329:   /* Now compute the adjoint state also. */
330:   VecCopy(user->state, user->tmp1);
331:   VecAXPY(user->tmp1, -1.0, user->data);
332:   MatMult(user->mass, user->tmp1, user->tmp2);
333:   VecDot(user->tmp1, user->tmp2, &inner);      /* misfit contribution to */
334:   *func += 0.5 * inner;                                             /* the functional         */

336:   VecSetValues(user->tmp2, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES);
337:   VecAssemblyBegin(user->tmp2);
338:   VecAssemblyEnd(user->tmp2);
339:   KSPSolve(user->ksp_laplace, user->tmp2, user->adjoint); /* adjoint solve */

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

345:   return 0;
346: }

348: int main(int argc, char **argv)
349: {
350:   DM             dm;
351:   Tao            tao;
352:   Vec            u, lb, ub;
353:   AppCtx         user;

355:   PetscInitialize(&argc, &argv, NULL,help);
356:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
357:   CreateCtx(dm, &user);

359:   DMCreateGlobalVector(dm, &u);
360:   VecSet(u, 0.0);
361:   VecDuplicate(u, &lb);
362:   VecDuplicate(u, &ub);
363:   VecSet(lb, 0.0); /* satisfied at the minimum anyway */
364:   VecSet(ub, 0.8); /* a nontrivial upper bound */

366:   TaoCreate(PETSC_COMM_WORLD, &tao);
367:   TaoSetSolution(tao, u);
368:   TaoSetObjectiveAndGradient(tao,NULL, ReducedFunctionGradient, &user);
369:   TaoSetVariableBounds(tao, lb, ub);
370:   TaoSetType(tao, TAOBLMVM);
371:   TaoSetFromOptions(tao);

373:   if (user.use_riesz) {
374:     TaoLMVMSetH0(tao, user.mass);       /* crucial for mesh independence */
375:     TaoSetGradientNorm(tao, user.mass);
376:   }

378:   TaoSolve(tao);

380:   DMViewFromOptions(dm, NULL, "-dm_view");
381:   VecViewFromOptions(u, NULL, "-sol_view");

383:   TaoDestroy(&tao);
384:   DMDestroy(&dm);
385:   VecDestroy(&u);
386:   VecDestroy(&lb);
387:   VecDestroy(&ub);
388:   DestroyCtx(&user);
389:   PetscFinalize();
390:   return 0;
391: }

393: /*TEST

395:     build:
396:       requires: !complex !single

398:     test:
399:       requires: hdf5 double datafilespath !defined(PETSC_USE_64BIT_INDICES) hypre !cuda
400:       args: -laplace_ksp_type cg -laplace_pc_type hypre -laplace_ksp_monitor_true_residual -laplace_ksp_max_it 5 -mat_lmvm_ksp_type cg -mat_lmvm_ksp_rtol 1e-5 -mat_lmvm_pc_type gamg -tao_monitor -petscspace_degree 1 -tao_converged_reason -tao_gatol 1.0e-6 -dm_view hdf5:solution.h5 -sol_view hdf5:solution.h5::append -use_riesz 1 -f $DATAFILESPATH/meshes/mesh-1.h5
401:       filter: sed -e "s/-nan/nan/g"

403:     test:
404:       suffix: guess_pod
405:       requires: double triangle
406:       args: -laplace_ksp_type cg -laplace_pc_type gamg -laplace_ksp_monitor_true_residual -laplace_ksp_max_it 8 -laplace_pc_gamg_esteig_ksp_type cg -laplace_pc_gamg_esteig_ksp_max_it 5 -laplace_mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.0 -laplace_ksp_converged_reason -mat_lmvm_ksp_type cg -mat_lmvm_ksp_rtol 1e-5 -mat_lmvm_pc_type gamg -mat_lmvm_pc_gamg_esteig_ksp_type cg -mat_lmvm_pc_gamg_esteig_ksp_max_it 3 -tao_monitor -petscspace_degree 1 -tao_converged_reason -dm_refine 0 -laplace_ksp_guess_type pod -tao_gatol 1e-6
407:       filter: sed -e "s/-nan/nan/g" -e "s/-NaN/nan/g" -e "s/NaN/nan/g" -e "s/CONVERGED_RTOL iterations 9/CONVERGED_RTOL iterations 8/g" -e "s/CONVERGED_RTOL iterations 4/CONVERGED_RTOL iterations 3/g"

409: TEST*/