Actual source code: dmlocalsnes.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/dmimpl.h>
  2:  #include <petsc/private/snesimpl.h>

  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;

 13: static PetscErrorCode DMSNESDestroy_DMLocal(DMSNES sdm)
 14: {

 18:   PetscFree(sdm->data);
 19:   sdm->data = NULL;
 20:   return(0);
 21: }

 23: static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm,DMSNES sdm)
 24: {

 28:   if (sdm->data != oldsdm->data) {
 29:     PetscFree(sdm->data);
 30:     PetscNewLog(sdm,(DMSNES_Local**)&sdm->data);
 31:     if (oldsdm->data) {PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_Local));}
 32:   }
 33:   return(0);
 34: }

 36: static PetscErrorCode DMLocalSNESGetContext(DM dm,DMSNES sdm,DMSNES_Local **dmlocalsnes)
 37: {

 41:   *dmlocalsnes = NULL;
 42:   if (!sdm->data) {
 43:     PetscNewLog(dm,(DMSNES_Local**)&sdm->data);

 45:     sdm->ops->destroy   = DMSNESDestroy_DMLocal;
 46:     sdm->ops->duplicate = DMSNESDuplicate_DMLocal;
 47:   }
 48:   *dmlocalsnes = (DMSNES_Local*)sdm->data;
 49:   return(0);
 50: }

 52: static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes,Vec X,Vec F,void *ctx)
 53: {
 55:   DM             dm;
 56:   DMSNES_Local   *dmlocalsnes = (DMSNES_Local*)ctx;
 57:   Vec            Xloc,Floc;

 63:   SNESGetDM(snes,&dm);
 64:   DMGetLocalVector(dm,&Xloc);
 65:   DMGetLocalVector(dm,&Floc);
 66:   VecZeroEntries(Xloc);
 67:   if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
 68:   VecZeroEntries(Floc);
 69:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 70:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 71:   CHKMEMQ;
 72:   (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);
 73:   CHKMEMQ;
 74:   VecZeroEntries(F);
 75:   DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
 76:   DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
 77:   DMRestoreLocalVector(dm,&Floc);
 78:   DMRestoreLocalVector(dm,&Xloc);
 79:   {
 80:     char        name[PETSC_MAX_PATH_LEN];
 81:     char        oldname[PETSC_MAX_PATH_LEN];
 82:     const char *tmp;
 83:     PetscInt    it;

 85:     SNESGetIterationNumber(snes, &it);
 86:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Solution, Iterate %d", (int) it);
 87:     PetscObjectGetName((PetscObject) X, &tmp);
 88:     PetscStrncpy(oldname, tmp, PETSC_MAX_PATH_LEN-1);
 89:     PetscObjectSetName((PetscObject) X, name);
 90:     VecViewFromOptions(X, (PetscObject) snes, "-dmsnes_solution_vec_view");
 91:     PetscObjectSetName((PetscObject) X, oldname);
 92:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Residual, Iterate %d", (int) it);
 93:     PetscObjectSetName((PetscObject) F, name);
 94:     VecViewFromOptions(F, (PetscObject) snes, "-dmsnes_residual_vec_view");
 95:   }
 96:   return(0);
 97: }

 99: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx)
100: {
102:   DM             dm;
103:   DMSNES_Local   *dmlocalsnes = (DMSNES_Local*)ctx;
104:   Vec            Xloc;

107:   SNESGetDM(snes,&dm);
108:   if (dmlocalsnes->jacobianlocal) {
109:     DMGetLocalVector(dm,&Xloc);
110:     VecZeroEntries(Xloc);
111:     if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
112:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
113:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
114:     CHKMEMQ;
115:     (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);
116:     CHKMEMQ;
117:     DMRestoreLocalVector(dm,&Xloc);
118:   } else {
119:     MatFDColoring fdcoloring;
120:     PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
121:     if (!fdcoloring) {
122:       ISColoring coloring;

124:       DMCreateColoring(dm,dm->coloringtype,&coloring);
125:       MatFDColoringCreate(B,coloring,&fdcoloring);
126:       ISColoringDestroy(&coloring);
127:       switch (dm->coloringtype) {
128:       case IS_COLORING_GLOBAL:
129:         MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);
130:         break;
131:       default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
132:       }
133:       PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
134:       MatFDColoringSetFromOptions(fdcoloring);
135:       MatFDColoringSetUp(B,coloring,fdcoloring);
136:       PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
137:       PetscObjectDereference((PetscObject)fdcoloring);

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:       PetscObjectDereference((PetscObject)dm);
146:     }
147:     MatFDColoringApply(B,fdcoloring,X,snes);
148:   }
149:   /* This will be redundant if the user called both, but it's too common to forget. */
150:   if (A != B) {
151:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
152:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
153:   }
154:   return(0);
155: }

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.

162:    Logically Collective

164:    Input Arguments:
165: +  dm - DM to associate callback with
166: .  func - local residual evaluation
167: -  ctx - optional context for local residual evaluation

169:    Level: beginner

171: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
172: @*/
173: PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx)
174: {
176:   DMSNES         sdm;
177:   DMSNES_Local   *dmlocalsnes;

181:   DMGetDMSNESWrite(dm,&sdm);
182:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

184:   dmlocalsnes->residuallocal    = func;
185:   dmlocalsnes->residuallocalctx = ctx;

187:   DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes);
188:   if (!sdm->ops->computejacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
189:     DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
190:   }
191:   return(0);
192: }

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.

199:    Logically Collective

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

206:    Level: intermediate

208: .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal()
209: @*/
210: PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx)
211: {
213:   DMSNES         sdm;
214:   DMSNES_Local   *dmlocalsnes;

218:   DMGetDMSNESWrite(dm,&sdm);
219:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

221:   dmlocalsnes->boundarylocal    = func;
222:   dmlocalsnes->boundarylocalctx = ctx;

224:   return(0);
225: }

227: /*@C
228:    DMSNESSetJacobianLocal - set a local Jacobian evaluation function

230:    Logically Collective

232:    Input Arguments:
233: +  dm - DM to associate callback with
234: .  func - local Jacobian evaluation
235: -  ctx - optional context for local Jacobian evaluation

237:    Level: beginner

239: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
240: @*/
241: PetscErrorCode DMSNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Mat,Mat,void*),void *ctx)
242: {
244:   DMSNES         sdm;
245:   DMSNES_Local   *dmlocalsnes;

249:   DMGetDMSNESWrite(dm,&sdm);
250:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

252:   dmlocalsnes->jacobianlocal    = func;
253:   dmlocalsnes->jacobianlocalctx = ctx;

255:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
256:   return(0);
257: }