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 MatSNESMFGetSNES - returns the SNES associated with a matrix created with MatCreateSNESMF() 55 56 Not collective 57 58 Input Parameter: 59 . J - the matrix 60 61 Output Parameter: 62 . snes - the SNES object 63 64 .seealso: MatCreateSNESMF() 65 @*/ 66 PetscErrorCode MatSNESMFGetSNES(Mat J,SNES *snes) 67 { 68 MatMFFD j = (MatMFFD)J->data; 69 70 PetscFunctionBegin; 71 *snes = (SNES)j->ctx; 72 PetscFunctionReturn(0); 73 } 74 75 /* 76 MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the 77 base from the SNES context 78 79 */ 80 static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt) 81 { 82 PetscErrorCode ierr; 83 MatMFFD j = (MatMFFD)J->data; 84 SNES snes = (SNES)j->ctx; 85 Vec u,f; 86 87 PetscFunctionBegin; 88 ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr); 89 90 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 91 if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) { 92 ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr); 93 ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr); 94 } else { 95 /* f value known by SNES is not correct for other differencing function */ 96 ierr = MatMFFDSetBase_MFFD(J,u,NULL);CHKERRQ(ierr); 97 } 98 PetscFunctionReturn(0); 99 } 100 101 /* 102 This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer 103 uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF(). 104 */ 105 static PetscErrorCode MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F) 106 { 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 ierr = MatMFFDSetBase_MFFD(J,U,F);CHKERRQ(ierr); 111 112 J->ops->assemblyend = MatAssemblyEnd_MFFD; 113 PetscFunctionReturn(0); 114 } 115 116 /*@ 117 MatCreateSNESMF - Creates a matrix-free matrix context for use with 118 a SNES solver. This matrix can be used as the Jacobian argument for 119 the routine SNESSetJacobian(). See MatCreateMFFD() for details on how 120 the finite difference computation is done. 121 122 Collective on SNES and Vec 123 124 Input Parameters: 125 . snes - the SNES context 126 127 Output Parameter: 128 . J - the matrix-free matrix 129 130 Level: advanced 131 132 133 Notes: 134 You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different 135 preconditioner matrix 136 137 If you wish to provide a different function to do differencing on to compute the matrix-free operator than 138 that provided to SNESSetFunction() then call MatMFFDSetFunction() with your function after this call. 139 140 The difference between this routine and MatCreateMFFD() is that this matrix 141 automatically gets the current base vector from the SNES object and not from an 142 explicit call to MatMFFDSetBase(). 143 144 Warning: 145 If MatMFFDSetBase() is ever called on jac then this routine will NO longer get 146 the x from the SNES object and MatMFFDSetBase() must from that point on be used to 147 change the base vector x. 148 149 Warning: 150 Using a different function for the differencing will not work if you are using non-linear left preconditioning. 151 152 153 .seealso: MatDestroy(), MatMFFDSetFunction(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin() 154 MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(), 155 MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian() 156 157 @*/ 158 PetscErrorCode MatCreateSNESMF(SNES snes,Mat *J) 159 { 160 PetscErrorCode ierr; 161 PetscInt n,N; 162 MatMFFD mf; 163 164 PetscFunctionBegin; 165 if (snes->vec_func) { 166 ierr = VecGetLocalSize(snes->vec_func,&n);CHKERRQ(ierr); 167 ierr = VecGetSize(snes->vec_func,&N);CHKERRQ(ierr); 168 } else if (snes->dm) { 169 Vec tmp; 170 ierr = DMGetGlobalVector(snes->dm,&tmp);CHKERRQ(ierr); 171 ierr = VecGetLocalSize(tmp,&n);CHKERRQ(ierr); 172 ierr = VecGetSize(tmp,&N);CHKERRQ(ierr); 173 ierr = DMRestoreGlobalVector(snes->dm,&tmp);CHKERRQ(ierr); 174 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first"); 175 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);CHKERRQ(ierr); 176 177 mf = (MatMFFD)(*J)->data; 178 mf->ctx = snes; 179 180 if (snes->npc && snes->npcside== PC_LEFT) { 181 ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunctionDefaultNPC,snes);CHKERRQ(ierr); 182 } else { 183 ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);CHKERRQ(ierr); 184 } 185 186 (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF; 187 188 ierr = PetscObjectComposeFunction((PetscObject)*J,"MatMFFDSetBase_C",MatMFFDSetBase_SNESMF);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191