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