Actual source code: ex26.c

      1: static char help[] = "'Good Cop' Helmholtz Problem in 2d and 3d with finite elements.\n\
      2: We solve the 'Good Cop' Helmholtz problem in a rectangular\n\
      3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
      4: This example supports automatic convergence estimation\n\
      5: and coarse space adaptivity.\n\n\n";

      7: /*
      8:    The model problem:
      9:       Solve "Good Cop" Helmholtz equation on the unit square: (0,1) x (0,1)
     10:           - \Delta u + u = f,
     11:            where \Delta = Laplace operator
     12:       Dirichlet b.c.'s on all sides

     14: */

     16: #include <petscdmplex.h>
     17: #include <petscsnes.h>
     18: #include <petscds.h>
     19: #include <petscconvest.h>

     21: typedef struct {
     22:   PetscBool trig; /* Use trig function as exact solution */
     23: } AppCtx;

     25: /* For Primal Problem */
     26: static void g0_uu(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[])
     27: {
     28:   PetscInt d;
     29:   for (d = 0; d < dim; ++d) g0[0] = 1.0;
     30: }

     32: static void g3_uu(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[])
     33: {
     34:   PetscInt d;
     35:   for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
     36: }

     38: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     39: {
     40:   PetscInt d;
     41:   *u = 0.0;
     42:   for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]);
     43:   return PETSC_SUCCESS;
     44: }

     46: static PetscErrorCode quad_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     47: {
     48:   PetscInt d;
     49:   *u = 1.0;
     50:   for (d = 0; d < dim; ++d) *u += (d + 1) * PetscSqr(x[d]);
     51:   return PETSC_SUCCESS;
     52: }

     54: static void f0_trig_u(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 f0[])
     55: {
     56:   PetscInt d;
     57:   f0[0] += u[0];
     58:   for (d = 0; d < dim; ++d) f0[0] -= 4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]) + PetscSinReal(2.0 * PETSC_PI * x[d]);
     59: }

     61: static void f0_quad_u(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 f0[])
     62: {
     63:   PetscInt d;
     64:   switch (dim) {
     65:   case 1:
     66:     f0[0] = 1.0;
     67:     break;
     68:   case 2:
     69:     f0[0] = 5.0;
     70:     break;
     71:   case 3:
     72:     f0[0] = 11.0;
     73:     break;
     74:   default:
     75:     f0[0] = 5.0;
     76:     break;
     77:   }
     78:   f0[0] += u[0];
     79:   for (d = 0; d < dim; ++d) f0[0] -= (d + 1) * PetscSqr(x[d]);
     80: }

     82: static void f1_u(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 f1[])
     83: {
     84:   PetscInt d;
     85:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
     86: }

     88: static PetscErrorCode ProcessOptions(DM dm, AppCtx *options)
     89: {
     90:   MPI_Comm comm;
     91:   PetscInt dim;

     93:   PetscFunctionBeginUser;
     94:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
     95:   PetscCall(DMGetDimension(dm, &dim));
     96:   options->trig = PETSC_FALSE;

     98:   PetscOptionsBegin(comm, "", "Helmholtz Problem Options", "DMPLEX");
     99:   PetscCall(PetscOptionsBool("-exact_trig", "Use trigonometric exact solution (better for more complex finite elements)", "ex26.c", options->trig, &options->trig, NULL));
    100:   PetscOptionsEnd();
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
    105: {
    106:   PetscFunctionBeginUser;
    107:   PetscCall(DMCreate(comm, dm));
    108:   PetscCall(DMSetType(*dm, DMPLEX));
    109:   PetscCall(DMSetFromOptions(*dm));

    111:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
    112:   PetscCall(DMSetApplicationContext(*dm, user));
    113:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    114:   PetscFunctionReturn(PETSC_SUCCESS);
    115: }

    117: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
    118: {
    119:   PetscDS        ds;
    120:   DMLabel        label;
    121:   const PetscInt id = 1;

    123:   PetscFunctionBeginUser;
    124:   PetscCall(DMGetDS(dm, &ds));
    125:   PetscCall(DMGetLabel(dm, "marker", &label));
    126:   if (user->trig) {
    127:     PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
    128:     PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
    129:     PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
    130:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
    131:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Trig Exact Solution\n"));
    132:   } else {
    133:     PetscCall(PetscDSSetResidual(ds, 0, f0_quad_u, f1_u));
    134:     PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
    135:     PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user));
    136:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quad_u, NULL, user, NULL));
    137:   }
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
    142: {
    143:   DM             cdm = dm;
    144:   PetscFE        fe;
    145:   DMPolytopeType ct;
    146:   PetscBool      simplex;
    147:   PetscInt       dim, cStart;
    148:   char           prefix[PETSC_MAX_PATH_LEN];

    150:   PetscFunctionBeginUser;
    151:   PetscCall(DMGetDimension(dm, &dim));

    153:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
    154:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
    155:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
    156:   /* Create finite element */
    157:   PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
    158:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
    159:   PetscCall(PetscObjectSetName((PetscObject)fe, name));
    160:   /* Set discretization and boundary conditions for each mesh */
    161:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
    162:   PetscCall(DMCreateDS(dm));
    163:   PetscCall((*setup)(dm, user));
    164:   while (cdm) {
    165:     PetscCall(DMCopyDisc(dm, cdm));
    166:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    167:   }
    168:   PetscCall(PetscFEDestroy(&fe));
    169:   PetscFunctionReturn(PETSC_SUCCESS);
    170: }

    172: int main(int argc, char **argv)
    173: {
    174:   DM      dm; /* Problem specification */
    175:   PetscDS ds;
    176:   SNES    snes; /* Nonlinear solver */
    177:   Vec     u;    /* Solutions */
    178:   AppCtx  user; /* User-defined work context */

    180:   PetscFunctionBeginUser;
    181:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    182:   /* Primal system */
    183:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    184:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    185:   PetscCall(ProcessOptions(dm, &user));
    186:   PetscCall(SNESSetDM(snes, dm));
    187:   PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
    188:   PetscCall(DMCreateGlobalVector(dm, &u));
    189:   PetscCall(VecSet(u, 0.0));
    190:   PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
    191:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
    192:   PetscCall(SNESSetFromOptions(snes));
    193:   PetscCall(DMSNESCheckFromOptions(snes, u));

    195:   /* Looking for field error */
    196:   PetscInt Nfields;
    197:   PetscCall(DMGetDS(dm, &ds));
    198:   PetscCall(PetscDSGetNumFields(ds, &Nfields));
    199:   PetscCall(SNESSolve(snes, NULL, u));
    200:   PetscCall(SNESGetSolution(snes, &u));
    201:   PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));

    203:   /* Cleanup */
    204:   PetscCall(VecDestroy(&u));
    205:   PetscCall(SNESDestroy(&snes));
    206:   PetscCall(DMDestroy(&dm));
    207:   PetscCall(PetscFinalize());
    208:   return 0;
    209: }

    211: /*TEST
    212: test:
    213:   # L_2 convergence rate: 1.9
    214:   suffix: 2d_p1_conv
    215:   requires: triangle
    216:   args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
    217: test:
    218:   # L_2 convergence rate: 1.9
    219:   suffix: 2d_q1_conv
    220:   args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
    221: test:
    222:   # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
    223:   suffix: 3d_p1_conv
    224:   requires: ctetgen
    225:   args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu
    226: test:
    227:   # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
    228:   suffix: 3d_q1_conv
    229:   args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu
    230: test:
    231:   # L_2 convergence rate: 1.9
    232:   suffix: 2d_p1_trig_conv
    233:   requires: triangle
    234:   args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
    235: test:
    236:   # L_2 convergence rate: 1.9
    237:   suffix: 2d_q1_trig_conv
    238:   args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
    239: test:
    240:   # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
    241:   suffix: 3d_p1_trig_conv
    242:   requires: ctetgen
    243:   args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig
    244: test:
    245:   # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
    246:   suffix: 3d_q1_trig_conv
    247:   args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig
    248: test:
    249:   suffix: 2d_p1_gmg_vcycle
    250:   requires: triangle
    251:   args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
    252:     -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
    253:     -mg_levels_ksp_max_it 1 \
    254:     -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
    255: test:
    256:   suffix: 2d_p1_gmg_fcycle
    257:   requires: triangle
    258:   args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
    259:     -ksp_type cg -ksp_rtol 1e-10 -pc_type mg -pc_mg_type full \
    260:     -mg_levels_ksp_max_it 2 \
    261:     -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
    262: test:
    263:   suffix: 2d_p1_gmg_vcycle_trig
    264:   requires: triangle
    265:   args: -exact_trig -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
    266:     -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
    267:     -mg_levels_ksp_max_it 1 \
    268:     -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
    269: test:
    270:   suffix: 2d_q3_cgns
    271:   requires: cgns
    272:   args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_box_faces 3,3 -snes_view_solution cgns:sol.cgns -potential_petscspace_degree 3 -dm_coord_petscspace_degree 3
    273:   output_file: output/empty.out
    274: TEST*/