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:     PetscInt       j;
     87:     DMLabel        label;

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

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

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

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

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

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

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

    125:       PetscCall(DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ));
    126:       if (detJ < 0) {
    127:         PetscCall(DMPlexOrientPoint(*dm, 0, -1));
    128:         PetscCall(DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ));
    129:         PetscCheck(detJ >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Something is wrong");
    130:       }
    131:     }
    132:     PetscCall(DMPlexOrient(*dm));
    133:     PetscCall(DMCreateLabel(*dm, "marker"));
    134:     PetscCall(DMGetLabel(*dm, "marker", &label));
    135:     PetscCall(DMPlexMarkBoundaryFaces(*dm, 1, label));
    136:     PetscCall(DMPlexLabelComplete(*dm, label));

    138:     PetscCall(PetscViewerDestroy(&viewer));
    139:     PetscCall(VecRestoreArray(coordinates, &coords));
    140:     PetscCall(VecRestoreArray(topology, &topo_f));
    141:     PetscCall(PetscFree(cells));
    142:     PetscCall(VecDestroy(&coordinates));
    143:     PetscCall(VecDestroy(&topology));
    144: #else
    145:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Reconfigure PETSc with --download-hdf5");
    146: #endif
    147:   }
    148:   PetscCall(DMSetFromOptions(*dm));
    149:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    150:   PetscFunctionReturn(PETSC_SUCCESS);
    151: }

    153: 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[])
    154: {
    155:   g0[0] = 1.0;
    156: }

    158: 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[])
    159: {
    160:   PetscInt d;
    161:   for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
    162: }

    164: /* data we seek to match */
    165: PetscErrorCode data_kernel(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *y, void *ctx)
    166: {
    167:   *y = 1.0 / (2 * PETSC_PI * PETSC_PI) * PetscSinReal(PETSC_PI * x[0]) * PetscSinReal(PETSC_PI * x[1]);
    168:   /* the associated control is sin(pi*x[0])*sin(pi*x[1]) */
    169:   return PETSC_SUCCESS;
    170: }
    171: PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
    172: {
    173:   *u = 0.0;
    174:   return PETSC_SUCCESS;
    175: }

    177: PetscErrorCode CreateCtx(DM dm, AppCtx *user)
    178: {
    179:   DM              dm_mass;
    180:   DM              dm_laplace;
    181:   PetscDS         prob_mass;
    182:   PetscDS         prob_laplace;
    183:   PetscFE         fe;
    184:   DMLabel         label;
    185:   PetscSection    section;
    186:   PetscInt        n, k, p, d;
    187:   PetscInt        dof, off;
    188:   IS              is;
    189:   const PetscInt *points;
    190:   const PetscInt  dim = 2;
    191:   PetscErrorCode (**wtf)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);

    193:   PetscFunctionBeginUser;
    194:   /* make the data we seek to match */
    195:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, PETSC_TRUE, NULL, 4, &fe));

    197:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
    198:   PetscCall(DMCreateDS(dm));
    199:   PetscCall(DMCreateGlobalVector(dm, &user->data));

    201:   /* ugh, this is hideous */
    202:   /* y_d = interpolate(Expression("sin(x[0]) + .."), V) */
    203:   PetscCall(PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &wtf));
    204:   wtf[0] = data_kernel;
    205:   PetscCall(DMProjectFunction(dm, 0.0, wtf, NULL, INSERT_VALUES, user->data));
    206:   PetscCall(PetscFree(wtf));

    208:   /* assemble(inner(u, v)*dx), almost */
    209:   PetscCall(DMClone(dm, &dm_mass));
    210:   PetscCall(DMCopyDisc(dm, dm_mass));
    211:   PetscCall(DMSetNumFields(dm_mass, 1));
    212:   PetscCall(DMPlexCopyCoordinates(dm, dm_mass)); /* why do I have to do this separately? */
    213:   PetscCall(DMGetDS(dm_mass, &prob_mass));
    214:   PetscCall(PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL, NULL, NULL));
    215:   PetscCall(PetscDSSetDiscretization(prob_mass, 0, (PetscObject)fe));
    216:   PetscCall(DMCreateMatrix(dm_mass, &user->mass));
    217:   PetscCall(DMPlexSNESComputeJacobianFEM(dm_mass, user->data, user->mass, user->mass, NULL));
    218:   PetscCall(MatSetOption(user->mass, MAT_SYMMETRIC, PETSC_TRUE));
    219:   PetscCall(DMDestroy(&dm_mass));

    221:   /* inner(grad(u), grad(v))*dx with homogeneous Dirichlet boundary conditions */
    222:   PetscCall(DMClone(dm, &dm_laplace));
    223:   PetscCall(DMCopyDisc(dm, dm_laplace));
    224:   PetscCall(DMSetNumFields(dm_laplace, 1));
    225:   PetscCall(DMPlexCopyCoordinates(dm, dm_laplace));
    226:   PetscCall(DMGetDS(dm_laplace, &prob_laplace));
    227:   PetscCall(PetscDSSetJacobian(prob_laplace, 0, 0, NULL, NULL, NULL, laplace_kernel));
    228:   PetscCall(PetscDSSetDiscretization(prob_laplace, 0, (PetscObject)fe));
    229:   PetscCall(DMCreateMatrix(dm_laplace, &user->laplace));
    230:   PetscCall(DMPlexSNESComputeJacobianFEM(dm_laplace, user->data, user->laplace, user->laplace, NULL));

    232:   PetscCall(DMGetLabel(dm_laplace, "marker", &label));
    233:   PetscCall(DMGetLocalSection(dm_laplace, &section));
    234:   PetscCall(DMLabelGetStratumSize(label, 1, &n));
    235:   PetscCall(DMLabelGetStratumIS(label, 1, &is));
    236:   PetscCall(ISGetIndices(is, &points));
    237:   user->num_bc_dofs = 0;
    238:   for (p = 0; p < n; ++p) {
    239:     PetscCall(PetscSectionGetDof(section, points[p], &dof));
    240:     user->num_bc_dofs += dof;
    241:   }
    242:   PetscCall(PetscMalloc1(user->num_bc_dofs, &user->bc_indices));
    243:   for (p = 0, k = 0; p < n; ++p) {
    244:     PetscCall(PetscSectionGetDof(section, points[p], &dof));
    245:     PetscCall(PetscSectionGetOffset(section, points[p], &off));
    246:     for (d = 0; d < dof; ++d) user->bc_indices[k++] = off + d;
    247:   }
    248:   PetscCall(ISRestoreIndices(is, &points));
    249:   PetscCall(ISDestroy(&is));
    250:   PetscCall(DMDestroy(&dm_laplace));

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

    258:   /* also create the KSP for solving the Laplace system */
    259:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->ksp_laplace));
    260:   PetscCall(KSPSetOperators(user->ksp_laplace, user->laplace, user->laplace));
    261:   PetscCall(KSPSetOptionsPrefix(user->ksp_laplace, "laplace_"));
    262:   PetscCall(KSPSetFromOptions(user->ksp_laplace));

    264:   /* A bit of setting up the user context */
    265:   user->dm = dm;
    266:   PetscCall(VecDuplicate(user->data, &user->state));
    267:   PetscCall(VecDuplicate(user->data, &user->adjoint));
    268:   PetscCall(VecDuplicate(user->data, &user->tmp1));
    269:   PetscCall(VecDuplicate(user->data, &user->tmp2));

    271:   PetscCall(PetscFEDestroy(&fe));
    272:   PetscFunctionReturn(PETSC_SUCCESS);
    273: }

    275: PetscErrorCode DestroyCtx(AppCtx *user)
    276: {
    277:   PetscFunctionBeginUser;
    278:   PetscCall(MatDestroy(&user->mass));
    279:   PetscCall(MatDestroy(&user->laplace));
    280:   PetscCall(KSPDestroy(&user->ksp_laplace));
    281:   PetscCall(VecDestroy(&user->data));
    282:   PetscCall(VecDestroy(&user->state));
    283:   PetscCall(VecDestroy(&user->adjoint));
    284:   PetscCall(VecDestroy(&user->tmp1));
    285:   PetscCall(VecDestroy(&user->tmp2));
    286:   PetscCall(PetscFree(user->bc_indices));
    287:   PetscCall(PetscFree(user->bc_values));
    288:   PetscFunctionReturn(PETSC_SUCCESS);
    289: }

    291: PetscErrorCode ReducedFunctionGradient(Tao tao, Vec u, PetscReal *func, Vec g, void *userv)
    292: {
    293:   AppCtx         *user  = (AppCtx *)userv;
    294:   const PetscReal alpha = 1.0e-6; /* regularisation parameter */
    295:   PetscReal       inner;

    297:   PetscFunctionBeginUser;
    298:   PetscCall(MatMult(user->mass, u, user->tmp1));
    299:   PetscCall(VecDot(u, user->tmp1, &inner)); /* regularisation contribution to */
    300:   *func = alpha * 0.5 * inner;              /* the functional                 */

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

    305:   /* Now compute the forward state. */
    306:   PetscCall(VecSetValues(user->tmp1, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES));
    307:   PetscCall(VecAssemblyBegin(user->tmp1));
    308:   PetscCall(VecAssemblyEnd(user->tmp1));
    309:   PetscCall(KSPSolve(user->ksp_laplace, user->tmp1, user->state)); /* forward solve */

    311:   /* Now compute the adjoint state also. */
    312:   PetscCall(VecCopy(user->state, user->tmp1));
    313:   PetscCall(VecAXPY(user->tmp1, -1.0, user->data));
    314:   PetscCall(MatMult(user->mass, user->tmp1, user->tmp2));
    315:   PetscCall(VecDot(user->tmp1, user->tmp2, &inner)); /* misfit contribution to */
    316:   *func += 0.5 * inner;                              /* the functional         */

    318:   PetscCall(VecSetValues(user->tmp2, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VALUES));
    319:   PetscCall(VecAssemblyBegin(user->tmp2));
    320:   PetscCall(VecAssemblyEnd(user->tmp2));
    321:   PetscCall(KSPSolve(user->ksp_laplace, user->tmp2, user->adjoint)); /* adjoint solve */

    323:   /* And bring it home with the gradient. */
    324:   PetscCall(MatMult(user->mass, user->adjoint, user->tmp1));
    325:   PetscCall(VecAXPY(g, 1.0, user->tmp1)); /* adjoint contribution to the gradient */
    326:   PetscFunctionReturn(PETSC_SUCCESS);
    327: }

    329: int main(int argc, char **argv)
    330: {
    331:   DM     dm;
    332:   Tao    tao;
    333:   Vec    u, lb, ub;
    334:   AppCtx user;

    336:   PetscFunctionBeginUser;
    337:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    338:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    339:   PetscCall(CreateCtx(dm, &user));

    341:   PetscCall(DMCreateGlobalVector(dm, &u));
    342:   PetscCall(VecSet(u, 0.0));
    343:   PetscCall(VecDuplicate(u, &lb));
    344:   PetscCall(VecDuplicate(u, &ub));
    345:   PetscCall(VecSet(lb, 0.0)); /* satisfied at the minimum anyway */
    346:   PetscCall(VecSet(ub, 0.8)); /* a nontrivial upper bound */

    348:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
    349:   PetscCall(TaoSetSolution(tao, u));
    350:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, ReducedFunctionGradient, &user));
    351:   PetscCall(TaoSetVariableBounds(tao, lb, ub));
    352:   PetscCall(TaoSetType(tao, TAOBLMVM));
    353:   PetscCall(TaoSetFromOptions(tao));

    355:   if (user.use_riesz) {
    356:     PetscCall(TaoLMVMSetH0(tao, user.mass)); /* crucial for mesh independence */
    357:     PetscCall(TaoSetGradientNorm(tao, user.mass));
    358:   }

    360:   PetscCall(TaoSolve(tao));

    362:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
    363:   PetscCall(VecViewFromOptions(u, NULL, "-sol_view"));

    365:   PetscCall(TaoDestroy(&tao));
    366:   PetscCall(DMDestroy(&dm));
    367:   PetscCall(VecDestroy(&u));
    368:   PetscCall(VecDestroy(&lb));
    369:   PetscCall(VecDestroy(&ub));
    370:   PetscCall(DestroyCtx(&user));
    371:   PetscCall(PetscFinalize());
    372:   return 0;
    373: }

    375: /*TEST

    377:     build:
    378:       requires: !complex !single

    380:     test:
    381:       requires: hdf5 double datafilespath !defined(PETSC_USE_64BIT_INDICES) hypre !cuda
    382:       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
    383:       filter: sed -e "s/-nan/nan/g"

    385:     test:
    386:       suffix: guess_pod
    387:       requires: double triangle
    388:       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
    389:       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"

    391: TEST*/