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 { 60 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 61 PetscReal normU, norma; 62 63 PetscFunctionBegin; 64 if (!(ctx->count % ctx->recomputeperiod)) { 65 if (hctx->computenormU || !ctx->ncurrenth) { 66 PetscCall(VecNorm(U, NORM_2, &normU)); 67 hctx->normUfact = PetscSqrtReal(1.0 + normU); 68 } 69 PetscCall(VecNorm(a, NORM_2, &norma)); 70 if (norma == 0.0) { 71 *zeroa = PETSC_TRUE; 72 PetscFunctionReturn(PETSC_SUCCESS); 73 } 74 *zeroa = PETSC_FALSE; 75 *h = ctx->error_rel * hctx->normUfact / norma; 76 } else { 77 *h = ctx->currenth; 78 } 79 PetscFunctionReturn(PETSC_SUCCESS); 80 } 81 82 /* 83 MatMFFDView_WP - Prints information about this particular 84 method for computing h. Note that this does not print the general 85 information about the matrix free, that is printed by the calling 86 routine. 87 88 Input Parameters: 89 + ctx - the matrix free context 90 - viewer - the PETSc viewer 91 92 */ 93 static PetscErrorCode MatMFFDView_WP(MatMFFD ctx, PetscViewer viewer) 94 { 95 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 96 PetscBool iascii; 97 98 PetscFunctionBegin; 99 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 100 if (iascii) { 101 if (hctx->computenormU) { 102 PetscCall(PetscViewerASCIIPrintf(viewer, " Computes normU\n")); 103 } else { 104 PetscCall(PetscViewerASCIIPrintf(viewer, " Does not compute normU\n")); 105 } 106 } 107 PetscFunctionReturn(PETSC_SUCCESS); 108 } 109 110 /* 111 MatMFFDSetFromOptions_WP - Looks in the options database for 112 any options appropriate for this method 113 114 Input Parameter: 115 . ctx - the matrix free context 116 117 */ 118 static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx, PetscOptionItems *PetscOptionsObject) 119 { 120 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 121 122 PetscFunctionBegin; 123 PetscOptionsHeadBegin(PetscOptionsObject, "Walker-Pernice options"); 124 PetscCall(PetscOptionsBool("-mat_mffd_compute_normu", "Compute the norm of u", "MatMFFDWPSetComputeNormU", hctx->computenormU, &hctx->computenormU, NULL)); 125 PetscOptionsHeadEnd(); 126 PetscFunctionReturn(PETSC_SUCCESS); 127 } 128 129 /* 130 MatMFFDDestroy_WP - Frees the space allocated by 131 MatCreateMFFD_WP(). 132 133 Input Parameter: 134 . ctx - the matrix free context 135 136 Notes: 137 does not free the ctx, that is handled by the calling routine 138 139 */ 140 static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx) 141 { 142 PetscFunctionBegin; 143 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL)); 144 PetscCall(PetscFree(ctx->hctx)); 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 148 PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag) 149 { 150 MatMFFD ctx = (MatMFFD)mat->data; 151 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 152 153 PetscFunctionBegin; 154 hctx->computenormU = flag; 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 /*@ 159 MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the Walker-Pernice 160 PETSc routine for computing h. With any Krylov solver this need only 161 be computed during the first iteration and kept for later. 162 163 Input Parameters: 164 + A - the `MATMFFD` matrix 165 - flag - `PETSC_TRUE` causes it to compute ||U||, `PETSC_FALSE` uses the previous value 166 167 Options Database Key: 168 . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you 169 must be sure that ||U|| has not changed in the mean time. 170 171 Level: advanced 172 173 Note: 174 See the manual page for `MATMFFD_WP` for a complete description of the 175 algorithm used to compute h. 176 177 .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()` 178 @*/ 179 PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag) 180 { 181 PetscFunctionBegin; 182 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 183 PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 /* 188 MatCreateMFFD_WP - Standard PETSc code for 189 computing h with matrix-free finite differences. 190 191 Input Parameter: 192 . ctx - the matrix free context created by MatCreateMFFD() 193 194 */ 195 PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx) 196 { 197 MatMFFD_WP *hctx; 198 199 PetscFunctionBegin; 200 /* allocate my own private data structure */ 201 PetscCall(PetscNew(&hctx)); 202 ctx->hctx = (void *)hctx; 203 hctx->computenormU = PETSC_FALSE; 204 205 /* set the functions I am providing */ 206 ctx->ops->compute = MatMFFDCompute_WP; 207 ctx->ops->destroy = MatMFFDDestroy_WP; 208 ctx->ops->view = MatMFFDView_WP; 209 ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP; 210 211 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P)); 212 PetscFunctionReturn(PETSC_SUCCESS); 213 } 214