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