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