1 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdm.h> /*I "petscdm.h" I*/ 4 #include <../src/mat/impls/mffd/mffdimpl.h> 5 #include <petsc/private/matimpl.h> 6 7 /*@C 8 MatMFFDComputeJacobian - Tells the matrix-free Jacobian object the new location at which 9 Jacobian matrix vector products will be computed at, i.e. J(x) * a. The x is obtained 10 from the SNES object (using SNESGetSolution()). 11 12 Logically Collective on SNES 13 14 Input Parameters: 15 + snes - the nonlinear solver context 16 . x - the point at which the Jacobian vector products will be performed 17 . jac - the matrix-free Jacobian object 18 . B - either the same as jac or another matrix type (ignored) 19 . flag - not relevent for matrix-free form 20 - dummy - the user context (ignored) 21 22 Level: developer 23 24 Warning: 25 If MatMFFDSetBase() is ever called on jac then this routine will NO longer get 26 the x from the SNES object and MatMFFDSetBase() must from that point on be used to 27 change the base vector x. 28 29 Notes: 30 This can be passed into SNESSetJacobian() as the Jacobian evaluation function argument 31 when using a completely matrix-free solver, 32 that is the B matrix is also the same matrix operator. This is used when you select 33 -snes_mf but rarely used directly by users. (All this routine does is call MatAssemblyBegin/End() on 34 the Mat jac.) 35 36 .seealso: MatMFFDGetH(), MatCreateSNESMF(), MatCreateMFFD(), MATMFFD, 37 MatMFFDSetHHistory(), MatMFFDSetFunctionError(), MatCreateMFFD(), SNESSetJacobian() 38 39 @*/ 40 PetscErrorCode MatMFFDComputeJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 41 { 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 46 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 47 PetscFunctionReturn(0); 48 } 49 50 PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat,MatAssemblyType); 51 PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat,Vec,Vec); 52 53 /* 54 MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the 55 base from the SNES context 56 57 */ 58 static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt) 59 { 60 PetscErrorCode ierr; 61 MatMFFD j = (MatMFFD)J->data; 62 SNES snes = (SNES)j->ctx; 63 Vec u,f; 64 65 PetscFunctionBegin; 66 ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr); 67 68 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 69 if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) { 70 ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr); 71 ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr); 72 } else { 73 /* f value known by SNES is not correct for other differencing function */ 74 ierr = MatMFFDSetBase_MFFD(J,u,NULL);CHKERRQ(ierr); 75 } 76 PetscFunctionReturn(0); 77 } 78 79 /* 80 MatAssemblyEnd_SNESMF_UseBase - Calls MatAssemblyEnd_MFFD() and then sets the 81 base from the SNES context. This version will cause the base to be used for differencing 82 even if the func is not SNESComputeFunction. See: MatSNESMFUseBase() 83 84 85 */ 86 static PetscErrorCode MatAssemblyEnd_SNESMF_UseBase(Mat J,MatAssemblyType mt) 87 { 88 PetscErrorCode ierr; 89 MatMFFD j = (MatMFFD)J->data; 90 SNES snes = (SNES)j->ctx; 91 Vec u,f; 92 93 PetscFunctionBegin; 94 ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr); 95 96 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 97 ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr); 98 ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 /* 103 This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer 104 uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF(). 105 */ 106 static PetscErrorCode MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F) 107 { 108 PetscErrorCode ierr; 109 110 PetscFunctionBegin; 111 ierr = MatMFFDSetBase_MFFD(J,U,F);CHKERRQ(ierr); 112 J->ops->assemblyend = MatAssemblyEnd_MFFD; 113 PetscFunctionReturn(0); 114 } 115 116 static PetscErrorCode MatSNESMFSetReuseBase_SNESMF(Mat J,PetscBool use) 117 { 118 PetscFunctionBegin; 119 if (use) { 120 J->ops->assemblyend = MatAssemblyEnd_SNESMF_UseBase; 121 } else { 122 J->ops->assemblyend = MatAssemblyEnd_SNESMF; 123 } 124 PetscFunctionReturn(0); 125 } 126 127 /*@ 128 MatSNESMFSetReuseBase - Causes the base vector to be used for differencing even if the function provided to SNESSetFunction() is not the 129 same as that provided to MatMFFDSetFunction(). 130 131 Logically Collective on Mat 132 133 Input Parameters: 134 + J - the MatMFFD matrix 135 - use - if true always reuse the base vector instead of recomputing f(u) even if the function in the MatSNESMF is 136 not SNESComputeFunction() 137 138 Notes: Care must be taken when using this routine to insure that the function provided to MatMFFDSetFunction(), call it F_MF() is compatible with 139 with that provided to SNESSetFunction(), call it F_SNES(). That is, (F_MF(u + h*d) - F_SNES(u))/h has to approximate the derivative 140 141 Developer Notes: This was provided for the MOOSE team who desired to have a SNESSetFunction() function that could change configurations due 142 to contacts while the function provided to MatMFFDSetFunction() cannot. Except for the possibility of changing the configuration 143 both functions compute the same mathematical function so the differencing makes sense. 144 145 Level: advanced 146 147 .seealso: MatCreateSNESMF(), MatSNESMFGetReuseBase() 148 @*/ 149 PetscErrorCode MatSNESMFSetReuseBase(Mat J,PetscBool use) 150 { 151 PetscErrorCode ierr; 152 153 PetscFunctionBegin; 154 PetscValidHeaderSpecific(J,MAT_CLASSID,1); 155 ierr = PetscTryMethod(J,"MatSNESMFSetReuseBase_C",(Mat,PetscBool),(J,use));CHKERRQ(ierr); 156 PetscFunctionReturn(0); 157 } 158 159 static PetscErrorCode MatSNESMFGetReuseBase_SNESMF(Mat J,PetscBool *use) 160 { 161 PetscFunctionBegin; 162 if (J->ops->assemblyend == MatAssemblyEnd_SNESMF_UseBase) *use = PETSC_TRUE; 163 else *use = PETSC_FALSE; 164 PetscFunctionReturn(0); 165 } 166 167 /*@ 168 MatSNESMFGetReuseBase - Determines if the base vector is to be used for differencing even if the function provided to SNESSetFunction() is not the 169 same as that provided to MatMFFDSetFunction(). 170 171 Logically Collective on Mat 172 173 Input Parameter: 174 . J - the MatMFFD matrix 175 176 Output Parameter: 177 . use - if true always reuse the base vector instead of recomputing f(u) even if the function in the MatSNESMF is 178 not SNESComputeFunction() 179 180 Notes: Care must be taken when using this routine to insure that the function provided to MatMFFDSetFunction(), call it F_MF() is compatible with 181 with that provided to SNESSetFunction(), call it F_SNES(). That is, (F_MF(u + h*d) - F_SNES(u))/h has to approximate the derivative 182 183 Developer Notes: This was provided for the MOOSE team who desired to have a SNESSetFunction() function that could change configurations due 184 to contacts while the function provided to MatMFFDSetFunction() cannot. Except for the possibility of changing the configuration 185 both functions compute the same mathematical function so the differencing makes sense. 186 187 Level: advanced 188 189 .seealso: MatCreateSNESMF(), MatSNESMFSetReuseBase() 190 @*/ 191 PetscErrorCode MatSNESMFGetReuseBase(Mat J,PetscBool *use) 192 { 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 PetscValidHeaderSpecific(J,MAT_CLASSID,1); 197 ierr = PetscUseMethod(J,"MatSNESMFGetReuseBase_C",(Mat,PetscBool*),(J,use));CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 /*@ 202 MatCreateSNESMF - Creates a matrix-free matrix context for use with 203 a SNES solver. This matrix can be used as the Jacobian argument for 204 the routine SNESSetJacobian(). See MatCreateMFFD() for details on how 205 the finite difference computation is done. 206 207 Collective on SNES and Vec 208 209 Input Parameters: 210 . snes - the SNES context 211 212 Output Parameter: 213 . J - the matrix-free matrix 214 215 Level: advanced 216 217 218 Notes: 219 You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different 220 preconditioner matrix 221 222 If you wish to provide a different function to do differencing on to compute the matrix-free operator than 223 that provided to SNESSetFunction() then call MatMFFDSetFunction() with your function after this call. 224 225 The difference between this routine and MatCreateMFFD() is that this matrix 226 automatically gets the current base vector from the SNES object and not from an 227 explicit call to MatMFFDSetBase(). 228 229 Warning: 230 If MatMFFDSetBase() is ever called on jac then this routine will NO longer get 231 the x from the SNES object and MatMFFDSetBase() must from that point on be used to 232 change the base vector x. 233 234 Warning: 235 Using a different function for the differencing will not work if you are using non-linear left preconditioning. 236 237 238 .seealso: MatDestroy(), MatMFFDSetFunction(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin() 239 MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(), 240 MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian(), MatSNESMFSetReuseBase(), MatSNESMFGetReuseBase() 241 242 @*/ 243 PetscErrorCode MatCreateSNESMF(SNES snes,Mat *J) 244 { 245 PetscErrorCode ierr; 246 PetscInt n,N; 247 MatMFFD mf; 248 249 PetscFunctionBegin; 250 if (snes->vec_func) { 251 ierr = VecGetLocalSize(snes->vec_func,&n);CHKERRQ(ierr); 252 ierr = VecGetSize(snes->vec_func,&N);CHKERRQ(ierr); 253 } else if (snes->dm) { 254 Vec tmp; 255 ierr = DMGetGlobalVector(snes->dm,&tmp);CHKERRQ(ierr); 256 ierr = VecGetLocalSize(tmp,&n);CHKERRQ(ierr); 257 ierr = VecGetSize(tmp,&N);CHKERRQ(ierr); 258 ierr = DMRestoreGlobalVector(snes->dm,&tmp);CHKERRQ(ierr); 259 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first"); 260 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);CHKERRQ(ierr); 261 262 mf = (MatMFFD)(*J)->data; 263 mf->ctx = snes; 264 265 if (snes->npc && snes->npcside== PC_LEFT) { 266 ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunctionDefaultNPC,snes);CHKERRQ(ierr); 267 } else { 268 ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);CHKERRQ(ierr); 269 } 270 271 (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF; 272 273 ierr = PetscObjectComposeFunction((PetscObject)*J,"MatMFFDSetBase_C",MatMFFDSetBase_SNESMF);CHKERRQ(ierr); 274 ierr = PetscObjectComposeFunction((PetscObject)*J,"MatSNESMFSetReuseBase_C",MatSNESMFSetReuseBase_SNESMF);CHKERRQ(ierr); 275 ierr = PetscObjectComposeFunction((PetscObject)*J,"MatSNESMFGetReuseBase_C",MatSNESMFGetReuseBase_SNESMF);CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278