Actual source code: ex3.c
1: static char help[] = "Reduced formulation of the mother problem of PDE-constrained optimisation.\n";
3: /*F
4: We solve the mother problem
6: min 1/2 || y - y_d ||^2_2 + \alpha/2 || u ||^2_2
8: subject to
10: - \laplace y = u on \Omega
11: y = 0 on \Gamma
13: where u is in L^2 and y is in H^1_0.
15: We formulate the reduced problem solely in terms of the control
16: by using the state equation to express y in terms of u, and then
17: apply LMVM/BLMVM to the resulting reduced problem.
19: Mesh independence is achieved by configuring the Riesz map for the control
20: space.
22: Example meshes where the Riesz map is crucial can be downloaded from the
23: http://people.maths.ox.ac.uk/~farrellp/meshes.tar.gz
25: Contributed by: Patrick Farrell <patrick.farrell@maths.ox.ac.uk>
27: Run with e.g.:
28: ./ex3 -laplace_ksp_type cg -laplace_pc_type hypre -tao_blmvm_mat_lmvm_J0_ksp_type cg -tao_blmvm_mat_lmvm_J0_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
30: and visualise in paraview with ../../../../petsc_gen_xdmf.py solution.h5.
32: Toggle the Riesz map (-use_riesz 0) to see the difference setting the Riesz maps makes.
34: TODO: broken for parallel runs
35: F*/
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];
62: PetscFunctionBeginUser;
63: filename[0] = '\0';
64: user->use_riesz = PETSC_TRUE;
66: PetscOptionsBegin(comm, "", "Poisson mother problem options", "DMPLEX");
67: PetscCall(PetscOptionsBool("-use_riesz", "Use the Riesz map to achieve mesh independence", "ex3.c", user->use_riesz, &user->use_riesz, NULL));
68: PetscCall(PetscOptionsString("-f", "filename to read", "ex3.c", filename, filename, sizeof(filename), &flg));
69: PetscOptionsEnd();
71: if (!flg) {
72: PetscCall(DMCreate(comm, dm));
73: PetscCall(DMSetType(*dm, DMPLEX));
74: } else {
75: /* TODO Eliminate this in favor of DMLoad() in new code */
76: #if defined(PETSC_HAVE_HDF5)
77: const PetscInt vertices_per_cell = 3;
78: PetscViewer viewer;
79: Vec coordinates;
80: Vec topology;
81: PetscInt dim = 2, numCells;
82: PetscInt numVertices;
83: PetscScalar *coords;
84: PetscScalar *topo_f;
85: PetscInt *cells;
86: DMLabel label;
88: /* Read in FEniCS HDF5 output */
89: PetscCall(PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer));
91: /* create Vecs to read in the data from H5 */
92: PetscCall(VecCreate(comm, &coordinates));
93: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
94: PetscCall(VecSetBlockSize(coordinates, dim));
95: PetscCall(VecCreate(comm, &topology));
96: PetscCall(PetscObjectSetName((PetscObject)topology, "topology"));
97: PetscCall(VecSetBlockSize(topology, vertices_per_cell));
99: /* navigate to the right group */
100: PetscCall(PetscViewerHDF5PushGroup(viewer, "/Mesh/mesh"));
102: /* Read the Vecs */
103: PetscCall(VecLoad(coordinates, viewer));
104: PetscCall(VecLoad(topology, viewer));
106: /* do some ugly calculations */
107: PetscCall(VecGetSize(topology, &numCells));
108: numCells = numCells / vertices_per_cell;
109: PetscCall(VecGetSize(coordinates, &numVertices));
110: numVertices = numVertices / dim;
112: PetscCall(VecGetArray(coordinates, &coords));
113: PetscCall(VecGetArray(topology, &topo_f));
114: /* and now we have to convert the double representation to integers to pass over, argh */
115: PetscCall(PetscMalloc1(numCells * vertices_per_cell, &cells));
116: for (PetscInt j = 0; j < numCells * vertices_per_cell; j++) cells[j] = (PetscInt)topo_f[j];
118: /* Now create the DM */
119: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, vertices_per_cell, PETSC_TRUE, cells, dim, coords, dm));
120: /* Check for flipped first cell */
121: {
122: PetscReal v0[3], J[9], invJ[9], detJ;
124: PetscCall(DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ));
125: if (detJ < 0) {
126: PetscCall(DMPlexOrientPoint(*dm, 0, -1));
127: PetscCall(DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ));
128: PetscCheck(detJ >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Something is wrong");
129: }
130: }
131: PetscCall(DMPlexOrient(*dm));
132: PetscCall(DMCreateLabel(*dm, "marker"));
133: PetscCall(DMGetLabel(*dm, "marker", &label));
134: PetscCall(DMPlexMarkBoundaryFaces(*dm, 1, label));
135: PetscCall(DMPlexLabelComplete(*dm, label));
137: PetscCall(PetscViewerDestroy(&viewer));
138: PetscCall(VecRestoreArray(coordinates, &coords));
139: PetscCall(VecRestoreArray(topology, &topo_f));
140: PetscCall(PetscFree(cells));
141: PetscCall(VecDestroy(&coordinates));
142: PetscCall(VecDestroy(&topology));
143: #else
144: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Reconfigure PETSc with --download-hdf5");
145: #endif
146: }
147: PetscCall(DMSetFromOptions(*dm));
148: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: void mass_kernel(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
153: {
154: g0[0] = 1.0;
155: }
157: void laplace_kernel(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
158: {
159: PetscInt d;
160: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
161: }
163: /* data we seek to match */
164: PetscErrorCode data_kernel(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *y, PetscCtx ctx)
165: {
166: *y = 1.0 / (2 * PETSC_PI * PETSC_PI) * PetscSinReal(PETSC_PI * x[0]) * PetscSinReal(PETSC_PI * x[1]);
167: /* the associated control is sin(pi*x[0])*sin(pi*x[1]) */
168: return PETSC_SUCCESS;
169: }
170: PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
171: {
172: *u = 0.0;
173: return PETSC_SUCCESS;
174: }
176: PetscErrorCode CreateCtx(DM dm, AppCtx *user)
177: {
178: DM dm_mass;
179: DM dm_laplace;
180: PetscDS prob_mass;
181: PetscDS prob_laplace;
182: PetscFE fe;
183: DMLabel label;
184: PetscSection section;
185: PetscInt n, k, p, d;
186: PetscInt dof, off;
187: IS is;
188: const PetscInt *points;
189: const PetscInt dim = 2;
190: PetscErrorCode (**wtf)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
192: PetscFunctionBeginUser;
193: /* make the data we seek to match */
194: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, PETSC_TRUE, NULL, 4, &fe));
196: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
197: PetscCall(DMCreateDS(dm));
198: PetscCall(DMCreateGlobalVector(dm, &user->data));
200: /* ugh, this is hideous */
201: /* y_d = interpolate(Expression("sin(x[0]) + .."), V) */
202: PetscCall(PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &wtf));
203: wtf[0] = data_kernel;
204: PetscCall(DMProjectFunction(dm, 0.0, wtf, NULL, INSERT_VALUES, user->data));
205: PetscCall(PetscFree(wtf));
207: /* assemble(inner(u, v)*dx), almost */
208: PetscCall(DMClone(dm, &dm_mass));
209: PetscCall(DMCopyDisc(dm, dm_mass));
210: PetscCall(DMSetNumFields(dm_mass, 1));
211: PetscCall(DMPlexCopyCoordinates(dm, dm_mass)); /* why do I have to do this separately? */
212: PetscCall(DMGetDS(dm_mass, &prob_mass));
213: PetscCall(PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL, NULL, NULL));
214: PetscCall(PetscDSSetDiscretization(prob_mass, 0, (PetscObject)fe));
215: PetscCall(DMCreateMatrix(dm_mass, &user->mass));
216: PetscCall(DMPlexSNESComputeJacobianFEM(dm_mass, user->data, user->mass, user->mass, NULL));
217: PetscCall(MatSetOption(user->mass, MAT_SYMMETRIC, PETSC_TRUE));
218: PetscCall(DMDestroy(&dm_mass));
220: /* inner(grad(u), grad(v))*dx with homogeneous Dirichlet boundary conditions */
221: PetscCall(DMClone(dm, &dm_laplace));
222: PetscCall(DMCopyDisc(dm, dm_laplace));
223: PetscCall(DMSetNumFields(dm_laplace, 1));
224: PetscCall(DMPlexCopyCoordinates(dm, dm_laplace));
225: PetscCall(DMGetDS(dm_laplace, &prob_laplace));
226: PetscCall(PetscDSSetJacobian(prob_laplace, 0, 0, NULL, NULL, NULL, laplace_kernel));
227: PetscCall(PetscDSSetDiscretization(prob_laplace, 0, (PetscObject)fe));
228: PetscCall(DMCreateMatrix(dm_laplace, &user->laplace));
229: PetscCall(DMPlexSNESComputeJacobianFEM(dm_laplace, user->data, user->laplace, user->laplace, NULL));
231: PetscCall(DMGetLabel(dm_laplace, "marker", &label));
232: PetscCall(DMGetLocalSection(dm_laplace, §ion));
233: PetscCall(DMLabelGetStratumSize(label, 1, &n));
234: PetscCall(DMLabelGetStratumIS(label, 1, &is));
235: PetscCall(ISGetIndices(is, &points));
236: user->num_bc_dofs = 0;
237: for (p = 0; p < n; ++p) {
238: PetscCall(PetscSectionGetDof(section, points[p], &dof));
239: user->num_bc_dofs += dof;
240: }
241: PetscCall(PetscMalloc1(user->num_bc_dofs, &user->bc_indices));
242: for (p = 0, k = 0; p < n; ++p) {
243: PetscCall(PetscSectionGetDof(section, points[p], &dof));
244: PetscCall(PetscSectionGetOffset(section, points[p], &off));
245: for (d = 0; d < dof; ++d) user->bc_indices[k++] = off + d;
246: }
247: PetscCall(ISRestoreIndices(is, &points));
248: PetscCall(ISDestroy(&is));
249: PetscCall(DMDestroy(&dm_laplace));
251: /* This is how I handle boundary conditions. I can't figure out how to get
252: plex to play with the way I want to impose the BCs. This loses symmetry,
253: but not in a disastrous way. If someone can improve it, please do! */
254: PetscCall(MatZeroRows(user->laplace, user->num_bc_dofs, user->bc_indices, 1.0, NULL, NULL));
255: PetscCall(PetscCalloc1(user->num_bc_dofs, &user->bc_values));
257: /* also create the KSP for solving the Laplace system */
258: PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->ksp_laplace));
259: PetscCall(KSPSetOperators(user->ksp_laplace, user->laplace, user->laplace));
260: PetscCall(KSPSetOptionsPrefix(user->ksp_laplace, "laplace_"));
261: PetscCall(KSPSetFromOptions(user->ksp_laplace));
263: /* A bit of setting up the user context */
264: user->dm = dm;
265: PetscCall(VecDuplicate(user->data, &user->state));
266: PetscCall(VecDuplicate(user->data, &user->adjoint));
267: PetscCall(VecDuplicate(user->data, &user->tmp1));
268: PetscCall(VecDuplicate(user->data, &user->tmp2));
270: PetscCall(PetscFEDestroy(&fe));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: PetscErrorCode DestroyCtx(AppCtx *user)
275: {
276: PetscFunctionBeginUser;
277: PetscCall(MatDestroy(&user->mass));
278: PetscCall(MatDestroy(&user->laplace));
279: PetscCall(KSPDestroy(&user->ksp_laplace));
280: PetscCall(VecDestroy(&user->data));
281: PetscCall(VecDestroy(&user->state));
282: PetscCall(VecDestroy(&user->adjoint));
283: PetscCall(VecDestroy(&user->tmp1));
284: PetscCall(VecDestroy(&user->tmp2));
285: PetscCall(PetscFree(user->bc_indices));
286: PetscCall(PetscFree(user->bc_values));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: PetscErrorCode ReducedFunctionGradient(Tao tao, Vec u, PetscReal *func, Vec g, void *userv)
291: {
292: AppCtx *user = (AppCtx *)userv;
293: const PetscReal alpha = 1.0e-6; /* regularisation parameter */
294: PetscReal inner;
296: PetscFunctionBeginUser;
297: PetscCall(MatMult(user->mass, u, user->tmp1));
298: PetscCall(VecDot(u, user->tmp1, &inner)); /* regularisation contribution to */
299: *func = alpha * 0.5 * inner; /* the functional */
301: PetscCall(VecSet(g, 0.0));
302: PetscCall(VecAXPY(g, alpha, user->tmp1)); /* regularisation contribution to the gradient */
304: /* Now compute the forward state. */
305: PetscCall(VecSetValues(user->tmp1, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES));
306: PetscCall(VecAssemblyBegin(user->tmp1));
307: PetscCall(VecAssemblyEnd(user->tmp1));
308: PetscCall(KSPSolve(user->ksp_laplace, user->tmp1, user->state)); /* forward solve */
310: /* Now compute the adjoint state also. */
311: PetscCall(VecCopy(user->state, user->tmp1));
312: PetscCall(VecAXPY(user->tmp1, -1.0, user->data));
313: PetscCall(MatMult(user->mass, user->tmp1, user->tmp2));
314: PetscCall(VecDot(user->tmp1, user->tmp2, &inner)); /* misfit contribution to */
315: *func += 0.5 * inner; /* the functional */
317: PetscCall(VecSetValues(user->tmp2, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES));
318: PetscCall(VecAssemblyBegin(user->tmp2));
319: PetscCall(VecAssemblyEnd(user->tmp2));
320: PetscCall(KSPSolve(user->ksp_laplace, user->tmp2, user->adjoint)); /* adjoint solve */
322: /* And bring it home with the gradient. */
323: PetscCall(MatMult(user->mass, user->adjoint, user->tmp1));
324: PetscCall(VecAXPY(g, 1.0, user->tmp1)); /* adjoint contribution to the gradient */
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: int main(int argc, char **argv)
329: {
330: DM dm;
331: Tao tao;
332: Vec u, lb, ub;
333: AppCtx user;
335: PetscFunctionBeginUser;
336: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
337: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
338: PetscCall(CreateCtx(dm, &user));
340: PetscCall(DMCreateGlobalVector(dm, &u));
341: PetscCall(VecDuplicate(u, &lb));
342: PetscCall(VecDuplicate(u, &ub));
343: PetscCall(VecSet(ub, 0.8)); /* a nontrivial upper bound */
345: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
346: PetscCall(TaoSetSolution(tao, u));
347: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, ReducedFunctionGradient, &user));
348: PetscCall(TaoSetVariableBounds(tao, lb, ub));
349: PetscCall(TaoSetType(tao, TAOBLMVM));
350: PetscCall(TaoSetFromOptions(tao));
352: if (user.use_riesz) {
353: PetscCall(TaoLMVMSetH0(tao, user.mass)); /* crucial for mesh independence */
354: PetscCall(TaoSetGradientNorm(tao, user.mass));
355: }
357: PetscCall(TaoSolve(tao));
359: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
360: PetscCall(VecViewFromOptions(u, NULL, "-sol_view"));
362: PetscCall(TaoDestroy(&tao));
363: PetscCall(DMDestroy(&dm));
364: PetscCall(VecDestroy(&u));
365: PetscCall(VecDestroy(&lb));
366: PetscCall(VecDestroy(&ub));
367: PetscCall(DestroyCtx(&user));
368: PetscCall(PetscFinalize());
369: return 0;
370: }
372: /*TEST
374: build:
375: requires: !complex !single
377: test:
378: requires: hdf5 double datafilespath !defined(PETSC_USE_64BIT_INDICES) hypre !cuda
379: args: -laplace_ksp_type cg -laplace_pc_type hypre -laplace_ksp_monitor_true_residual -laplace_ksp_max_it 5 -tao_blmvm_mat_lmvm_J0_ksp_type cg -tao_blmvm_mat_lmvm_J0_ksp_rtol 1e-5 -tao_blmvm_mat_lmvm_J0_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
380: filter: sed -e "s/-nan/nan/g"
382: test:
383: suffix: guess_pod
384: requires: double triangle
385: args: -laplace_ksp_type cg -laplace_pc_type gamg -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 -tao_blmvm_mat_lmvm_J0_ksp_type cg -tao_blmvm_mat_lmvm_J0_ksp_rtol 1e-5 -tao_blmvm_mat_lmvm_J0_pc_type gamg -tao_blmvm_mat_lmvm_J0_pc_gamg_esteig_ksp_type cg -tao_blmvm_mat_lmvm_J0_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 -laplace_pc_gamg_aggressive_coarsening 0
386: 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"
388: TEST*/