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 #undef __FUNCT__ 8 #define __FUNCT__ "MatMFFDComputeJacobian" 9 /*@C 10 MatMFFDComputeJacobian - Tells the matrix-free Jacobian object the new location at which 11 Jacobian matrix vector products will be computed at, i.e. J(x) * a. The x is obtained 12 from the SNES object (using SNESGetSolution()). 13 14 Logically Collective on SNES 15 16 Input Parameters: 17 + snes - the nonlinear solver context 18 . x - the point at which the Jacobian vector products will be performed 19 . jac - the matrix-free Jacobian object 20 . B - either the same as jac or another matrix type (ignored) 21 . flag - not relevent for matrix-free form 22 - dummy - the user context (ignored) 23 24 Level: developer 25 26 Warning: 27 If MatMFFDSetBase() is ever called on jac then this routine will NO longer get 28 the x from the SNES object and MatMFFDSetBase() must from that point on be used to 29 change the base vector x. 30 31 Notes: 32 This can be passed into SNESSetJacobian() as the Jacobian evaluation function argument 33 when using a completely matrix-free solver, 34 that is the B matrix is also the same matrix operator. This is used when you select 35 -snes_mf but rarely used directly by users. (All this routine does is call MatAssemblyBegin/End() on 36 the Mat jac.) 37 38 .seealso: MatMFFDGetH(), MatCreateSNESMF(), MatCreateMFFD(), MATMFFD, 39 MatMFFDSetHHistory(), MatMFFDSetFunctionError(), MatCreateMFFD(), SNESSetJacobian() 40 41 @*/ 42 PetscErrorCode MatMFFDComputeJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 43 { 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 48 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 49 PetscFunctionReturn(0); 50 } 51 52 PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat,MatAssemblyType); 53 PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat,Vec,Vec); 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "MatAssemblyEnd_SNESMF" 57 /* 58 MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the 59 base from the SNES context 60 61 */ 62 PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt) 63 { 64 PetscErrorCode ierr; 65 MatMFFD j = (MatMFFD)J->data; 66 SNES snes = (SNES)j->ctx; 67 Vec u,f; 68 69 PetscFunctionBegin; 70 ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr); 71 72 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 73 if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) { 74 ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr); 75 ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr); 76 } else { 77 /* f value known by SNES is not correct for other differencing function */ 78 ierr = MatMFFDSetBase_MFFD(J,u,NULL);CHKERRQ(ierr); 79 } 80 PetscFunctionReturn(0); 81 } 82 83 /* 84 This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer 85 uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF(). 86 */ 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatMFFDSetBase_SNESMF" 89 PetscErrorCode MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F) 90 { 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 ierr = MatMFFDSetBase_MFFD(J,U,F);CHKERRQ(ierr); 95 96 J->ops->assemblyend = MatAssemblyEnd_MFFD; 97 PetscFunctionReturn(0); 98 } 99 100 #undef __FUNCT__ 101 #define __FUNCT__ "MatCreateSNESMF" 102 /*@ 103 MatCreateSNESMF - Creates a matrix-free matrix context for use with 104 a SNES solver. This matrix can be used as the Jacobian argument for 105 the routine SNESSetJacobian(). See MatCreateMFFD() for details on how 106 the finite difference computation is done. 107 108 Collective on SNES and Vec 109 110 Input Parameters: 111 . snes - the SNES context 112 113 Output Parameter: 114 . J - the matrix-free matrix 115 116 Level: advanced 117 118 119 Notes: 120 You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different 121 preconditioner matrix 122 123 If you wish to provide a different function to do differencing on to compute the matrix-free operator than 124 that provided to SNESSetFunction() then call MatMFFDSetFunction() with your function after this call. 125 126 The difference between this routine and MatCreateMFFD() is that this matrix 127 automatically gets the current base vector from the SNES object and not from an 128 explicit call to MatMFFDSetBase(). 129 130 Warning: 131 If MatMFFDSetBase() is ever called on jac then this routine will NO longer get 132 the x from the SNES object and MatMFFDSetBase() must from that point on be used to 133 change the base vector x. 134 135 Warning: 136 Using a different function for the differencing will not work if you are using non-linear left preconditioning. 137 138 139 .seealso: MatDestroy(), MatMFFDSetFunction(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin() 140 MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(), 141 MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian() 142 143 @*/ 144 PetscErrorCode MatCreateSNESMF(SNES snes,Mat *J) 145 { 146 PetscErrorCode ierr; 147 PetscInt n,N; 148 MatMFFD mf; 149 150 PetscFunctionBegin; 151 if (snes->vec_func) { 152 ierr = VecGetLocalSize(snes->vec_func,&n);CHKERRQ(ierr); 153 ierr = VecGetSize(snes->vec_func,&N);CHKERRQ(ierr); 154 } else if (snes->dm) { 155 Vec tmp; 156 ierr = DMGetGlobalVector(snes->dm,&tmp);CHKERRQ(ierr); 157 ierr = VecGetLocalSize(tmp,&n);CHKERRQ(ierr); 158 ierr = VecGetSize(tmp,&N);CHKERRQ(ierr); 159 ierr = DMRestoreGlobalVector(snes->dm,&tmp);CHKERRQ(ierr); 160 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first"); 161 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);CHKERRQ(ierr); 162 mf = (MatMFFD)(*J)->data; 163 mf->ctx = snes; 164 165 if (snes->pc && snes->pcside == PC_LEFT) { 166 ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunctionDefaultNPC,snes);CHKERRQ(ierr); 167 } else { 168 ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);CHKERRQ(ierr); 169 } 170 171 (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF; 172 173 ierr = PetscObjectComposeFunction((PetscObject)*J,"MatMFFDSetBase_C",MatMFFDSetBase_SNESMF);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176 177 178 179 180 181 182