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