1 /*MC 2 MATMFFD_WP - Implements an approach for computing the differencing parameter 3 h used with the finite difference based matrix-free Jacobian. 4 5 h = error_rel * sqrt(1 + ||U||) / ||a|| 6 7 Options Database Key: 8 . -mat_mffd_compute_normu -Compute the norm of u every time see `MatMFFDWPSetComputeNormU()` 9 10 Level: intermediate 11 12 Notes: 13 || U || does not change between linear iterations so is reused 14 15 In `KSPGMRES` || a || == 1 and so does not need to ever be computed except at restart 16 when it is recomputed. Thus requires no global collectives when used with `KSPGMRES` 17 18 Formula used: 19 F'(u)*a = [F(u+h*a) - F(u)]/h where 20 21 Reference: 22 . * - M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative 23 Solver for Nonlinear Systems", SIAM J. Sci. Stat. Comput.", 1998, 24 vol 19, pp. 302--318. 25 26 .seealso: `MATMFFD`, `MATMFFD_DS`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_DS` 27 M*/ 28 29 /* 30 This include file defines the data structure MatMFFD that 31 includes information about the computation of h. It is shared by 32 all implementations that people provide. 33 34 See snesmfjdef.c for a full set of comments on the routines below. 35 */ 36 #include <petsc/private/matimpl.h> 37 #include <../src/mat/impls/mffd/mffdimpl.h> /*I "petscmat.h" I*/ 38 39 typedef struct { 40 PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */ 41 PetscBool computenormU; 42 } MatMFFD_WP; 43 44 /* 45 MatMFFDCompute_WP - code for 46 computing h with matrix-free finite differences. 47 48 Input Parameters: 49 + ctx - the matrix-free context 50 . U - the location at which you want the Jacobian 51 - a - the direction you want the derivative 52 53 Output Parameter: 54 . h - the scale computed 55 56 */ 57 static PetscErrorCode MatMFFDCompute_WP(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa) 58 { 59 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 60 PetscReal normU, norma; 61 62 PetscFunctionBegin; 63 if (!(ctx->count % ctx->recomputeperiod)) { 64 if (hctx->computenormU || !ctx->ncurrenth) { 65 PetscCall(VecNorm(U, NORM_2, &normU)); 66 hctx->normUfact = PetscSqrtReal(1.0 + normU); 67 } 68 PetscCall(VecNorm(a, NORM_2, &norma)); 69 if (norma == 0.0) { 70 *zeroa = PETSC_TRUE; 71 PetscFunctionReturn(PETSC_SUCCESS); 72 } 73 *zeroa = PETSC_FALSE; 74 *h = ctx->error_rel * hctx->normUfact / norma; 75 } else { 76 *h = ctx->currenth; 77 } 78 PetscFunctionReturn(PETSC_SUCCESS); 79 } 80 81 /* 82 MatMFFDView_WP - Prints information about this particular 83 method for computing h. Note that this does not print the general 84 information about the matrix-free, that is printed by the calling 85 routine. 86 87 Input Parameters: 88 + ctx - the matrix-free context 89 - viewer - the PETSc viewer 90 91 */ 92 static PetscErrorCode MatMFFDView_WP(MatMFFD ctx, PetscViewer viewer) 93 { 94 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 95 PetscBool iascii; 96 97 PetscFunctionBegin; 98 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 99 if (iascii) { 100 if (hctx->computenormU) { 101 PetscCall(PetscViewerASCIIPrintf(viewer, " Computes normU\n")); 102 } else { 103 PetscCall(PetscViewerASCIIPrintf(viewer, " Does not compute normU\n")); 104 } 105 } 106 PetscFunctionReturn(PETSC_SUCCESS); 107 } 108 109 /* 110 MatMFFDSetFromOptions_WP - Looks in the options database for 111 any options appropriate for this method 112 113 Input Parameter: 114 . ctx - the matrix-free context 115 116 */ 117 static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx, PetscOptionItems *PetscOptionsObject) 118 { 119 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 120 121 PetscFunctionBegin; 122 PetscOptionsHeadBegin(PetscOptionsObject, "Walker-Pernice options"); 123 PetscCall(PetscOptionsBool("-mat_mffd_compute_normu", "Compute the norm of u", "MatMFFDWPSetComputeNormU", hctx->computenormU, &hctx->computenormU, NULL)); 124 PetscOptionsHeadEnd(); 125 PetscFunctionReturn(PETSC_SUCCESS); 126 } 127 128 static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx) 129 { 130 PetscFunctionBegin; 131 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL)); 132 PetscCall(PetscFree(ctx->hctx)); 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 136 static PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag) 137 { 138 MatMFFD ctx = (MatMFFD)mat->data; 139 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 140 141 PetscFunctionBegin; 142 hctx->computenormU = flag; 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /*@ 147 MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the Walker-Pernice 148 PETSc routine for computing h. With any Krylov solver this need only 149 be computed during the first iteration and kept for later. 150 151 Input Parameters: 152 + A - the `MATMFFD` matrix 153 - flag - `PETSC_TRUE` causes it to compute ||U||, `PETSC_FALSE` uses the previous value 154 155 Options Database Key: 156 . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you 157 must be sure that ||U|| has not changed in the mean time. 158 159 Level: advanced 160 161 Note: 162 See the manual page for `MATMFFD_WP` for a complete description of the 163 algorithm used to compute h. 164 165 .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()` 166 @*/ 167 PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag) 168 { 169 PetscFunctionBegin; 170 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 171 PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag)); 172 PetscFunctionReturn(PETSC_SUCCESS); 173 } 174 175 /* 176 MatCreateMFFD_WP - Standard PETSc code for 177 computing h with matrix-free finite differences. 178 179 Input Parameter: 180 . ctx - the matrix-free context created by MatCreateMFFD() 181 182 */ 183 PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx) 184 { 185 MatMFFD_WP *hctx; 186 187 PetscFunctionBegin; 188 /* allocate my own private data structure */ 189 PetscCall(PetscNew(&hctx)); 190 ctx->hctx = (void *)hctx; 191 hctx->computenormU = PETSC_FALSE; 192 193 /* set the functions I am providing */ 194 ctx->ops->compute = MatMFFDCompute_WP; 195 ctx->ops->destroy = MatMFFDDestroy_WP; 196 ctx->ops->view = MatMFFDView_WP; 197 ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP; 198 199 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P)); 200 PetscFunctionReturn(PETSC_SUCCESS); 201 } 202