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, iN;
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: PetscCall(PetscMPIIntCast(red->N, &iN));
115: PetscCallMPI(MPI_Reduce(source, gv, iN, MPIU_SCALAR, (imode == ADD_VALUES) ? MPIU_SUM : MPIU_MAX, red->rank, PetscObjectComm((PetscObject)dm)));
116: } break;
117: case INSERT_VALUES:
118: PetscCall(PetscArraycpy(gv, lv, red->n));
119: break;
120: default:
121: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "InsertMode not supported");
122: }
123: PetscCall(VecRestoreArrayRead(l, &lv));
124: PetscCall(VecRestoreArray(g, &gv));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: static PetscErrorCode DMLocalToGlobalEnd_Redundant(DM dm, Vec l, InsertMode imode, Vec g)
129: {
130: PetscFunctionBegin;
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: static PetscErrorCode DMGlobalToLocalBegin_Redundant(DM dm, Vec g, InsertMode imode, Vec l)
135: {
136: DM_Redundant *red = (DM_Redundant *)dm->data;
137: const PetscScalar *gv;
138: PetscScalar *lv;
139: PetscMPIInt iN;
141: PetscFunctionBegin;
142: PetscCall(VecGetArrayRead(g, &gv));
143: PetscCall(VecGetArray(l, &lv));
144: switch (imode) {
145: case INSERT_VALUES:
146: if (red->n) PetscCall(PetscArraycpy(lv, gv, red->n));
147: PetscCall(PetscMPIIntCast(red->N, &iN));
148: PetscCallMPI(MPI_Bcast(lv, iN, MPIU_SCALAR, red->rank, PetscObjectComm((PetscObject)dm)));
149: break;
150: default:
151: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "InsertMode not supported");
152: }
153: PetscCall(VecRestoreArrayRead(g, &gv));
154: PetscCall(VecRestoreArray(l, &lv));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: static PetscErrorCode DMGlobalToLocalEnd_Redundant(DM dm, Vec g, InsertMode imode, Vec l)
159: {
160: PetscFunctionBegin;
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: static PetscErrorCode DMSetUp_Redundant(DM dm)
165: {
166: PetscFunctionBegin;
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode DMView_Redundant(DM dm, PetscViewer viewer)
171: {
172: DM_Redundant *red = (DM_Redundant *)dm->data;
173: PetscBool iascii;
175: PetscFunctionBegin;
176: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
177: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "redundant: rank=%d N=%" PetscInt_FMT "\n", red->rank, red->N));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode DMCreateColoring_Redundant(DM dm, ISColoringType ctype, ISColoring *coloring)
182: {
183: DM_Redundant *red = (DM_Redundant *)dm->data;
184: PetscInt i, nloc;
185: ISColoringValue *colors;
187: PetscFunctionBegin;
188: switch (ctype) {
189: case IS_COLORING_GLOBAL:
190: nloc = red->n;
191: break;
192: case IS_COLORING_LOCAL:
193: nloc = red->N;
194: break;
195: default:
196: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown ISColoringType %d", (int)ctype);
197: }
198: PetscCall(PetscMalloc1(nloc, &colors));
199: for (i = 0; i < nloc; i++) PetscCall(ISColoringValueCast(i, colors + i));
200: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), red->N, nloc, colors, PETSC_OWN_POINTER, coloring));
201: PetscCall(ISColoringSetType(*coloring, ctype));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode DMRefine_Redundant(DM dmc, MPI_Comm comm, DM *dmf)
206: {
207: PetscMPIInt flag;
208: DM_Redundant *redc = (DM_Redundant *)dmc->data;
210: PetscFunctionBegin;
211: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmc, &comm));
212: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)dmc), comm, &flag));
213: PetscCheck(flag == MPI_CONGRUENT || flag == MPI_IDENT, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "cannot change communicators");
214: PetscCall(DMRedundantCreate(comm, redc->rank, redc->N, dmf));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: static PetscErrorCode DMCoarsen_Redundant(DM dmf, MPI_Comm comm, DM *dmc)
219: {
220: PetscMPIInt flag;
221: DM_Redundant *redf = (DM_Redundant *)dmf->data;
223: PetscFunctionBegin;
224: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmf, &comm));
225: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)dmf), comm, &flag));
226: PetscCheck(flag == MPI_CONGRUENT || flag == MPI_IDENT, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "cannot change communicators");
227: PetscCall(DMRedundantCreate(comm, redf->rank, redf->N, dmc));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: static PetscErrorCode DMCreateInterpolation_Redundant(DM dmc, DM dmf, Mat *P, Vec *scale)
232: {
233: DM_Redundant *redc = (DM_Redundant *)dmc->data;
234: DM_Redundant *redf = (DM_Redundant *)dmf->data;
235: PetscMPIInt flag;
236: PetscInt i, rstart, rend;
238: PetscFunctionBegin;
239: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)dmc), PetscObjectComm((PetscObject)dmf), &flag));
240: PetscCheck(flag == MPI_CONGRUENT || flag == MPI_IDENT, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "cannot change communicators");
241: PetscCheck(redc->rank == redf->rank, PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_INCOMP, "Owning rank does not match");
242: PetscCheck(redc->N == redf->N, PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_INCOMP, "Global size does not match");
243: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), P));
244: PetscCall(MatSetSizes(*P, redc->n, redc->n, redc->N, redc->N));
245: PetscCall(MatSetType(*P, MATAIJ));
246: PetscCall(MatSeqAIJSetPreallocation(*P, 1, NULL));
247: PetscCall(MatMPIAIJSetPreallocation(*P, 1, NULL, 0, NULL));
248: PetscCall(MatGetOwnershipRange(*P, &rstart, &rend));
249: for (i = rstart; i < rend; i++) PetscCall(MatSetValue(*P, i, i, 1.0, INSERT_VALUES));
250: PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
251: PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
252: if (scale) PetscCall(DMCreateInterpolationScale(dmc, dmf, *P, scale));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: DMRedundantSetSize - Sets the size of a densely coupled redundant object
259: Collective
261: Input Parameters:
262: + dm - `DM` object of type `DMREDUNDANT`
263: . rank - rank of process to own the redundant degrees of freedom
264: - N - total number of redundant degrees of freedom
266: Level: advanced
268: .seealso: `DM`, `DMREDUNDANT`, `DMDestroy()`, `DMCreateGlobalVector()`, `DMRedundantCreate()`, `DMRedundantGetSize()`
269: @*/
270: PetscErrorCode DMRedundantSetSize(DM dm, PetscMPIInt rank, PetscInt N)
271: {
272: PetscFunctionBegin;
277: PetscTryMethod(dm, "DMRedundantSetSize_C", (DM, PetscMPIInt, PetscInt), (dm, rank, N));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*@
282: DMRedundantGetSize - Gets the size of a densely coupled redundant object
284: Not Collective
286: Input Parameter:
287: . dm - `DM` object of type `DMREDUNDANT`
289: Output Parameters:
290: + rank - rank of process to own the redundant degrees of freedom (or `NULL`)
291: - N - total number of redundant degrees of freedom (or `NULL`)
293: Level: advanced
295: .seealso: `DM`, `DMREDUNDANT`, `DMDestroy()`, `DMCreateGlobalVector()`, `DMRedundantCreate()`, `DMRedundantSetSize()`
296: @*/
297: PetscErrorCode DMRedundantGetSize(DM dm, PetscMPIInt *rank, PetscInt *N)
298: {
299: PetscFunctionBegin;
302: PetscUseMethod(dm, "DMRedundantGetSize_C", (DM, PetscMPIInt *, PetscInt *), (dm, rank, N));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: static PetscErrorCode DMRedundantSetSize_Redundant(DM dm, PetscMPIInt rank, PetscInt N)
307: {
308: DM_Redundant *red = (DM_Redundant *)dm->data;
309: PetscMPIInt myrank;
310: PetscInt i, *globals;
312: PetscFunctionBegin;
313: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &myrank));
314: red->rank = rank;
315: red->N = N;
316: red->n = (myrank == rank) ? N : 0;
318: /* mapping is setup here */
319: PetscCall(PetscMalloc1(red->N, &globals));
320: for (i = 0; i < red->N; i++) globals[i] = i;
321: PetscCall(ISLocalToGlobalMappingDestroy(&dm->ltogmap));
322: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, red->N, globals, PETSC_OWN_POINTER, &dm->ltogmap));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode DMRedundantGetSize_Redundant(DM dm, PetscInt *rank, PetscInt *N)
327: {
328: DM_Redundant *red = (DM_Redundant *)dm->data;
330: PetscFunctionBegin;
331: if (rank) *rank = red->rank;
332: if (N) *N = red->N;
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: static PetscErrorCode DMSetUpGLVisViewer_Redundant(PetscObject odm, PetscViewer viewer)
337: {
338: PetscFunctionBegin;
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*MC
343: DMREDUNDANT = "redundant" - A DM object that is used to manage data for a small set of dense globally coupled variables.
344: In the global representation of the vector the variables are all stored on a single MPI process (all the other MPI processes have
345: no variables) in the local representation all the variables are stored on ALL the MPI processes (because they are all needed for each
346: processes local computations).
348: This DM is generally used inside a DMCOMPOSITE object. For example, it may be used to store continuation parameters for a bifurcation problem.
350: Level: intermediate
352: .seealso: `DMType`, `DMCOMPOSITE`, `DMCreate()`, `DMRedundantSetSize()`, `DMRedundantGetSize()`
353: M*/
355: PETSC_EXTERN PetscErrorCode DMCreate_Redundant(DM dm)
356: {
357: DM_Redundant *red;
359: PetscFunctionBegin;
360: PetscCall(PetscNew(&red));
361: dm->data = red;
363: dm->ops->setup = DMSetUp_Redundant;
364: dm->ops->view = DMView_Redundant;
365: dm->ops->createglobalvector = DMCreateGlobalVector_Redundant;
366: dm->ops->createlocalvector = DMCreateLocalVector_Redundant;
367: dm->ops->creatematrix = DMCreateMatrix_Redundant;
368: dm->ops->destroy = DMDestroy_Redundant;
369: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Redundant;
370: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Redundant;
371: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Redundant;
372: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Redundant;
373: dm->ops->refine = DMRefine_Redundant;
374: dm->ops->coarsen = DMCoarsen_Redundant;
375: dm->ops->createinterpolation = DMCreateInterpolation_Redundant;
376: dm->ops->getcoloring = DMCreateColoring_Redundant;
378: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMRedundantSetSize_C", DMRedundantSetSize_Redundant));
379: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMRedundantGetSize_C", DMRedundantGetSize_Redundant));
380: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Redundant));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@
385: DMRedundantCreate - Creates a `DM` object, used to manage data for dense globally coupled variables
387: Collective
389: Input Parameters:
390: + comm - the processors that will share the global vector
391: . rank - the MPI rank to own the redundant values
392: - N - total number of degrees of freedom
394: Output Parameter:
395: . dm - the `DM` object of type `DMREDUNDANT`
397: Level: advanced
399: .seealso: `DM`, `DMREDUNDANT`, `DMDestroy()`, `DMCreateGlobalVector()`, `DMCreateMatrix()`, `DMCompositeAddDM()`, `DMSetType()`, `DMRedundantSetSize()`, `DMRedundantGetSize()`
400: @*/
401: PetscErrorCode DMRedundantCreate(MPI_Comm comm, PetscMPIInt rank, PetscInt N, DM *dm)
402: {
403: PetscFunctionBegin;
404: PetscAssertPointer(dm, 4);
405: PetscCall(DMCreate(comm, dm));
406: PetscCall(DMSetType(*dm, DMREDUNDANT));
407: PetscCall(DMRedundantSetSize(*dm, rank, N));
408: PetscCall(DMSetUp(*dm));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }