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 static PetscErrorCode DMSNESDestroy_DMLocal(DMSNES sdm) 14 { 15 PetscFunctionBegin; 16 PetscCall(PetscFree(sdm->data)); 17 sdm->data = NULL; 18 PetscFunctionReturn(0); 19 } 20 21 static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm,DMSNES sdm) 22 { 23 PetscFunctionBegin; 24 if (sdm->data != oldsdm->data) { 25 PetscCall(PetscFree(sdm->data)); 26 PetscCall(PetscNewLog(sdm,(DMSNES_Local**)&sdm->data)); 27 if (oldsdm->data) PetscCall(PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_Local))); 28 } 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode DMLocalSNESGetContext(DM dm,DMSNES sdm,DMSNES_Local **dmlocalsnes) 33 { 34 PetscFunctionBegin; 35 *dmlocalsnes = NULL; 36 if (!sdm->data) { 37 PetscCall(PetscNewLog(dm,(DMSNES_Local**)&sdm->data)); 38 39 sdm->ops->destroy = DMSNESDestroy_DMLocal; 40 sdm->ops->duplicate = DMSNESDuplicate_DMLocal; 41 } 42 *dmlocalsnes = (DMSNES_Local*)sdm->data; 43 PetscFunctionReturn(0); 44 } 45 46 static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes,Vec X,Vec F,void *ctx) 47 { 48 DMSNES_Local *dmlocalsnes = (DMSNES_Local *) ctx; 49 DM dm; 50 Vec Xloc,Floc; 51 PetscBool transform; 52 53 PetscFunctionBegin; 54 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 55 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 56 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 57 PetscCall(SNESGetDM(snes,&dm)); 58 PetscCall(DMGetLocalVector(dm,&Xloc)); 59 PetscCall(DMGetLocalVector(dm,&Floc)); 60 PetscCall(VecZeroEntries(Xloc)); 61 PetscCall(VecZeroEntries(Floc)); 62 /* Non-conforming routines needs boundary values before G2L */ 63 if (dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx)); 64 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 65 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 66 /* Need to reset boundary values if we transformed */ 67 PetscCall(DMHasBasisTransform(dm, &transform)); 68 if (transform && dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx)); 69 CHKMEMQ; 70 PetscCall((*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx)); 71 CHKMEMQ; 72 PetscCall(VecZeroEntries(F)); 73 PetscCall(DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F)); 74 PetscCall(DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F)); 75 PetscCall(DMRestoreLocalVector(dm,&Floc)); 76 PetscCall(DMRestoreLocalVector(dm,&Xloc)); 77 { 78 char name[PETSC_MAX_PATH_LEN]; 79 char oldname[PETSC_MAX_PATH_LEN]; 80 const char *tmp; 81 PetscInt it; 82 83 PetscCall(SNESGetIterationNumber(snes, &it)); 84 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Solution, Iterate %d", (int) it)); 85 PetscCall(PetscObjectGetName((PetscObject) X, &tmp)); 86 PetscCall(PetscStrncpy(oldname, tmp, PETSC_MAX_PATH_LEN-1)); 87 PetscCall(PetscObjectSetName((PetscObject) X, name)); 88 PetscCall(VecViewFromOptions(X, (PetscObject) snes, "-dmsnes_solution_vec_view")); 89 PetscCall(PetscObjectSetName((PetscObject) X, oldname)); 90 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Residual, Iterate %d", (int) it)); 91 PetscCall(PetscObjectSetName((PetscObject) F, name)); 92 PetscCall(VecViewFromOptions(F, (PetscObject) snes, "-dmsnes_residual_vec_view")); 93 } 94 PetscFunctionReturn(0); 95 } 96 97 static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx) 98 { 99 DMSNES_Local *dmlocalsnes = (DMSNES_Local *) ctx; 100 DM dm; 101 Vec Xloc; 102 PetscBool transform; 103 104 PetscFunctionBegin; 105 PetscCall(SNESGetDM(snes,&dm)); 106 if (dmlocalsnes->jacobianlocal) { 107 PetscCall(DMGetLocalVector(dm,&Xloc)); 108 PetscCall(VecZeroEntries(Xloc)); 109 /* Non-conforming routines needs boundary values before G2L */ 110 if (dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx)); 111 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 112 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 113 /* Need to reset boundary values if we transformed */ 114 PetscCall(DMHasBasisTransform(dm, &transform)); 115 if (transform && dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx)); 116 CHKMEMQ; 117 PetscCall((*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx)); 118 CHKMEMQ; 119 PetscCall(DMRestoreLocalVector(dm,&Xloc)); 120 } else { 121 MatFDColoring fdcoloring; 122 PetscCall(PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring)); 123 if (!fdcoloring) { 124 ISColoring coloring; 125 126 PetscCall(DMCreateColoring(dm,dm->coloringtype,&coloring)); 127 PetscCall(MatFDColoringCreate(B,coloring,&fdcoloring)); 128 PetscCall(ISColoringDestroy(&coloring)); 129 switch (dm->coloringtype) { 130 case IS_COLORING_GLOBAL: 131 PetscCall(MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes)); 132 break; 133 default: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]); 134 } 135 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix)); 136 PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 137 PetscCall(MatFDColoringSetUp(B,coloring,fdcoloring)); 138 PetscCall(PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring)); 139 PetscCall(PetscObjectDereference((PetscObject)fdcoloring)); 140 141 /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call 142 * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the 143 * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually 144 * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be 145 * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed. 146 */ 147 PetscCall(PetscObjectDereference((PetscObject)dm)); 148 } 149 PetscCall(MatFDColoringApply(B,fdcoloring,X,snes)); 150 } 151 /* This will be redundant if the user called both, but it's too common to forget. */ 152 if (A != B) { 153 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 154 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 155 } 156 PetscFunctionReturn(0); 157 } 158 159 /*@C 160 DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector 161 containing the local vector information PLUS ghost point information. It should compute a result for all local 162 elements and DMSNES will automatically accumulate the overlapping values. 163 164 Logically Collective 165 166 Input Parameters: 167 + dm - DM to associate callback with 168 . func - local residual evaluation 169 - ctx - optional context for local residual evaluation 170 171 Level: beginner 172 173 .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 174 @*/ 175 PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx) 176 { 177 DMSNES sdm; 178 DMSNES_Local *dmlocalsnes; 179 180 PetscFunctionBegin; 181 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 182 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 183 PetscCall(DMLocalSNESGetContext(dm,sdm,&dmlocalsnes)); 184 185 dmlocalsnes->residuallocal = func; 186 dmlocalsnes->residuallocalctx = ctx; 187 188 PetscCall(DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes)); 189 if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 190 PetscCall(DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes)); 191 } 192 PetscFunctionReturn(0); 193 } 194 195 /*@C 196 DMSNESSetBoundaryLocal - set a local boundary value function. This function is called with local vector 197 containing the local vector information PLUS ghost point information. It should insert values into the local 198 vector that do not come from the global vector, such as essential boundary condition data. 199 200 Logically Collective 201 202 Input Parameters: 203 + dm - DM to associate callback with 204 . func - local boundary value evaluation 205 - ctx - optional context for local boundary value evaluation 206 207 Level: intermediate 208 209 .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal() 210 @*/ 211 PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx) 212 { 213 DMSNES sdm; 214 DMSNES_Local *dmlocalsnes; 215 216 PetscFunctionBegin; 217 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 218 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 219 PetscCall(DMLocalSNESGetContext(dm,sdm,&dmlocalsnes)); 220 221 dmlocalsnes->boundarylocal = func; 222 dmlocalsnes->boundarylocalctx = ctx; 223 224 PetscFunctionReturn(0); 225 } 226 227 /*@C 228 DMSNESSetJacobianLocal - set a local Jacobian evaluation function 229 230 Logically Collective 231 232 Input Parameters: 233 + dm - DM to associate callback with 234 . func - local Jacobian evaluation 235 - ctx - optional context for local Jacobian evaluation 236 237 Level: beginner 238 239 .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 240 @*/ 241 PetscErrorCode DMSNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Mat,Mat,void*),void *ctx) 242 { 243 DMSNES sdm; 244 DMSNES_Local *dmlocalsnes; 245 246 PetscFunctionBegin; 247 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 248 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 249 PetscCall(DMLocalSNESGetContext(dm,sdm,&dmlocalsnes)); 250 251 dmlocalsnes->jacobianlocal = func; 252 dmlocalsnes->jacobianlocalctx = ctx; 253 254 PetscCall(DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes)); 255 PetscFunctionReturn(0); 256 } 257 258 /*@C 259 DMSNESGetFunctionLocal - get the local residual evaluation function information set with DMSNESSetFunctionLocal. 260 261 Not Collective 262 263 Input Parameter: 264 . dm - DM with the associated callback 265 266 Output Parameters: 267 + func - local residual evaluation 268 - ctx - context for local residual evaluation 269 270 Level: beginner 271 272 .seealso: DMSNESSetFunction(), DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 273 @*/ 274 PetscErrorCode DMSNESGetFunctionLocal(DM dm,PetscErrorCode (**func)(DM,Vec,Vec,void*),void **ctx) 275 { 276 DMSNES sdm; 277 DMSNES_Local *dmlocalsnes; 278 279 PetscFunctionBegin; 280 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 281 PetscCall(DMGetDMSNES(dm,&sdm)); 282 PetscCall(DMLocalSNESGetContext(dm,sdm,&dmlocalsnes)); 283 if (func) *func = dmlocalsnes->residuallocal; 284 if (ctx) *ctx = dmlocalsnes->residuallocalctx; 285 PetscFunctionReturn(0); 286 } 287 288 /*@C 289 DMSNESGetBoundaryLocal - get the local boundary value function set with DMSNESSetBoundaryLocal. 290 291 Not Collective 292 293 Input Parameter: 294 . dm - DM with the associated callback 295 296 Output Parameters: 297 + func - local boundary value evaluation 298 - ctx - context for local boundary value evaluation 299 300 Level: intermediate 301 302 .seealso: DMSNESSetFunctionLocal(), DMSNESSetBoundaryLocal(), DMDASNESSetJacobianLocal() 303 @*/ 304 PetscErrorCode DMSNESGetBoundaryLocal(DM dm,PetscErrorCode (**func)(DM,Vec,void*),void **ctx) 305 { 306 DMSNES sdm; 307 DMSNES_Local *dmlocalsnes; 308 309 PetscFunctionBegin; 310 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 311 PetscCall(DMGetDMSNES(dm,&sdm)); 312 PetscCall(DMLocalSNESGetContext(dm,sdm,&dmlocalsnes)); 313 if (func) *func = dmlocalsnes->boundarylocal; 314 if (ctx) *ctx = dmlocalsnes->boundarylocalctx; 315 PetscFunctionReturn(0); 316 } 317 318 /*@C 319 DMSNESGetJacobianLocal - the local Jacobian evaluation function set with DMSNESSetJacobianLocal. 320 321 Logically Collective 322 323 Input Parameter: 324 . dm - DM with the associated callback 325 326 Output Parameters: 327 + func - local Jacobian evaluation 328 - ctx - context for local Jacobian evaluation 329 330 Level: beginner 331 332 .seealso: DMSNESSetJacobianLocal(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 333 @*/ 334 PetscErrorCode DMSNESGetJacobianLocal(DM dm,PetscErrorCode (**func)(DM,Vec,Mat,Mat,void*),void **ctx) 335 { 336 DMSNES sdm; 337 DMSNES_Local *dmlocalsnes; 338 339 PetscFunctionBegin; 340 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 341 PetscCall(DMGetDMSNES(dm,&sdm)); 342 PetscCall(DMLocalSNESGetContext(dm,sdm,&dmlocalsnes)); 343 if (func) *func = dmlocalsnes->jacobianlocal; 344 if (ctx) *ctx = dmlocalsnes->jacobianlocalctx; 345 PetscFunctionReturn(0); 346 } 347