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