Actual source code: dmdasnes.c
petsc-3.14.6 2021-03-30
1: #include <petscdmda.h>
2: #include <petsc/private/dmimpl.h>
3: #include <petsc/private/snesimpl.h>
5: /* This structure holds the user-provided DMDA callbacks */
6: typedef struct {
7: PetscErrorCode (*residuallocal)(DMDALocalInfo*,void*,void*,void*);
8: PetscErrorCode (*jacobianlocal)(DMDALocalInfo*,void*,Mat,Mat,void*);
9: PetscErrorCode (*objectivelocal)(DMDALocalInfo*,void*,PetscReal*,void*);
10: void *residuallocalctx;
11: void *jacobianlocalctx;
12: void *objectivelocalctx;
13: InsertMode residuallocalimode;
15: /* For Picard iteration defined locally */
16: PetscErrorCode (*rhsplocal)(DMDALocalInfo*,void*,void*,void*);
17: PetscErrorCode (*jacobianplocal)(DMDALocalInfo*,void*,Mat,Mat,void*);
18: void *picardlocalctx;
19: } DMSNES_DA;
21: static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
22: {
26: PetscFree(sdm->data);
27: return(0);
28: }
30: static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm,DMSNES sdm)
31: {
35: PetscNewLog(sdm,(DMSNES_DA**)&sdm->data);
36: if (oldsdm->data) {
37: PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_DA));
38: }
39: return(0);
40: }
43: static PetscErrorCode DMDASNESGetContext(DM dm,DMSNES sdm,DMSNES_DA **dmdasnes)
44: {
48: *dmdasnes = NULL;
49: if (!sdm->data) {
50: PetscNewLog(dm,(DMSNES_DA**)&sdm->data);
51: sdm->ops->destroy = DMSNESDestroy_DMDA;
52: sdm->ops->duplicate = DMSNESDuplicate_DMDA;
53: }
54: *dmdasnes = (DMSNES_DA*)sdm->data;
55: return(0);
56: }
58: static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx)
59: {
61: DM dm;
62: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
63: DMDALocalInfo info;
64: Vec Xloc;
65: void *x,*f;
71: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
72: SNESGetDM(snes,&dm);
73: DMGetLocalVector(dm,&Xloc);
74: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
75: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
76: DMDAGetLocalInfo(dm,&info);
77: DMDAVecGetArray(dm,Xloc,&x);
78: switch (dmdasnes->residuallocalimode) {
79: case INSERT_VALUES: {
80: DMDAVecGetArray(dm,F,&f);
81: PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
82: CHKMEMQ;
83: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
84: CHKMEMQ;
85: PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
86: DMDAVecRestoreArray(dm,F,&f);
87: } break;
88: case ADD_VALUES: {
89: Vec Floc;
90: DMGetLocalVector(dm,&Floc);
91: VecZeroEntries(Floc);
92: DMDAVecGetArray(dm,Floc,&f);
93: PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
94: CHKMEMQ;
95: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
96: CHKMEMQ;
97: PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
98: DMDAVecRestoreArray(dm,Floc,&f);
99: VecZeroEntries(F);
100: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
101: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
102: DMRestoreLocalVector(dm,&Floc);
103: } break;
104: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
105: }
106: DMDAVecRestoreArray(dm,Xloc,&x);
107: DMRestoreLocalVector(dm,&Xloc);
108: if (snes->domainerror) {
109: VecSetInf(F);
110: }
111: return(0);
112: }
114: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
115: {
117: DM dm;
118: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
119: DMDALocalInfo info;
120: Vec Xloc;
121: void *x;
127: if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
128: SNESGetDM(snes,&dm);
129: DMGetLocalVector(dm,&Xloc);
130: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
131: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
132: DMDAGetLocalInfo(dm,&info);
133: DMDAVecGetArray(dm,Xloc,&x);
134: CHKMEMQ;
135: (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);
136: CHKMEMQ;
137: DMDAVecRestoreArray(dm,Xloc,&x);
138: DMRestoreLocalVector(dm,&Xloc);
139: return(0);
140: }
143: /* Routine is called by example, hence must be labeled PETSC_EXTERN */
144: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
145: {
147: DM dm;
148: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
149: DMDALocalInfo info;
150: Vec Xloc;
151: void *x;
154: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
155: SNESGetDM(snes,&dm);
157: if (dmdasnes->jacobianlocal) {
158: DMGetLocalVector(dm,&Xloc);
159: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
160: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
161: DMDAGetLocalInfo(dm,&info);
162: DMDAVecGetArray(dm,Xloc,&x);
163: CHKMEMQ;
164: (*dmdasnes->jacobianlocal)(&info,x,A,B,dmdasnes->jacobianlocalctx);
165: CHKMEMQ;
166: DMDAVecRestoreArray(dm,Xloc,&x);
167: DMRestoreLocalVector(dm,&Xloc);
168: } else {
169: MatFDColoring fdcoloring;
170: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
171: if (!fdcoloring) {
172: ISColoring coloring;
174: DMCreateColoring(dm,dm->coloringtype,&coloring);
175: MatFDColoringCreate(B,coloring,&fdcoloring);
176: switch (dm->coloringtype) {
177: case IS_COLORING_GLOBAL:
178: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
179: break;
180: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
181: }
182: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
183: MatFDColoringSetFromOptions(fdcoloring);
184: MatFDColoringSetUp(B,coloring,fdcoloring);
185: ISColoringDestroy(&coloring);
186: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
187: PetscObjectDereference((PetscObject)fdcoloring);
189: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
190: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
191: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
192: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
193: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
194: */
195: PetscObjectDereference((PetscObject)dm);
196: }
197: MatFDColoringApply(B,fdcoloring,X,snes);
198: }
199: /* This will be redundant if the user called both, but it's too common to forget. */
200: if (A != B) {
201: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
202: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
203: }
204: return(0);
205: }
207: /*@C
208: DMDASNESSetFunctionLocal - set a local residual evaluation function
210: Logically Collective
212: Input Arguments:
213: + dm - DM to associate callback with
214: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
215: . func - local residual evaluation
216: - ctx - optional context for local residual evaluation
218: Calling sequence:
219: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx),
220: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
221: . x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
222: . f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
223: - ctx - optional context passed above
225: Level: beginner
227: .seealso: DMDASNESSetJacobianLocal(), DMSNESSetFunction(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
228: @*/
229: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
230: {
232: DMSNES sdm;
233: DMSNES_DA *dmdasnes;
237: DMGetDMSNESWrite(dm,&sdm);
238: DMDASNESGetContext(dm,sdm,&dmdasnes);
240: dmdasnes->residuallocalimode = imode;
241: dmdasnes->residuallocal = func;
242: dmdasnes->residuallocalctx = ctx;
244: DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
245: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
246: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
247: }
248: return(0);
249: }
251: /*@C
252: DMDASNESSetJacobianLocal - set a local Jacobian evaluation function
254: Logically Collective
256: Input Arguments:
257: + dm - DM to associate callback with
258: . func - local Jacobian evaluation
259: - ctx - optional context for local Jacobian evaluation
261: Calling sequence:
262: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x,Mat J,Mat M,void *ctx),
263: + info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at
264: . x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
265: . J - Mat object for the Jacobian
266: . M - Mat object for the Jacobian preconditioner matrix
267: - ctx - optional context passed above
269: Level: beginner
271: .seealso: DMDASNESSetFunctionLocal(), DMSNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
272: @*/
273: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
274: {
276: DMSNES sdm;
277: DMSNES_DA *dmdasnes;
281: DMGetDMSNESWrite(dm,&sdm);
282: DMDASNESGetContext(dm,sdm,&dmdasnes);
284: dmdasnes->jacobianlocal = func;
285: dmdasnes->jacobianlocalctx = ctx;
287: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
288: return(0);
289: }
292: /*@C
293: DMDASNESSetObjectiveLocal - set a local residual evaluation function
295: Logically Collective
297: Input Arguments:
298: + dm - DM to associate callback with
299: . func - local objective evaluation
300: - ctx - optional context for local residual evaluation
302: Calling sequence for func:
303: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
304: . x - dimensional pointer to state at which to evaluate residual
305: . ob - eventual objective value
306: - ctx - optional context passed above
308: Level: beginner
310: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
311: @*/
312: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
313: {
315: DMSNES sdm;
316: DMSNES_DA *dmdasnes;
320: DMGetDMSNESWrite(dm,&sdm);
321: DMDASNESGetContext(dm,sdm,&dmdasnes);
323: dmdasnes->objectivelocal = func;
324: dmdasnes->objectivelocalctx = ctx;
326: DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
327: return(0);
328: }
330: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
331: {
333: DM dm;
334: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
335: DMDALocalInfo info;
336: Vec Xloc;
337: void *x,*f;
343: if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
344: SNESGetDM(snes,&dm);
345: DMGetLocalVector(dm,&Xloc);
346: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
347: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
348: DMDAGetLocalInfo(dm,&info);
349: DMDAVecGetArray(dm,Xloc,&x);
350: switch (dmdasnes->residuallocalimode) {
351: case INSERT_VALUES: {
352: DMDAVecGetArray(dm,F,&f);
353: CHKMEMQ;
354: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
355: CHKMEMQ;
356: DMDAVecRestoreArray(dm,F,&f);
357: } break;
358: case ADD_VALUES: {
359: Vec Floc;
360: DMGetLocalVector(dm,&Floc);
361: VecZeroEntries(Floc);
362: DMDAVecGetArray(dm,Floc,&f);
363: CHKMEMQ;
364: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
365: CHKMEMQ;
366: DMDAVecRestoreArray(dm,Floc,&f);
367: VecZeroEntries(F);
368: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
369: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
370: DMRestoreLocalVector(dm,&Floc);
371: } break;
372: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
373: }
374: DMDAVecRestoreArray(dm,Xloc,&x);
375: DMRestoreLocalVector(dm,&Xloc);
376: return(0);
377: }
379: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
380: {
382: DM dm;
383: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
384: DMDALocalInfo info;
385: Vec Xloc;
386: void *x;
389: if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
390: SNESGetDM(snes,&dm);
392: DMGetLocalVector(dm,&Xloc);
393: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
394: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
395: DMDAGetLocalInfo(dm,&info);
396: DMDAVecGetArray(dm,Xloc,&x);
397: CHKMEMQ;
398: (*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx);
399: CHKMEMQ;
400: DMDAVecRestoreArray(dm,Xloc,&x);
401: DMRestoreLocalVector(dm,&Xloc);
402: if (A != B) {
403: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
404: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
405: }
406: return(0);
407: }
409: /*@C
410: DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
412: Logically Collective
414: Input Arguments:
415: + dm - DM to associate callback with
416: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
417: . func - local residual evaluation
418: - ctx - optional context for local residual evaluation
420: Calling sequence for func:
421: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
422: . x - dimensional pointer to state at which to evaluate residual
423: . f - dimensional pointer to residual, write the residual here
424: - ctx - optional context passed above
426: Notes:
427: The user must use
428: SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
429: in their code before calling this routine.
432: Level: beginner
434: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
435: @*/
436: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
437: PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
438: {
440: DMSNES sdm;
441: DMSNES_DA *dmdasnes;
445: DMGetDMSNESWrite(dm,&sdm);
446: DMDASNESGetContext(dm,sdm,&dmdasnes);
448: dmdasnes->residuallocalimode = imode;
449: dmdasnes->rhsplocal = func;
450: dmdasnes->jacobianplocal = jac;
451: dmdasnes->picardlocalctx = ctx;
453: DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
454: return(0);
455: }