1 #include <petsc/private/dmimpl.h> 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 4 typedef struct { 5 PetscErrorCode (*residuallocal)(DM,Vec,Vec,void*); 6 PetscErrorCode (*jacobianlocal)(DM,Vec,Mat,Mat,void*); 7 PetscErrorCode (*boundarylocal)(DM,Vec,void*); 8 void *residuallocalctx; 9 void *jacobianlocalctx; 10 void *boundarylocalctx; 11 } DMSNES_Local; 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "DMSNESDestroy_DMLocal" 15 static PetscErrorCode DMSNESDestroy_DMLocal(DMSNES 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__ "DMSNESDuplicate_DMLocal" 26 static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm,DMSNES sdm) 27 { 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 ierr = PetscNewLog(sdm,(DMSNES_Local**)&sdm->data);CHKERRQ(ierr); 32 if (oldsdm->data) { 33 ierr = PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_Local));CHKERRQ(ierr); 34 } 35 PetscFunctionReturn(0); 36 } 37 38 #undef __FUNCT__ 39 #define __FUNCT__ "DMLocalSNESGetContext" 40 static PetscErrorCode DMLocalSNESGetContext(DM dm,DMSNES sdm,DMSNES_Local **dmlocalsnes) 41 { 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 *dmlocalsnes = NULL; 46 if (!sdm->data) { 47 ierr = PetscNewLog(dm,(DMSNES_Local**)&sdm->data);CHKERRQ(ierr); 48 49 sdm->ops->destroy = DMSNESDestroy_DMLocal; 50 sdm->ops->duplicate = DMSNESDuplicate_DMLocal; 51 } 52 *dmlocalsnes = (DMSNES_Local*)sdm->data; 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "SNESComputeFunction_DMLocal" 58 static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes,Vec X,Vec F,void *ctx) 59 { 60 PetscErrorCode ierr; 61 DM dm; 62 DMSNES_Local *dmlocalsnes = (DMSNES_Local*)ctx; 63 Vec Xloc,Floc; 64 65 PetscFunctionBegin; 66 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 67 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 68 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 69 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 70 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 71 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr); 72 ierr = VecZeroEntries(Xloc);CHKERRQ(ierr); 73 if (dmlocalsnes->boundarylocal) {ierr = (*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);CHKERRQ(ierr);} 74 ierr = VecZeroEntries(Floc);CHKERRQ(ierr); 75 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 76 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 77 CHKMEMQ; 78 ierr = (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);CHKERRQ(ierr); 79 CHKMEMQ; 80 ierr = VecZeroEntries(F);CHKERRQ(ierr); 81 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 82 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 83 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr); 84 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "SNESComputeJacobian_DMLocal" 90 static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx) 91 { 92 PetscErrorCode ierr; 93 DM dm; 94 DMSNES_Local *dmlocalsnes = (DMSNES_Local*)ctx; 95 Vec Xloc; 96 97 PetscFunctionBegin; 98 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 99 if (dmlocalsnes->jacobianlocal) { 100 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 101 ierr = VecZeroEntries(Xloc);CHKERRQ(ierr); 102 if (dmlocalsnes->boundarylocal) {ierr = (*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);CHKERRQ(ierr);} 103 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 104 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 105 CHKMEMQ; 106 ierr = (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);CHKERRQ(ierr); 107 CHKMEMQ; 108 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 109 } else { 110 MatFDColoring fdcoloring; 111 ierr = PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);CHKERRQ(ierr); 112 if (!fdcoloring) { 113 ISColoring coloring; 114 115 ierr = DMCreateColoring(dm,dm->coloringtype,&coloring);CHKERRQ(ierr); 116 ierr = MatFDColoringCreate(B,coloring,&fdcoloring);CHKERRQ(ierr); 117 ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr); 118 switch (dm->coloringtype) { 119 case IS_COLORING_GLOBAL: 120 ierr = MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);CHKERRQ(ierr); 121 break; 122 default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]); 123 } 124 ierr = PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);CHKERRQ(ierr); 125 ierr = MatFDColoringSetFromOptions(fdcoloring);CHKERRQ(ierr); 126 ierr = MatFDColoringSetUp(B,coloring,fdcoloring);CHKERRQ(ierr); 127 ierr = PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);CHKERRQ(ierr); 128 ierr = PetscObjectDereference((PetscObject)fdcoloring);CHKERRQ(ierr); 129 130 /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call 131 * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the 132 * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually 133 * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be 134 * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed. 135 */ 136 ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr); 137 } 138 ierr = MatFDColoringApply(B,fdcoloring,X,snes);CHKERRQ(ierr); 139 } 140 /* This will be redundant if the user called both, but it's too common to forget. */ 141 if (A != B) { 142 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 143 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 144 } 145 PetscFunctionReturn(0); 146 } 147 148 #undef __FUNCT__ 149 #define __FUNCT__ "DMSNESSetFunctionLocal" 150 /*@C 151 DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector 152 containing the local vector information PLUS ghost point information. It should compute a result for all local 153 elements and DMSNES will automatically accumulate the overlapping values. 154 155 Logically Collective 156 157 Input Arguments: 158 + dm - DM to associate callback with 159 . func - local residual evaluation 160 - ctx - optional context for local residual evaluation 161 162 Level: beginner 163 164 .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 165 @*/ 166 PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx) 167 { 168 PetscErrorCode ierr; 169 DMSNES sdm; 170 DMSNES_Local *dmlocalsnes; 171 172 PetscFunctionBegin; 173 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 174 ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr); 175 ierr = DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);CHKERRQ(ierr); 176 177 dmlocalsnes->residuallocal = func; 178 dmlocalsnes->residuallocalctx = ctx; 179 180 ierr = DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes);CHKERRQ(ierr); 181 if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 182 ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);CHKERRQ(ierr); 183 } 184 PetscFunctionReturn(0); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "DMSNESSetBoundaryLocal" 189 /*@C 190 DMSNESSetBoundaryLocal - set a local boundary value function. This function is called with local vector 191 containing the local vector information PLUS ghost point information. It should insert values into the local 192 vector that do not come from the global vector, such as essential boundary condition data. 193 194 Logically Collective 195 196 Input Arguments: 197 + dm - DM to associate callback with 198 . func - local boundary value evaluation 199 - ctx - optional context for local boundary value evaluation 200 201 Level: intermediate 202 203 .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal() 204 @*/ 205 PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx) 206 { 207 PetscErrorCode ierr; 208 DMSNES sdm; 209 DMSNES_Local *dmlocalsnes; 210 211 PetscFunctionBegin; 212 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 213 ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr); 214 ierr = DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);CHKERRQ(ierr); 215 216 dmlocalsnes->boundarylocal = func; 217 dmlocalsnes->boundarylocalctx = ctx; 218 219 PetscFunctionReturn(0); 220 } 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "DMSNESSetJacobianLocal" 224 /*@C 225 DMSNESSetJacobianLocal - set a local Jacobian evaluation function 226 227 Logically Collective 228 229 Input Arguments: 230 + dm - DM to associate callback with 231 . func - local Jacobian evaluation 232 - ctx - optional context for local Jacobian evaluation 233 234 Level: beginner 235 236 .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 237 @*/ 238 PetscErrorCode DMSNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Mat,Mat,void*),void *ctx) 239 { 240 PetscErrorCode ierr; 241 DMSNES sdm; 242 DMSNES_Local *dmlocalsnes; 243 244 PetscFunctionBegin; 245 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 246 ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr); 247 ierr = DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);CHKERRQ(ierr); 248 249 dmlocalsnes->jacobianlocal = func; 250 dmlocalsnes->jacobianlocalctx = ctx; 251 252 ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256