1 static char help[] = "Reduced formulation of the mother problem of PDE-constrained optimisation.\n";
2
3 /*F
4 We solve the mother problem
5
6 min 1/2 || y - y_d ||^2_2 + \alpha/2 || u ||^2_2
7
8 subject to
9
10 - \laplace y = u on \Omega
11 y = 0 on \Gamma
12
13 where u is in L^2 and y is in H^1_0.
14
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.
18
19 Mesh independence is achieved by configuring the Riesz map for the control
20 space.
21
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
24
25 Contributed by: Patrick Farrell <patrick.farrell@maths.ox.ac.uk>
26
27 Run with e.g.:
28 ./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
29
30 and visualise in paraview with ../../../../petsc_gen_xdmf.py solution.h5.
31
32 Toggle the Riesz map (-use_riesz 0) to see the difference setting the Riesz maps makes.
33
34 TODO: broken for parallel runs
35 F*/
36
37 #include <petsc.h>
38 #include <petscfe.h>
39 #include <petscviewerhdf5.h>
40
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;
56
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)57 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
58 {
59 DM distributedMesh = NULL;
60 PetscErrorCode ierr;
61 const PetscInt dim = 2;
62 char filename[2048];
63 PetscBool flg;
64
65 PetscFunctionBeginUser;
66 ierr = PetscOptionsBegin(comm, "", "Poisson mother problem options", "DMPLEX");CHKERRQ(ierr);
67 filename[0] = '\0';
68 user->use_riesz = PETSC_TRUE;
69
70 ierr = PetscOptionsBool("-use_riesz", "Use the Riesz map to achieve mesh independence", "ex2.c", user->use_riesz, &user->use_riesz, NULL);CHKERRQ(ierr);
71 ierr = PetscOptionsString("-f", "filename to read", "ex2.c", filename, filename, sizeof(filename), &flg);CHKERRQ(ierr);
72 ierr = PetscOptionsEnd();CHKERRQ(ierr);
73
74 if (!flg) {
75 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
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;
89
90 /* Read in FEniCS HDF5 output */
91 ierr = PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer);CHKERRQ(ierr);
92
93 /* create Vecs to read in the data from H5 */
94 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr);
95 ierr = PetscObjectSetName((PetscObject)coordinates, "coordinates");CHKERRQ(ierr);
96 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
97 ierr = VecCreate(comm, &topology);CHKERRQ(ierr);
98 ierr = PetscObjectSetName((PetscObject)topology, "topology");CHKERRQ(ierr);
99 ierr = VecSetBlockSize(topology, vertices_per_cell);CHKERRQ(ierr);
100
101 /* navigate to the right group */
102 ierr = PetscViewerHDF5PushGroup(viewer, "/Mesh/mesh");CHKERRQ(ierr);
103
104 /* Read the Vecs */
105 ierr = VecLoad(coordinates, viewer);CHKERRQ(ierr);
106 ierr = VecLoad(topology, viewer);CHKERRQ(ierr);
107
108 /* do some ugly calculations */
109 ierr = VecGetSize(topology, &numCells);CHKERRQ(ierr);
110 numCells = numCells / vertices_per_cell;
111 ierr = VecGetSize(coordinates, &numVertices);CHKERRQ(ierr);
112 numVertices = numVertices / dim;
113
114 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
115 ierr = VecGetArray(topology, &topo_f);CHKERRQ(ierr);
116 /* and now we have to convert the double representation to integers to pass over, argh */
117 ierr = PetscMalloc1(numCells*vertices_per_cell, &cells);CHKERRQ(ierr);
118 for (j = 0; j < numCells*vertices_per_cell; j++) {
119 cells[j] = (PetscInt) topo_f[j];
120 }
121
122 /* Now create the DM */
123 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, vertices_per_cell, PETSC_TRUE, cells, dim, coords, dm);CHKERRQ(ierr);
124 /* Check for flipped first cell */
125 {
126 PetscReal v0[3], J[9], invJ[9], detJ;
127
128 ierr = DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
129 if (detJ < 0) {
130 ierr = DMPlexReverseCell(*dm, 0);CHKERRQ(ierr);
131 ierr = DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
132 if (detJ < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Something is wrong");
133 }
134 }
135 ierr = DMPlexOrient(*dm);CHKERRQ(ierr);
136 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr);
137 ierr = DMGetLabel(*dm, "marker", &label);CHKERRQ(ierr);
138 ierr = DMPlexMarkBoundaryFaces(*dm, 1, label);CHKERRQ(ierr);
139 ierr = DMPlexLabelComplete(*dm, label);CHKERRQ(ierr);
140
141 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
142 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
143 ierr = VecRestoreArray(topology, &topo_f);CHKERRQ(ierr);
144 ierr = PetscFree(cells);CHKERRQ(ierr);
145 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
146 ierr = VecDestroy(&topology);CHKERRQ(ierr);
147 #else
148 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Reconfigure PETSc with --download-hdf5");
149 #endif
150 }
151
152 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
153 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);
154 if (distributedMesh) {
155 ierr = DMDestroy(dm);CHKERRQ(ierr);
156 *dm = distributedMesh;
157 }
158
159 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
160 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
161 PetscFunctionReturn(0);
162 }
163
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[])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 }
171
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[])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 }
180
181 /* data we seek to match */
data_kernel(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * y,void * ctx)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 }
zero(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)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 }
193
CreateCtx(DM dm,AppCtx * user)194 PetscErrorCode CreateCtx(DM dm, AppCtx* user)
195 {
196 PetscErrorCode ierr;
197
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);
212
213 PetscFunctionBeginUser;
214
215 /* make the data we seek to match */
216 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, PETSC_TRUE, NULL, 4, &fe);CHKERRQ(ierr);
217
218 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr);
219 ierr = DMCreateDS(dm);CHKERRQ(ierr);
220 ierr = DMCreateGlobalVector(dm, &user->data);CHKERRQ(ierr);
221
222 /* ugh, this is hideous */
223 /* y_d = interpolate(Expression("sin(x[0]) + .."), V) */
224 ierr = PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &wtf);CHKERRQ(ierr);
225 wtf[0] = data_kernel;
226 ierr = DMProjectFunction(dm, 0.0, wtf, NULL, INSERT_VALUES, user->data);CHKERRQ(ierr);
227 ierr = PetscFree(wtf);CHKERRQ(ierr);
228
229 /* assemble(inner(u, v)*dx), almost */
230 ierr = DMClone(dm, &dm_mass);CHKERRQ(ierr);
231 ierr = DMCopyDisc(dm, dm_mass);CHKERRQ(ierr);
232 ierr = DMSetNumFields(dm_mass, 1);CHKERRQ(ierr);
233 ierr = DMPlexCopyCoordinates(dm, dm_mass);CHKERRQ(ierr); /* why do I have to do this separately? */
234 ierr = DMGetDS(dm_mass, &prob_mass);CHKERRQ(ierr);
235 ierr = PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL, NULL, NULL);CHKERRQ(ierr);
236 ierr = PetscDSSetDiscretization(prob_mass, 0, (PetscObject) fe);CHKERRQ(ierr);
237 ierr = DMCreateMatrix(dm_mass, &user->mass);CHKERRQ(ierr);
238 ierr = DMPlexSNESComputeJacobianFEM(dm_mass, user->data, user->mass, user->mass, NULL);CHKERRQ(ierr);
239 ierr = MatSetOption(user->mass, MAT_SYMMETRIC, PETSC_TRUE);CHKERRQ(ierr);
240 ierr = DMDestroy(&dm_mass);CHKERRQ(ierr);
241
242 /* inner(grad(u), grad(v))*dx with homogeneous Dirichlet boundary conditions */
243 ierr = DMClone(dm, &dm_laplace);CHKERRQ(ierr);
244 ierr = DMCopyDisc(dm, dm_laplace);CHKERRQ(ierr);
245 ierr = DMSetNumFields(dm_laplace, 1);CHKERRQ(ierr);
246 ierr = DMPlexCopyCoordinates(dm, dm_laplace);CHKERRQ(ierr);
247 ierr = DMGetDS(dm_laplace, &prob_laplace);CHKERRQ(ierr);
248 ierr = PetscDSSetJacobian(prob_laplace, 0, 0, NULL, NULL, NULL, laplace_kernel);CHKERRQ(ierr);
249 ierr = PetscDSSetDiscretization(prob_laplace, 0, (PetscObject) fe);CHKERRQ(ierr);
250 ierr = DMCreateMatrix(dm_laplace, &user->laplace);CHKERRQ(ierr);
251 ierr = DMPlexSNESComputeJacobianFEM(dm_laplace, user->data, user->laplace, user->laplace, NULL);CHKERRQ(ierr);
252
253 /* Code from Matt to get the indices associated with the boundary dofs */
254 ierr = DMAddBoundary(dm_laplace, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) zero, NULL, 1, &id, NULL);CHKERRQ(ierr);
255 ierr = DMGetLocalSection(dm_laplace, §ion);CHKERRQ(ierr);
256 ierr = DMGetLabel(dm_laplace, "marker", &label);CHKERRQ(ierr);
257 ierr = DMLabelGetStratumSize(label, 1, &n);CHKERRQ(ierr);
258 ierr = DMLabelGetStratumIS(label, 1, &is);CHKERRQ(ierr);
259 ierr = ISGetIndices(is, &points);CHKERRQ(ierr);
260 user->num_bc_dofs = 0;
261 for (p = 0; p < n; ++p) {
262 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr);
263 user->num_bc_dofs += dof;
264 }
265 ierr = PetscMalloc1(user->num_bc_dofs, &user->bc_indices);CHKERRQ(ierr);
266 for (p = 0, k = 0; p < n; ++p) {
267 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr);
268 ierr = PetscSectionGetOffset(section, points[p], &off);CHKERRQ(ierr);
269 for (d = 0; d < dof; ++d) user->bc_indices[k++] = off+d;
270 }
271 ierr = ISRestoreIndices(is, &points);CHKERRQ(ierr);
272 ierr = ISDestroy(&is);CHKERRQ(ierr);
273 ierr = DMDestroy(&dm_laplace);CHKERRQ(ierr);
274
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 ierr = MatZeroRows(user->laplace, user->num_bc_dofs, user->bc_indices, 1.0, NULL, NULL);CHKERRQ(ierr);
279 ierr = PetscCalloc1(user->num_bc_dofs, &user->bc_values);CHKERRQ(ierr);
280
281 /* also create the KSP for solving the Laplace system */
282 ierr = KSPCreate(PETSC_COMM_WORLD, &user->ksp_laplace);CHKERRQ(ierr);
283 ierr = KSPSetOperators(user->ksp_laplace, user->laplace, user->laplace);CHKERRQ(ierr);
284 ierr = KSPSetOptionsPrefix(user->ksp_laplace, "laplace_");CHKERRQ(ierr);
285 ierr = KSPSetFromOptions(user->ksp_laplace);CHKERRQ(ierr);
286
287 /* A bit of setting up the user context */
288 user->dm = dm;
289 ierr = VecDuplicate(user->data, &user->state);CHKERRQ(ierr);
290 ierr = VecDuplicate(user->data, &user->adjoint);CHKERRQ(ierr);
291 ierr = VecDuplicate(user->data, &user->tmp1);CHKERRQ(ierr);
292 ierr = VecDuplicate(user->data, &user->tmp2);CHKERRQ(ierr);
293
294 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
295
296 PetscFunctionReturn(0);
297 }
298
DestroyCtx(AppCtx * user)299 PetscErrorCode DestroyCtx(AppCtx* user)
300 {
301 PetscErrorCode ierr;
302
303 PetscFunctionBeginUser;
304
305 ierr = MatDestroy(&user->mass);CHKERRQ(ierr);
306 ierr = MatDestroy(&user->laplace);CHKERRQ(ierr);
307 ierr = KSPDestroy(&user->ksp_laplace);CHKERRQ(ierr);
308 ierr = VecDestroy(&user->data);CHKERRQ(ierr);
309 ierr = VecDestroy(&user->state);CHKERRQ(ierr);
310 ierr = VecDestroy(&user->adjoint);CHKERRQ(ierr);
311 ierr = VecDestroy(&user->tmp1);CHKERRQ(ierr);
312 ierr = VecDestroy(&user->tmp2);CHKERRQ(ierr);
313 ierr = PetscFree(user->bc_indices);CHKERRQ(ierr);
314 ierr = PetscFree(user->bc_values);CHKERRQ(ierr);
315
316 PetscFunctionReturn(0);
317 }
318
ReducedFunctionGradient(Tao tao,Vec u,PetscReal * func,Vec g,void * userv)319 PetscErrorCode ReducedFunctionGradient(Tao tao, Vec u, PetscReal* func, Vec g, void* userv)
320 {
321 PetscErrorCode ierr;
322 AppCtx* user = (AppCtx*) userv;
323 const PetscReal alpha = 1.0e-6; /* regularisation parameter */
324 PetscReal inner;
325
326 PetscFunctionBeginUser;
327
328 ierr = MatMult(user->mass, u, user->tmp1);CHKERRQ(ierr);
329 ierr = VecDot(u, user->tmp1, &inner);CHKERRQ(ierr); /* regularisation contribution to */
330 *func = alpha * 0.5 * inner; /* the functional */
331
332 ierr = VecSet(g, 0.0);CHKERRQ(ierr);
333 ierr = VecAXPY(g, alpha, user->tmp1);CHKERRQ(ierr); /* regularisation contribution to the gradient */
334
335 /* Now compute the forward state. */
336 ierr = VecSetValues(user->tmp1, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES);CHKERRQ(ierr);
337 ierr = VecAssemblyBegin(user->tmp1);CHKERRQ(ierr);
338 ierr = VecAssemblyEnd(user->tmp1);CHKERRQ(ierr);
339 ierr = KSPSolve(user->ksp_laplace, user->tmp1, user->state);CHKERRQ(ierr); /* forward solve */
340
341 /* Now compute the adjoint state also. */
342 ierr = VecCopy(user->state, user->tmp1);CHKERRQ(ierr);
343 ierr = VecAXPY(user->tmp1, -1.0, user->data);CHKERRQ(ierr);
344 ierr = MatMult(user->mass, user->tmp1, user->tmp2);CHKERRQ(ierr);
345 ierr = VecDot(user->tmp1, user->tmp2, &inner);CHKERRQ(ierr); /* misfit contribution to */
346 *func += 0.5 * inner; /* the functional */
347
348 ierr = VecSetValues(user->tmp2, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES);CHKERRQ(ierr);
349 ierr = VecAssemblyBegin(user->tmp2);CHKERRQ(ierr);
350 ierr = VecAssemblyEnd(user->tmp2);CHKERRQ(ierr);
351 ierr = KSPSolve(user->ksp_laplace, user->tmp2, user->adjoint);CHKERRQ(ierr); /* adjoint solve */
352
353 /* And bring it home with the gradient. */
354 ierr = MatMult(user->mass, user->adjoint, user->tmp1);CHKERRQ(ierr);
355 ierr = VecAXPY(g, 1.0, user->tmp1);CHKERRQ(ierr); /* adjoint contribution to the gradient */
356
357 PetscFunctionReturn(0);
358 }
359
main(int argc,char ** argv)360 int main(int argc, char **argv)
361 {
362 DM dm;
363 Tao tao;
364 Vec u, lb, ub;
365 AppCtx user;
366 PetscErrorCode ierr;
367
368 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
369 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
370 ierr = CreateCtx(dm, &user);CHKERRQ(ierr);
371
372 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
373 ierr = VecSet(u, 0.0);CHKERRQ(ierr);
374 ierr = VecDuplicate(u, &lb);CHKERRQ(ierr);
375 ierr = VecDuplicate(u, &ub);CHKERRQ(ierr);
376 ierr = VecSet(lb, 0.0);CHKERRQ(ierr); /* satisfied at the minimum anyway */
377 ierr = VecSet(ub, 0.8);CHKERRQ(ierr); /* a nontrivial upper bound */
378
379 ierr = TaoCreate(PETSC_COMM_WORLD, &tao);CHKERRQ(ierr);
380 ierr = TaoSetInitialVector(tao, u);CHKERRQ(ierr);
381 ierr = TaoSetObjectiveAndGradientRoutine(tao, ReducedFunctionGradient, &user);CHKERRQ(ierr);
382 ierr = TaoSetVariableBounds(tao, lb, ub);CHKERRQ(ierr);
383 ierr = TaoSetType(tao, TAOBLMVM);CHKERRQ(ierr);
384 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
385
386 if (user.use_riesz) {
387 ierr = TaoLMVMSetH0(tao, user.mass);CHKERRQ(ierr); /* crucial for mesh independence */
388 ierr = TaoSetGradientNorm(tao, user.mass);CHKERRQ(ierr);
389 }
390
391 ierr = TaoSolve(tao);CHKERRQ(ierr);
392
393 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
394 ierr = VecViewFromOptions(u, NULL, "-sol_view");CHKERRQ(ierr);
395
396 ierr = TaoDestroy(&tao);CHKERRQ(ierr);
397 ierr = DMDestroy(&dm);CHKERRQ(ierr);
398 ierr = VecDestroy(&u);CHKERRQ(ierr);
399 ierr = VecDestroy(&lb);CHKERRQ(ierr);
400 ierr = VecDestroy(&ub);CHKERRQ(ierr);
401 ierr = DestroyCtx(&user);CHKERRQ(ierr);
402 ierr = PetscFinalize();
403 return ierr;
404 }
405
406 /*TEST
407
408 build:
409 requires: !complex !single
410
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"
415
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"
421
422 TEST*/
423