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,void*); 9 PetscErrorCode (*objectivelocal)(DMDALocalInfo*,void*,PetscReal*,void*); 10 11 PetscErrorCode (*residuallocalvec)(DMDALocalInfo*,Vec,Vec,void*); 12 PetscErrorCode (*jacobianlocalvec)(DMDALocalInfo*,Vec,Mat,Mat,void*); 13 PetscErrorCode (*objectivelocalvec)(DMDALocalInfo*,Vec,PetscReal*,void*); 14 void *residuallocalctx; 15 void *jacobianlocalctx; 16 void *objectivelocalctx; 17 InsertMode residuallocalimode; 18 19 /* For Picard iteration defined locally */ 20 PetscErrorCode (*rhsplocal)(DMDALocalInfo*,void*,void*,void*); 21 PetscErrorCode (*jacobianplocal)(DMDALocalInfo*,void*,Mat,Mat,void*); 22 void *picardlocalctx; 23 } DMSNES_DA; 24 25 static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm) 26 { 27 PetscFunctionBegin; 28 PetscCall(PetscFree(sdm->data)); 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm,DMSNES sdm) 33 { 34 PetscFunctionBegin; 35 PetscCall(PetscNewLog(sdm,(DMSNES_DA**)&sdm->data)); 36 if (oldsdm->data) { 37 PetscCall(PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_DA))); 38 } 39 PetscFunctionReturn(0); 40 } 41 42 static PetscErrorCode DMDASNESGetContext(DM dm,DMSNES sdm,DMSNES_DA **dmdasnes) 43 { 44 PetscFunctionBegin; 45 *dmdasnes = NULL; 46 if (!sdm->data) { 47 PetscCall(PetscNewLog(dm,(DMSNES_DA**)&sdm->data)); 48 sdm->ops->destroy = DMSNESDestroy_DMDA; 49 sdm->ops->duplicate = DMSNESDuplicate_DMDA; 50 } 51 *dmdasnes = (DMSNES_DA*)sdm->data; 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx) 56 { 57 DM dm; 58 DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx; 59 DMDALocalInfo info; 60 Vec Xloc; 61 void *x,*f; 62 63 PetscFunctionBegin; 64 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 65 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 66 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 67 PetscCheck(dmdasnes->residuallocal || dmdasnes->residuallocalvec,PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context"); 68 PetscCall(SNESGetDM(snes,&dm)); 69 PetscCall(DMGetLocalVector(dm,&Xloc)); 70 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 71 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 72 PetscCall(DMDAGetLocalInfo(dm,&info)); 73 switch (dmdasnes->residuallocalimode) { 74 case INSERT_VALUES: { 75 PetscCall(PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0)); 76 CHKMEMQ; 77 if (dmdasnes->residuallocalvec) { 78 PetscCall((*dmdasnes->residuallocalvec)(&info,Xloc,F,dmdasnes->residuallocalctx)); 79 } else { 80 PetscCall(DMDAVecGetArray(dm,Xloc,&x)); 81 PetscCall(DMDAVecGetArray(dm,F,&f)); 82 PetscCall((*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx)); 83 PetscCall(DMDAVecRestoreArray(dm,Xloc,&x)); 84 PetscCall(DMDAVecRestoreArray(dm,F,&f)); 85 } 86 CHKMEMQ; 87 PetscCall(PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0)); 88 } break; 89 case ADD_VALUES: { 90 Vec Floc; 91 PetscCall(DMGetLocalVector(dm,&Floc)); 92 PetscCall(VecZeroEntries(Floc)); 93 PetscCall(PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0)); 94 CHKMEMQ; 95 if (dmdasnes->residuallocalvec) { 96 PetscCall((*dmdasnes->residuallocalvec)(&info,Xloc,Floc,dmdasnes->residuallocalctx)); 97 } else { 98 PetscCall(DMDAVecGetArray(dm,Xloc,&x)); 99 PetscCall(DMDAVecGetArray(dm,Floc,&f)); 100 PetscCall((*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx)); 101 PetscCall(DMDAVecRestoreArray(dm,Xloc,&x)); 102 PetscCall(DMDAVecRestoreArray(dm,Floc,&f)); 103 } 104 CHKMEMQ; 105 PetscCall(PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0)); 106 PetscCall(VecZeroEntries(F)); 107 PetscCall(DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F)); 108 PetscCall(DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F)); 109 PetscCall(DMRestoreLocalVector(dm,&Floc)); 110 } break; 111 default: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode); 112 } 113 PetscCall(DMRestoreLocalVector(dm,&Xloc)); 114 if (snes->domainerror) { 115 PetscCall(VecSetInf(F)); 116 } 117 PetscFunctionReturn(0); 118 } 119 120 static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx) 121 { 122 DM dm; 123 DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx; 124 DMDALocalInfo info; 125 Vec Xloc; 126 void *x; 127 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 130 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 131 PetscValidRealPointer(ob,3); 132 PetscCheck(dmdasnes->objectivelocal || dmdasnes->objectivelocalvec,PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context"); 133 PetscCall(SNESGetDM(snes,&dm)); 134 PetscCall(DMGetLocalVector(dm,&Xloc)); 135 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 136 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 137 PetscCall(DMDAGetLocalInfo(dm,&info)); 138 CHKMEMQ; 139 if (dmdasnes->objectivelocalvec) { 140 PetscCall((*dmdasnes->objectivelocalvec)(&info,Xloc,ob,dmdasnes->objectivelocalctx)); 141 } else { 142 PetscCall(DMDAVecGetArray(dm,Xloc,&x)); 143 PetscCall((*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx)); 144 PetscCall(DMDAVecRestoreArray(dm,Xloc,&x)); 145 } 146 CHKMEMQ; 147 PetscCall(DMRestoreLocalVector(dm,&Xloc)); 148 PetscFunctionReturn(0); 149 } 150 151 /* Routine is called by example, hence must be labeled PETSC_EXTERN */ 152 PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx) 153 { 154 DM dm; 155 DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx; 156 DMDALocalInfo info; 157 Vec Xloc; 158 void *x; 159 160 PetscFunctionBegin; 161 PetscCheck(dmdasnes->residuallocal || dmdasnes->residuallocalvec,PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context"); 162 PetscCall(SNESGetDM(snes,&dm)); 163 164 if (dmdasnes->jacobianlocal || dmdasnes->jacobianlocalctx) { 165 PetscCall(DMGetLocalVector(dm,&Xloc)); 166 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 167 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 168 PetscCall(DMDAGetLocalInfo(dm,&info)); 169 CHKMEMQ; 170 if (dmdasnes->jacobianlocalvec) { 171 PetscCall((*dmdasnes->jacobianlocalvec)(&info,Xloc,A,B,dmdasnes->jacobianlocalctx)); 172 } else { 173 PetscCall(DMDAVecGetArray(dm,Xloc,&x)); 174 PetscCall((*dmdasnes->jacobianlocal)(&info,x,A,B,dmdasnes->jacobianlocalctx)); 175 PetscCall(DMDAVecRestoreArray(dm,Xloc,&x)); 176 } 177 CHKMEMQ; 178 PetscCall(DMRestoreLocalVector(dm,&Xloc)); 179 } else { 180 MatFDColoring fdcoloring; 181 PetscCall(PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring)); 182 if (!fdcoloring) { 183 ISColoring coloring; 184 185 PetscCall(DMCreateColoring(dm,dm->coloringtype,&coloring)); 186 PetscCall(MatFDColoringCreate(B,coloring,&fdcoloring)); 187 switch (dm->coloringtype) { 188 case IS_COLORING_GLOBAL: 189 PetscCall(MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes)); 190 break; 191 default: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]); 192 } 193 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix)); 194 PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 195 PetscCall(MatFDColoringSetUp(B,coloring,fdcoloring)); 196 PetscCall(ISColoringDestroy(&coloring)); 197 PetscCall(PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring)); 198 PetscCall(PetscObjectDereference((PetscObject)fdcoloring)); 199 200 /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call 201 * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the 202 * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually 203 * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be 204 * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed. 205 */ 206 PetscCall(PetscObjectDereference((PetscObject)dm)); 207 } 208 PetscCall(MatFDColoringApply(B,fdcoloring,X,snes)); 209 } 210 /* This will be redundant if the user called both, but it's too common to forget. */ 211 if (A != B) { 212 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 213 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 214 } 215 PetscFunctionReturn(0); 216 } 217 218 /*@C 219 DMDASNESSetFunctionLocal - set a local residual evaluation function 220 221 Logically Collective 222 223 Input Parameters: 224 + dm - DM to associate callback with 225 . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part 226 . func - local residual evaluation 227 - ctx - optional context for local residual evaluation 228 229 Calling sequence: 230 For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx), 231 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 232 . x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x) 233 . f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f) 234 - ctx - optional context passed above 235 236 Level: beginner 237 238 .seealso: `DMDASNESSetJacobianLocal()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 239 @*/ 240 PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx) 241 { 242 DMSNES sdm; 243 DMSNES_DA *dmdasnes; 244 245 PetscFunctionBegin; 246 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 247 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 248 PetscCall(DMDASNESGetContext(dm,sdm,&dmdasnes)); 249 250 dmdasnes->residuallocalimode = imode; 251 dmdasnes->residuallocal = func; 252 dmdasnes->residuallocalctx = ctx; 253 254 PetscCall(DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes)); 255 if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 256 PetscCall(DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes)); 257 } 258 PetscFunctionReturn(0); 259 } 260 261 /*@C 262 DMDASNESSetFunctionLocalVec - set a local residual evaluation function that operates on a local vector 263 264 Logically Collective 265 266 Input Parameters: 267 + dm - DM to associate callback with 268 . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part 269 . func - local residual evaluation 270 - ctx - optional context for local residual evaluation 271 272 Calling sequence: 273 For PetscErrorCode (*func)(DMDALocalInfo *info,Vec x, Vec f, void *ctx), 274 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 275 . x - state vector at which to evaluate residual 276 . f - residual vector 277 - ctx - optional context passed above 278 279 Level: beginner 280 281 .seealso: `DMDASNESSetFunctionLocal()`, `DMDASNESSetJacobianLocalVec()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 282 @*/ 283 PetscErrorCode DMDASNESSetFunctionLocalVec(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,Vec,Vec,void*),void *ctx) 284 { 285 DMSNES sdm; 286 DMSNES_DA *dmdasnes; 287 288 PetscFunctionBegin; 289 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 290 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 291 PetscCall(DMDASNESGetContext(dm,sdm,&dmdasnes)); 292 293 dmdasnes->residuallocalimode = imode; 294 dmdasnes->residuallocalvec = func; 295 dmdasnes->residuallocalctx = ctx; 296 297 PetscCall(DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes)); 298 if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 299 PetscCall(DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes)); 300 } 301 PetscFunctionReturn(0); 302 } 303 304 /*@C 305 DMDASNESSetJacobianLocal - set a local Jacobian evaluation function 306 307 Logically Collective 308 309 Input Parameters: 310 + dm - DM to associate callback with 311 . func - local Jacobian evaluation 312 - ctx - optional context for local Jacobian evaluation 313 314 Calling sequence: 315 For PetscErrorCode (*func)(DMDALocalInfo *info,void *x,Mat J,Mat M,void *ctx), 316 + info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at 317 . x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x) 318 . J - Mat object for the Jacobian 319 . M - Mat object for the Jacobian preconditioner matrix 320 - ctx - optional context passed above 321 322 Level: beginner 323 324 .seealso: `DMDASNESSetFunctionLocal()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 325 @*/ 326 PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx) 327 { 328 DMSNES sdm; 329 DMSNES_DA *dmdasnes; 330 331 PetscFunctionBegin; 332 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 333 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 334 PetscCall(DMDASNESGetContext(dm,sdm,&dmdasnes)); 335 336 dmdasnes->jacobianlocal = func; 337 dmdasnes->jacobianlocalctx = ctx; 338 339 PetscCall(DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes)); 340 PetscFunctionReturn(0); 341 } 342 343 /*@C 344 DMDASNESSetJacobianLocalVec - set a local Jacobian evaluation function that operates on a local vector 345 346 Logically Collective 347 348 Input Parameters: 349 + dm - DM to associate callback with 350 . func - local Jacobian evaluation 351 - ctx - optional context for local Jacobian evaluation 352 353 Calling sequence: 354 For PetscErrorCode (*func)(DMDALocalInfo *info,Vec x,Mat J,Mat M,void *ctx), 355 + info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at 356 . x - state vector at which to evaluate Jacobian 357 . J - Mat object for the Jacobian 358 . M - Mat object for the Jacobian preconditioner matrix 359 - ctx - optional context passed above 360 361 Level: beginner 362 363 .seealso: `DMDASNESSetJacobianLocal()`, `DMDASNESSetFunctionLocalVec()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 364 @*/ 365 PetscErrorCode DMDASNESSetJacobianLocalVec(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,Vec,Mat,Mat,void*),void *ctx) 366 { 367 DMSNES sdm; 368 DMSNES_DA *dmdasnes; 369 370 PetscFunctionBegin; 371 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 372 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 373 PetscCall(DMDASNESGetContext(dm,sdm,&dmdasnes)); 374 375 dmdasnes->jacobianlocalvec = func; 376 dmdasnes->jacobianlocalctx = ctx; 377 378 PetscCall(DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes)); 379 PetscFunctionReturn(0); 380 } 381 382 /*@C 383 DMDASNESSetObjectiveLocal - set a local residual evaluation function 384 385 Logically Collective 386 387 Input Parameters: 388 + dm - DM to associate callback with 389 . func - local objective evaluation 390 - ctx - optional context for local residual evaluation 391 392 Calling sequence for func: 393 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 394 . x - dimensional pointer to state at which to evaluate residual 395 . ob - eventual objective value 396 - ctx - optional context passed above 397 398 Level: beginner 399 400 .seealso: `DMSNESSetFunction()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 401 @*/ 402 PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx) 403 { 404 DMSNES sdm; 405 DMSNES_DA *dmdasnes; 406 407 PetscFunctionBegin; 408 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 409 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 410 PetscCall(DMDASNESGetContext(dm,sdm,&dmdasnes)); 411 412 dmdasnes->objectivelocal = func; 413 dmdasnes->objectivelocalctx = ctx; 414 415 PetscCall(DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes)); 416 PetscFunctionReturn(0); 417 } 418 419 /*@C 420 DMDASNESSetObjectiveLocal - set a local residual evaluation function that operates on a local vector 421 422 Logically Collective 423 424 Input Parameters: 425 + dm - DM to associate callback with 426 . func - local objective evaluation 427 - ctx - optional context for local residual evaluation 428 429 Calling sequence 430 For PetscErrorCode (*func)(DMDALocalInfo *info,Vec x,PetscReal *ob,void *ctx); 431 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 432 . x - state vector at which to evaluate residual 433 . ob - eventual objective value 434 - ctx - optional context passed above 435 436 Level: beginner 437 438 .seealso: `DMDASNESSetObjectiveLocal()`, `DMSNESSetFunction()`, `DMDASNESSetJacobianLocalVec()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 439 @*/ 440 PetscErrorCode DMDASNESSetObjectiveLocalVec(DM dm,DMDASNESObjectiveVec func,void *ctx) 441 { 442 DMSNES sdm; 443 DMSNES_DA *dmdasnes; 444 445 PetscFunctionBegin; 446 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 447 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 448 PetscCall(DMDASNESGetContext(dm,sdm,&dmdasnes)); 449 450 dmdasnes->objectivelocalvec = func; 451 dmdasnes->objectivelocalctx = ctx; 452 453 PetscCall(DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes)); 454 PetscFunctionReturn(0); 455 } 456 457 static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx) 458 { 459 DM dm; 460 DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx; 461 DMDALocalInfo info; 462 Vec Xloc; 463 void *x,*f; 464 465 PetscFunctionBegin; 466 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 467 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 468 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 469 PetscCheck(dmdasnes->rhsplocal,PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context"); 470 PetscCall(SNESGetDM(snes,&dm)); 471 PetscCall(DMGetLocalVector(dm,&Xloc)); 472 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 473 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 474 PetscCall(DMDAGetLocalInfo(dm,&info)); 475 PetscCall(DMDAVecGetArray(dm,Xloc,&x)); 476 switch (dmdasnes->residuallocalimode) { 477 case INSERT_VALUES: { 478 PetscCall(DMDAVecGetArray(dm,F,&f)); 479 CHKMEMQ; 480 PetscCall((*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx)); 481 CHKMEMQ; 482 PetscCall(DMDAVecRestoreArray(dm,F,&f)); 483 } break; 484 case ADD_VALUES: { 485 Vec Floc; 486 PetscCall(DMGetLocalVector(dm,&Floc)); 487 PetscCall(VecZeroEntries(Floc)); 488 PetscCall(DMDAVecGetArray(dm,Floc,&f)); 489 CHKMEMQ; 490 PetscCall((*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx)); 491 CHKMEMQ; 492 PetscCall(DMDAVecRestoreArray(dm,Floc,&f)); 493 PetscCall(VecZeroEntries(F)); 494 PetscCall(DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F)); 495 PetscCall(DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F)); 496 PetscCall(DMRestoreLocalVector(dm,&Floc)); 497 } break; 498 default: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode); 499 } 500 PetscCall(DMDAVecRestoreArray(dm,Xloc,&x)); 501 PetscCall(DMRestoreLocalVector(dm,&Xloc)); 502 PetscFunctionReturn(0); 503 } 504 505 static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx) 506 { 507 DM dm; 508 DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx; 509 DMDALocalInfo info; 510 Vec Xloc; 511 void *x; 512 513 PetscFunctionBegin; 514 PetscCheck(dmdasnes->jacobianplocal,PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context"); 515 PetscCall(SNESGetDM(snes,&dm)); 516 517 PetscCall(DMGetLocalVector(dm,&Xloc)); 518 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 519 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 520 PetscCall(DMDAGetLocalInfo(dm,&info)); 521 PetscCall(DMDAVecGetArray(dm,Xloc,&x)); 522 CHKMEMQ; 523 PetscCall((*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx)); 524 CHKMEMQ; 525 PetscCall(DMDAVecRestoreArray(dm,Xloc,&x)); 526 PetscCall(DMRestoreLocalVector(dm,&Xloc)); 527 if (A != B) { 528 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 529 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 530 } 531 PetscFunctionReturn(0); 532 } 533 534 /*@C 535 DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration 536 537 Logically Collective 538 539 Input Parameters: 540 + dm - DM to associate callback with 541 . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part 542 . func - local residual evaluation 543 - ctx - optional context for local residual evaluation 544 545 Calling sequence for func: 546 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 547 . x - dimensional pointer to state at which to evaluate residual 548 . f - dimensional pointer to residual, write the residual here 549 - ctx - optional context passed above 550 551 Notes: 552 The user must use 553 PetscCall(SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user)); 554 in their code before calling this routine. 555 556 Level: beginner 557 558 .seealso: `DMSNESSetFunction()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()` 559 @*/ 560 PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*), 561 PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx) 562 { 563 DMSNES sdm; 564 DMSNES_DA *dmdasnes; 565 566 PetscFunctionBegin; 567 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 568 PetscCall(DMGetDMSNESWrite(dm,&sdm)); 569 PetscCall(DMDASNESGetContext(dm,sdm,&dmdasnes)); 570 571 dmdasnes->residuallocalimode = imode; 572 dmdasnes->rhsplocal = func; 573 dmdasnes->jacobianplocal = jac; 574 dmdasnes->picardlocalctx = ctx; 575 576 PetscCall(DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes)); 577 PetscCall(DMSNESSetMFFunction(dm,SNESComputeFunction_DMDA,dmdasnes)); 578 PetscFunctionReturn(0); 579 } 580