Actual source code: ex47.c

      1: static char help[] = "Pure advection with finite elements.\n\
      2: We solve the hyperbolic problem in a rectangular\n\
      3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";

      5: /*
      6: The continuity equation (https://en.wikipedia.org/wiki/Continuity_equation) for advection
      7: (https://en.wikipedia.org/wiki/Advection) of a conserved scalar quantity phi, with source q,

      9:   phi_t + div (phi u) = q

     11: if used with a solenoidal velocity field u (div u = 0) is given by

     13:   phi_t + u . grad phi = q

     15: For a vector quantity a, we likewise have

     17:   a_t + u . grad a = q
     18: */

     20: /*
     21:   r1: 8 SOR
     22:   r2: 1128 SOR
     23:   r3: > 10000 SOR

     25:   SOR is completely unreliable as a smoother, use Jacobi
     26:   r1: 8 MG
     27:   r2:
     28: */

     30: #include <petscdmplex.h>
     31: #include <petscts.h>
     32: #include <petscds.h>

     34: typedef enum {
     35:   PRIMITIVE,
     36:   INT_BY_PARTS
     37: } WeakFormType;

     39: typedef struct {
     40:   WeakFormType formType;
     41: } AppCtx;

     43: /* MMS1:

     45:   2D:
     46:   u   = <1, 1>
     47:   phi = x + y - 2t

     49:   phi_t + u . grad phi = -2 + <1, 1> . <1, 1> = 0

     51:   3D:
     52:   u   = <1, 1, 1>
     53:   phi = x + y + z - 3t

     55:   phi_t + u . grad phi = -3 + <1, 1, 1> . <1, 1, 1> = 0
     56: */

     58: static PetscErrorCode analytic_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
     59: {
     60:   PetscInt d;

     62:   *u = -dim * time;
     63:   for (d = 0; d < dim; ++d) *u += x[d];
     64:   return PETSC_SUCCESS;
     65: }

     67: static PetscErrorCode velocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
     68: {
     69:   PetscInt d;
     70:   for (d = 0; d < dim; ++d) u[d] = 1.0;
     71:   return PETSC_SUCCESS;
     72: }

     74: /* <psi, phi_t> + <psi, u . grad phi> */
     75: static void f0_prim_phi(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[])
     76: {
     77:   PetscInt d;

     79:   f0[0] = u_t[0];
     80:   for (d = 0; d < dim; ++d) f0[0] += a[d] * u_x[d];
     81: }

     83: /* <psi, phi_t> */
     84: static void f0_ibp_phi(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[])
     85: {
     86:   f0[0] = u_t[0];
     87: }

     89: /* <grad psi, u phi> */
     90: static void f1_ibp_phi(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[])
     91: {
     92:   PetscInt d;
     93:   for (d = 0; d < dim; ++d) f1[d] = a[d] * u[0];
     94: }

     96: /* <psi, phi_t> */
     97: static void g0_prim_phi(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[])
     98: {
     99:   g0[0] = u_tShift * 1.0;
    100: }

    102: /* <psi, u . grad phi> */
    103: static void g1_prim_phi(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 g1[])
    104: {
    105:   PetscInt d;
    106:   for (d = 0; d < dim; ++d) g1[d] = a[d];
    107: }

    109: /* <grad psi, u phi> */
    110: static void g2_ibp_phi(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 g2[])
    111: {
    112:   PetscInt d;
    113:   for (d = 0; d < dim; ++d) g2[d] = a[d];
    114: }

    116: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
    117: {
    118:   const char *formTypes[2] = {"primitive", "int_by_parts"};
    119:   PetscInt    form;

    121:   PetscFunctionBeginUser;
    122:   options->formType = PRIMITIVE;
    123:   PetscOptionsBegin(comm, "", "Advection Equation Options", "DMPLEX");
    124:   form = options->formType;
    125:   PetscCall(PetscOptionsEList("-form_type", "The weak form type", "ex47.c", formTypes, 2, formTypes[options->formType], &form, NULL));
    126:   options->formType = (WeakFormType)form;
    127:   PetscOptionsEnd();
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

    131: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
    132: {
    133:   PetscFunctionBeginUser;
    134:   PetscCall(DMCreate(comm, dm));
    135:   PetscCall(DMSetType(*dm, DMPLEX));
    136:   PetscCall(DMSetFromOptions(*dm));
    137:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
    142: {
    143:   PetscDS        ds;
    144:   DMLabel        label;
    145:   const PetscInt id = 1;

    147:   PetscFunctionBeginUser;
    148:   PetscCall(DMGetDS(dm, &ds));
    149:   switch (ctx->formType) {
    150:   case PRIMITIVE:
    151:     PetscCall(PetscDSSetResidual(ds, 0, f0_prim_phi, NULL));
    152:     PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_prim_phi, g1_prim_phi, NULL, NULL));
    153:     break;
    154:   case INT_BY_PARTS:
    155:     PetscCall(PetscDSSetResidual(ds, 0, f0_ibp_phi, f1_ibp_phi));
    156:     PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_prim_phi, NULL, g2_ibp_phi, NULL));
    157:     break;
    158:   }
    159:   PetscCall(PetscDSSetExactSolution(ds, 0, analytic_phi, ctx));
    160:   PetscCall(DMGetLabel(dm, "marker", &label));
    161:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)analytic_phi, NULL, ctx, NULL));
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }

    165: static PetscErrorCode SetupVelocity(DM dm, DM dmAux, AppCtx *user)
    166: {
    167:   PetscSimplePointFn *funcs[1] = {velocity};
    168:   Vec                 v;

    170:   PetscFunctionBeginUser;
    171:   PetscCall(DMCreateLocalVector(dmAux, &v));
    172:   PetscCall(DMProjectFunctionLocal(dmAux, 0.0, funcs, NULL, INSERT_ALL_VALUES, v));
    173:   PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, v));
    174:   PetscCall(VecDestroy(&v));
    175:   PetscFunctionReturn(PETSC_SUCCESS);
    176: }

    178: static PetscErrorCode SetupAuxDM(DM dm, PetscFE feAux, AppCtx *user)
    179: {
    180:   DM dmAux, coordDM;

    182:   PetscFunctionBeginUser;
    183:   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
    184:   PetscCall(DMGetCoordinateDM(dm, &coordDM));
    185:   if (!feAux) PetscFunctionReturn(PETSC_SUCCESS);
    186:   PetscCall(DMClone(dm, &dmAux));
    187:   PetscCall(DMSetCoordinateDM(dmAux, coordDM));
    188:   PetscCall(DMSetField(dmAux, 0, NULL, (PetscObject)feAux));
    189:   PetscCall(DMCreateDS(dmAux));
    190:   PetscCall(SetupVelocity(dm, dmAux, user));
    191:   PetscCall(DMDestroy(&dmAux));
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

    195: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
    196: {
    197:   DM        cdm = dm;
    198:   PetscFE   fe, feAux;
    199:   MPI_Comm  comm;
    200:   PetscInt  dim;
    201:   PetscBool simplex;

    203:   PetscFunctionBeginUser;
    204:   PetscCall(DMGetDimension(dm, &dim));
    205:   PetscCall(DMPlexIsSimplex(dm, &simplex));
    206:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    207:   PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "phi_", -1, &fe));
    208:   PetscCall(PetscObjectSetName((PetscObject)fe, "phi"));
    209:   PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "vel_", -1, &feAux));
    210:   PetscCall(PetscFECopyQuadrature(fe, feAux));
    211:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
    212:   PetscCall(DMCreateDS(dm));
    213:   PetscCall(SetupProblem(dm, ctx));
    214:   while (cdm) {
    215:     PetscCall(SetupAuxDM(cdm, feAux, ctx));
    216:     PetscCall(DMCopyDisc(dm, cdm));
    217:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    218:   }
    219:   PetscCall(PetscFEDestroy(&fe));
    220:   PetscCall(PetscFEDestroy(&feAux));
    221:   PetscFunctionReturn(PETSC_SUCCESS);
    222: }

    224: static PetscErrorCode MonitorError(KSP ksp, PetscInt it, PetscReal rnorm, PetscCtx ctx)
    225: {
    226:   DM                  dm;
    227:   PetscDS             ds;
    228:   PetscSimplePointFn *func[1];
    229:   void               *ctxs[1];
    230:   Vec                 u, r, error;
    231:   PetscReal           time = 0.5, res;

    233:   PetscFunctionBeginUser;
    234:   PetscCall(KSPGetDM(ksp, &dm));
    235:   PetscCall(DMSetOutputSequenceNumber(dm, it, time));
    236:   /* Calculate residual */
    237:   PetscCall(KSPBuildResidual(ksp, NULL, NULL, &r));
    238:   PetscCall(VecNorm(r, NORM_2, &res));
    239:   PetscCall(DMSetOutputSequenceNumber(dm, it, res));
    240:   PetscCall(PetscObjectSetName((PetscObject)r, "residual"));
    241:   PetscCall(VecViewFromOptions(r, NULL, "-res_vec_view"));
    242:   PetscCall(VecDestroy(&r));
    243:   /* Calculate error */
    244:   PetscCall(DMGetDS(dm, &ds));
    245:   PetscCall(PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]));
    246:   PetscCall(KSPBuildSolution(ksp, NULL, &u));
    247:   PetscCall(DMGetGlobalVector(dm, &error));
    248:   PetscCall(DMProjectFunction(dm, time, func, ctxs, INSERT_ALL_VALUES, error));
    249:   PetscCall(VecAXPY(error, -1.0, u));
    250:   PetscCall(PetscObjectSetName((PetscObject)error, "error"));
    251:   PetscCall(VecViewFromOptions(error, NULL, "-err_vec_view"));
    252:   PetscCall(DMRestoreGlobalVector(dm, &error));
    253:   PetscFunctionReturn(PETSC_SUCCESS);
    254: }

    256: static PetscErrorCode MyTSMonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, PetscCtx ctx)
    257: {
    258:   DM                  dm;
    259:   PetscDS             ds;
    260:   PetscSimplePointFn *func[1];
    261:   void               *ctxs[1];
    262:   PetscReal           error;

    264:   PetscFunctionBeginUser;
    265:   PetscCall(TSGetDM(ts, &dm));
    266:   PetscCall(DMGetDS(dm, &ds));
    267:   PetscCall(PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]));
    268:   PetscCall(DMComputeL2Diff(dm, crtime, func, ctxs, u, &error));
    269:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: %2.5g\n", (int)step, (double)crtime, (double)error));
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: int main(int argc, char **argv)
    274: {
    275:   AppCtx    ctx;
    276:   DM        dm;
    277:   TS        ts;
    278:   Vec       u, r;
    279:   PetscReal t = 0.0;

    281:   PetscFunctionBeginUser;
    282:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    283:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
    284:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &ctx));
    285:   PetscCall(DMSetApplicationContext(dm, &ctx));
    286:   PetscCall(SetupDiscretization(dm, &ctx));

    288:   PetscCall(DMCreateGlobalVector(dm, &u));
    289:   PetscCall(PetscObjectSetName((PetscObject)u, "phi"));
    290:   PetscCall(VecDuplicate(u, &r));

    292:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    293:   PetscCall(TSMonitorSet(ts, MyTSMonitorError, &ctx, NULL));
    294:   PetscCall(TSSetDM(ts, dm));
    295:   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
    296:   PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
    297:   PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
    298:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    299:   PetscCall(TSSetFromOptions(ts));

    301:   {
    302:     PetscDS             ds;
    303:     PetscSimplePointFn *func[1];
    304:     void               *ctxs[1];

    306:     PetscCall(DMGetDS(dm, &ds));
    307:     PetscCall(PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]));
    308:     PetscCall(DMProjectFunction(dm, t, func, ctxs, INSERT_ALL_VALUES, u));
    309:   }
    310:   {
    311:     SNES snes;
    312:     KSP  ksp;

    314:     PetscCall(TSGetSNES(ts, &snes));
    315:     PetscCall(SNESGetKSP(snes, &ksp));
    316:     PetscCall(KSPMonitorSet(ksp, MonitorError, &ctx, NULL));
    317:   }
    318:   PetscCall(TSSolve(ts, u));
    319:   PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));

    321:   PetscCall(VecDestroy(&u));
    322:   PetscCall(VecDestroy(&r));
    323:   PetscCall(TSDestroy(&ts));
    324:   PetscCall(DMDestroy(&dm));
    325:   PetscCall(PetscFinalize());
    326:   return 0;
    327: }

    329: /*TEST

    331:   # Full solves
    332:   test:
    333:     suffix: 2d_p1p1_r1
    334:     requires: triangle
    335:     args: -dm_refine 1 -phi_petscspace_degree 1 -vel_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_time_step 0.1 -pc_type lu -snes_monitor_short -snes_converged_reason -ts_monitor

    337:   test:
    338:     suffix: 2d_p1p1_sor_r1
    339:     requires: triangle !single
    340:     args: -dm_refine 1 -phi_petscspace_degree 1 -vel_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_time_step 0.1 -ksp_rtol 1.0e-9 -pc_type sor -snes_monitor_short -snes_converged_reason -ksp_monitor_short -ts_monitor

    342:   test:
    343:     suffix: 2d_p1p1_mg_r1
    344:     requires: triangle !single
    345:     args: -dm_refine_hierarchy 1 -phi_petscspace_degree 1 -vel_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_time_step 0.1 -ksp_type fgmres -ksp_rtol 1.0e-9 -pc_type mg -pc_mg_levels 2 -snes_monitor_short -snes_converged_reason -snes_view -ksp_monitor_true_residual -ts_monitor

    347: TEST*/