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 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 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 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 */ 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 } 193 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 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 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 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