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:   PetscScalar diffusion;
     36:   PetscReal   u, v;
     37:   PetscScalar delta_x, delta_y;
     38: } AppCtx;

     40: /* Options for the scenario */
     41: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
     42: {
     43:   PetscFunctionBeginUser;
     44:   options->u         = 2.5;
     45:   options->v         = 0.0;
     46:   options->diffusion = 0.0;
     47:   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
     48:   PetscCall(PetscOptionsReal("-u", "The x component of the convective coefficient", "advection_DMPLEX.c", options->u, &options->u, NULL));
     49:   PetscCall(PetscOptionsReal("-v", "The y component of the convective coefficient", "advection_DMPLEX.c", options->v, &options->v, NULL));
     50:   PetscCall(PetscOptionsScalar("-diffus", "The diffusive coefficient", "advection_DMPLEX.c", options->diffusion, &options->diffusion, NULL));
     51:   PetscOptionsEnd();
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: /*
     56:   User can provide the file containing the mesh.
     57:   Or can generate the mesh using DMPlexCreateBoxMesh with the specified options.
     58: */
     59: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
     60: {
     61:   PetscFunctionBeginUser;
     62:   PetscCall(DMCreate(comm, dm));
     63:   PetscCall(DMSetType(*dm, DMPLEX));
     64:   PetscCall(DMSetFromOptions(*dm));
     65:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
     66:   {
     67:     DMLabel label;
     68:     PetscCall(DMGetLabel(*dm, "boundary", &label));
     69:     PetscCall(DMPlexLabelComplete(*dm, label));
     70:   }
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: /* This routine is responsible for defining the local solution vector x
     75:     with a given initial solution.
     76:     The initial solution can be modified accordingly inside the loops.
     77:     No need for a local vector because there is exchange of information
     78:     across the processors. Unlike for FormFunction which depends on the neighbours */
     79: PetscErrorCode FormInitialSolution(DM da, Vec U)
     80: {
     81:   PetscScalar *u;
     82:   PetscInt     cell, cStart, cEnd;
     83:   PetscReal    cellvol, centroid[3], normal[3];

     85:   PetscFunctionBeginUser;
     86:   /* Get pointers to vector data */
     87:   PetscCall(VecGetArray(U, &u));
     88:   /* Get local grid boundaries */
     89:   PetscCall(DMPlexGetHeightStratum(da, 0, &cStart, &cEnd));
     90:   /* Assigning the values at the cell centers based on x and y directions */
     91:   PetscCall(DMGetCoordinatesLocalSetUp(da));
     92:   for (cell = cStart; cell < cEnd; cell++) {
     93:     PetscCall(DMPlexComputeCellGeometryFVM(da, cell, &cellvol, centroid, normal));
     94:     if (centroid[0] > 0.9 && centroid[0] < 0.95) {
     95:       if (centroid[1] > 0.9 && centroid[1] < 0.95) u[cell] = 2.0;
     96:     } else u[cell] = 0;
     97:   }
     98:   PetscCall(VecRestoreArray(U, &u));
     99:   PetscFunctionReturn(PETSC_SUCCESS);
    100: }

    102: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx)
    103: {
    104:   PetscReal norm;
    105:   MPI_Comm  comm;

    107:   PetscFunctionBeginUser;
    108:   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* step of -1 indicates an interpolated solution */
    109:   PetscCall(VecNorm(v, NORM_2, &norm));
    110:   PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
    111:   PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm));
    112:   PetscFunctionReturn(PETSC_SUCCESS);
    113: }

    115: /*
    116:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
    117:    Input Parameters:
    118:      snes - the SNES context
    119:      its - iteration number
    120:      fnorm - 2-norm function value (may be estimated)
    121:      ctx - optional user-defined context for private data for the
    122:          monitor routine, as set by SNESMonitorSet()
    123: */
    124: PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
    125: {
    126:   PetscFunctionBeginUser;
    127:   PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf));
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

    131: /*
    132:    FormFunction - Evaluates nonlinear function, F(x).

    134:    Input Parameters:
    135: .  ts - the TS context
    136: .  X - input vector
    137: .  ctx - optional user-defined context, as set by SNESSetFunction()

    139:    Output Parameter:
    140: .  F - function vector
    141:  */
    142: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ctx)
    143: {
    144:   AppCtx            *user = (AppCtx *)ctx;
    145:   DM                 da;
    146:   PetscScalar       *x, *f;
    147:   Vec                localX;
    148:   PetscInt           fStart, fEnd, nF;
    149:   PetscInt           cell, cStart, cEnd, nC;
    150:   DM                 dmFace;             /* DMPLEX for face geometry */
    151:   PetscFV            fvm;                /* specify type of FVM discretization */
    152:   Vec                cellGeom, faceGeom; /* vector of structs related to cell/face geometry*/
    153:   const PetscScalar *fgeom;              /* values stored in the vector facegeom */
    154:   PetscFVFaceGeom   *fgA;                /* struct with face geometry information */
    155:   const PetscInt    *cellcone, *cellsupport;
    156:   PetscScalar        flux_east, flux_west, flux_north, flux_south, flux_centre;
    157:   PetscScalar        centroid_x[2], centroid_y[2], boundary = 0.0;
    158:   PetscScalar        boundary_left = 0.0;
    159:   PetscReal          u_plus, u_minus, v_plus, v_minus, zero = 0.0;
    160:   PetscScalar        delta_x, delta_y;

    162:   /* Get the local vector from the DM object. */
    163:   PetscFunctionBeginUser;
    164:   PetscCall(TSGetDM(ts, &da));
    165:   PetscCall(DMGetLocalVector(da, &localX));

    167:   /* Scatter ghost points to local vector,using the 2-step process
    168:        DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). */
    169:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
    170:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
    171:   /* Get pointers to vector data. */
    172:   PetscCall(VecGetArray(localX, &x));
    173:   PetscCall(VecGetArray(F, &f));

    175:   /* Obtaining local cell and face ownership */
    176:   PetscCall(DMPlexGetHeightStratum(da, 0, &cStart, &cEnd));
    177:   PetscCall(DMPlexGetHeightStratum(da, 1, &fStart, &fEnd));

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

    182:   PetscCall(PetscFVCreate(PETSC_COMM_WORLD, &fvm));
    183:   PetscCall(PetscFVSetType(fvm, PETSCFVUPWIND));
    184:   /*....Retrieve precomputed cell geometry....*/
    185:   PetscCall(DMPlexGetDataFVM(da, fvm, &cellGeom, &faceGeom, NULL));
    186:   PetscCall(VecGetDM(faceGeom, &dmFace));
    187:   PetscCall(VecGetArrayRead(faceGeom, &fgeom));

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

    192:   u_plus  = PetscMax(user->u, zero);
    193:   u_minus = PetscMin(user->u, zero);
    194:   v_plus  = PetscMax(user->v, zero);
    195:   v_minus = PetscMin(user->v, zero);

    197:   for (cell = cStart; cell < cEnd; cell++) {
    198:     /* Obtaining the faces of the cell */
    199:     PetscCall(DMPlexGetConeSize(da, cell, &nF));
    200:     PetscCall(DMPlexGetCone(da, cell, &cellcone));

    202:     /* south */
    203:     PetscCall(DMPlexPointLocalRead(dmFace, cellcone[0], fgeom, &fgA));
    204:     centroid_y[0] = fgA->centroid[1];
    205:     /* North */
    206:     PetscCall(DMPlexPointLocalRead(dmFace, cellcone[2], fgeom, &fgA));
    207:     centroid_y[1] = fgA->centroid[1];
    208:     /* West */
    209:     PetscCall(DMPlexPointLocalRead(dmFace, cellcone[3], fgeom, &fgA));
    210:     centroid_x[0] = fgA->centroid[0];
    211:     /* East */
    212:     PetscCall(DMPlexPointLocalRead(dmFace, cellcone[1], fgeom, &fgA));
    213:     centroid_x[1] = fgA->centroid[0];

    215:     /* Computing the cell widths in the x and y direction */
    216:     delta_x = centroid_x[1] - centroid_x[0];
    217:     delta_y = centroid_y[1] - centroid_y[0];

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

    222:     /* cellcone[0] - south */
    223:     PetscCall(DMPlexGetSupportSize(da, cellcone[0], &nC));
    224:     PetscCall(DMPlexGetSupport(da, cellcone[0], &cellsupport));
    225:     if (nC == 2) flux_south = (x[cellsupport[0]] * (-v_plus - user->diffusion * delta_x)) / delta_y;
    226:     else flux_south = (boundary * (-v_plus - user->diffusion * delta_x)) / delta_y;

    228:     /* cellcone[1] - east */
    229:     PetscCall(DMPlexGetSupportSize(da, cellcone[1], &nC));
    230:     PetscCall(DMPlexGetSupport(da, cellcone[1], &cellsupport));
    231:     if (nC == 2) flux_east = (x[cellsupport[1]] * (u_minus - user->diffusion * delta_y)) / delta_x;
    232:     else flux_east = (boundary * (u_minus - user->diffusion * delta_y)) / delta_x;

    234:     /* cellcone[2] - north */
    235:     PetscCall(DMPlexGetSupportSize(da, cellcone[2], &nC));
    236:     PetscCall(DMPlexGetSupport(da, cellcone[2], &cellsupport));
    237:     if (nC == 2) flux_north = (x[cellsupport[1]] * (v_minus - user->diffusion * delta_x)) / delta_y;
    238:     else flux_north = (boundary * (v_minus - user->diffusion * delta_x)) / delta_y;

    240:     /* cellcone[3] - west */
    241:     PetscCall(DMPlexGetSupportSize(da, cellcone[3], &nC));
    242:     PetscCall(DMPlexGetSupport(da, cellcone[3], &cellsupport));
    243:     if (nC == 2) flux_west = (x[cellsupport[0]] * (-u_plus - user->diffusion * delta_y)) / delta_x;
    244:     else flux_west = (boundary_left * (-u_plus - user->diffusion * delta_y)) / delta_x;

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

    249:     /* Calculating the net flux for each cell
    250:            and computing the RHS time derivative f[.] */
    251:     f[cell] = -(flux_centre + flux_east + flux_west + flux_north + flux_south);
    252:   }
    253:   PetscCall(PetscFVDestroy(&fvm));
    254:   PetscCall(VecRestoreArray(localX, &x));
    255:   PetscCall(VecRestoreArray(F, &f));
    256:   PetscCall(DMRestoreLocalVector(da, &localX));
    257:   PetscFunctionReturn(PETSC_SUCCESS);
    258: }

    260: int main(int argc, char **argv)
    261: {
    262:   TS                    ts; /* time integrator */
    263:   SNES                  snes;
    264:   Vec                   x, r; /* solution, residual vectors */
    265:   DM                    da;
    266:   PetscMPIInt           rank;
    267:   PetscViewerAndFormat *vf;
    268:   AppCtx                user; /* mesh context */
    269:   PetscInt              dim, numFields = 1, numBC, i;
    270:   PetscInt              numComp[1];
    271:   PetscInt              numDof[12];
    272:   PetscInt              bcField[1];
    273:   PetscSection          section;
    274:   IS                    bcPointIS[1];

    276:   /* Initialize program */
    277:   PetscFunctionBeginUser;
    278:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    279:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    280:   /* Create distributed array (DMPLEX) to manage parallel grid and vectors */
    281:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    282:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &da));
    283:   PetscCall(DMGetDimension(da, &dim));

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

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

    293:   numDof[0 * (dim + 1)]           = 1;
    294:   numDof[0 * (dim + 1) + dim - 1] = 1;
    295:   numDof[0 * (dim + 1) + dim]     = 1;

    297:   /* Setup boundary conditions */
    298:   numBC = 1;
    299:   /* Prescribe a Dirichlet condition on u on the boundary
    300:        Label "marker" is made by the mesh creation routine  */
    301:   bcField[0] = 0;
    302:   PetscCall(DMGetStratumIS(da, "marker", 1, &bcPointIS[0]));

    304:   /* Create a PetscSection with this data layout */
    305:   PetscCall(DMSetNumFields(da, numFields));
    306:   PetscCall(DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section));

    308:   /* Name the Field variables */
    309:   PetscCall(PetscSectionSetFieldName(section, 0, "u"));

    311:   /* Tell the DM to use this section (with the specified fields and dof) */
    312:   PetscCall(DMSetLocalSection(da, section));

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

    317:   /* Create a Vec with this layout and view it */
    318:   PetscCall(DMGetGlobalVector(da, &x));
    319:   PetscCall(VecDuplicate(x, &r));

    321:   /* Create timestepping solver context */
    322:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    323:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    324:   PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &user));

    326:   PetscCall(TSSetMaxTime(ts, 1.0));
    327:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    328:   PetscCall(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL));
    329:   PetscCall(TSSetDM(ts, da));

    331:   /* Customize nonlinear solver */
    332:   PetscCall(TSSetType(ts, TSEULER));
    333:   PetscCall(TSGetSNES(ts, &snes));
    334:   PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
    335:   PetscCall(SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));

    337:   /* Set initial conditions */
    338:   PetscCall(FormInitialSolution(da, x));
    339:   PetscCall(TSSetTimeStep(ts, .0001));
    340:   PetscCall(TSSetSolution(ts, x));
    341:   /* Set runtime options */
    342:   PetscCall(TSSetFromOptions(ts));
    343:   /* Solve nonlinear system */
    344:   PetscCall(TSSolve(ts, x));

    346:   /* Clean up routine */
    347:   PetscCall(DMRestoreGlobalVector(da, &x));
    348:   PetscCall(ISDestroy(&bcPointIS[0]));
    349:   PetscCall(PetscSectionDestroy(&section));
    350:   PetscCall(VecDestroy(&r));
    351:   PetscCall(TSDestroy(&ts));
    352:   PetscCall(DMDestroy(&da));
    353:   PetscCall(PetscFinalize());
    354:   return 0;
    355: }

    357: /*TEST

    359:     test:
    360:       suffix: 0
    361:       args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -dm_plex_boundary_label boundary -ts_max_steps 5 -ts_type rk
    362:       requires: !single !complex triangle ctetgen

    364: TEST*/