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(PETSC_SUCCESS); 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(PetscNew((DMSNES_Local **)&sdm->data)); 27 if (oldsdm->data) PetscCall(PetscMemcpy(sdm->data, oldsdm->data, sizeof(DMSNES_Local))); 28 } 29 PetscFunctionReturn(PETSC_SUCCESS); 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(PetscNew((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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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: 134 SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]); 135 } 136 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix)); 137 PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 138 PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring)); 139 PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring)); 140 PetscCall(PetscObjectDereference((PetscObject)fdcoloring)); 141 142 /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call 143 * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the 144 * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually 145 * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be 146 * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed. 147 */ 148 PetscCall(PetscObjectDereference((PetscObject)dm)); 149 } 150 PetscCall(MatFDColoringApply(B, fdcoloring, X, snes)); 151 } 152 /* This will be redundant if the user called both, but it's too common to forget. */ 153 if (A != B) { 154 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 155 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 156 } 157 PetscFunctionReturn(PETSC_SUCCESS); 158 } 159 160 /*@C 161 DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector 162 containing the local vector information PLUS ghost point information. It should compute a result for all local 163 elements and `DMSNES` will automatically accumulate the overlapping values. 164 165 Logically Collective 166 167 Input Parameters: 168 + dm - `DM` to associate callback with 169 . func - local residual evaluation 170 - ctx - optional context for local residual evaluation 171 172 Level: advanced 173 174 .seealso: `DMSNESSetFunction()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 175 @*/ 176 PetscErrorCode DMSNESSetFunctionLocal(DM dm, PetscErrorCode (*func)(DM, Vec, Vec, void *), void *ctx) 177 { 178 DMSNES sdm; 179 DMSNES_Local *dmlocalsnes; 180 181 PetscFunctionBegin; 182 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 183 PetscCall(DMGetDMSNESWrite(dm, &sdm)); 184 PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes)); 185 186 dmlocalsnes->residuallocal = func; 187 dmlocalsnes->residuallocalctx = ctx; 188 189 PetscCall(DMSNESSetFunction(dm, SNESComputeFunction_DMLocal, dmlocalsnes)); 190 if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 191 PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes)); 192 } 193 PetscFunctionReturn(PETSC_SUCCESS); 194 } 195 196 /*@C 197 DMSNESSetBoundaryLocal - set a function to insert, for example, essential boundary conditions into a ghosted solution vector 198 199 Logically Collective 200 201 Input Parameters: 202 + dm - `DM` to associate callback with 203 . func - local boundary value evaluation 204 - ctx - optional context for local boundary value evaluation 205 206 Calling sequence of `func`: 207 + dm - the `DM` context 208 . X - ghosted solution vector, appropriate locations (such as essential boundary condition nodes) should be filled 209 - ctx - a user provided context 210 211 Level: advanced 212 213 .seealso: `DMSNESSetFunctionLocal()`, `DMDASNESSetJacobianLocal()` 214 @*/ 215 PetscErrorCode DMSNESSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec X, void *ctx), void *ctx) 216 { 217 DMSNES sdm; 218 DMSNES_Local *dmlocalsnes; 219 220 PetscFunctionBegin; 221 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 222 PetscCall(DMGetDMSNESWrite(dm, &sdm)); 223 PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes)); 224 225 dmlocalsnes->boundarylocal = func; 226 dmlocalsnes->boundarylocalctx = ctx; 227 228 PetscFunctionReturn(PETSC_SUCCESS); 229 } 230 231 /*@C 232 DMSNESSetJacobianLocal - set a local Jacobian evaluation function 233 234 Logically Collective 235 236 Input Parameters: 237 + dm - DM to associate callback with 238 . func - local Jacobian evaluation 239 - ctx - optional context for local Jacobian evaluation 240 241 Calling sequence of `func`: 242 + dm - the `DM` context 243 . X - current solution vector (ghosted or not?) 244 . J - the Jacobian 245 . Jp - approximate Jacobian used to compute the preconditioner, often `J` 246 - ctx - a user provided context 247 248 Level: advanced 249 250 .seealso: `DMSNESSetJacobian()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 251 @*/ 252 PetscErrorCode DMSNESSetJacobianLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec X, Mat J, Mat Jp, void *ctx), void *ctx) 253 { 254 DMSNES sdm; 255 DMSNES_Local *dmlocalsnes; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 259 PetscCall(DMGetDMSNESWrite(dm, &sdm)); 260 PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes)); 261 262 dmlocalsnes->jacobianlocal = func; 263 dmlocalsnes->jacobianlocalctx = ctx; 264 265 PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes)); 266 PetscFunctionReturn(PETSC_SUCCESS); 267 } 268 269 /*@C 270 DMSNESGetFunctionLocal - get the local residual evaluation function information set with `DMSNESSetFunctionLocal()`. 271 272 Not Collective 273 274 Input Parameter: 275 . dm - `DM` with the associated callback 276 277 Output Parameters: 278 + func - local residual evaluation 279 - ctx - context for local residual evaluation 280 281 Level: beginner 282 283 .seealso: `DMSNESSetFunction()`, `DMSNESSetFunctionLocal()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 284 @*/ 285 PetscErrorCode DMSNESGetFunctionLocal(DM dm, PetscErrorCode (**func)(DM, Vec, Vec, void *), void **ctx) 286 { 287 DMSNES sdm; 288 DMSNES_Local *dmlocalsnes; 289 290 PetscFunctionBegin; 291 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 292 PetscCall(DMGetDMSNES(dm, &sdm)); 293 PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes)); 294 if (func) *func = dmlocalsnes->residuallocal; 295 if (ctx) *ctx = dmlocalsnes->residuallocalctx; 296 PetscFunctionReturn(PETSC_SUCCESS); 297 } 298 299 /*@C 300 DMSNESGetBoundaryLocal - get the local boundary value function set with `DMSNESSetBoundaryLocal()`. 301 302 Not Collective 303 304 Input Parameter: 305 . dm - `DM` with the associated callback 306 307 Output Parameters: 308 + func - local boundary value evaluation 309 - ctx - context for local boundary value evaluation 310 311 Level: intermediate 312 313 .seealso: `DMSNESSetFunctionLocal()`, `DMSNESSetBoundaryLocal()`, `DMDASNESSetJacobianLocal()` 314 @*/ 315 PetscErrorCode DMSNESGetBoundaryLocal(DM dm, PetscErrorCode (**func)(DM, Vec, void *), void **ctx) 316 { 317 DMSNES sdm; 318 DMSNES_Local *dmlocalsnes; 319 320 PetscFunctionBegin; 321 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 322 PetscCall(DMGetDMSNES(dm, &sdm)); 323 PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes)); 324 if (func) *func = dmlocalsnes->boundarylocal; 325 if (ctx) *ctx = dmlocalsnes->boundarylocalctx; 326 PetscFunctionReturn(PETSC_SUCCESS); 327 } 328 329 /*@C 330 DMSNESGetJacobianLocal - the local Jacobian evaluation function set with `DMSNESSetJacobianLocal()`. 331 332 Logically Collective 333 334 Input Parameter: 335 . dm - `DM` with the associated callback 336 337 Output Parameters: 338 + func - local Jacobian evaluation 339 - ctx - context for local Jacobian evaluation 340 341 Level: beginner 342 343 .seealso: `DMSNESSetJacobianLocal()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 344 @*/ 345 PetscErrorCode DMSNESGetJacobianLocal(DM dm, PetscErrorCode (**func)(DM, Vec, Mat, Mat, void *), void **ctx) 346 { 347 DMSNES sdm; 348 DMSNES_Local *dmlocalsnes; 349 350 PetscFunctionBegin; 351 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 352 PetscCall(DMGetDMSNES(dm, &sdm)); 353 PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes)); 354 if (func) *func = dmlocalsnes->jacobianlocal; 355 if (ctx) *ctx = dmlocalsnes->jacobianlocalctx; 356 PetscFunctionReturn(PETSC_SUCCESS); 357 } 358