Actual source code: dmlocalsnes.c

petsc-master 2016-08-28
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;

 15: static PetscErrorCode DMSNESDestroy_DMLocal(DMSNES sdm)
 16: {

 20:   PetscFree(sdm->data);
 21:   return(0);
 22: }

 26: static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm,DMSNES sdm)
 27: {

 31:   PetscNewLog(sdm,(DMSNES_Local**)&sdm->data);
 32:   if (oldsdm->data) {
 33:     PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_Local));
 34:   }
 35:   return(0);
 36: }

 40: static PetscErrorCode DMLocalSNESGetContext(DM dm,DMSNES sdm,DMSNES_Local **dmlocalsnes)
 41: {

 45:   *dmlocalsnes = NULL;
 46:   if (!sdm->data) {
 47:     PetscNewLog(dm,(DMSNES_Local**)&sdm->data);

 49:     sdm->ops->destroy   = DMSNESDestroy_DMLocal;
 50:     sdm->ops->duplicate = DMSNESDuplicate_DMLocal;
 51:   }
 52:   *dmlocalsnes = (DMSNES_Local*)sdm->data;
 53:   return(0);
 54: }

 58: static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes,Vec X,Vec F,void *ctx)
 59: {
 61:   DM             dm;
 62:   DMSNES_Local   *dmlocalsnes = (DMSNES_Local*)ctx;
 63:   Vec            Xloc,Floc;

 69:   SNESGetDM(snes,&dm);
 70:   DMGetLocalVector(dm,&Xloc);
 71:   DMGetLocalVector(dm,&Floc);
 72:   VecZeroEntries(Xloc);
 73:   if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
 74:   VecZeroEntries(Floc);
 75:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 76:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 77:   CHKMEMQ;
 78:   (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);
 79:   CHKMEMQ;
 80:   VecZeroEntries(F);
 81:   DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
 82:   DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
 83:   DMRestoreLocalVector(dm,&Floc);
 84:   DMRestoreLocalVector(dm,&Xloc);
 85:   return(0);
 86: }

 90: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx)
 91: {
 93:   DM             dm;
 94:   DMSNES_Local   *dmlocalsnes = (DMSNES_Local*)ctx;
 95:   Vec            Xloc;

 98:   SNESGetDM(snes,&dm);
 99:   if (dmlocalsnes->jacobianlocal) {
100:     DMGetLocalVector(dm,&Xloc);
101:     VecZeroEntries(Xloc);
102:     if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
103:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
104:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
105:     CHKMEMQ;
106:     (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);
107:     CHKMEMQ;
108:     DMRestoreLocalVector(dm,&Xloc);
109:   } else {
110:     MatFDColoring fdcoloring;
111:     PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
112:     if (!fdcoloring) {
113:       ISColoring coloring;

115:       DMCreateColoring(dm,dm->coloringtype,&coloring);
116:       MatFDColoringCreate(B,coloring,&fdcoloring);
117:       ISColoringDestroy(&coloring);
118:       switch (dm->coloringtype) {
119:       case IS_COLORING_GLOBAL:
120:         MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);
121:         break;
122:       default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
123:       }
124:       PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
125:       MatFDColoringSetFromOptions(fdcoloring);
126:       MatFDColoringSetUp(B,coloring,fdcoloring);
127:       PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
128:       PetscObjectDereference((PetscObject)fdcoloring);

130:       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
131:        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
132:        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
133:        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
134:        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
135:        */
136:       PetscObjectDereference((PetscObject)dm);
137:     }
138:     MatFDColoringApply(B,fdcoloring,X,snes);
139:   }
140:   /* This will be redundant if the user called both, but it's too common to forget. */
141:   if (A != B) {
142:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
143:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
144:   }
145:   return(0);
146: }

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.

155:    Logically Collective

157:    Input Arguments:
158: +  dm - DM to associate callback with
159: .  func - local residual evaluation
160: -  ctx - optional context for local residual evaluation

162:    Level: beginner

164: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
165: @*/
166: PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx)
167: {
169:   DMSNES         sdm;
170:   DMSNES_Local   *dmlocalsnes;

174:   DMGetDMSNESWrite(dm,&sdm);
175:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

177:   dmlocalsnes->residuallocal    = func;
178:   dmlocalsnes->residuallocalctx = ctx;

180:   DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes);
181:   if (!sdm->ops->computejacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
182:     DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
183:   }
184:   return(0);
185: }

189: /*@C
190:    DMSNESSetBoundaryLocal - set a local boundary value function. This function is called with local vector
191:       containing the local vector information PLUS ghost point information. It should insert values into the local
192:       vector that do not come from the global vector, such as essential boundary condition data.

194:    Logically Collective

196:    Input Arguments:
197: +  dm - DM to associate callback with
198: .  func - local boundary value evaluation
199: -  ctx - optional context for local boundary value evaluation

201:    Level: intermediate

203: .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal()
204: @*/
205: PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx)
206: {
208:   DMSNES         sdm;
209:   DMSNES_Local   *dmlocalsnes;

213:   DMGetDMSNESWrite(dm,&sdm);
214:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

216:   dmlocalsnes->boundarylocal    = func;
217:   dmlocalsnes->boundarylocalctx = ctx;

219:   return(0);
220: }

224: /*@C
225:    DMSNESSetJacobianLocal - set a local Jacobian evaluation function

227:    Logically Collective

229:    Input Arguments:
230: +  dm - DM to associate callback with
231: .  func - local Jacobian evaluation
232: -  ctx - optional context for local Jacobian evaluation

234:    Level: beginner

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

246:   DMGetDMSNESWrite(dm,&sdm);
247:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

249:   dmlocalsnes->jacobianlocal    = func;
250:   dmlocalsnes->jacobianlocalctx = ctx;

252:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
253:   return(0);
254: }