Actual source code: ex52.c

      1: static char help[] = "Simple Advection-diffusion equation solved using FVM in DMPLEX\n";

      3: /*
      4:    Solves the simple advection equation given by

      6:    q_t + u (q_x) + v (q_y) - D (q_xx + q_yy) = 0 using FVM and First Order Upwind discretization.

      8:    with a user defined initial condition.

     10:    with dirichlet/neumann conditions on the four boundaries of the domain.

     12:    User can define the mesh parameters either in the command line or inside
     13:    the ProcessOptions() routine.

     15:    Contributed by: Mukkund Sunjii, Domenico Lahaye
     16: */

     18: #include <petscdmplex.h>
     19: #include <petscts.h>
     20: #include <petscblaslapack.h>

     22: #if defined(PETSC_HAVE_CGNS)
     23: #undef I
     24: #include <cgnslib.h>
     25: #endif
     26: /*
     27:    User-defined routines
     28: */
     29: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec);
     30: extern PetscErrorCode MyTSMonitor(TS, PetscInt, PetscReal, Vec, void *);
     31: extern PetscErrorCode MySNESMonitor(SNES, PetscInt, PetscReal, PetscViewerAndFormat *);

     33: /* Defining the usr defined context */
     34: typedef struct {
     35:     DM          da;
     36:     PetscBool   interpolate;                  /* Generate intermediate mesh elements */
     37:     char        filename[PETSC_MAX_PATH_LEN]; /* Mesh filename */
     38:     PetscInt    dim;
     39:     PetscScalar diffusion;
     40:     PetscReal   u, v;
     41:     PetscScalar delta_x, delta_y;
     42:     PetscInt    cells[2];
     43: } AppCtx;

     45: /* Options for the scenario */
     46: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
     47: {

     51:     options->interpolate = PETSC_TRUE;
     52:     options->filename[0] = '\0';
     53:     options->dim = 2;
     54:     options->u = 2.5;
     55:     options->v = 0.0;
     56:     options->cells[0] = 20;
     57:     options->cells[1] = 20;
     58:     options->diffusion = 0.0;

     60:     PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
     61:     PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "advection_DMPLEX.c",options->interpolate, &options->interpolate, NULL);
     62:     PetscOptionsString("-filename", "The mesh file", "advection_DMPLEX.c", options->filename, options->filename, sizeof(options->filename), NULL);
     63:     PetscOptionsRangeInt("-dim", "Problem dimension used for the non-file mesh.", "ex7.c", options->dim, &options->dim, NULL,1,3);
     64:     PetscOptionsReal("-u", "The x component of the convective coefficient", "advection_DMPLEX.c", options->u, &options->u, NULL);
     65:     PetscOptionsReal("-v", "The y component of the convective coefficient", "advection_DMPLEX.c", options->v, &options->v, NULL);
     66:     PetscOptionsScalar("-diffus", "The diffusive coefficient", "advection_DMPLEX.c", options->diffusion, &options->diffusion, NULL);
     67:     PetscOptionsEnd();
     68:     return(0);
     69: }

     71: /*
     72:   User can provide the file containing the mesh.
     73:   Or can generate the mesh using DMPlexCreateBoxMesh with the specified options.
     74: */
     75: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
     76: {
     77:     size_t len;
     80:     PetscStrlen(user->filename, &len);
     81:     if (!len) {
     82:         DMLabel label;
     83:         DMPlexCreateBoxMesh(comm, user->dim, PETSC_FALSE, user->cells, NULL, NULL, NULL, user->interpolate, dm);
     84:         /* Mark boundary and set BC */
     85:         DMCreateLabel(*dm, "boundary");
     86:         DMGetLabel(*dm, "boundary", &label);
     87:         DMPlexMarkBoundaryFaces(*dm, 1, label);
     88:         DMPlexLabelComplete(*dm, label);
     89:     } else {
     90:         DMPlexCreateFromFile(comm, user->filename, user->interpolate, dm);
     91:     }
     92:     PetscObjectSetName((PetscObject) * dm, "Mesh");
     93:     DMSetFromOptions(*dm);
     94:     DMViewFromOptions(*dm, NULL, "-dm_view");
     95:     return(0);
     96: }

     98:     /* This routine is responsible for defining the local solution vector x
     99:     with a given initial solution.
    100:     The initial solution can be modified accordingly inside the loops.
    101:     No need for a local vector because there is exchange of information
    102:     across the processors. Unlike for FormFunction which depends on the neighbours */
    103: PetscErrorCode FormInitialSolution(DM da, Vec U)
    104: {
    106:     PetscScalar    *u;
    107:     PetscInt       cell, cStart, cEnd;
    108:     PetscReal      cellvol, centroid[3], normal[3];

    111:     /* Get pointers to vector data */
    112:     VecGetArray(U, &u);
    113:     /* Get local grid boundaries */
    114:     DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
    115:     /* Assigning the values at the cell centers based on x and y directions */
    116:     for (cell = cStart; cell < cEnd; cell++) {
    117:         DMPlexComputeCellGeometryFVM(da, cell, &cellvol, centroid, normal);
    118:         if (centroid[0] > 0.9 && centroid[0] < 0.95) {
    119:             if (centroid[1] > 0.9 && centroid[1] < 0.95) u[cell] = 2.0;
    120:         }
    121:         else u[cell] = 0;
    122:     }
    123:     VecRestoreArray(U, &u);
    124:     return(0);
    125: }

    127: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx)
    128: {
    130:     PetscReal      norm;
    131:     MPI_Comm       comm;

    134:     if (step < 0) return(0); /* step of -1 indicates an interpolated solution */
    135:     VecNorm(v, NORM_2, &norm);
    136:     PetscObjectGetComm((PetscObject) ts, &comm);
    137:     PetscPrintf(comm, "timestep %D time %g norm %g\n", step, (double) ptime, (double) norm);
    138:     return(0);
    139: }

    141: /*
    142:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
    143:    Input Parameters:
    144:      snes - the SNES context
    145:      its - iteration number
    146:      fnorm - 2-norm function value (may be estimated)
    147:      ctx - optional user-defined context for private data for the
    148:          monitor routine, as set by SNESMonitorSet()
    149: */
    150: PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
    151: {

    155:     SNESMonitorDefaultShort(snes, its, fnorm, vf);
    156:     return(0);
    157: }

    159: /*
    160:    FormFunction - Evaluates nonlinear function, F(x).

    162:    Input Parameters:
    163: .  ts - the TS context
    164: .  X - input vector
    165: .  ctx - optional user-defined context, as set by SNESSetFunction()

    167:    Output Parameter:
    168: .  F - function vector
    169:  */
    170: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ctx)
    171: {
    172:     AppCtx *user = (AppCtx *) ctx;
    173:     DM da = (DM) user->da;
    175:     PetscScalar *x, *f;
    176:     Vec localX;
    177:     PetscInt fStart, fEnd, nF;
    178:     PetscInt cell, cStart, cEnd, nC;
    179:     DM dmFace;      /* DMPLEX for face geometry */
    180:     PetscFV fvm;                /* specify type of FVM discretization */
    181:     Vec cellGeom, faceGeom; /* vector of structs related to cell/face geometry*/
    182:     const PetscScalar *fgeom;             /* values stored in the vector facegeom */
    183:     PetscFVFaceGeom *fgA;               /* struct with face geometry information */
    184:     const PetscInt *cellcone, *cellsupport;
    185:     PetscScalar flux_east, flux_west, flux_north, flux_south, flux_centre;
    186:     PetscScalar centroid_x[2], centroid_y[2], boundary = 0.0;
    187:     PetscScalar boundary_left = 0.0;
    188:     PetscReal u_plus, u_minus, v_plus, v_minus, zero = 0.0;
    189:     PetscScalar delta_x, delta_y;

    191:     /* Get the local vector from the DM object. */
    193:     DMGetLocalVector(da, &localX);

    195:     /* Scatter ghost points to local vector,using the 2-step process
    196:        DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). */
    197:     DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
    198:     DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
    199:     /* Get pointers to vector data. */
    200:     VecGetArray(localX, &x);
    201:     VecGetArray(F, &f);

    203:     /* Obtaining local cell and face ownership */
    204:     DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
    205:     DMPlexGetHeightStratum(da, 1, &fStart, &fEnd);

    207:     /* Creating the PetscFV object to obtain face and cell geometry.
    208:     Later to be used to compute face centroid to find cell widths. */

    210:     PetscFVCreate(PETSC_COMM_WORLD, &fvm);
    211:     PetscFVSetType(fvm, PETSCFVUPWIND);
    212:     /*....Retrieve precomputed cell geometry....*/
    213:     DMPlexGetDataFVM(da, fvm, &cellGeom, &faceGeom, NULL);
    214:     VecGetDM(faceGeom, &dmFace);
    215:     VecGetArrayRead(faceGeom, &fgeom);

    217:     /* Spanning through all the cells and an inner loop through the faces. Find the
    218:     face neighbors and pick the upwinded cell value for flux. */

    220:     u_plus = PetscMax(user->u, zero);
    221:     u_minus = PetscMin(user->u, zero);
    222:     v_plus = PetscMax(user->v, zero);
    223:     v_minus = PetscMin(user->v, zero);

    225:     for (cell = cStart; cell < cEnd; cell++) {
    226:         /* Obtaining the faces of the cell */
    227:         DMPlexGetConeSize(da, cell, &nF);
    228:         DMPlexGetCone(da, cell, &cellcone);

    230:         /* south */
    231:         DMPlexPointLocalRead(dmFace, cellcone[0], fgeom, &fgA);
    232:         centroid_y[0] = fgA->centroid[1];
    233:         /* North */
    234:         DMPlexPointLocalRead(dmFace, cellcone[2], fgeom, &fgA);
    235:         centroid_y[1] = fgA->centroid[1];
    236:         /* West */
    237:         DMPlexPointLocalRead(dmFace, cellcone[3], fgeom, &fgA);
    238:         centroid_x[0] = fgA->centroid[0];
    239:         /* East */
    240:         DMPlexPointLocalRead(dmFace, cellcone[1], fgeom, &fgA);
    241:         centroid_x[1] = fgA->centroid[0];

    243:         /* Computing the cell widths in the x and y direction */
    244:         delta_x = centroid_x[1] - centroid_x[0];
    245:         delta_y = centroid_y[1] - centroid_y[0];

    247:         /* Getting the neighbors of each face
    248:            Going through the faces by the order (cellcone) */

    250:         /* cellcone[0] - south */
    251:         DMPlexGetSupportSize(da, cellcone[0], &nC);
    252:         DMPlexGetSupport(da, cellcone[0], &cellsupport);
    253:         if (nC == 2) flux_south = (x[cellsupport[0]] * (-v_plus - user->diffusion * delta_x)) / delta_y;
    254:         else flux_south = (boundary * (-v_plus - user->diffusion * delta_x)) / delta_y;

    256:         /* cellcone[1] - east */
    257:         DMPlexGetSupportSize(da, cellcone[1], &nC);
    258:         DMPlexGetSupport(da, cellcone[1], &cellsupport);
    259:         if (nC == 2) flux_east = (x[cellsupport[1]] * (u_minus - user->diffusion * delta_y)) / delta_x;
    260:         else flux_east = (boundary * (u_minus - user->diffusion * delta_y)) / delta_x;

    262:         /* cellcone[2] - north */
    263:         DMPlexGetSupportSize(da, cellcone[2], &nC);
    264:         DMPlexGetSupport(da, cellcone[2], &cellsupport);
    265:         if (nC == 2) flux_north = (x[cellsupport[1]] * (v_minus - user->diffusion * delta_x)) / delta_y;
    266:         else flux_north = (boundary * (v_minus - user->diffusion * delta_x)) / delta_y;

    268:         /* cellcone[3] - west */
    269:         DMPlexGetSupportSize(da, cellcone[3], &nC);
    270:         DMPlexGetSupport(da, cellcone[3], &cellsupport);
    271:         if (nC == 2) flux_west = (x[cellsupport[0]] * (-u_plus - user->diffusion * delta_y)) / delta_x;
    272:         else flux_west = (boundary_left * (-u_plus - user->diffusion * delta_y)) / delta_x;

    274:         /* Contribution by the cell to the fluxes */
    275:         flux_centre = x[cell] * ((u_plus - u_minus + 2 * user->diffusion * delta_y) / delta_x +
    276:                                  (v_plus - v_minus + 2 * user->diffusion * delta_x) / delta_y);

    278:         /* Calculating the net flux for each cell
    279:            and computing the RHS time derivative f[.] */
    280:         f[cell] = -(flux_centre + flux_east + flux_west + flux_north + flux_south);
    281:     }
    282:     PetscFVDestroy(&fvm);
    283:     VecRestoreArray(localX, &x);
    284:     VecRestoreArray(F, &f);
    285:     DMRestoreLocalVector(da, &localX);
    286:     return(0);
    287: }

    289: int main(int argc, char **argv)
    290: {
    291:     TS                   ts;                         /* time integrator */
    292:     SNES                 snes;
    293:     Vec                  x, r;                        /* solution, residual vectors */
    294:     PetscErrorCode       ierr;
    295:     DM                   da;
    296:     PetscMPIInt          rank;
    297:     PetscViewerAndFormat *vf;
    298:     AppCtx               user;                             /* mesh context */
    299:     PetscInt             numFields = 1, numBC, i;
    300:     PetscInt             numComp[1];
    301:     PetscInt             numDof[12];
    302:     PetscInt             bcField[1];
    303:     PetscSection         section;
    304:     IS                   bcPointIS[1];

    306:     /* Initialize program */
    307:     PetscInitialize(&argc, &argv, (char *) 0, help);if (ierr) return ierr;
    308:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
    309:     /* Create distributed array (DMPLEX) to manage parallel grid and vectors */
    310:     ProcessOptions(PETSC_COMM_WORLD, &user);
    311:     CreateMesh(PETSC_COMM_WORLD, &user, &da);

    313:     /* Specifying the fields and dof for the formula through PETSc Section
    314:     Create a scalar field u with 1 component on cells, faces and edges.
    315:     Alternatively, the field information could be added through a PETSCFV object
    316:     using DMAddField(...).*/
    317:     numComp[0] = 1;

    319:     for (i = 0; i < numFields * (user.dim + 1); ++i) numDof[i] = 0;

    321:     numDof[0 * (user.dim + 1)] = 1;
    322:     numDof[0 * (user.dim + 1) + user.dim - 1] = 1;
    323:     numDof[0 * (user.dim + 1) + user.dim] = 1;

    325:     /* Setup boundary conditions */
    326:     numBC = 1;
    327:     /* Prescribe a Dirichlet condition on u on the boundary
    328:        Label "marker" is made by the mesh creation routine  */
    329:     bcField[0] = 0;
    330:     DMGetStratumIS(da, "marker", 1, &bcPointIS[0]);

    332:     /* Create a PetscSection with this data layout */
    333:     DMSetNumFields(da, numFields);
    334:     DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section);

    336:     /* Name the Field variables */
    337:     PetscSectionSetFieldName(section, 0, "u");

    339:     /* Tell the DM to use this section (with the specified fields and dof) */
    340:     DMSetLocalSection(da, section);
    341:     user.da = da;

    343:     /* Extract global vectors from DMDA; then duplicate for remaining
    344:        vectors that are the same types */

    346:     /* Create a Vec with this layout and view it */
    347:     DMGetGlobalVector(da, &x);
    348:     VecDuplicate(x, &r);

    350:     /* Create timestepping solver context */
    351:     TSCreate(PETSC_COMM_WORLD, &ts);
    352:     TSSetProblemType(ts, TS_NONLINEAR);
    353:     TSSetRHSFunction(ts, NULL, FormFunction, &user);

    355:     TSSetMaxTime(ts, 1.0);
    356:     TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
    357:     TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL);
    358:     TSSetDM(ts, da);

    360:     /* Customize nonlinear solver */
    361:     TSSetType(ts, TSEULER);
    362:     TSGetSNES(ts, &snes);
    363:     PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf);
    364:     SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *)) MySNESMonitor, vf,(PetscErrorCode (*)(void **)) PetscViewerAndFormatDestroy);

    366:      /* Set initial conditions */
    367:     FormInitialSolution(da, x);
    368:     TSSetTimeStep(ts, .0001);
    369:     TSSetSolution(ts, x);
    370:     /* Set runtime options */
    371:     TSSetFromOptions(ts);
    372:     /* Solve nonlinear system */
    373:     TSSolve(ts, x);

    375:     /* Clean up routine */
    376:     DMRestoreGlobalVector(da, &x);
    377:     ISDestroy(&bcPointIS[0]);
    378:     PetscSectionDestroy(&section);
    379:     VecDestroy(&r);
    380:     TSDestroy(&ts);
    381:     DMDestroy(&da);
    382:     PetscFinalize();
    383:     return ierr;
    384: }

    386: /*TEST

    388:     test:
    389:       suffix: 0
    390:       args: -ts_max_steps 5 -ts_type rk
    391:       requires: !single !complex triangle ctetgen

    393: TEST*/