Actual source code: dmredundant.c
1: #include <petsc/private/dmimpl.h>
2: #include <petscdmredundant.h>
4: typedef struct {
5: PetscMPIInt rank; /* owner */
6: PetscInt N; /* total number of dofs */
7: PetscInt n; /* owned number of dofs, n=N on owner, n=0 on non-owners */
8: } DM_Redundant;
10: static PetscErrorCode DMCreateMatrix_Redundant(DM dm, Mat *J)
11: {
12: DM_Redundant *red = (DM_Redundant *)dm->data;
13: ISLocalToGlobalMapping ltog;
14: PetscInt i, rstart, rend, *cols;
15: PetscScalar *vals;
17: PetscFunctionBegin;
18: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
19: PetscCall(MatSetSizes(*J, red->n, red->n, red->N, red->N));
20: PetscCall(MatSetType(*J, dm->mattype));
21: PetscCall(MatSeqAIJSetPreallocation(*J, red->n, NULL));
22: PetscCall(MatSeqBAIJSetPreallocation(*J, 1, red->n, NULL));
23: PetscCall(MatMPIAIJSetPreallocation(*J, red->n, NULL, red->N - red->n, NULL));
24: PetscCall(MatMPIBAIJSetPreallocation(*J, 1, red->n, NULL, red->N - red->n, NULL));
26: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
27: PetscCall(MatSetLocalToGlobalMapping(*J, ltog, ltog));
28: PetscCall(MatSetDM(*J, dm));
30: PetscCall(PetscMalloc2(red->N, &cols, red->N, &vals));
31: for (i = 0; i < red->N; i++) {
32: cols[i] = i;
33: vals[i] = 0.0;
34: }
35: PetscCall(MatGetOwnershipRange(*J, &rstart, &rend));
36: for (i = rstart; i < rend; i++) PetscCall(MatSetValues(*J, 1, &i, red->N, cols, vals, INSERT_VALUES));
37: PetscCall(PetscFree2(cols, vals));
38: PetscCall(MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY));
39: PetscCall(MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode DMDestroy_Redundant(DM dm)
44: {
45: PetscFunctionBegin;
46: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMRedundantSetSize_C", NULL));
47: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMRedundantGetSize_C", NULL));
48: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL));
49: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
50: PetscCall(PetscFree(dm->data));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode DMCreateGlobalVector_Redundant(DM dm, Vec *gvec)
55: {
56: DM_Redundant *red = (DM_Redundant *)dm->data;
57: ISLocalToGlobalMapping ltog;
59: PetscFunctionBegin;
61: PetscAssertPointer(gvec, 2);
62: *gvec = NULL;
63: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
64: PetscCall(VecSetSizes(*gvec, red->n, red->N));
65: PetscCall(VecSetType(*gvec, dm->vectype));
66: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
67: PetscCall(VecSetLocalToGlobalMapping(*gvec, ltog));
68: PetscCall(VecSetDM(*gvec, dm));
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: static PetscErrorCode DMCreateLocalVector_Redundant(DM dm, Vec *lvec)
73: {
74: DM_Redundant *red = (DM_Redundant *)dm->data;
76: PetscFunctionBegin;
78: PetscAssertPointer(lvec, 2);
79: *lvec = NULL;
80: PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
81: PetscCall(VecSetSizes(*lvec, red->N, red->N));
82: PetscCall(VecSetType(*lvec, dm->vectype));
83: PetscCall(VecSetDM(*lvec, dm));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode DMLocalToGlobalBegin_Redundant(DM dm, Vec l, InsertMode imode, Vec g)
88: {
89: DM_Redundant *red = (DM_Redundant *)dm->data;
90: const PetscScalar *lv;
91: PetscScalar *gv;
92: PetscMPIInt rank;
94: PetscFunctionBegin;
95: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
96: PetscCall(VecGetArrayRead(l, &lv));
97: PetscCall(VecGetArray(g, &gv));
98: switch (imode) {
99: case ADD_VALUES:
100: case MAX_VALUES: {
101: void *source;
102: PetscScalar *buffer;
103: PetscInt i;
104: if (rank == red->rank) {
105: buffer = gv;
106: source = MPI_IN_PLACE;
107: if (imode == ADD_VALUES)
108: for (i = 0; i < red->N; i++) buffer[i] = gv[i] + lv[i];
109: #if !defined(PETSC_USE_COMPLEX)
110: if (imode == MAX_VALUES)
111: for (i = 0; i < red->N; i++) buffer[i] = PetscMax(gv[i], lv[i]);
112: #endif
113: } else source = (void *)lv;
114: PetscCallMPI(MPI_Reduce(source, gv, red->N, MPIU_SCALAR, (imode == ADD_VALUES) ? MPIU_SUM : MPIU_MAX, red->rank, PetscObjectComm((PetscObject)dm)));
115: } break;
116: case INSERT_VALUES:
117: PetscCall(PetscArraycpy(gv, lv, red->n));
118: break;
119: default:
120: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "InsertMode not supported");
121: }
122: PetscCall(VecRestoreArrayRead(l, &lv));
123: PetscCall(VecRestoreArray(g, &gv));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode DMLocalToGlobalEnd_Redundant(DM dm, Vec l, InsertMode imode, Vec g)
128: {
129: PetscFunctionBegin;
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode DMGlobalToLocalBegin_Redundant(DM dm, Vec g, InsertMode imode, Vec l)
134: {
135: DM_Redundant *red = (DM_Redundant *)dm->data;
136: const PetscScalar *gv;
137: PetscScalar *lv;
139: PetscFunctionBegin;
140: PetscCall(VecGetArrayRead(g, &gv));
141: PetscCall(VecGetArray(l, &lv));
142: switch (imode) {
143: case INSERT_VALUES:
144: if (red->n) PetscCall(PetscArraycpy(lv, gv, red->n));
145: PetscCallMPI(MPI_Bcast(lv, red->N, MPIU_SCALAR, red->rank, PetscObjectComm((PetscObject)dm)));
146: break;
147: default:
148: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "InsertMode not supported");
149: }
150: PetscCall(VecRestoreArrayRead(g, &gv));
151: PetscCall(VecRestoreArray(l, &lv));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode DMGlobalToLocalEnd_Redundant(DM dm, Vec g, InsertMode imode, Vec l)
156: {
157: PetscFunctionBegin;
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: static PetscErrorCode DMSetUp_Redundant(DM dm)
162: {
163: PetscFunctionBegin;
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode DMView_Redundant(DM dm, PetscViewer viewer)
168: {
169: DM_Redundant *red = (DM_Redundant *)dm->data;
170: PetscBool iascii;
172: PetscFunctionBegin;
173: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
174: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "redundant: rank=%d N=%" PetscInt_FMT "\n", red->rank, red->N));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode DMCreateColoring_Redundant(DM dm, ISColoringType ctype, ISColoring *coloring)
179: {
180: DM_Redundant *red = (DM_Redundant *)dm->data;
181: PetscInt i, nloc;
182: ISColoringValue *colors;
184: PetscFunctionBegin;
185: switch (ctype) {
186: case IS_COLORING_GLOBAL:
187: nloc = red->n;
188: break;
189: case IS_COLORING_LOCAL:
190: nloc = red->N;
191: break;
192: default:
193: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown ISColoringType %d", (int)ctype);
194: }
195: PetscCall(PetscMalloc1(nloc, &colors));
196: for (i = 0; i < nloc; i++) colors[i] = i;
197: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), red->N, nloc, colors, PETSC_OWN_POINTER, coloring));
198: PetscCall(ISColoringSetType(*coloring, ctype));
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: static PetscErrorCode DMRefine_Redundant(DM dmc, MPI_Comm comm, DM *dmf)
203: {
204: PetscMPIInt flag;
205: DM_Redundant *redc = (DM_Redundant *)dmc->data;
207: PetscFunctionBegin;
208: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmc, &comm));
209: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)dmc), comm, &flag));
210: PetscCheck(flag == MPI_CONGRUENT || flag == MPI_IDENT, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "cannot change communicators");
211: PetscCall(DMRedundantCreate(comm, redc->rank, redc->N, dmf));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: static PetscErrorCode DMCoarsen_Redundant(DM dmf, MPI_Comm comm, DM *dmc)
216: {
217: PetscMPIInt flag;
218: DM_Redundant *redf = (DM_Redundant *)dmf->data;
220: PetscFunctionBegin;
221: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmf, &comm));
222: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)dmf), comm, &flag));
223: PetscCheck(flag == MPI_CONGRUENT || flag == MPI_IDENT, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "cannot change communicators");
224: PetscCall(DMRedundantCreate(comm, redf->rank, redf->N, dmc));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode DMCreateInterpolation_Redundant(DM dmc, DM dmf, Mat *P, Vec *scale)
229: {
230: DM_Redundant *redc = (DM_Redundant *)dmc->data;
231: DM_Redundant *redf = (DM_Redundant *)dmf->data;
232: PetscMPIInt flag;
233: PetscInt i, rstart, rend;
235: PetscFunctionBegin;
236: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)dmc), PetscObjectComm((PetscObject)dmf), &flag));
237: PetscCheck(flag == MPI_CONGRUENT || flag == MPI_IDENT, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "cannot change communicators");
238: PetscCheck(redc->rank == redf->rank, PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_INCOMP, "Owning rank does not match");
239: PetscCheck(redc->N == redf->N, PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_INCOMP, "Global size does not match");
240: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), P));
241: PetscCall(MatSetSizes(*P, redc->n, redc->n, redc->N, redc->N));
242: PetscCall(MatSetType(*P, MATAIJ));
243: PetscCall(MatSeqAIJSetPreallocation(*P, 1, NULL));
244: PetscCall(MatMPIAIJSetPreallocation(*P, 1, NULL, 0, NULL));
245: PetscCall(MatGetOwnershipRange(*P, &rstart, &rend));
246: for (i = rstart; i < rend; i++) PetscCall(MatSetValue(*P, i, i, 1.0, INSERT_VALUES));
247: PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
248: PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
249: if (scale) PetscCall(DMCreateInterpolationScale(dmc, dmf, *P, scale));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: /*@
254: DMRedundantSetSize - Sets the size of a densely coupled redundant object
256: Collective
258: Input Parameters:
259: + dm - `DM` object of type `DMREDUNDANT`
260: . rank - rank of process to own the redundant degrees of freedom
261: - N - total number of redundant degrees of freedom
263: Level: advanced
265: .seealso: `DM`, `DMREDUNDANT`, `DMDestroy()`, `DMCreateGlobalVector()`, `DMRedundantCreate()`, `DMRedundantGetSize()`
266: @*/
267: PetscErrorCode DMRedundantSetSize(DM dm, PetscMPIInt rank, PetscInt N)
268: {
269: PetscFunctionBegin;
274: PetscTryMethod(dm, "DMRedundantSetSize_C", (DM, PetscMPIInt, PetscInt), (dm, rank, N));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: /*@
279: DMRedundantGetSize - Gets the size of a densely coupled redundant object
281: Not Collective
283: Input Parameter:
284: . dm - `DM` object of type `DMREDUNDANT`
286: Output Parameters:
287: + rank - rank of process to own the redundant degrees of freedom (or `NULL`)
288: - N - total number of redundant degrees of freedom (or `NULL`)
290: Level: advanced
292: .seealso: `DM`, `DMREDUNDANT`, `DMDestroy()`, `DMCreateGlobalVector()`, `DMRedundantCreate()`, `DMRedundantSetSize()`
293: @*/
294: PetscErrorCode DMRedundantGetSize(DM dm, PetscMPIInt *rank, PetscInt *N)
295: {
296: PetscFunctionBegin;
299: PetscUseMethod(dm, "DMRedundantGetSize_C", (DM, PetscMPIInt *, PetscInt *), (dm, rank, N));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode DMRedundantSetSize_Redundant(DM dm, PetscMPIInt rank, PetscInt N)
304: {
305: DM_Redundant *red = (DM_Redundant *)dm->data;
306: PetscMPIInt myrank;
307: PetscInt i, *globals;
309: PetscFunctionBegin;
310: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &myrank));
311: red->rank = rank;
312: red->N = N;
313: red->n = (myrank == rank) ? N : 0;
315: /* mapping is setup here */
316: PetscCall(PetscMalloc1(red->N, &globals));
317: for (i = 0; i < red->N; i++) globals[i] = i;
318: PetscCall(ISLocalToGlobalMappingDestroy(&dm->ltogmap));
319: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, red->N, globals, PETSC_OWN_POINTER, &dm->ltogmap));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: static PetscErrorCode DMRedundantGetSize_Redundant(DM dm, PetscInt *rank, PetscInt *N)
324: {
325: DM_Redundant *red = (DM_Redundant *)dm->data;
327: PetscFunctionBegin;
328: if (rank) *rank = red->rank;
329: if (N) *N = red->N;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode DMSetUpGLVisViewer_Redundant(PetscObject odm, PetscViewer viewer)
334: {
335: PetscFunctionBegin;
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: /*MC
340: DMREDUNDANT = "redundant" - A DM object that is used to manage data for a small set of dense globally coupled variables.
341: In the global representation of the vector the variables are all stored on a single MPI process (all the other MPI processes have
342: no variables) in the local representation all the variables are stored on ALL the MPI processes (because they are all needed for each
343: processes local computations).
345: This DM is generally used inside a DMCOMPOSITE object. For example, it may be used to store continuation parameters for a bifurcation problem.
347: Level: intermediate
349: .seealso: `DMType`, `DMCOMPOSITE`, `DMCreate()`, `DMRedundantSetSize()`, `DMRedundantGetSize()`
350: M*/
352: PETSC_EXTERN PetscErrorCode DMCreate_Redundant(DM dm)
353: {
354: DM_Redundant *red;
356: PetscFunctionBegin;
357: PetscCall(PetscNew(&red));
358: dm->data = red;
360: dm->ops->setup = DMSetUp_Redundant;
361: dm->ops->view = DMView_Redundant;
362: dm->ops->createglobalvector = DMCreateGlobalVector_Redundant;
363: dm->ops->createlocalvector = DMCreateLocalVector_Redundant;
364: dm->ops->creatematrix = DMCreateMatrix_Redundant;
365: dm->ops->destroy = DMDestroy_Redundant;
366: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Redundant;
367: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Redundant;
368: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Redundant;
369: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Redundant;
370: dm->ops->refine = DMRefine_Redundant;
371: dm->ops->coarsen = DMCoarsen_Redundant;
372: dm->ops->createinterpolation = DMCreateInterpolation_Redundant;
373: dm->ops->getcoloring = DMCreateColoring_Redundant;
375: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMRedundantSetSize_C", DMRedundantSetSize_Redundant));
376: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMRedundantGetSize_C", DMRedundantGetSize_Redundant));
377: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Redundant));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*@C
382: DMRedundantCreate - Creates a `DM` object, used to manage data for dense globally coupled variables
384: Collective
386: Input Parameters:
387: + comm - the processors that will share the global vector
388: . rank - the MPI rank to own the redundant values
389: - N - total number of degrees of freedom
391: Output Parameter:
392: . dm - the `DM` object of type `DMREDUNDANT`
394: Level: advanced
396: .seealso: `DM`, `DMREDUNDANT`, `DMDestroy()`, `DMCreateGlobalVector()`, `DMCreateMatrix()`, `DMCompositeAddDM()`, `DMSetType()`, `DMRedundantSetSize()`, `DMRedundantGetSize()`
397: @*/
398: PetscErrorCode DMRedundantCreate(MPI_Comm comm, PetscMPIInt rank, PetscInt N, DM *dm)
399: {
400: PetscFunctionBegin;
401: PetscAssertPointer(dm, 4);
402: PetscCall(DMCreate(comm, dm));
403: PetscCall(DMSetType(*dm, DMREDUNDANT));
404: PetscCall(DMRedundantSetSize(*dm, rank, N));
405: PetscCall(DMSetUp(*dm));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }