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