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