Actual source code: ex8.c

      1: static const char help[] = "\
      2: Minimum surface area problem in 2D using DMPLEX.\n\
      3: It solves the unconstrained minimization problem \n\
      4: \n\
      5: argmin \\int_\\Omega (1 + ||u||^2), u = g on \\partial\\Omega.\n\
      6: \n\
      7: This example shows how to setup an optimization problem using DMPLEX FEM routines.\n\
      8: It supports nonlinear domain decomposition and multilevel solvers.\n\
      9: \n\n";

     11: #include <petscdmplex.h>
     12: #include <petscsnes.h>
     13: #include <petscds.h>

     15: /* The pointwise evaluation routine for the objective function */
     16: static void objective_2d(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
     17: {
     18:   const PetscScalar ux2 = PetscSqr(u_x[0]);
     19:   const PetscScalar uy2 = PetscSqr(u_x[1]);

     21:   obj[0] = PetscSqrtReal(PetscAbsScalar(1 + ux2 + uy2));
     22: }

     24: /* The pointwise evaluation routine for the gradient wrt the gradients of the FEM basis */
     25: static void gradient_1_2d(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
     26: {
     27:   const PetscScalar ux2 = PetscSqr(u_x[0]);
     28:   const PetscScalar uy2 = PetscSqr(u_x[1]);
     29:   const PetscScalar v   = 1. / PetscSqrtReal(PetscAbsScalar(1 + ux2 + uy2));

     31:   f[0] = v * u_x[0];
     32:   f[1] = v * u_x[1];
     33: }

     35: /* The pointwise evaluation routine for the hessian wrt the gradients of the FEM basis */
     36: static void hessian_11_2d(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 jac[])
     37: {
     38:   const PetscScalar ux2 = PetscSqr(u_x[0]);
     39:   const PetscScalar uy2 = PetscSqr(u_x[1]);
     40:   const PetscScalar uxy = u_x[0] * u_x[1];
     41:   const PetscScalar v1  = 1. / PetscSqrtReal(PetscAbsScalar(1 + ux2 + uy2));
     42:   const PetscScalar v2  = v1 / (1 + ux2 + uy2);

     44:   jac[0] = v1 - v2 * ux2;
     45:   jac[1] = -v2 * uxy;
     46:   jac[2] = -v2 * uxy;
     47:   jac[3] = v1 - v2 * uy2;
     48: }

     50: /* The boundary conditions we impose */
     51: static PetscErrorCode sins_2d(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
     52: {
     53:   PetscFunctionBeginUser;
     54:   const PetscReal pi2 = PETSC_PI * 2;
     55:   const PetscReal x   = xx[0];
     56:   const PetscReal y   = xx[1];

     58:   *u = (y - 0.5) * PetscSinReal(pi2 * x) + (x - 0.5) * PetscSinReal(pi2 * y);
     59:   PetscFunctionReturn(PETSC_SUCCESS);
     60: }

     62: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
     63: {
     64:   DM      plex;
     65:   DMLabel label;

     67:   PetscFunctionBeginUser;
     68:   PetscCall(DMCreateLabel(dm, name));
     69:   PetscCall(DMGetLabel(dm, name, &label));
     70:   PetscCall(DMConvert(dm, DMPLEX, &plex));
     71:   PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label));
     72:   PetscCall(DMDestroy(&plex));
     73:   PetscFunctionReturn(PETSC_SUCCESS);
     74: }

     76: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
     77: {
     78:   PetscFunctionBeginUser;
     79:   PetscCall(DMCreate(comm, dm));
     80:   PetscCall(DMSetType(*dm, DMPLEX));
     81:   PetscCall(DMSetFromOptions(*dm));
     82:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
     83:   PetscFunctionReturn(PETSC_SUCCESS);
     84: }

     86: static PetscErrorCode SetupProblem(DM dm)
     87: {
     88:   PetscDS        ds;
     89:   DMLabel        label;
     90:   PetscInt       dim;
     91:   const PetscInt id = 1;

     93:   PetscFunctionBeginUser;
     94:   PetscCall(DMGetDS(dm, &ds));
     95:   PetscCall(DMGetDimension(dm, &dim));
     96:   PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only for 2-D");
     97:   PetscCall(PetscDSSetObjective(ds, 0, objective_2d));
     98:   PetscCall(PetscDSSetResidual(ds, 0, NULL, gradient_1_2d));
     99:   PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, hessian_11_2d));
    100:   PetscCall(DMGetLabel(dm, "marker", &label));
    101:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "data", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sins_2d, NULL, NULL, NULL));
    102:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_TRUE, NULL));
    103:   PetscFunctionReturn(PETSC_SUCCESS);
    104: }

    106: static PetscErrorCode SetupDiscretization(DM dm)
    107: {
    108:   DM        plex, cdm = dm;
    109:   PetscFE   fe;
    110:   PetscBool simplex;
    111:   PetscInt  dim;

    113:   PetscFunctionBeginUser;
    114:   PetscCall(DMGetDimension(dm, &dim));
    115:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    116:   PetscCall(DMPlexIsSimplex(plex, &simplex));
    117:   PetscCall(DMDestroy(&plex));
    118:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe));
    119:   PetscCall(PetscObjectSetName((PetscObject)fe, "potential"));
    120:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
    121:   PetscCall(DMCreateDS(dm));
    122:   PetscCall(SetupProblem(dm));
    123:   while (cdm) {
    124:     PetscBool hasLabel;

    126:     PetscCall(DMHasLabel(cdm, "marker", &hasLabel));
    127:     if (!hasLabel) PetscCall(CreateBCLabel(cdm, "marker"));
    128:     PetscCall(DMCopyDisc(dm, cdm));
    129:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    130:   }
    131:   PetscCall(PetscFEDestroy(&fe));
    132:   PetscFunctionReturn(PETSC_SUCCESS);
    133: }

    135: int main(int argc, char **argv)
    136: {
    137:   DM   dm;   /* Problem specification */
    138:   SNES snes; /* nonlinear solver */
    139:   Vec  u;    /* solution vector */

    141:   PetscFunctionBeginUser;
    142:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    143:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    144:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
    145:   PetscCall(SNESSetDM(snes, dm));

    147:   PetscCall(SetupDiscretization(dm));

    149:   PetscCall(SNESSetFromOptions(snes));

    151:   PetscCall(DMCreateGlobalVector(dm, &u));
    152:   PetscCall(SNESSolve(snes, NULL, u));

    154:   PetscCall(VecDestroy(&u));
    155:   PetscCall(SNESDestroy(&snes));
    156:   PetscCall(DMDestroy(&dm));
    157:   PetscCall(PetscFinalize());
    158:   return 0;
    159: }

    161: /*TEST

    163:   test:
    164:     requires: !single
    165:     suffix: qn_nasm
    166:     filter: sed -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"
    167:     nsize: 4
    168:     args: -petscspace_degree 1 -dm_refine 2 -snes_type qn -snes_npc_side {{left right}separate output} -npc_snes_type nasm -npc_snes_nasm_type restrict -npc_sub_snes_linesearch_order 1 -npc_sub_snes_linesearch_type bt -dm_plex_dd_overlap 1 -snes_linesearch_type bt -snes_linesearch_order 1 -npc_sub_pc_type lu -npc_sub_ksp_type preonly -snes_converged_reason -snes_monitor_short -petscpartitioner_type simple -npc_sub_snes_max_it 1 -dm_plex_simplex 0 -snes_rtol 1.e-6

    170: TEST*/