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 static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx) 130 { 131 PetscFunctionBegin; 132 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL)); 133 PetscCall(PetscFree(ctx->hctx)); 134 PetscFunctionReturn(PETSC_SUCCESS); 135 } 136 137 PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag) 138 { 139 MatMFFD ctx = (MatMFFD)mat->data; 140 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx; 141 142 PetscFunctionBegin; 143 hctx->computenormU = flag; 144 PetscFunctionReturn(PETSC_SUCCESS); 145 } 146 147 /*@ 148 MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the Walker-Pernice 149 PETSc routine for computing h. With any Krylov solver this need only 150 be computed during the first iteration and kept for later. 151 152 Input Parameters: 153 + A - the `MATMFFD` matrix 154 - flag - `PETSC_TRUE` causes it to compute ||U||, `PETSC_FALSE` uses the previous value 155 156 Options Database Key: 157 . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you 158 must be sure that ||U|| has not changed in the mean time. 159 160 Level: advanced 161 162 Note: 163 See the manual page for `MATMFFD_WP` for a complete description of the 164 algorithm used to compute h. 165 166 .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()` 167 @*/ 168 PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag) 169 { 170 PetscFunctionBegin; 171 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 172 PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag)); 173 PetscFunctionReturn(PETSC_SUCCESS); 174 } 175 176 /* 177 MatCreateMFFD_WP - Standard PETSc code for 178 computing h with matrix-free finite differences. 179 180 Input Parameter: 181 . ctx - the matrix-free context created by MatCreateMFFD() 182 183 */ 184 PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx) 185 { 186 MatMFFD_WP *hctx; 187 188 PetscFunctionBegin; 189 /* allocate my own private data structure */ 190 PetscCall(PetscNew(&hctx)); 191 ctx->hctx = (void *)hctx; 192 hctx->computenormU = PETSC_FALSE; 193 194 /* set the functions I am providing */ 195 ctx->ops->compute = MatMFFDCompute_WP; 196 ctx->ops->destroy = MatMFFDDestroy_WP; 197 ctx->ops->view = MatMFFDView_WP; 198 ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP; 199 200 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P)); 201 PetscFunctionReturn(PETSC_SUCCESS); 202 } 203