Actual source code: characteristic.c
petsc-3.9.4 2018-09-11
2: #include <petsc/private/characteristicimpl.h>
3: #include <petscdmda.h>
4: #include <petscviewer.h>
6: PetscClassId CHARACTERISTIC_CLASSID;
7: PetscLogEvent CHARACTERISTIC_SetUp, CHARACTERISTIC_Solve, CHARACTERISTIC_QueueSetup, CHARACTERISTIC_DAUpdate;
8: PetscLogEvent CHARACTERISTIC_HalfTimeLocal, CHARACTERISTIC_HalfTimeRemote, CHARACTERISTIC_HalfTimeExchange;
9: PetscLogEvent CHARACTERISTIC_FullTimeLocal, CHARACTERISTIC_FullTimeRemote, CHARACTERISTIC_FullTimeExchange;
10: /*
11: Contains the list of registered characteristic routines
12: */
13: PetscFunctionList CharacteristicList = NULL;
14: PetscBool CharacteristicRegisterAllCalled = PETSC_FALSE;
16: PetscErrorCode DMDAGetNeighborsRank(DM, PetscMPIInt []);
17: PetscInt DMDAGetNeighborRelative(DM, PetscReal, PetscReal);
18: PetscErrorCode DMDAMapToPeriodicDomain(DM, PetscScalar []);
20: PetscErrorCode CharacteristicHeapSort(Characteristic, Queue, PetscInt);
21: PetscErrorCode CharacteristicSiftDown(Characteristic, Queue, PetscInt, PetscInt);
23: PetscErrorCode CharacteristicView(Characteristic c, PetscViewer viewer)
24: {
25: PetscBool iascii;
30: if (!viewer) {
31: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
32: }
36: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
37: if (!iascii) {
38: if (c->ops->view) {
39: (*c->ops->view)(c, viewer);
40: }
41: }
42: return(0);
43: }
45: PetscErrorCode CharacteristicDestroy(Characteristic *c)
46: {
50: if (!*c) return(0);
52: if (--((PetscObject)(*c))->refct > 0) return(0);
54: if ((*c)->ops->destroy) {
55: (*(*c)->ops->destroy)((*c));
56: }
57: MPI_Type_free(&(*c)->itemType);
58: PetscFree((*c)->queue);
59: PetscFree((*c)->queueLocal);
60: PetscFree((*c)->queueRemote);
61: PetscFree((*c)->neighbors);
62: PetscFree((*c)->needCount);
63: PetscFree((*c)->localOffsets);
64: PetscFree((*c)->fillCount);
65: PetscFree((*c)->remoteOffsets);
66: PetscFree((*c)->request);
67: PetscFree((*c)->status);
68: PetscHeaderDestroy(c);
69: return(0);
70: }
72: PetscErrorCode CharacteristicCreate(MPI_Comm comm, Characteristic *c)
73: {
74: Characteristic newC;
79: *c = NULL;
80: CharacteristicInitializePackage();
82: PetscHeaderCreate(newC, CHARACTERISTIC_CLASSID, "Characteristic", "Characteristic", "Characteristic", comm, CharacteristicDestroy, CharacteristicView);
83: *c = newC;
85: newC->structured = PETSC_TRUE;
86: newC->numIds = 0;
87: newC->velocityDA = NULL;
88: newC->velocity = NULL;
89: newC->velocityOld = NULL;
90: newC->numVelocityComp = 0;
91: newC->velocityComp = NULL;
92: newC->velocityInterp = NULL;
93: newC->velocityInterpLocal = NULL;
94: newC->velocityCtx = NULL;
95: newC->fieldDA = NULL;
96: newC->field = NULL;
97: newC->numFieldComp = 0;
98: newC->fieldComp = NULL;
99: newC->fieldInterp = NULL;
100: newC->fieldInterpLocal = NULL;
101: newC->fieldCtx = NULL;
102: newC->itemType = 0;
103: newC->queue = NULL;
104: newC->queueSize = 0;
105: newC->queueMax = 0;
106: newC->queueLocal = NULL;
107: newC->queueLocalSize = 0;
108: newC->queueLocalMax = 0;
109: newC->queueRemote = NULL;
110: newC->queueRemoteSize = 0;
111: newC->queueRemoteMax = 0;
112: newC->numNeighbors = 0;
113: newC->neighbors = NULL;
114: newC->needCount = NULL;
115: newC->localOffsets = NULL;
116: newC->fillCount = NULL;
117: newC->remoteOffsets = NULL;
118: newC->request = NULL;
119: newC->status = NULL;
120: return(0);
121: }
123: /*@C
124: CharacteristicSetType - Builds Characteristic for a particular solver.
126: Logically Collective on Characteristic
128: Input Parameters:
129: + c - the method of characteristics context
130: - type - a known method
132: Options Database Key:
133: . -characteristic_type <method> - Sets the method; use -help for a list
134: of available methods
136: Notes:
137: See "include/petsccharacteristic.h" for available methods
139: Normally, it is best to use the CharacteristicSetFromOptions() command and
140: then set the Characteristic type from the options database rather than by using
141: this routine. Using the options database provides the user with
142: maximum flexibility in evaluating the many different Krylov methods.
143: The CharacteristicSetType() routine is provided for those situations where it
144: is necessary to set the iterative solver independently of the command
145: line or options database. This might be the case, for example, when
146: the choice of iterative solver changes during the execution of the
147: program, and the user's application is taking responsibility for
148: choosing the appropriate method. In other words, this routine is
149: not for beginners.
151: Level: intermediate
153: .keywords: Characteristic, set, method
155: .seealso: CharacteristicType
157: @*/
158: PetscErrorCode CharacteristicSetType(Characteristic c, CharacteristicType type)
159: {
160: PetscErrorCode ierr, (*r)(Characteristic);
161: PetscBool match;
167: PetscObjectTypeCompare((PetscObject) c, type, &match);
168: if (match) return(0);
170: if (c->data) {
171: /* destroy the old private Characteristic context */
172: (*c->ops->destroy)(c);
173: c->ops->destroy = NULL;
174: c->data = 0;
175: }
177: PetscFunctionListFind(CharacteristicList,type,&r);
178: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Characteristic type given: %s", type);
179: c->setupcalled = 0;
180: (*r)(c);
181: PetscObjectChangeTypeName((PetscObject) c, type);
182: return(0);
183: }
185: /*@
186: CharacteristicSetUp - Sets up the internal data structures for the
187: later use of an iterative solver.
189: Collective on Characteristic
191: Input Parameter:
192: . ksp - iterative context obtained from CharacteristicCreate()
194: Level: developer
196: .keywords: Characteristic, setup
198: .seealso: CharacteristicCreate(), CharacteristicSolve(), CharacteristicDestroy()
199: @*/
200: PetscErrorCode CharacteristicSetUp(Characteristic c)
201: {
207: if (!((PetscObject)c)->type_name) {
208: CharacteristicSetType(c, CHARACTERISTICDA);
209: }
211: if (c->setupcalled == 2) return(0);
213: PetscLogEventBegin(CHARACTERISTIC_SetUp,c,0,0,0);
214: if (!c->setupcalled) {
215: (*c->ops->setup)(c);
216: }
217: PetscLogEventEnd(CHARACTERISTIC_SetUp,c,0,0,0);
218: c->setupcalled = 2;
219: return(0);
220: }
222: /*@C
223: CharacteristicRegister - Adds a solver to the method of characteristics package.
225: Not Collective
227: Input Parameters:
228: + name_solver - name of a new user-defined solver
229: - routine_create - routine to create method context
231: Sample usage:
232: .vb
233: CharacteristicRegister("my_char", MyCharCreate);
234: .ve
236: Then, your Characteristic type can be chosen with the procedural interface via
237: .vb
238: CharacteristicCreate(MPI_Comm, Characteristic* &char);
239: CharacteristicSetType(char,"my_char");
240: .ve
241: or at runtime via the option
242: .vb
243: -characteristic_type my_char
244: .ve
246: Notes:
247: CharacteristicRegister() may be called multiple times to add several user-defined solvers.
249: .keywords: Characteristic, register
251: .seealso: CharacteristicRegisterAll(), CharacteristicRegisterDestroy()
253: Level: advanced
254: @*/
255: PetscErrorCode CharacteristicRegister(const char sname[],PetscErrorCode (*function)(Characteristic))
256: {
260: PetscFunctionListAdd(&CharacteristicList,sname,function);
261: return(0);
262: }
264: PetscErrorCode CharacteristicSetVelocityInterpolation(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal[], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx)
265: {
267: c->velocityDA = da;
268: c->velocity = v;
269: c->velocityOld = vOld;
270: c->numVelocityComp = numComponents;
271: c->velocityComp = components;
272: c->velocityInterp = interp;
273: c->velocityCtx = ctx;
274: return(0);
275: }
277: PetscErrorCode CharacteristicSetVelocityInterpolationLocal(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void*, PetscReal [], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx)
278: {
280: c->velocityDA = da;
281: c->velocity = v;
282: c->velocityOld = vOld;
283: c->numVelocityComp = numComponents;
284: c->velocityComp = components;
285: c->velocityInterpLocal = interp;
286: c->velocityCtx = ctx;
287: return(0);
288: }
290: PetscErrorCode CharacteristicSetFieldInterpolation(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal[], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx)
291: {
293: #if 0
294: if (numComponents > 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Fields with more than 2 components are not supported. Send mail to petsc-maint@mcs.anl.gov.");
295: #endif
296: c->fieldDA = da;
297: c->field = v;
298: c->numFieldComp = numComponents;
299: c->fieldComp = components;
300: c->fieldInterp = interp;
301: c->fieldCtx = ctx;
302: return(0);
303: }
305: PetscErrorCode CharacteristicSetFieldInterpolationLocal(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void*, PetscReal[], PetscInt, PetscInt[], PetscScalar [], void*), void *ctx)
306: {
308: #if 0
309: if (numComponents > 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Fields with more than 2 components are not supported. Send mail to petsc-maint@mcs.anl.gov.");
310: #endif
311: c->fieldDA = da;
312: c->field = v;
313: c->numFieldComp = numComponents;
314: c->fieldComp = components;
315: c->fieldInterpLocal = interp;
316: c->fieldCtx = ctx;
317: return(0);
318: }
320: PetscErrorCode CharacteristicSolve(Characteristic c, PetscReal dt, Vec solution)
321: {
322: CharacteristicPointDA2D Qi;
323: DM da = c->velocityDA;
324: Vec velocityLocal, velocityLocalOld;
325: Vec fieldLocal;
326: DMDALocalInfo info;
327: PetscScalar **solArray;
328: void *velocityArray;
329: void *velocityArrayOld;
330: void *fieldArray;
331: PetscScalar *interpIndices;
332: PetscScalar *velocityValues, *velocityValuesOld;
333: PetscScalar *fieldValues;
334: PetscMPIInt rank;
335: PetscInt dim;
336: PetscMPIInt neighbors[9];
337: PetscInt dof;
338: PetscInt gx, gy;
339: PetscInt n, is, ie, js, je, comp;
340: PetscErrorCode ierr;
343: c->queueSize = 0;
344: MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);
345: DMDAGetNeighborsRank(da, neighbors);
346: CharacteristicSetNeighbors(c, 9, neighbors);
347: CharacteristicSetUp(c);
348: /* global and local grid info */
349: DMDAGetInfo(da, &dim, &gx, &gy, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
350: DMDAGetLocalInfo(da, &info);
351: is = info.xs; ie = info.xs+info.xm;
352: js = info.ys; je = info.ys+info.ym;
353: /* Allocation */
354: PetscMalloc1(dim, &interpIndices);
355: PetscMalloc1(c->numVelocityComp, &velocityValues);
356: PetscMalloc1(c->numVelocityComp, &velocityValuesOld);
357: PetscMalloc1(c->numFieldComp, &fieldValues);
358: PetscLogEventBegin(CHARACTERISTIC_Solve,0,0,0,0);
360: /* -----------------------------------------------------------------------
361: PART 1, AT t-dt/2
362: -----------------------------------------------------------------------*/
363: PetscLogEventBegin(CHARACTERISTIC_QueueSetup,0,0,0,0);
364: /* GET POSITION AT HALF TIME IN THE PAST */
365: if (c->velocityInterpLocal) {
366: DMGetLocalVector(c->velocityDA, &velocityLocal);
367: DMGetLocalVector(c->velocityDA, &velocityLocalOld);
368: DMGlobalToLocalBegin(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);
369: DMGlobalToLocalEnd(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);
370: DMGlobalToLocalBegin(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);
371: DMGlobalToLocalEnd(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);
372: DMDAVecGetArray(c->velocityDA, velocityLocal, &velocityArray);
373: DMDAVecGetArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);
374: }
375: PetscInfo(NULL, "Calculating position at t_{n - 1/2}\n");
376: for (Qi.j = js; Qi.j < je; Qi.j++) {
377: for (Qi.i = is; Qi.i < ie; Qi.i++) {
378: interpIndices[0] = Qi.i;
379: interpIndices[1] = Qi.j;
380: if (c->velocityInterpLocal) {c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);}
381: else {c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);}
382: Qi.x = Qi.i - velocityValues[0]*dt/2.0;
383: Qi.y = Qi.j - velocityValues[1]*dt/2.0;
385: /* Determine whether the position at t - dt/2 is local */
386: Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y);
388: /* Check for Periodic boundaries and move all periodic points back onto the domain */
389: DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));
390: CharacteristicAddPoint(c, &Qi);
391: }
392: }
393: PetscLogEventEnd(CHARACTERISTIC_QueueSetup,0,0,0,0);
395: PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);
396: CharacteristicSendCoordinatesBegin(c);
397: PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);
399: PetscLogEventBegin(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);
400: /* Calculate velocity at t_n+1/2 (local values) */
401: PetscInfo(NULL, "Calculating local velocities at t_{n - 1/2}\n");
402: for (n = 0; n < c->queueSize; n++) {
403: Qi = c->queue[n];
404: if (c->neighbors[Qi.proc] == rank) {
405: interpIndices[0] = Qi.x;
406: interpIndices[1] = Qi.y;
407: if (c->velocityInterpLocal) {
408: c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);
409: c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);
410: } else {
411: c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);
412: c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);
413: }
414: Qi.x = 0.5*(velocityValues[0] + velocityValuesOld[0]);
415: Qi.y = 0.5*(velocityValues[1] + velocityValuesOld[1]);
416: }
417: c->queue[n] = Qi;
418: }
419: PetscLogEventEnd(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);
421: PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);
422: CharacteristicSendCoordinatesEnd(c);
423: PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);
426: /* Calculate velocity at t_n+1/2 (fill remote requests) */
427: PetscLogEventBegin(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);
428: PetscInfo1(NULL, "Calculating %d remote velocities at t_{n - 1/2}\n", c->queueRemoteSize);
429: for (n = 0; n < c->queueRemoteSize; n++) {
430: Qi = c->queueRemote[n];
431: interpIndices[0] = Qi.x;
432: interpIndices[1] = Qi.y;
433: if (c->velocityInterpLocal) {
434: c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);
435: c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);
436: } else {
437: c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);
438: c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);
439: }
440: Qi.x = 0.5*(velocityValues[0] + velocityValuesOld[0]);
441: Qi.y = 0.5*(velocityValues[1] + velocityValuesOld[1]);
442: c->queueRemote[n] = Qi;
443: }
444: PetscLogEventEnd(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);
445: PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);
446: CharacteristicGetValuesBegin(c);
447: CharacteristicGetValuesEnd(c);
448: if (c->velocityInterpLocal) {
449: DMDAVecRestoreArray(c->velocityDA, velocityLocal, &velocityArray);
450: DMDAVecRestoreArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);
451: DMRestoreLocalVector(c->velocityDA, &velocityLocal);
452: DMRestoreLocalVector(c->velocityDA, &velocityLocalOld);
453: }
454: PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);
456: /* -----------------------------------------------------------------------
457: PART 2, AT t-dt
458: -----------------------------------------------------------------------*/
460: /* GET POSITION AT t_n (local values) */
461: PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);
462: PetscInfo(NULL, "Calculating position at t_{n}\n");
463: for (n = 0; n < c->queueSize; n++) {
464: Qi = c->queue[n];
465: Qi.x = Qi.i - Qi.x*dt;
466: Qi.y = Qi.j - Qi.y*dt;
468: /* Determine whether the position at t-dt is local */
469: Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y);
471: /* Check for Periodic boundaries and move all periodic points back onto the domain */
472: DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));
474: c->queue[n] = Qi;
475: }
476: PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);
478: PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
479: CharacteristicSendCoordinatesBegin(c);
480: PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
482: /* GET VALUE AT FULL TIME IN THE PAST (LOCAL REQUESTS) */
483: PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);
484: if (c->fieldInterpLocal) {
485: DMGetLocalVector(c->fieldDA, &fieldLocal);
486: DMGlobalToLocalBegin(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);
487: DMGlobalToLocalEnd(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);
488: DMDAVecGetArray(c->fieldDA, fieldLocal, &fieldArray);
489: }
490: PetscInfo(NULL, "Calculating local field at t_{n}\n");
491: for (n = 0; n < c->queueSize; n++) {
492: if (c->neighbors[c->queue[n].proc] == rank) {
493: interpIndices[0] = c->queue[n].x;
494: interpIndices[1] = c->queue[n].y;
495: if (c->fieldInterpLocal) {c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);}
496: else {c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);}
497: for (comp = 0; comp < c->numFieldComp; comp++) c->queue[n].field[comp] = fieldValues[comp];
498: }
499: }
500: PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);
502: PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
503: CharacteristicSendCoordinatesEnd(c);
504: PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
506: /* GET VALUE AT FULL TIME IN THE PAST (REMOTE REQUESTS) */
507: PetscLogEventBegin(CHARACTERISTIC_FullTimeRemote,0,0,0,0);
508: PetscInfo1(NULL, "Calculating %d remote field points at t_{n}\n", c->queueRemoteSize);
509: for (n = 0; n < c->queueRemoteSize; n++) {
510: interpIndices[0] = c->queueRemote[n].x;
511: interpIndices[1] = c->queueRemote[n].y;
513: /* for debugging purposes */
514: if (1) { /* hacked bounds test...let's do better */
515: PetscScalar im = interpIndices[0]; PetscScalar jm = interpIndices[1];
517: if ((im < (PetscScalar) is - 1.) || (im > (PetscScalar) ie) || (jm < (PetscScalar) js - 1.) || (jm > (PetscScalar) je)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "Nonlocal point: (%g,%g)", im, jm);
518: }
520: if (c->fieldInterpLocal) {c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);}
521: else {c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);}
522: for (comp = 0; comp < c->numFieldComp; comp++) c->queueRemote[n].field[comp] = fieldValues[comp];
523: }
524: PetscLogEventEnd(CHARACTERISTIC_FullTimeRemote,0,0,0,0);
526: PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
527: CharacteristicGetValuesBegin(c);
528: CharacteristicGetValuesEnd(c);
529: if (c->fieldInterpLocal) {
530: DMDAVecRestoreArray(c->fieldDA, fieldLocal, &fieldArray);
531: DMRestoreLocalVector(c->fieldDA, &fieldLocal);
532: }
533: PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
535: /* Return field of characteristics at t_n-1 */
536: PetscLogEventBegin(CHARACTERISTIC_DAUpdate,0,0,0,0);
537: DMDAGetInfo(c->fieldDA,0,0,0,0,0,0,0,&dof,0,0,0,0,0);
538: DMDAVecGetArray(c->fieldDA, solution, &solArray);
539: for (n = 0; n < c->queueSize; n++) {
540: Qi = c->queue[n];
541: for (comp = 0; comp < c->numFieldComp; comp++) solArray[Qi.j][Qi.i*dof+c->fieldComp[comp]] = Qi.field[comp];
542: }
543: DMDAVecRestoreArray(c->fieldDA, solution, &solArray);
544: PetscLogEventEnd(CHARACTERISTIC_DAUpdate,0,0,0,0);
545: PetscLogEventEnd(CHARACTERISTIC_Solve,0,0,0,0);
547: /* Cleanup */
548: PetscFree(interpIndices);
549: PetscFree(velocityValues);
550: PetscFree(velocityValuesOld);
551: PetscFree(fieldValues);
552: return(0);
553: }
555: PetscErrorCode CharacteristicSetNeighbors(Characteristic c, PetscInt numNeighbors, PetscMPIInt neighbors[])
556: {
560: c->numNeighbors = numNeighbors;
561: PetscFree(c->neighbors);
562: PetscMalloc1(numNeighbors, &c->neighbors);
563: PetscMemcpy(c->neighbors, neighbors, numNeighbors * sizeof(PetscMPIInt));
564: return(0);
565: }
567: PetscErrorCode CharacteristicAddPoint(Characteristic c, CharacteristicPointDA2D *point)
568: {
570: if (c->queueSize >= c->queueMax) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Exceeeded maximum queue size %d", c->queueMax);
571: c->queue[c->queueSize++] = *point;
572: return(0);
573: }
575: int CharacteristicSendCoordinatesBegin(Characteristic c)
576: {
577: PetscMPIInt rank, tag = 121;
578: PetscInt i, n;
582: MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);
583: CharacteristicHeapSort(c, c->queue, c->queueSize);
584: PetscMemzero(c->needCount, c->numNeighbors * sizeof(PetscInt));
585: for (i = 0; i < c->queueSize; i++) c->needCount[c->queue[i].proc]++;
586: c->fillCount[0] = 0;
587: for (n = 1; n < c->numNeighbors; n++) {
588: MPI_Irecv(&(c->fillCount[n]), 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &(c->request[n-1]));
589: }
590: for (n = 1; n < c->numNeighbors; n++) {
591: MPI_Send(&(c->needCount[n]), 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));
592: }
593: MPI_Waitall(c->numNeighbors-1, c->request, c->status);
594: /* Initialize the remote queue */
595: c->queueLocalMax = c->localOffsets[0] = 0;
596: c->queueRemoteMax = c->remoteOffsets[0] = 0;
597: for (n = 1; n < c->numNeighbors; n++) {
598: c->remoteOffsets[n] = c->queueRemoteMax;
599: c->queueRemoteMax += c->fillCount[n];
600: c->localOffsets[n] = c->queueLocalMax;
601: c->queueLocalMax += c->needCount[n];
602: }
603: /* HACK BEGIN */
604: for (n = 1; n < c->numNeighbors; n++) c->localOffsets[n] += c->needCount[0];
605: c->needCount[0] = 0;
606: /* HACK END */
607: if (c->queueRemoteMax) {
608: PetscMalloc1(c->queueRemoteMax, &c->queueRemote);
609: } else c->queueRemote = NULL;
610: c->queueRemoteSize = c->queueRemoteMax;
612: /* Send and Receive requests for values at t_n+1/2, giving the coordinates for interpolation */
613: for (n = 1; n < c->numNeighbors; n++) {
614: PetscInfo2(NULL, "Receiving %d requests for values from proc %d\n", c->fillCount[n], c->neighbors[n]);
615: MPI_Irecv(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &(c->request[n-1]));
616: }
617: for (n = 1; n < c->numNeighbors; n++) {
618: PetscInfo2(NULL, "Sending %d requests for values from proc %d\n", c->needCount[n], c->neighbors[n]);
619: MPI_Send(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));
620: }
621: return(0);
622: }
624: PetscErrorCode CharacteristicSendCoordinatesEnd(Characteristic c)
625: {
626: #if 0
627: PetscMPIInt rank;
628: PetscInt n;
629: #endif
633: MPI_Waitall(c->numNeighbors-1, c->request, c->status);
634: #if 0
635: MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);
636: for (n = 0; n < c->queueRemoteSize; n++) {
637: if (c->neighbors[c->queueRemote[n].proc] == rank) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "This is messed up, n = %d proc = %d", n, c->queueRemote[n].proc);
638: }
639: #endif
640: return(0);
641: }
643: PetscErrorCode CharacteristicGetValuesBegin(Characteristic c)
644: {
645: PetscMPIInt tag = 121;
646: PetscInt n;
650: /* SEND AND RECIEVE FILLED REQUESTS for velocities at t_n+1/2 */
651: for (n = 1; n < c->numNeighbors; n++) {
652: MPI_Irecv(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &(c->request[n-1]));
653: }
654: for (n = 1; n < c->numNeighbors; n++) {
655: MPI_Send(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));
656: }
657: return(0);
658: }
660: PetscErrorCode CharacteristicGetValuesEnd(Characteristic c)
661: {
665: MPI_Waitall(c->numNeighbors-1, c->request, c->status);
666: /* Free queue of requests from other procs */
667: PetscFree(c->queueRemote);
668: return(0);
669: }
671: /*---------------------------------------------------------------------*/
672: /*
673: Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html
674: */
675: PetscErrorCode CharacteristicHeapSort(Characteristic c, Queue queue, PetscInt size)
676: /*---------------------------------------------------------------------*/
677: {
678: PetscErrorCode ierr;
679: CharacteristicPointDA2D temp;
680: PetscInt n;
683: if (0) { /* Check the order of the queue before sorting */
684: PetscInfo(NULL, "Before Heap sort\n");
685: for (n=0; n<size; n++) {
686: PetscInfo2(NULL,"%d %d\n",n,queue[n].proc);
687: }
688: }
690: /* SORTING PHASE */
691: for (n = (size / 2)-1; n >= 0; n--) {
692: CharacteristicSiftDown(c, queue, n, size-1); /* Rich had size-1 here, Matt had size*/
693: }
694: for (n = size-1; n >= 1; n--) {
695: temp = queue[0];
696: queue[0] = queue[n];
697: queue[n] = temp;
698: CharacteristicSiftDown(c, queue, 0, n-1);
699: }
700: if (0) { /* Check the order of the queue after sorting */
701: PetscInfo(NULL, "Avter Heap sort\n");
702: for (n=0; n<size; n++) {
703: PetscInfo2(NULL,"%d %d\n",n,queue[n].proc);
704: }
705: }
706: return(0);
707: }
709: /*---------------------------------------------------------------------*/
710: /*
711: Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html
712: */
713: PetscErrorCode CharacteristicSiftDown(Characteristic c, Queue queue, PetscInt root, PetscInt bottom)
714: /*---------------------------------------------------------------------*/
715: {
716: PetscBool done = PETSC_FALSE;
717: PetscInt maxChild;
718: CharacteristicPointDA2D temp;
721: while ((root*2 <= bottom) && (!done)) {
722: if (root*2 == bottom) maxChild = root * 2;
723: else if (queue[root*2].proc > queue[root*2+1].proc) maxChild = root * 2;
724: else maxChild = root * 2 + 1;
726: if (queue[root].proc < queue[maxChild].proc) {
727: temp = queue[root];
728: queue[root] = queue[maxChild];
729: queue[maxChild] = temp;
730: root = maxChild;
731: } else done = PETSC_TRUE;
732: }
733: return(0);
734: }
736: /* [center, left, top-left, top, top-right, right, bottom-right, bottom, bottom-left] */
737: PetscErrorCode DMDAGetNeighborsRank(DM da, PetscMPIInt neighbors[])
738: {
739: DMBoundaryType bx, by;
740: PetscBool IPeriodic = PETSC_FALSE, JPeriodic = PETSC_FALSE;
741: MPI_Comm comm;
742: PetscMPIInt rank;
743: PetscInt **procs,pi,pj,pim,pip,pjm,pjp,PI,PJ;
744: PetscErrorCode ierr;
747: PetscObjectGetComm((PetscObject) da, &comm);
748: MPI_Comm_rank(comm, &rank);
749: DMDAGetInfo(da, 0, 0, 0, 0, &PI,&PJ, 0, 0, 0, &bx, &by,0, 0);
751: if (bx == DM_BOUNDARY_PERIODIC) IPeriodic = PETSC_TRUE;
752: if (by == DM_BOUNDARY_PERIODIC) JPeriodic = PETSC_TRUE;
754: neighbors[0] = rank;
755: rank = 0;
756: PetscMalloc1(PJ,&procs);
757: for (pj=0; pj<PJ; pj++) {
758: PetscMalloc1(PI,&(procs[pj]));
759: for (pi=0; pi<PI; pi++) {
760: procs[pj][pi] = rank;
761: rank++;
762: }
763: }
765: pi = neighbors[0] % PI;
766: pj = neighbors[0] / PI;
767: pim = pi-1; if (pim<0) pim=PI-1;
768: pip = (pi+1)%PI;
769: pjm = pj-1; if (pjm<0) pjm=PJ-1;
770: pjp = (pj+1)%PJ;
772: neighbors[1] = procs[pj] [pim];
773: neighbors[2] = procs[pjp][pim];
774: neighbors[3] = procs[pjp][pi];
775: neighbors[4] = procs[pjp][pip];
776: neighbors[5] = procs[pj] [pip];
777: neighbors[6] = procs[pjm][pip];
778: neighbors[7] = procs[pjm][pi];
779: neighbors[8] = procs[pjm][pim];
781: if (!IPeriodic) {
782: if (pi==0) neighbors[1]=neighbors[2]=neighbors[8]=neighbors[0];
783: if (pi==PI-1) neighbors[4]=neighbors[5]=neighbors[6]=neighbors[0];
784: }
786: if (!JPeriodic) {
787: if (pj==0) neighbors[6]=neighbors[7]=neighbors[8]=neighbors[0];
788: if (pj==PJ-1) neighbors[2]=neighbors[3]=neighbors[4]=neighbors[0];
789: }
791: for (pj = 0; pj < PJ; pj++) {
792: PetscFree(procs[pj]);
793: }
794: PetscFree(procs);
795: return(0);
796: }
798: /*
799: SUBDOMAIN NEIGHBORHOOD PROCESS MAP:
800: 2 | 3 | 4
801: __|___|__
802: 1 | 0 | 5
803: __|___|__
804: 8 | 7 | 6
805: | |
806: */
807: PetscInt DMDAGetNeighborRelative(DM da, PetscReal ir, PetscReal jr)
808: {
809: DMDALocalInfo info;
810: PetscReal is,ie,js,je;
813: DMDAGetLocalInfo(da, &info);
814: is = (PetscReal) info.xs - 0.5; ie = (PetscReal) info.xs + info.xm - 0.5;
815: js = (PetscReal) info.ys - 0.5; je = (PetscReal) info.ys + info.ym - 0.5;
817: if (ir >= is && ir <= ie) { /* center column */
818: if (jr >= js && jr <= je) return 0;
819: else if (jr < js) return 7;
820: else return 3;
821: } else if (ir < is) { /* left column */
822: if (jr >= js && jr <= je) return 1;
823: else if (jr < js) return 8;
824: else return 2;
825: } else { /* right column */
826: if (jr >= js && jr <= je) return 5;
827: else if (jr < js) return 6;
828: else return 4;
829: }
830: }