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