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