1 #include <petscdmda.h> /*I "petscdmda.h" I*/ 2 #include <petsc-private/dmimpl.h> 3 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 4 5 /* This structure holds the user-provided DMDA callbacks */ 6 typedef struct { 7 PetscErrorCode (*residuallocal)(DMDALocalInfo*,void*,void*,void*); 8 PetscErrorCode (*jacobianlocal)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*); 9 void *residuallocalctx; 10 void *jacobianlocalctx; 11 InsertMode residuallocalimode; 12 } DM_DA_SNES; 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "SNESDMDestroy_DMDA" 16 static PetscErrorCode SNESDMDestroy_DMDA(SNESDM sdm) 17 { 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 ierr = PetscFree(sdm->data);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "DMDASNESGetContext" 27 static PetscErrorCode DMDASNESGetContext(DM dm,SNESDM sdm,DM_DA_SNES **dmdasnes) 28 { 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 *dmdasnes = PETSC_NULL; 33 if (!sdm->data) { 34 ierr = PetscNewLog(dm,DM_DA_SNES,&sdm->data);CHKERRQ(ierr); 35 sdm->destroy = SNESDMDestroy_DMDA; 36 } 37 *dmdasnes = (DM_DA_SNES*)sdm->data; 38 PetscFunctionReturn(0); 39 } 40 41 #undef __FUNCT__ 42 #define __FUNCT__ "SNESComputeFunction_DMDA" 43 /* 44 This function should eventually replace: 45 DMDAComputeFunction() and DMDAComputeFunction1() 46 */ 47 static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx) 48 { 49 PetscErrorCode ierr; 50 DM dm; 51 DM_DA_SNES *dmdasnes = (DM_DA_SNES*)ctx; 52 DMDALocalInfo info; 53 Vec Xloc; 54 void *x,*f; 55 56 PetscFunctionBegin; 57 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 58 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 59 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 60 if (!dmdasnes->residuallocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context"); 61 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 62 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 63 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 64 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 65 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 66 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 67 switch (dmdasnes->residuallocalimode) { 68 case INSERT_VALUES: { 69 ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr); 70 CHKMEMQ; 71 ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr); 72 CHKMEMQ; 73 ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr); 74 } break; 75 case ADD_VALUES: { 76 Vec Floc; 77 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr); 78 ierr = VecZeroEntries(Floc);CHKERRQ(ierr); 79 ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr); 80 CHKMEMQ; 81 ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr); 82 CHKMEMQ; 83 ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr); 84 ierr = VecZeroEntries(F);CHKERRQ(ierr); 85 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 86 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 87 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr); 88 } break; 89 default: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode); 90 } 91 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 92 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "SNESComputeJacobian_DMDA" 98 /* 99 This function should eventually replace: 100 DMComputeJacobian() and DMDAComputeJacobian1() 101 */ 102 static PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx) 103 { 104 PetscErrorCode ierr; 105 DM dm; 106 DM_DA_SNES *dmdasnes = (DM_DA_SNES*)ctx; 107 DMDALocalInfo info; 108 Vec Xloc; 109 void *x; 110 111 PetscFunctionBegin; 112 if (!dmdasnes->residuallocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context"); 113 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 114 115 if (dmdasnes->jacobianlocal) { 116 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 117 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 118 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 119 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 120 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 121 CHKMEMQ; 122 ierr = (*dmdasnes->jacobianlocal)(&info,x,*A,*B,mstr,dmdasnes->jacobianlocalctx);CHKERRQ(ierr); 123 CHKMEMQ; 124 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 125 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 126 } else { 127 MatFDColoring fdcoloring; 128 ierr = PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);CHKERRQ(ierr); 129 if (!fdcoloring) { 130 ISColoring coloring; 131 132 ierr = DMCreateColoring(dm,dm->coloringtype,MATAIJ,&coloring);CHKERRQ(ierr); 133 ierr = MatFDColoringCreate(*B,coloring,&fdcoloring);CHKERRQ(ierr); 134 ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr); 135 switch (dm->coloringtype) { 136 case IS_COLORING_GLOBAL: 137 ierr = MatFDColoringSetFunction(fdcoloring,(PetscErrorCode(*)(void))SNESComputeFunction_DMDA,dmdasnes);CHKERRQ(ierr); 138 break; 139 default: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]); 140 } 141 ierr = PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);CHKERRQ(ierr); 142 ierr = MatFDColoringSetFromOptions(fdcoloring);CHKERRQ(ierr); 143 ierr = PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);CHKERRQ(ierr); 144 ierr = PetscObjectDereference((PetscObject)fdcoloring);CHKERRQ(ierr); 145 146 /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call 147 * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the 148 * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually 149 * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be 150 * taken when the PetscOList for the Vec inside MatFDColoring is destroyed. 151 */ 152 ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr); 153 } 154 *mstr = SAME_NONZERO_PATTERN; 155 ierr = MatFDColoringApply(*B,fdcoloring,X,mstr,snes);CHKERRQ(ierr); 156 } 157 /* This will be redundant if the user called both, but it's too common to forget. */ 158 if (*A != *B) { 159 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 160 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 161 } 162 PetscFunctionReturn(0); 163 } 164 165 #undef __FUNCT__ 166 #define __FUNCT__ "DMDASNESSetFunctionLocal" 167 /*@C 168 DMDASNESSetFunctionLocal - set a local residual evaluation function 169 170 Logically Collective 171 172 Input Arguments: 173 + dm - DM to associate callback with 174 . func - local residual evaluation 175 - ctx - optional context for local residual evaluation 176 177 Calling sequence for func: 178 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 179 . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part 180 . x - dimensional pointer to state at which to evaluate residual 181 . f - dimensional pointer to residual, write the residual here 182 - ctx - optional context passed above 183 184 Level: beginner 185 186 .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 187 @*/ 188 PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx) 189 { 190 PetscErrorCode ierr; 191 SNESDM sdm; 192 DM_DA_SNES *dmdasnes; 193 194 PetscFunctionBegin; 195 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 196 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 197 ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr); 198 dmdasnes->residuallocalimode = imode; 199 dmdasnes->residuallocal = func; 200 dmdasnes->residuallocalctx = ctx; 201 ierr = DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);CHKERRQ(ierr); 202 if (!sdm->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 203 ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);CHKERRQ(ierr); 204 } 205 PetscFunctionReturn(0); 206 } 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "DMDASNESSetJacobianLocal" 210 /*@C 211 DMDASNESSetJacobianLocal - set a local residual evaluation function 212 213 Logically Collective 214 215 Input Arguments: 216 + dm - DM to associate callback with 217 . func - local residual evaluation 218 - ctx - optional context for local residual evaluation 219 220 Calling sequence for func: 221 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 222 . x - dimensional pointer to state at which to evaluate residual 223 . f - dimensional pointer to residual, write the residual here 224 - ctx - optional context passed above 225 226 Level: beginner 227 228 .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 229 @*/ 230 PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx) 231 { 232 PetscErrorCode ierr; 233 SNESDM sdm; 234 DM_DA_SNES *dmdasnes; 235 236 PetscFunctionBegin; 237 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 238 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 239 ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr); 240 dmdasnes->jacobianlocal = func; 241 dmdasnes->jacobianlocalctx = ctx; 242 ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245