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