Actual source code: wp.c

      1: /*MC
      2:    MATMFFD_WP - Implements an approach for computing the differencing parameter
      3:    h used with the finite difference based matrix-free Jacobian. From Walker-Pernice {cite}`pw98`

      5:    $$
      6:    h = error_rel * sqrt(1 + ||u||) / ||a||
      7:    $$

      9:    Options Database Key:
     10: .   -mat_mffd_compute_normu -Compute the norm of u every time see `MatMFFDWPSetComputeNormU()`

     12:    Level: intermediate

     14:    Notes:
     15:    $ || U || $ does not change between linear iterations so is reused

     17:    In `KSPGMRES` $ || a || == 1 $ and so does not need to ever be computed except at restart
     18:     when it is recomputed.  Thus requires no global collectives when used with `KSPGMRES`

     20: .seealso: `MATMFFD`, `MATMFFD_DS`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_DS`
     21: M*/

     23: /*
     24:     This include file defines the data structure  MatMFFD that
     25:    includes information about the computation of h. It is shared by
     26:    all implementations that people provide.

     28:    See snesmfjdef.c for  a full set of comments on the routines below.
     29: */
     30: #include <petsc/private/matimpl.h>
     31: #include <../src/mat/impls/mffd/mffdimpl.h>

     33: typedef struct {
     34:   PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
     35:   PetscBool computenormU;
     36: } MatMFFD_WP;

     38: /*
     39:      MatMFFDCompute_WP - code for
     40:    computing h with matrix-free finite differences.

     42:   Input Parameters:
     43: +   ctx - the matrix-free context
     44: .   U - the location at which you want the Jacobian
     45: -   a - the direction you want the derivative

     47:   Output Parameter:
     48: .   h - the scale computed
     49: */
     50: static PetscErrorCode MatMFFDCompute_WP(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa)
     51: {
     52:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
     53:   PetscReal   normU, norma;

     55:   PetscFunctionBegin;
     56:   if (!(ctx->count % ctx->recomputeperiod)) {
     57:     if (hctx->computenormU || !ctx->ncurrenth) {
     58:       PetscCall(VecNorm(U, NORM_2, &normU));
     59:       hctx->normUfact = PetscSqrtReal(1.0 + normU);
     60:     }
     61:     PetscCall(VecNorm(a, NORM_2, &norma));
     62:     if (norma == 0.0) {
     63:       *zeroa = PETSC_TRUE;
     64:       PetscFunctionReturn(PETSC_SUCCESS);
     65:     }
     66:     *zeroa = PETSC_FALSE;
     67:     *h     = ctx->error_rel * hctx->normUfact / norma;
     68:   } else {
     69:     *h = ctx->currenth;
     70:   }
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: /*
     75:    MatMFFDView_WP - Prints information about this particular
     76:      method for computing h. Note that this does not print the general
     77:      information about the matrix-free, that is printed by the calling
     78:      routine.

     80:   Input Parameters:
     81: +   ctx - the matrix-free context
     82: -   viewer - the PETSc viewer

     84: */
     85: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx, PetscViewer viewer)
     86: {
     87:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
     88:   PetscBool   iascii;

     90:   PetscFunctionBegin;
     91:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
     92:   if (iascii) {
     93:     if (hctx->computenormU) {
     94:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Computes normU\n"));
     95:     } else {
     96:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Does not compute normU\n"));
     97:     }
     98:   }
     99:   PetscFunctionReturn(PETSC_SUCCESS);
    100: }

    102: /*
    103:    MatMFFDSetFromOptions_WP - Looks in the options database for
    104:      any options appropriate for this method

    106:   Input Parameter:
    107: .  ctx - the matrix-free context

    109: */
    110: static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
    111: {
    112:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;

    114:   PetscFunctionBegin;
    115:   PetscOptionsHeadBegin(PetscOptionsObject, "Walker-Pernice options");
    116:   PetscCall(PetscOptionsBool("-mat_mffd_compute_normu", "Compute the norm of u", "MatMFFDWPSetComputeNormU", hctx->computenormU, &hctx->computenormU, NULL));
    117:   PetscOptionsHeadEnd();
    118:   PetscFunctionReturn(PETSC_SUCCESS);
    119: }

    121: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
    122: {
    123:   PetscFunctionBegin;
    124:   PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL));
    125:   PetscCall(PetscFree(ctx->hctx));
    126:   PetscFunctionReturn(PETSC_SUCCESS);
    127: }

    129: static PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag)
    130: {
    131:   MatMFFD     ctx  = (MatMFFD)mat->data;
    132:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;

    134:   PetscFunctionBegin;
    135:   hctx->computenormU = flag;
    136:   PetscFunctionReturn(PETSC_SUCCESS);
    137: }

    139: /*@
    140:   MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the Walker-Pernice {cite}`pw98`
    141:   PETSc routine for computing h. With any Krylov solver this need only
    142:   be computed during the first iteration and kept for later.

    144:   Input Parameters:
    145: + A    - the `MATMFFD` matrix
    146: - flag - `PETSC_TRUE` causes it to compute $||U||$, `PETSC_FALSE` uses the previous value

    148:   Options Database Key:
    149: . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you
    150:               must be sure that $||U||$ has not changed in the mean time.

    152:   Level: advanced

    154:   Note:
    155:   See the manual page for `MATMFFD_WP` for a complete description of the
    156:   algorithm used to compute h.

    158: .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
    159: @*/
    160: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag)
    161: {
    162:   PetscFunctionBegin;
    164:   PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag));
    165:   PetscFunctionReturn(PETSC_SUCCESS);
    166: }

    168: /*
    169:      MatCreateMFFD_WP - Standard PETSc code for
    170:    computing h with matrix-free finite differences.

    172:    Input Parameter:
    173: .  ctx - the matrix-free context created by MatCreateMFFD()

    175: */
    176: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
    177: {
    178:   MatMFFD_WP *hctx;

    180:   PetscFunctionBegin;
    181:   /* allocate my own private data structure */
    182:   PetscCall(PetscNew(&hctx));
    183:   ctx->hctx          = (void *)hctx;
    184:   hctx->computenormU = PETSC_FALSE;

    186:   /* set the functions I am providing */
    187:   ctx->ops->compute        = MatMFFDCompute_WP;
    188:   ctx->ops->destroy        = MatMFFDDestroy_WP;
    189:   ctx->ops->view           = MatMFFDView_WP;
    190:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;

    192:   PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P));
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }