Actual source code: characteristic.c
1: #include <petsc/private/characteristicimpl.h>
2: #include <petscdmda.h>
3: #include <petscviewer.h>
5: PetscClassId CHARACTERISTIC_CLASSID;
6: PetscLogEvent CHARACTERISTIC_SetUp, CHARACTERISTIC_Solve, CHARACTERISTIC_QueueSetup, CHARACTERISTIC_DAUpdate;
7: PetscLogEvent CHARACTERISTIC_HalfTimeLocal, CHARACTERISTIC_HalfTimeRemote, CHARACTERISTIC_HalfTimeExchange;
8: PetscLogEvent CHARACTERISTIC_FullTimeLocal, CHARACTERISTIC_FullTimeRemote, CHARACTERISTIC_FullTimeExchange;
9: /*
10: Contains the list of registered characteristic routines
11: */
12: PetscFunctionList CharacteristicList = NULL;
13: PetscBool CharacteristicRegisterAllCalled = PETSC_FALSE;
15: static PetscErrorCode DMDAGetNeighborsRank(DM, PetscMPIInt[]);
16: static PetscMPIInt DMDAGetNeighborRelative(DM, PetscReal, PetscReal);
18: static PetscErrorCode CharacteristicHeapSort(Characteristic, Queue, PetscInt);
19: static PetscErrorCode CharacteristicSiftDown(Characteristic, Queue, PetscInt, PetscInt);
21: static PetscErrorCode CharacteristicView(Characteristic c, PetscViewer viewer)
22: {
23: PetscBool isascii;
25: PetscFunctionBegin;
27: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c), &viewer));
29: PetscCheckSameComm(c, 1, viewer, 2);
31: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
32: if (!isascii) PetscTryTypeMethod(c, view, viewer);
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: /*@
37: CharacteristicDestroy - Destroys a `Characteristic` context created with `CharacteristicCreate()`
39: Collective
41: Input Parameter:
42: . c - the `Characteristic` context
44: Level: beginner
46: .seealso: `Characteristic`, `CharacteristicCreate()`
47: @*/
48: PetscErrorCode CharacteristicDestroy(Characteristic *c)
49: {
50: PetscFunctionBegin;
51: if (!*c) PetscFunctionReturn(PETSC_SUCCESS);
53: if (--((PetscObject)*c)->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
55: PetscTryTypeMethod(*c, destroy);
56: PetscCallMPI(MPI_Type_free(&(*c)->itemType));
57: PetscCall(PetscFree((*c)->queue));
58: PetscCall(PetscFree((*c)->queueLocal));
59: PetscCall(PetscFree((*c)->queueRemote));
60: PetscCall(PetscFree((*c)->neighbors));
61: PetscCall(PetscFree((*c)->needCount));
62: PetscCall(PetscFree((*c)->localOffsets));
63: PetscCall(PetscFree((*c)->fillCount));
64: PetscCall(PetscFree((*c)->remoteOffsets));
65: PetscCall(PetscFree((*c)->request));
66: PetscCall(PetscFree((*c)->status));
67: PetscCall(PetscHeaderDestroy(c));
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: /*@
72: CharacteristicCreate - Creates a `Characteristic` context for use with the Method of Characteristics
74: Collective
76: Input Parameter:
77: . comm - MPI communicator
79: Output Parameter:
80: . c - the `Characteristic` context
82: Level: beginner
84: .seealso: `Characteristic`, `CharacteristicDestroy()`
85: @*/
86: PetscErrorCode CharacteristicCreate(MPI_Comm comm, Characteristic *c)
87: {
88: Characteristic newC;
90: PetscFunctionBegin;
91: PetscAssertPointer(c, 2);
92: *c = NULL;
93: PetscCall(CharacteristicInitializePackage());
95: PetscCall(PetscHeaderCreate(newC, CHARACTERISTIC_CLASSID, "Characteristic", "Characteristic", "Characteristic", comm, CharacteristicDestroy, CharacteristicView));
96: *c = newC;
98: newC->structured = PETSC_TRUE;
99: newC->numIds = 0;
100: newC->velocityDA = NULL;
101: newC->velocity = NULL;
102: newC->velocityOld = NULL;
103: newC->numVelocityComp = 0;
104: newC->velocityComp = NULL;
105: newC->velocityInterp = NULL;
106: newC->velocityInterpLocal = NULL;
107: newC->velocityCtx = NULL;
108: newC->fieldDA = NULL;
109: newC->field = NULL;
110: newC->numFieldComp = 0;
111: newC->fieldComp = NULL;
112: newC->fieldInterp = NULL;
113: newC->fieldInterpLocal = NULL;
114: newC->fieldCtx = NULL;
115: newC->itemType = 0;
116: newC->queue = NULL;
117: newC->queueSize = 0;
118: newC->queueMax = 0;
119: newC->queueLocal = NULL;
120: newC->queueLocalSize = 0;
121: newC->queueLocalMax = 0;
122: newC->queueRemote = NULL;
123: newC->queueRemoteSize = 0;
124: newC->queueRemoteMax = 0;
125: newC->numNeighbors = 0;
126: newC->neighbors = NULL;
127: newC->needCount = NULL;
128: newC->localOffsets = NULL;
129: newC->fillCount = NULL;
130: newC->remoteOffsets = NULL;
131: newC->request = NULL;
132: newC->status = NULL;
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /*@
137: CharacteristicSetType - Builds Characteristic for a particular solver.
139: Logically Collective
141: Input Parameters:
142: + c - the method of characteristics context
143: - type - a known method
145: Options Database Key:
146: . -characteristic_type method - Sets the method; use -help for a list
147: of available methods
149: Level: intermediate
151: Note:
152: See "include/petsccharacteristic.h" for available methods
154: .seealso: [](ch_ts), `CharacteristicType`
155: @*/
156: PetscErrorCode CharacteristicSetType(Characteristic c, CharacteristicType type)
157: {
158: PetscBool match;
159: PetscErrorCode (*r)(Characteristic);
161: PetscFunctionBegin;
163: PetscAssertPointer(type, 2);
165: PetscCall(PetscObjectTypeCompare((PetscObject)c, type, &match));
166: if (match) PetscFunctionReturn(PETSC_SUCCESS);
168: if (c->data) {
169: /* destroy the old private Characteristic context */
170: PetscUseTypeMethod(c, destroy);
171: c->ops->destroy = NULL;
172: c->data = NULL;
173: }
175: PetscCall(PetscFunctionListFind(CharacteristicList, type, &r));
176: PetscCheck(r, PetscObjectComm((PetscObject)c), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Characteristic type given: %s", type);
177: c->setupcalled = PETSC_FALSE;
178: PetscCall((*r)(c));
179: PetscCall(PetscObjectChangeTypeName((PetscObject)c, type));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@
184: CharacteristicSetUp - Sets up the internal data structures for the
185: later use of a `Charactoristic` .
187: Collective
189: Input Parameter:
190: . c - context obtained from CharacteristicCreate()
192: Level: developer
194: .seealso: [](ch_ts), `Characteristic`, `CharacteristicCreate()`, `CharacteristicSolve()`, `CharacteristicDestroy()`
195: @*/
196: PetscErrorCode CharacteristicSetUp(Characteristic c)
197: {
198: PetscFunctionBegin;
201: if (!((PetscObject)c)->type_name) PetscCall(CharacteristicSetType(c, CHARACTERISTICDA));
203: if (c->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
205: PetscCall(PetscLogEventBegin(CHARACTERISTIC_SetUp, c, NULL, NULL, NULL));
206: if (!c->setupcalled) PetscUseTypeMethod(c, setup);
207: PetscCall(PetscLogEventEnd(CHARACTERISTIC_SetUp, c, NULL, NULL, NULL));
208: c->setupcalled = PETSC_TRUE;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*@C
213: CharacteristicRegister - Adds an approarch to the method of characteristics package.
215: Not Collective, No Fortran Support
217: Input Parameters:
218: + sname - name of a new approach
219: - function - routine to create method context
221: Level: advanced
223: Example Usage:
224: .vb
225: CharacteristicRegister("my_char", MyCharCreate);
226: .ve
228: Then, your Characteristic type can be chosen with the procedural interface via
229: .vb
230: CharacteristicCreate(MPI_Comm, Characteristic* &char);
231: CharacteristicSetType(char,"my_char");
232: .ve
233: or at runtime via the option
234: .vb
235: -characteristic_type my_char
236: .ve
238: Notes:
239: `CharacteristicRegister()` may be called multiple times to add several approaches.
241: .seealso: [](ch_ts), `CharacteristicRegisterAll()`, `CharacteristicRegisterDestroy()`
242: @*/
243: PetscErrorCode CharacteristicRegister(const char sname[], PetscErrorCode (*function)(Characteristic))
244: {
245: PetscFunctionBegin;
246: PetscCall(CharacteristicInitializePackage());
247: PetscCall(PetscFunctionListAdd(&CharacteristicList, sname, function));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: /*@C
252: CharacteristicSetVelocityInterpolation - Sets the routine used to interpolate the velocity field at points along a characteristic
254: Not Collective
256: Input Parameters:
257: + c - the `Characteristic` context
258: . da - the `DM` describing the layout of the velocity vectors
259: . v - the current velocity vector
260: . vOld - the previous-time-step velocity vector
261: . numComponents - the number of velocity components to interpolate
262: . components - the indices of the velocity components in `v` and `vOld`
263: . interp - the interpolation routine, called with the global vector
264: - ctx - context passed to the interpolation routine
266: Calling sequence of `interp`:
267: + v - the velocity `Vec` from which to interpolate
268: . interpIndices - the coordinates at which to interpolate
269: . numComponents - the number of components to interpolate
270: . components - the indices of the components in `v`
271: . values - the interpolated values, one per component per point
272: - ctx - the application context
274: Level: developer
276: .seealso: [](ch_ts), `Characteristic`, `CharacteristicSetVelocityInterpolationLocal()`, `CharacteristicSetFieldInterpolation()`
277: @*/
278: PetscErrorCode CharacteristicSetVelocityInterpolation(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec v, PetscReal interpIndices[], PetscInt numComponents, PetscInt components[], PetscScalar values[], PetscCtx ctx), PetscCtx ctx)
279: {
280: PetscFunctionBegin;
281: c->velocityDA = da;
282: c->velocity = v;
283: c->velocityOld = vOld;
284: c->numVelocityComp = numComponents;
285: c->velocityComp = components;
286: c->velocityInterp = interp;
287: c->velocityCtx = ctx;
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /*@C
292: CharacteristicSetVelocityInterpolationLocal - Sets the routine used to interpolate the velocity field along a characteristic using a locally-accessible array
294: Not Collective
296: Input Parameters:
297: + c - the `Characteristic` context
298: . da - the `DM` describing the layout of the velocity vectors
299: . v - the current velocity vector
300: . vOld - the previous-time-step velocity vector
301: . numComponents - the number of velocity components to interpolate
302: . components - the indices of the velocity components in `v` and `vOld`
303: . interp - the interpolation routine, called with a local array pointer rather than a `Vec`
304: - ctx - context passed to the interpolation routine
306: Calling sequence of `interp`:
307: + array - the locally-accessible array of the velocity vector obtained from the `DM`
308: . interpIndices - the coordinates at which to interpolate
309: . numComponents - the number of components to interpolate
310: . components - the indices of the components in the array
311: . values - the interpolated values, one per component per point
312: - ctx - the application context
314: Level: developer
316: .seealso: [](ch_ts), `Characteristic`, `CharacteristicSetVelocityInterpolation()`, `CharacteristicSetFieldInterpolationLocal()`
317: @*/
318: PetscErrorCode CharacteristicSetVelocityInterpolationLocal(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void *array, PetscReal interpIndices[], PetscInt numComponents, PetscInt components[], PetscScalar values[], PetscCtx ctx), PetscCtx ctx)
319: {
320: PetscFunctionBegin;
321: c->velocityDA = da;
322: c->velocity = v;
323: c->velocityOld = vOld;
324: c->numVelocityComp = numComponents;
325: c->velocityComp = components;
326: c->velocityInterpLocal = interp;
327: c->velocityCtx = ctx;
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /*@C
332: CharacteristicSetFieldInterpolation - Sets the routine used to interpolate the field being advected at the foot of a characteristic
334: Not Collective
336: Input Parameters:
337: + c - the `Characteristic` context
338: . da - the `DM` describing the layout of the field vector
339: . v - the field vector to be advected
340: . numComponents - the number of field components to interpolate
341: . components - the indices of the field components in `v`
342: . interp - the interpolation routine, called with the global vector
343: - ctx - context passed to the interpolation routine
345: Calling sequence of `interp`:
346: + v - the field `Vec` from which to interpolate
347: . interpIndices - the coordinates at which to interpolate
348: . numComponents - the number of components to interpolate
349: . components - the indices of the components in `v`
350: . values - the interpolated values, one per component per point
351: - ctx - the application context
353: Level: developer
355: .seealso: [](ch_ts), `Characteristic`, `CharacteristicSetFieldInterpolationLocal()`, `CharacteristicSetVelocityInterpolation()`
356: @*/
357: PetscErrorCode CharacteristicSetFieldInterpolation(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec v, PetscReal interpIndices[], PetscInt numComponents, PetscInt components[], PetscScalar values[], PetscCtx ctx), PetscCtx ctx)
358: {
359: PetscFunctionBegin;
360: #if 0
361: PetscCheck(numComponents <= 2,PETSC_COMM_SELF,PETSC_ERR_SUP, "Fields with more than 2 components are not supported. Send mail to petsc-maint@mcs.anl.gov.");
362: #endif
363: c->fieldDA = da;
364: c->field = v;
365: c->numFieldComp = numComponents;
366: c->fieldComp = components;
367: c->fieldInterp = interp;
368: c->fieldCtx = ctx;
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: /*@C
373: CharacteristicSetFieldInterpolationLocal - Sets the routine used to interpolate the field being advected at the foot of a characteristic using a locally-accessible array
375: Not Collective
377: Input Parameters:
378: + c - the `Characteristic` context
379: . da - the `DM` describing the layout of the field vector
380: . v - the field vector to be advected
381: . numComponents - the number of field components to interpolate
382: . components - the indices of the field components in `v`
383: . interp - the interpolation routine, called with a local array pointer rather than a `Vec`
384: - ctx - context passed to the interpolation routine
386: Calling sequence of `interp`:
387: + array - the locally-accessible array of the field vector obtained from the `DM`
388: . interpIndices - the coordinates at which to interpolate
389: . numComponents - the number of components to interpolate
390: . components - the indices of the components in the array
391: . values - the interpolated values, one per component per point
392: - ctx - the application context
394: Level: developer
396: .seealso: [](ch_ts), `Characteristic`, `CharacteristicSetFieldInterpolation()`, `CharacteristicSetVelocityInterpolationLocal()`
397: @*/
398: PetscErrorCode CharacteristicSetFieldInterpolationLocal(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void *array, PetscReal interpIndices[], PetscInt numComponents, PetscInt components[], PetscScalar values[], PetscCtx ctx), PetscCtx ctx)
399: {
400: PetscFunctionBegin;
401: #if 0
402: PetscCheck(numComponents <= 2,PETSC_COMM_SELF,PETSC_ERR_SUP, "Fields with more than 2 components are not supported. Send mail to petsc-maint@mcs.anl.gov.");
403: #endif
404: c->fieldDA = da;
405: c->field = v;
406: c->numFieldComp = numComponents;
407: c->fieldComp = components;
408: c->fieldInterpLocal = interp;
409: c->fieldCtx = ctx;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@
414: CharacteristicSolve - Apply the Method of Characteristics solver
416: Collective
418: Input Parameters:
419: + c - context obtained from `CharacteristicCreate()`
420: . dt - the time-step
421: - solution - vector holding the solution
423: Level: developer
425: .seealso: [](ch_ts), `Characteristic`, `CharacteristicCreate()`, `CharacteristicDestroy()`
426: @*/
427: PetscErrorCode CharacteristicSolve(Characteristic c, PetscReal dt, Vec solution)
428: {
429: CharacteristicPointDA2D Qi;
430: DM da = c->velocityDA;
431: Vec velocityLocal, velocityLocalOld;
432: Vec fieldLocal;
433: DMDALocalInfo info;
434: PetscScalar **solArray;
435: void *velocityArray;
436: void *velocityArrayOld;
437: void *fieldArray;
438: PetscScalar *interpIndices;
439: PetscScalar *velocityValues, *velocityValuesOld;
440: PetscScalar *fieldValues;
441: PetscMPIInt rank;
442: PetscInt dim;
443: PetscMPIInt neighbors[9];
444: PetscInt dof;
445: PetscInt gx, gy;
446: PetscInt n, is, ie, js, je, comp;
448: PetscFunctionBegin;
449: c->queueSize = 0;
450: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank));
451: PetscCall(DMDAGetNeighborsRank(da, neighbors));
452: PetscCall(CharacteristicSetNeighbors(c, 9, neighbors));
453: PetscCall(CharacteristicSetUp(c));
454: /* global and local grid info */
455: PetscCall(DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
456: PetscCall(DMDAGetLocalInfo(da, &info));
457: is = info.xs;
458: ie = info.xs + info.xm;
459: js = info.ys;
460: je = info.ys + info.ym;
461: /* Allocation */
462: PetscCall(PetscMalloc1(dim, &interpIndices));
463: PetscCall(PetscMalloc1(c->numVelocityComp, &velocityValues));
464: PetscCall(PetscMalloc1(c->numVelocityComp, &velocityValuesOld));
465: PetscCall(PetscMalloc1(c->numFieldComp, &fieldValues));
466: PetscCall(PetscLogEventBegin(CHARACTERISTIC_Solve, NULL, NULL, NULL, NULL));
468: /*
469: PART 1, AT t-dt/2
470: */
471: PetscCall(PetscLogEventBegin(CHARACTERISTIC_QueueSetup, NULL, NULL, NULL, NULL));
472: /* GET POSITION AT HALF TIME IN THE PAST */
473: if (c->velocityInterpLocal) {
474: PetscCall(DMGetLocalVector(c->velocityDA, &velocityLocal));
475: PetscCall(DMGetLocalVector(c->velocityDA, &velocityLocalOld));
476: PetscCall(DMGlobalToLocalBegin(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal));
477: PetscCall(DMGlobalToLocalEnd(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal));
478: PetscCall(DMGlobalToLocalBegin(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld));
479: PetscCall(DMGlobalToLocalEnd(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld));
480: PetscCall(DMDAVecGetArray(c->velocityDA, velocityLocal, &velocityArray));
481: PetscCall(DMDAVecGetArray(c->velocityDA, velocityLocalOld, &velocityArrayOld));
482: }
483: PetscCall(PetscInfo(NULL, "Calculating position at t_{n - 1/2}\n"));
484: for (Qi.j = js; Qi.j < je; Qi.j++) {
485: for (Qi.i = is; Qi.i < ie; Qi.i++) {
486: interpIndices[0] = Qi.i;
487: interpIndices[1] = Qi.j;
488: if (c->velocityInterpLocal) PetscCall(c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx));
489: else PetscCall(c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx));
490: Qi.x = Qi.i - velocityValues[0] * dt / 2.0;
491: Qi.y = Qi.j - velocityValues[1] * dt / 2.0;
493: /* Determine whether the position at t - dt/2 is local */
494: Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y);
496: /* Check for Periodic boundaries and move all periodic points back onto the domain */
497: PetscCall(DMDAMapCoordsToPeriodicDomain(da, &Qi.x, &Qi.y));
498: PetscCall(CharacteristicAddPoint(c, &Qi));
499: }
500: }
501: PetscCall(PetscLogEventEnd(CHARACTERISTIC_QueueSetup, NULL, NULL, NULL, NULL));
503: PetscCall(PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange, NULL, NULL, NULL, NULL));
504: PetscCall(CharacteristicSendCoordinatesBegin(c));
505: PetscCall(PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange, NULL, NULL, NULL, NULL));
507: PetscCall(PetscLogEventBegin(CHARACTERISTIC_HalfTimeLocal, NULL, NULL, NULL, NULL));
508: /* Calculate velocity at t_n+1/2 (local values) */
509: PetscCall(PetscInfo(NULL, "Calculating local velocities at t_{n - 1/2}\n"));
510: for (n = 0; n < c->queueSize; n++) {
511: Qi = c->queue[n];
512: if (c->neighbors[Qi.proc] == rank) {
513: interpIndices[0] = Qi.x;
514: interpIndices[1] = Qi.y;
515: if (c->velocityInterpLocal) {
516: PetscCall(c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx));
517: PetscCall(c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx));
518: } else {
519: PetscCall(c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx));
520: PetscCall(c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx));
521: }
522: Qi.x = 0.5 * (velocityValues[0] + velocityValuesOld[0]);
523: Qi.y = 0.5 * (velocityValues[1] + velocityValuesOld[1]);
524: }
525: c->queue[n] = Qi;
526: }
527: PetscCall(PetscLogEventEnd(CHARACTERISTIC_HalfTimeLocal, NULL, NULL, NULL, NULL));
529: PetscCall(PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange, NULL, NULL, NULL, NULL));
530: PetscCall(CharacteristicSendCoordinatesEnd(c));
531: PetscCall(PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange, NULL, NULL, NULL, NULL));
533: /* Calculate velocity at t_n+1/2 (fill remote requests) */
534: PetscCall(PetscLogEventBegin(CHARACTERISTIC_HalfTimeRemote, NULL, NULL, NULL, NULL));
535: PetscCall(PetscInfo(NULL, "Calculating %" PetscInt_FMT " remote velocities at t_{n - 1/2}\n", c->queueRemoteSize));
536: for (n = 0; n < c->queueRemoteSize; n++) {
537: Qi = c->queueRemote[n];
538: interpIndices[0] = Qi.x;
539: interpIndices[1] = Qi.y;
540: if (c->velocityInterpLocal) {
541: PetscCall(c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx));
542: PetscCall(c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx));
543: } else {
544: PetscCall(c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx));
545: PetscCall(c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx));
546: }
547: Qi.x = 0.5 * (velocityValues[0] + velocityValuesOld[0]);
548: Qi.y = 0.5 * (velocityValues[1] + velocityValuesOld[1]);
549: c->queueRemote[n] = Qi;
550: }
551: PetscCall(PetscLogEventEnd(CHARACTERISTIC_HalfTimeRemote, NULL, NULL, NULL, NULL));
552: PetscCall(PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange, NULL, NULL, NULL, NULL));
553: PetscCall(CharacteristicGetValuesBegin(c));
554: PetscCall(CharacteristicGetValuesEnd(c));
555: if (c->velocityInterpLocal) {
556: PetscCall(DMDAVecRestoreArray(c->velocityDA, velocityLocal, &velocityArray));
557: PetscCall(DMDAVecRestoreArray(c->velocityDA, velocityLocalOld, &velocityArrayOld));
558: PetscCall(DMRestoreLocalVector(c->velocityDA, &velocityLocal));
559: PetscCall(DMRestoreLocalVector(c->velocityDA, &velocityLocalOld));
560: }
561: PetscCall(PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange, NULL, NULL, NULL, NULL));
563: /*
564: PART 2, AT t-dt
565: */
567: /* GET POSITION AT t_n (local values) */
568: PetscCall(PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal, NULL, NULL, NULL, NULL));
569: PetscCall(PetscInfo(NULL, "Calculating position at t_{n}\n"));
570: for (n = 0; n < c->queueSize; n++) {
571: Qi = c->queue[n];
572: Qi.x = Qi.i - Qi.x * dt;
573: Qi.y = Qi.j - Qi.y * dt;
575: /* Determine whether the position at t-dt is local */
576: Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y);
578: /* Check for Periodic boundaries and move all periodic points back onto the domain */
579: PetscCall(DMDAMapCoordsToPeriodicDomain(da, &Qi.x, &Qi.y));
581: c->queue[n] = Qi;
582: }
583: PetscCall(PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal, NULL, NULL, NULL, NULL));
585: PetscCall(PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange, NULL, NULL, NULL, NULL));
586: PetscCall(CharacteristicSendCoordinatesBegin(c));
587: PetscCall(PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange, NULL, NULL, NULL, NULL));
589: /* GET VALUE AT FULL TIME IN THE PAST (LOCAL REQUESTS) */
590: PetscCall(PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal, NULL, NULL, NULL, NULL));
591: if (c->fieldInterpLocal) {
592: PetscCall(DMGetLocalVector(c->fieldDA, &fieldLocal));
593: PetscCall(DMGlobalToLocalBegin(c->fieldDA, c->field, INSERT_VALUES, fieldLocal));
594: PetscCall(DMGlobalToLocalEnd(c->fieldDA, c->field, INSERT_VALUES, fieldLocal));
595: PetscCall(DMDAVecGetArray(c->fieldDA, fieldLocal, &fieldArray));
596: }
597: PetscCall(PetscInfo(NULL, "Calculating local field at t_{n}\n"));
598: for (n = 0; n < c->queueSize; n++) {
599: if (c->neighbors[c->queue[n].proc] == rank) {
600: interpIndices[0] = c->queue[n].x;
601: interpIndices[1] = c->queue[n].y;
602: if (c->fieldInterpLocal) PetscCall(c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx));
603: else PetscCall(c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx));
604: for (comp = 0; comp < c->numFieldComp; comp++) c->queue[n].field[comp] = fieldValues[comp];
605: }
606: }
607: PetscCall(PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal, NULL, NULL, NULL, NULL));
609: PetscCall(PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange, NULL, NULL, NULL, NULL));
610: PetscCall(CharacteristicSendCoordinatesEnd(c));
611: PetscCall(PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange, NULL, NULL, NULL, NULL));
613: /* GET VALUE AT FULL TIME IN THE PAST (REMOTE REQUESTS) */
614: PetscCall(PetscLogEventBegin(CHARACTERISTIC_FullTimeRemote, NULL, NULL, NULL, NULL));
615: PetscCall(PetscInfo(NULL, "Calculating %" PetscInt_FMT " remote field points at t_{n}\n", c->queueRemoteSize));
616: for (n = 0; n < c->queueRemoteSize; n++) {
617: interpIndices[0] = c->queueRemote[n].x;
618: interpIndices[1] = c->queueRemote[n].y;
620: /* for debugging purposes */
621: if (1) { /* hacked bounds test...let's do better */
622: PetscScalar im = interpIndices[0];
623: PetscScalar jm = interpIndices[1];
625: PetscCheck((im >= (PetscScalar)is - 1.) && (im <= (PetscScalar)ie) && (jm >= (PetscScalar)js - 1.) && (jm <= (PetscScalar)je), PETSC_COMM_SELF, PETSC_ERR_LIB, "Nonlocal point: (%g,%g)", (double)PetscAbsScalar(im), (double)PetscAbsScalar(jm));
626: }
628: if (c->fieldInterpLocal) PetscCall(c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx));
629: else PetscCall(c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx));
630: for (comp = 0; comp < c->numFieldComp; comp++) c->queueRemote[n].field[comp] = fieldValues[comp];
631: }
632: PetscCall(PetscLogEventEnd(CHARACTERISTIC_FullTimeRemote, NULL, NULL, NULL, NULL));
634: PetscCall(PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange, NULL, NULL, NULL, NULL));
635: PetscCall(CharacteristicGetValuesBegin(c));
636: PetscCall(CharacteristicGetValuesEnd(c));
637: if (c->fieldInterpLocal) {
638: PetscCall(DMDAVecRestoreArray(c->fieldDA, fieldLocal, &fieldArray));
639: PetscCall(DMRestoreLocalVector(c->fieldDA, &fieldLocal));
640: }
641: PetscCall(PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange, NULL, NULL, NULL, NULL));
643: /* Return field of characteristics at t_n-1 */
644: PetscCall(PetscLogEventBegin(CHARACTERISTIC_DAUpdate, NULL, NULL, NULL, NULL));
645: PetscCall(DMDAGetInfo(c->fieldDA, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
646: PetscCall(DMDAVecGetArray(c->fieldDA, solution, &solArray));
647: for (n = 0; n < c->queueSize; n++) {
648: Qi = c->queue[n];
649: for (comp = 0; comp < c->numFieldComp; comp++) solArray[Qi.j][Qi.i * dof + c->fieldComp[comp]] = Qi.field[comp];
650: }
651: PetscCall(DMDAVecRestoreArray(c->fieldDA, solution, &solArray));
652: PetscCall(PetscLogEventEnd(CHARACTERISTIC_DAUpdate, NULL, NULL, NULL, NULL));
653: PetscCall(PetscLogEventEnd(CHARACTERISTIC_Solve, NULL, NULL, NULL, NULL));
655: /* Cleanup */
656: PetscCall(PetscFree(interpIndices));
657: PetscCall(PetscFree(velocityValues));
658: PetscCall(PetscFree(velocityValuesOld));
659: PetscCall(PetscFree(fieldValues));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: PetscErrorCode CharacteristicSetNeighbors(Characteristic c, PetscInt numNeighbors, PetscMPIInt neighbors[])
664: {
665: PetscFunctionBegin;
666: PetscCall(PetscMPIIntCast(numNeighbors, &c->numNeighbors));
667: PetscCall(PetscFree(c->neighbors));
668: PetscCall(PetscMalloc1(numNeighbors, &c->neighbors));
669: PetscCall(PetscArraycpy(c->neighbors, neighbors, numNeighbors));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: PetscErrorCode CharacteristicAddPoint(Characteristic c, CharacteristicPointDA2D *point)
674: {
675: PetscFunctionBegin;
676: PetscCheck(c->queueSize < c->queueMax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Exceeded maximum queue size %" PetscInt_FMT, c->queueMax);
677: c->queue[c->queueSize++] = *point;
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: PetscErrorCode CharacteristicSendCoordinatesBegin(Characteristic c)
682: {
683: PetscMPIInt rank, tag = 121;
684: PetscInt n;
686: PetscFunctionBegin;
687: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank));
688: PetscCall(CharacteristicHeapSort(c, c->queue, c->queueSize));
689: PetscCall(PetscArrayzero(c->needCount, c->numNeighbors));
690: for (PetscInt i = 0; i < c->queueSize; i++) c->needCount[c->queue[i].proc]++;
691: c->fillCount[0] = 0;
692: for (n = 1; n < c->numNeighbors; n++) PetscCallMPI(MPIU_Irecv(&c->fillCount[n], 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &c->request[n - 1]));
693: for (n = 1; n < c->numNeighbors; n++) PetscCallMPI(MPIU_Send(&c->needCount[n], 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c)));
694: PetscCallMPI(MPI_Waitall(c->numNeighbors - 1, c->request, c->status));
695: /* Initialize the remote queue */
696: c->queueLocalMax = c->localOffsets[0] = 0;
697: c->queueRemoteMax = c->remoteOffsets[0] = 0;
698: for (n = 1; n < c->numNeighbors; n++) {
699: c->remoteOffsets[n] = c->queueRemoteMax;
700: c->queueRemoteMax += c->fillCount[n];
701: c->localOffsets[n] = c->queueLocalMax;
702: c->queueLocalMax += c->needCount[n];
703: }
704: /* HACK BEGIN */
705: for (n = 1; n < c->numNeighbors; n++) c->localOffsets[n] += c->needCount[0];
706: c->needCount[0] = 0;
707: /* HACK END */
708: if (c->queueRemoteMax) PetscCall(PetscMalloc1(c->queueRemoteMax, &c->queueRemote));
709: else c->queueRemote = NULL;
710: c->queueRemoteSize = c->queueRemoteMax;
712: /* Send and Receive requests for values at t_n+1/2, giving the coordinates for interpolation */
713: for (n = 1; n < c->numNeighbors; n++) {
714: PetscCall(PetscInfo(NULL, "Receiving %" PetscInt_FMT " requests for values from proc %d\n", c->fillCount[n], c->neighbors[n]));
715: PetscCallMPI(MPIU_Irecv(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &c->request[n - 1]));
716: }
717: for (n = 1; n < c->numNeighbors; n++) {
718: PetscCall(PetscInfo(NULL, "Sending %" PetscInt_FMT " requests for values from proc %d\n", c->needCount[n], c->neighbors[n]));
719: PetscCallMPI(MPIU_Send(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c)));
720: }
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: PetscErrorCode CharacteristicSendCoordinatesEnd(Characteristic c)
725: {
726: #if 0
727: PetscMPIInt rank;
728: PetscInt n;
729: #endif
731: PetscFunctionBegin;
732: PetscCallMPI(MPI_Waitall(c->numNeighbors - 1, c->request, c->status));
733: #if 0
734: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank));
735: for (n = 0; n < c->queueRemoteSize; n++) PetscCheck(c->neighbors[c->queueRemote[n].proc] != rank,PETSC_COMM_SELF,PETSC_ERR_PLIB, "This is messed up, n = %d proc = %d", n, c->queueRemote[n].proc);
736: #endif
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: PetscErrorCode CharacteristicGetValuesBegin(Characteristic c)
741: {
742: PetscMPIInt tag = 121;
743: PetscInt n;
745: PetscFunctionBegin;
746: /* SEND AND RECEIVE FILLED REQUESTS for velocities at t_n+1/2 */
747: for (n = 1; n < c->numNeighbors; n++) PetscCallMPI(MPIU_Irecv(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &c->request[n - 1]));
748: for (n = 1; n < c->numNeighbors; n++) PetscCallMPI(MPIU_Send(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c)));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: PetscErrorCode CharacteristicGetValuesEnd(Characteristic c)
753: {
754: PetscFunctionBegin;
755: PetscCallMPI(MPI_Waitall(c->numNeighbors - 1, c->request, c->status));
756: /* Free queue of requests from other procs */
757: PetscCall(PetscFree(c->queueRemote));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: /*
762: Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html
763: */
764: static PetscErrorCode CharacteristicHeapSort(Characteristic c, Queue queue, PetscInt size)
765: {
766: CharacteristicPointDA2D temp;
767: PetscInt n;
769: PetscFunctionBegin;
770: if (0) { /* Check the order of the queue before sorting */
771: PetscCall(PetscInfo(NULL, "Before Heap sort\n"));
772: for (n = 0; n < size; n++) PetscCall(PetscInfo(NULL, "%" PetscInt_FMT " %d\n", n, queue[n].proc));
773: }
775: /* SORTING PHASE */
776: for (n = (size / 2) - 1; n >= 0; n--) PetscCall(CharacteristicSiftDown(c, queue, n, size - 1)); /* Rich had size-1 here, Matt had size*/
777: for (n = size - 1; n >= 1; n--) {
778: temp = queue[0];
779: queue[0] = queue[n];
780: queue[n] = temp;
781: PetscCall(CharacteristicSiftDown(c, queue, 0, n - 1));
782: }
783: if (0) { /* Check the order of the queue after sorting */
784: PetscCall(PetscInfo(NULL, "Avter Heap sort\n"));
785: for (n = 0; n < size; n++) PetscCall(PetscInfo(NULL, "%" PetscInt_FMT " %d\n", n, queue[n].proc));
786: }
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /*
791: Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html
792: */
793: static PetscErrorCode CharacteristicSiftDown(Characteristic c, Queue queue, PetscInt root, PetscInt bottom)
794: {
795: PetscBool done = PETSC_FALSE;
796: PetscInt maxChild;
797: CharacteristicPointDA2D temp;
799: PetscFunctionBegin;
800: while ((root * 2 <= bottom) && (!done)) {
801: if (root * 2 == bottom) maxChild = root * 2;
802: else if (queue[root * 2].proc > queue[root * 2 + 1].proc) maxChild = root * 2;
803: else maxChild = root * 2 + 1;
805: if (queue[root].proc < queue[maxChild].proc) {
806: temp = queue[root];
807: queue[root] = queue[maxChild];
808: queue[maxChild] = temp;
809: root = maxChild;
810: } else done = PETSC_TRUE;
811: }
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: /* [center, left, top-left, top, top-right, right, bottom-right, bottom, bottom-left] */
816: static PetscErrorCode DMDAGetNeighborsRank(DM da, PetscMPIInt neighbors[])
817: {
818: DMBoundaryType bx, by;
819: PetscBool IPeriodic = PETSC_FALSE, JPeriodic = PETSC_FALSE;
820: MPI_Comm comm;
821: PetscMPIInt rank;
822: PetscMPIInt **procs, pi, pj, pim, pip, pjm, pjp, PIi, PJi;
823: PetscInt PI, PJ;
825: PetscFunctionBegin;
826: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
827: PetscCallMPI(MPI_Comm_rank(comm, &rank));
828: PetscCall(DMDAGetInfo(da, NULL, NULL, NULL, NULL, &PI, &PJ, NULL, NULL, NULL, &bx, &by, NULL, NULL));
829: PetscCall(PetscMPIIntCast(PI, &PIi));
830: PetscCall(PetscMPIIntCast(PJ, &PJi));
831: if (bx == DM_BOUNDARY_PERIODIC) IPeriodic = PETSC_TRUE;
832: if (by == DM_BOUNDARY_PERIODIC) JPeriodic = PETSC_TRUE;
834: neighbors[0] = rank;
835: rank = 0;
836: PetscCall(PetscMalloc1(PJ, &procs));
837: for (pj = 0; pj < PJ; pj++) {
838: PetscCall(PetscMalloc1(PI, &procs[pj]));
839: for (pi = 0; pi < PI; pi++) {
840: procs[pj][pi] = rank;
841: rank++;
842: }
843: }
845: pi = neighbors[0] % PI;
846: pj = neighbors[0] / PI;
847: pim = pi - 1;
848: if (pim < 0) pim = PIi - 1;
849: pip = (pi + 1) % PIi;
850: pjm = pj - 1;
851: if (pjm < 0) pjm = PJi - 1;
852: pjp = (pj + 1) % PJi;
854: neighbors[1] = procs[pj][pim];
855: neighbors[2] = procs[pjp][pim];
856: neighbors[3] = procs[pjp][pi];
857: neighbors[4] = procs[pjp][pip];
858: neighbors[5] = procs[pj][pip];
859: neighbors[6] = procs[pjm][pip];
860: neighbors[7] = procs[pjm][pi];
861: neighbors[8] = procs[pjm][pim];
863: if (!IPeriodic) {
864: if (pi == 0) neighbors[1] = neighbors[2] = neighbors[8] = neighbors[0];
865: if (pi == PI - 1) neighbors[4] = neighbors[5] = neighbors[6] = neighbors[0];
866: }
868: if (!JPeriodic) {
869: if (pj == 0) neighbors[6] = neighbors[7] = neighbors[8] = neighbors[0];
870: if (pj == PJ - 1) neighbors[2] = neighbors[3] = neighbors[4] = neighbors[0];
871: }
873: for (pj = 0; pj < PJ; pj++) PetscCall(PetscFree(procs[pj]));
874: PetscCall(PetscFree(procs));
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: /*
879: SUBDOMAIN NEIGHBORHOOD PROCESS MAP:
880: 2 | 3 | 4
881: __|___|__
882: 1 | 0 | 5
883: __|___|__
884: 8 | 7 | 6
885: | |
886: */
887: static PetscMPIInt DMDAGetNeighborRelative(DM da, PetscReal ir, PetscReal jr)
888: {
889: DMDALocalInfo info;
890: PetscReal is, ie, js, je;
892: PetscCallAbort(PETSC_COMM_SELF, DMDAGetLocalInfo(da, &info));
893: is = (PetscReal)info.xs - 0.5;
894: ie = (PetscReal)info.xs + info.xm - 0.5;
895: js = (PetscReal)info.ys - 0.5;
896: je = (PetscReal)info.ys + info.ym - 0.5;
898: if (ir >= is && ir <= ie) { /* center column */
899: if (jr >= js && jr <= je) return 0;
900: else if (jr < js) return 7;
901: else return 3;
902: } else if (ir < is) { /* left column */
903: if (jr >= js && jr <= je) return 1;
904: else if (jr < js) return 8;
905: else return 2;
906: } else { /* right column */
907: if (jr >= js && jr <= je) return 5;
908: else if (jr < js) return 6;
909: else return 4;
910: }
911: }