Actual source code: dmplexts.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/tsimpl.h>
      3: #include <petsc/private/snesimpl.h>
      4: #include <petscds.h>
      5: #include <petscfv.h>

      7: static PetscErrorCode DMTSConvertPlex(DM dm, DM *plex, PetscBool copy)
      8: {
      9:   PetscBool      isPlex;

     13:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
     14:   if (isPlex) {
     15:     *plex = dm;
     16:     PetscObjectReference((PetscObject) dm);
     17:   } else {
     18:     PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
     19:     if (!*plex) {
     20:       DMConvert(dm,DMPLEX,plex);
     21:       PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
     22:       if (copy) {
     23:         PetscInt    i;
     24:         PetscObject obj;
     25:         const char *comps[3] = {"A","dmAux","dmCh"};

     27:         DMCopyDMTS(dm, *plex);
     28:         DMCopyDMSNES(dm, *plex);
     29:         for (i = 0; i < 3; i++) {
     30:           PetscObjectQuery((PetscObject) dm, comps[i], &obj);
     31:           PetscObjectCompose((PetscObject) *plex, comps[i], obj);
     32:         }
     33:       }
     34:     } else {
     35:       PetscObjectReference((PetscObject) *plex);
     36:     }
     37:   }
     38:   return(0);
     39: }

     41: /*@
     42:   DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user

     44:   Input Parameters:
     45: + dm - The mesh
     46: . t - The time
     47: . locX  - Local solution
     48: - user - The user context

     50:   Output Parameter:
     51: . F  - Global output vector

     53:   Level: developer

     55: .seealso: DMPlexComputeJacobianActionFEM()
     56: @*/
     57: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
     58: {
     59:   Vec            locF;
     60:   IS             cellIS;
     61:   DM             plex;
     62:   PetscInt       depth;
     63:   PetscHashFormKey key = {NULL, 0, 0};

     67:   DMTSConvertPlex(dm,&plex,PETSC_TRUE);
     68:   DMPlexGetDepth(plex, &depth);
     69:   DMGetStratumIS(plex, "dim", depth, &cellIS);
     70:   if (!cellIS) {
     71:     DMGetStratumIS(plex, "depth", depth, &cellIS);
     72:   }
     73:   DMGetLocalVector(plex, &locF);
     74:   VecZeroEntries(locF);
     75:   DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locF, user);
     76:   DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
     77:   DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
     78:   DMRestoreLocalVector(plex, &locF);
     79:   ISDestroy(&cellIS);
     80:   DMDestroy(&plex);
     81:   return(0);
     82: }

     84: /*@
     85:   DMPlexTSComputeBoundary - Insert the essential boundary values for the local input X and/or its time derivative X_t using pointwise functions specified by the user

     87:   Input Parameters:
     88: + dm - The mesh
     89: . t - The time
     90: . locX  - Local solution
     91: . locX_t - Local solution time derivative, or NULL
     92: - user - The user context

     94:   Level: developer

     96: .seealso: DMPlexComputeJacobianActionFEM()
     97: @*/
     98: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
     99: {
    100:   DM             plex;
    101:   Vec            faceGeometryFVM = NULL;
    102:   PetscInt       Nf, f;

    106:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
    107:   DMGetNumFields(plex, &Nf);
    108:   if (!locX_t) {
    109:     /* This is the RHS part */
    110:     for (f = 0; f < Nf; f++) {
    111:       PetscObject  obj;
    112:       PetscClassId id;

    114:       DMGetField(plex, f, NULL, &obj);
    115:       PetscObjectGetClassId(obj, &id);
    116:       if (id == PETSCFV_CLASSID) {
    117:         DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
    118:         break;
    119:       }
    120:     }
    121:   }
    122:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
    123:   DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);
    124:   DMDestroy(&plex);
    125:   return(0);
    126: }

    128: /*@
    129:   DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user

    131:   Input Parameters:
    132: + dm - The mesh
    133: . t - The time
    134: . locX  - Local solution
    135: . locX_t - Local solution time derivative, or NULL
    136: - user - The user context

    138:   Output Parameter:
    139: . locF  - Local output vector

    141:   Level: developer

    143: .seealso: DMPlexComputeJacobianActionFEM()
    144: @*/
    145: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
    146: {
    147:   DM             plex;
    148:   IS             allcellIS;
    149:   PetscInt       Nds, s;

    153:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
    154:   DMPlexGetAllCells_Internal(plex, &allcellIS);
    155:   DMGetNumDS(dm, &Nds);
    156:   for (s = 0; s < Nds; ++s) {
    157:     PetscDS          ds;
    158:     IS               cellIS;
    159:     PetscHashFormKey key;

    161:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
    162:     key.value = 0;
    163:     key.field = 0;
    164:     if (!key.label) {
    165:       PetscObjectReference((PetscObject) allcellIS);
    166:       cellIS = allcellIS;
    167:     } else {
    168:       IS pointIS;

    170:       key.value = 1;
    171:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
    172:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
    173:       ISDestroy(&pointIS);
    174:     }
    175:     DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user);
    176:     ISDestroy(&cellIS);
    177:   }
    178:   ISDestroy(&allcellIS);
    179:   DMDestroy(&plex);
    180:   return(0);
    181: }

    183: /*@
    184:   DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user

    186:   Input Parameters:
    187: + dm - The mesh
    188: . t - The time
    189: . locX  - Local solution
    190: . locX_t - Local solution time derivative, or NULL
    191: . X_tshift - The multiplicative parameter for dF/du_t
    192: - user - The user context

    194:   Output Parameter:
    195: . locF  - Local output vector

    197:   Level: developer

    199: .seealso: DMPlexComputeJacobianActionFEM()
    200: @*/
    201: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
    202: {
    203:   DM             plex;
    204:   IS             allcellIS;
    205:   PetscBool      hasJac, hasPrec;
    206:   PetscInt       Nds, s;

    210:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
    211:   DMPlexGetAllCells_Internal(plex, &allcellIS);
    212:   DMGetNumDS(dm, &Nds);
    213:   for (s = 0; s < Nds; ++s) {
    214:     PetscDS          ds;
    215:     IS               cellIS;
    216:     PetscHashFormKey key;

    218:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
    219:     key.value = 0;
    220:     key.field = 0;
    221:     if (!key.label) {
    222:       PetscObjectReference((PetscObject) allcellIS);
    223:       cellIS = allcellIS;
    224:     } else {
    225:       IS pointIS;

    227:       key.value = 1;
    228:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
    229:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
    230:       ISDestroy(&pointIS);
    231:     }
    232:     if (!s) {
    233:       PetscDSHasJacobian(ds, &hasJac);
    234:       PetscDSHasJacobianPreconditioner(ds, &hasPrec);
    235:       if (hasJac && hasPrec) {MatZeroEntries(Jac);}
    236:       MatZeroEntries(JacP);
    237:     }
    238:     DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
    239:     ISDestroy(&cellIS);
    240:   }
    241:   ISDestroy(&allcellIS);
    242:   DMDestroy(&plex);
    243:   return(0);
    244: }

    246: /*@C
    247:   DMTSCheckResidual - Check the residual of the exact solution

    249:   Input Parameters:
    250: + ts  - the TS object
    251: . dm  - the DM
    252: . t   - the time
    253: . u   - a DM vector
    254: . u_t - a DM vector
    255: - tol - A tolerance for the check, or -1 to print the results instead

    257:   Output Parameters:
    258: . residual - The residual norm of the exact solution, or NULL

    260:   Level: developer

    262: .seealso: DNTSCheckFromOptions(), DMTSCheckJacobian(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
    263: @*/
    264: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
    265: {
    266:   MPI_Comm       comm;
    267:   Vec            r;
    268:   PetscReal      res;

    276:   PetscObjectGetComm((PetscObject) ts, &comm);
    277:   DMComputeExactSolution(dm, t, u, u_t);
    278:   VecDuplicate(u, &r);
    279:   TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
    280:   VecNorm(r, NORM_2, &res);
    281:   if (tol >= 0.0) {
    282:     if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
    283:   } else if (residual) {
    284:     *residual = res;
    285:   } else {
    286:     PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
    287:     VecChop(r, 1.0e-10);
    288:     PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) dm);
    289:     PetscObjectSetName((PetscObject) r, "Initial Residual");
    290:     PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
    291:     VecViewFromOptions(r, NULL, "-vec_view");
    292:     PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
    293:   }
    294:   VecDestroy(&r);
    295:   return(0);
    296: }

    298: /*@C
    299:   DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test

    301:   Input Parameters:
    302: + ts  - the TS object
    303: . dm  - the DM
    304: . t   - the time
    305: . u   - a DM vector
    306: . u_t - a DM vector
    307: - tol - A tolerance for the check, or -1 to print the results instead

    309:   Output Parameters:
    310: + isLinear - Flag indicaing that the function looks linear, or NULL
    311: - convRate - The rate of convergence of the linear model, or NULL

    313:   Level: developer

    315: .seealso: DNTSCheckFromOptions(), DMTSCheckResidual(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
    316: @*/
    317: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
    318: {
    319:   MPI_Comm       comm;
    320:   PetscDS        ds;
    321:   Mat            J, M;
    322:   MatNullSpace   nullspace;
    323:   PetscReal      dt, shift, slope, intercept;
    324:   PetscBool      hasJac, hasPrec, isLin = PETSC_FALSE;

    333:   PetscObjectGetComm((PetscObject) ts, &comm);
    334:   DMComputeExactSolution(dm, t, u, u_t);
    335:   /* Create and view matrices */
    336:   TSGetTimeStep(ts, &dt);
    337:   shift = 1.0/dt;
    338:   DMCreateMatrix(dm, &J);
    339:   DMGetDS(dm, &ds);
    340:   PetscDSHasJacobian(ds, &hasJac);
    341:   PetscDSHasJacobianPreconditioner(ds, &hasPrec);
    342:   if (hasJac && hasPrec) {
    343:     DMCreateMatrix(dm, &M);
    344:     TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);
    345:     PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
    346:     PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
    347:     MatViewFromOptions(M, NULL, "-mat_view");
    348:     MatDestroy(&M);
    349:   } else {
    350:     TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);
    351:   }
    352:   PetscObjectSetName((PetscObject) J, "Jacobian");
    353:   PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
    354:   MatViewFromOptions(J, NULL, "-mat_view");
    355:   /* Check nullspace */
    356:   MatGetNullSpace(J, &nullspace);
    357:   if (nullspace) {
    358:     PetscBool isNull;
    359:     MatNullSpaceTest(nullspace, J, &isNull);
    360:     if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
    361:   }
    362:   /* Taylor test */
    363:   {
    364:     PetscRandom rand;
    365:     Vec         du, uhat, uhat_t, r, rhat, df;
    366:     PetscReal   h;
    367:     PetscReal  *es, *hs, *errors;
    368:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
    369:     PetscInt    Nv, v;

    371:     /* Choose a perturbation direction */
    372:     PetscRandomCreate(comm, &rand);
    373:     VecDuplicate(u, &du);
    374:     VecSetRandom(du, rand);
    375:     PetscRandomDestroy(&rand);
    376:     VecDuplicate(u, &df);
    377:     MatMult(J, du, df);
    378:     /* Evaluate residual at u, F(u), save in vector r */
    379:     VecDuplicate(u, &r);
    380:     TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
    381:     /* Look at the convergence of our Taylor approximation as we approach u */
    382:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
    383:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
    384:     VecDuplicate(u, &uhat);
    385:     VecDuplicate(u, &uhat_t);
    386:     VecDuplicate(u, &rhat);
    387:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
    388:       VecWAXPY(uhat, h, du, u);
    389:       VecWAXPY(uhat_t, h*shift, du, u_t);
    390:       /* F(\hat u, \hat u_t) \approx F(u, u_t) + J(u, u_t) (uhat - u) + J_t(u, u_t) (uhat_t - u_t) = F(u) + h * J(u) du + h * shift * J_t(u) du = F(u) + h F' du */
    391:       TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
    392:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
    393:       VecNorm(rhat, NORM_2, &errors[Nv]);

    395:       es[Nv] = PetscLog10Real(errors[Nv]);
    396:       hs[Nv] = PetscLog10Real(h);
    397:     }
    398:     VecDestroy(&uhat);
    399:     VecDestroy(&uhat_t);
    400:     VecDestroy(&rhat);
    401:     VecDestroy(&df);
    402:     VecDestroy(&r);
    403:     VecDestroy(&du);
    404:     for (v = 0; v < Nv; ++v) {
    405:       if ((tol >= 0) && (errors[v] > tol)) break;
    406:       else if (errors[v] > PETSC_SMALL)    break;
    407:     }
    408:     if (v == Nv) isLin = PETSC_TRUE;
    409:     PetscLinearRegression(Nv, hs, es, &slope, &intercept);
    410:     PetscFree3(es, hs, errors);
    411:     /* Slope should be about 2 */
    412:     if (tol >= 0) {
    413:       if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
    414:     } else if (isLinear || convRate) {
    415:       if (isLinear) *isLinear = isLin;
    416:       if (convRate) *convRate = slope;
    417:     } else {
    418:       if (!isLin) {PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);}
    419:       else        {PetscPrintf(comm, "Function appears to be linear\n");}
    420:     }
    421:   }
    422:   MatDestroy(&J);
    423:   return(0);
    424: }

    426: /*@C
    427:   DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information

    429:   Input Parameters:
    430: + ts - the TS object
    431: - u  - representative TS vector

    433:   Note: The user must call PetscDSSetExactSolution() beforehand

    435:   Level: developer
    436: @*/
    437: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
    438: {
    439:   DM             dm;
    440:   SNES           snes;
    441:   Vec            sol, u_t;
    442:   PetscReal      t;
    443:   PetscBool      check;

    447:   PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
    448:   if (!check) return(0);
    449:   VecDuplicate(u, &sol);
    450:   VecCopy(u, sol);
    451:   TSSetSolution(ts, u);
    452:   TSGetDM(ts, &dm);
    453:   TSSetUp(ts);
    454:   TSGetSNES(ts, &snes);
    455:   SNESSetSolution(snes, u);

    457:   TSGetTime(ts, &t);
    458:   DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);
    459:   DMGetGlobalVector(dm, &u_t);
    460:   DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);
    461:   DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);
    462:   DMRestoreGlobalVector(dm, &u_t);

    464:   VecDestroy(&sol);
    465:   return(0);
    466: }