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