Actual source code: dmlocalsnes.c

petsc-master 2017-05-22
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:   return(0);
 20: }

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

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

 34: static PetscErrorCode DMLocalSNESGetContext(DM dm,DMSNES sdm,DMSNES_Local **dmlocalsnes)
 35: {

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

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

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

 61:   SNESGetDM(snes,&dm);
 62:   DMGetLocalVector(dm,&Xloc);
 63:   DMGetLocalVector(dm,&Floc);
 64:   VecZeroEntries(Xloc);
 65:   if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
 66:   VecZeroEntries(Floc);
 67:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 68:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 69:   CHKMEMQ;
 70:   (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);
 71:   CHKMEMQ;
 72:   VecZeroEntries(F);
 73:   DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
 74:   DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
 75:   DMRestoreLocalVector(dm,&Floc);
 76:   DMRestoreLocalVector(dm,&Xloc);
 77:   return(0);
 78: }

 80: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx)
 81: {
 83:   DM             dm;
 84:   DMSNES_Local   *dmlocalsnes = (DMSNES_Local*)ctx;
 85:   Vec            Xloc;

 88:   SNESGetDM(snes,&dm);
 89:   if (dmlocalsnes->jacobianlocal) {
 90:     DMGetLocalVector(dm,&Xloc);
 91:     VecZeroEntries(Xloc);
 92:     if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
 93:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 94:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 95:     CHKMEMQ;
 96:     (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);
 97:     CHKMEMQ;
 98:     DMRestoreLocalVector(dm,&Xloc);
 99:   } else {
100:     MatFDColoring fdcoloring;
101:     PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
102:     if (!fdcoloring) {
103:       ISColoring coloring;

105:       DMCreateColoring(dm,dm->coloringtype,&coloring);
106:       MatFDColoringCreate(B,coloring,&fdcoloring);
107:       ISColoringDestroy(&coloring);
108:       switch (dm->coloringtype) {
109:       case IS_COLORING_GLOBAL:
110:         MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);
111:         break;
112:       default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
113:       }
114:       PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
115:       MatFDColoringSetFromOptions(fdcoloring);
116:       MatFDColoringSetUp(B,coloring,fdcoloring);
117:       PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
118:       PetscObjectDereference((PetscObject)fdcoloring);

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

138: /*@C
139:    DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector
140:       containing the local vector information PLUS ghost point information. It should compute a result for all local
141:       elements and DMSNES will automatically accumulate the overlapping values.

143:    Logically Collective

145:    Input Arguments:
146: +  dm - DM to associate callback with
147: .  func - local residual evaluation
148: -  ctx - optional context for local residual evaluation

150:    Level: beginner

152: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
153: @*/
154: PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx)
155: {
157:   DMSNES         sdm;
158:   DMSNES_Local   *dmlocalsnes;

162:   DMGetDMSNESWrite(dm,&sdm);
163:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

165:   dmlocalsnes->residuallocal    = func;
166:   dmlocalsnes->residuallocalctx = ctx;

168:   DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes);
169:   if (!sdm->ops->computejacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
170:     DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
171:   }
172:   return(0);
173: }

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

180:    Logically Collective

182:    Input Arguments:
183: +  dm - DM to associate callback with
184: .  func - local boundary value evaluation
185: -  ctx - optional context for local boundary value evaluation

187:    Level: intermediate

189: .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal()
190: @*/
191: PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx)
192: {
194:   DMSNES         sdm;
195:   DMSNES_Local   *dmlocalsnes;

199:   DMGetDMSNESWrite(dm,&sdm);
200:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

202:   dmlocalsnes->boundarylocal    = func;
203:   dmlocalsnes->boundarylocalctx = ctx;

205:   return(0);
206: }

208: /*@C
209:    DMSNESSetJacobianLocal - set a local Jacobian evaluation function

211:    Logically Collective

213:    Input Arguments:
214: +  dm - DM to associate callback with
215: .  func - local Jacobian evaluation
216: -  ctx - optional context for local Jacobian evaluation

218:    Level: beginner

220: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
221: @*/
222: PetscErrorCode DMSNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Mat,Mat,void*),void *ctx)
223: {
225:   DMSNES         sdm;
226:   DMSNES_Local   *dmlocalsnes;

230:   DMGetDMSNESWrite(dm,&sdm);
231:   DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);

233:   dmlocalsnes->jacobianlocal    = func;
234:   dmlocalsnes->jacobianlocalctx = ctx;

236:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
237:   return(0);
238: }