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