Actual source code: characteristic.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/characteristicimpl.h> /*I "petsccharacteristic.h" I*/
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: PetscBool CharacteristicRegisterAllCalled = PETSC_FALSE;
9: /*
10: Contains the list of registered characteristic routines
11: */
12: PetscFList CharacteristicList = PETSC_NULL;
14: PetscErrorCode DMDAGetNeighborsRank(DM, PetscMPIInt []);
15: PetscInt DMDAGetNeighborRelative(DM, PassiveReal, PassiveReal);
16: PetscErrorCode DMDAMapToPeriodicDomain(DM, PetscScalar [] );
18: PetscErrorCode CharacteristicHeapSort(Characteristic, Queue, PetscInt);
19: PetscErrorCode CharacteristicSiftDown(Characteristic, Queue, PetscInt, PetscInt);
23: PetscErrorCode CharacteristicView(Characteristic c, PetscViewer viewer)
24: {
25: PetscBool iascii;
30: if (!viewer) {
31: PetscViewerASCIIGetStdout(((PetscObject)c)->comm,&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: }
47: PetscErrorCode CharacteristicDestroy(Characteristic *c)
48: {
52: if (!*c) return(0);
54: if (--((PetscObject)(*c))->refct > 0) return(0);
56: if ((*c)->ops->destroy) {
57: (*(*c)->ops->destroy)((*c));
58: }
59: MPI_Type_free(&(*c)->itemType);
60: PetscFree((*c)->queue);
61: PetscFree((*c)->queueLocal);
62: PetscFree((*c)->queueRemote);
63: PetscFree((*c)->neighbors);
64: PetscFree((*c)->needCount);
65: PetscFree((*c)->localOffsets);
66: PetscFree((*c)->fillCount);
67: PetscFree((*c)->remoteOffsets);
68: PetscFree((*c)->request);
69: PetscFree((*c)->status);
70: PetscHeaderDestroy(c);
71: return(0);
72: }
76: PetscErrorCode CharacteristicCreate(MPI_Comm comm, Characteristic *c)
77: {
78: Characteristic newC;
83: *c = PETSC_NULL;
84: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
85: CharacteristicInitializePackage(PETSC_NULL);
86: #endif
88: PetscHeaderCreate(newC, _p_Characteristic, struct _CharacteristicOps, CHARACTERISTIC_CLASSID, -1, "Characteristic", "Characteristic", "SemiLagrange", comm, CharacteristicDestroy, CharacteristicView);
89: PetscLogObjectCreate(newC);
90: *c = newC;
92: newC->structured = PETSC_TRUE;
93: newC->numIds = 0;
94: newC->velocityDA = PETSC_NULL;
95: newC->velocity = PETSC_NULL;
96: newC->velocityOld = PETSC_NULL;
97: newC->numVelocityComp = 0;
98: newC->velocityComp = PETSC_NULL;
99: newC->velocityInterp = PETSC_NULL;
100: newC->velocityInterpLocal = PETSC_NULL;
101: newC->velocityCtx = PETSC_NULL;
102: newC->fieldDA = PETSC_NULL;
103: newC->field = PETSC_NULL;
104: newC->numFieldComp = 0;
105: newC->fieldComp = PETSC_NULL;
106: newC->fieldInterp = PETSC_NULL;
107: newC->fieldInterpLocal = PETSC_NULL;
108: newC->fieldCtx = PETSC_NULL;
109: newC->itemType = PETSC_NULL;
110: newC->queue = PETSC_NULL;
111: newC->queueSize = 0;
112: newC->queueMax = 0;
113: newC->queueLocal = PETSC_NULL;
114: newC->queueLocalSize = 0;
115: newC->queueLocalMax = 0;
116: newC->queueRemote = PETSC_NULL;
117: newC->queueRemoteSize = 0;
118: newC->queueRemoteMax = 0;
119: newC->numNeighbors = 0;
120: newC->neighbors = PETSC_NULL;
121: newC->needCount = PETSC_NULL;
122: newC->localOffsets = PETSC_NULL;
123: newC->fillCount = PETSC_NULL;
124: newC->remoteOffsets = PETSC_NULL;
125: newC->request = PETSC_NULL;
126: newC->status = PETSC_NULL;
127: return(0);
128: }
132: /*@C
133: CharacteristicSetType - Builds Characteristic for a particular solver.
135: Logically Collective on Characteristic
137: Input Parameters:
138: + c - the method of characteristics context
139: - type - a known method
141: Options Database Key:
142: . -characteristic_type <method> - Sets the method; use -help for a list
143: of available methods
145: Notes:
146: See "include/petsccharacteristic.h" for available methods
148: Normally, it is best to use the CharacteristicSetFromOptions() command and
149: then set the Characteristic type from the options database rather than by using
150: this routine. Using the options database provides the user with
151: maximum flexibility in evaluating the many different Krylov methods.
152: The CharacteristicSetType() routine is provided for those situations where it
153: is necessary to set the iterative solver independently of the command
154: line or options database. This might be the case, for example, when
155: the choice of iterative solver changes during the execution of the
156: program, and the user's application is taking responsibility for
157: choosing the appropriate method. In other words, this routine is
158: not for beginners.
160: Level: intermediate
162: .keywords: Characteristic, set, method
164: .seealso: CharacteristicType
166: @*/
167: PetscErrorCode CharacteristicSetType(Characteristic c, const CharacteristicType type)
168: {
169: PetscErrorCode ierr, (*r)(Characteristic);
170: PetscBool match;
176: PetscObjectTypeCompare((PetscObject) c, type, &match);
177: if (match) return(0);
179: if (c->data) {
180: /* destroy the old private Characteristic context */
181: (*c->ops->destroy)(c);
182: c->ops->destroy = PETSC_NULL;
183: c->data = 0;
184: }
186: PetscFListFind(CharacteristicList, ((PetscObject)c)->comm,type,PETSC_TRUE, (void (**)(void)) &r);
187: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Characteristic type given: %s", type);
188: c->setupcalled = 0;
189: (*r)(c);
190: PetscObjectChangeTypeName((PetscObject) c, type);
191: return(0);
192: }
196: /*@
197: CharacteristicSetUp - Sets up the internal data structures for the
198: later use of an iterative solver.
200: Collective on Characteristic
202: Input Parameter:
203: . ksp - iterative context obtained from CharacteristicCreate()
205: Level: developer
207: .keywords: Characteristic, setup
209: .seealso: CharacteristicCreate(), CharacteristicSolve(), CharacteristicDestroy()
210: @*/
211: PetscErrorCode CharacteristicSetUp(Characteristic c)
212: {
218: if (!((PetscObject)c)->type_name){
219: CharacteristicSetType(c, CHARACTERISTICDA);
220: }
222: if (c->setupcalled == 2) return(0);
224: PetscLogEventBegin(CHARACTERISTIC_SetUp,c,0,0,0);
225: if (!c->setupcalled) {
226: (*c->ops->setup)(c);
227: }
228: PetscLogEventEnd(CHARACTERISTIC_SetUp,c,0,0,0);
229: c->setupcalled = 2;
230: return(0);
231: }
235: /*@C
236: CharacteristicRegister - See CharacteristicRegisterDynamic()
238: Level: advanced
239: @*/
240: PetscErrorCode CharacteristicRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Characteristic))
241: {
243: char fullname[PETSC_MAX_PATH_LEN];
246: PetscFListConcat(path,name,fullname);
247: PetscFListAdd(&CharacteristicList,sname,fullname,(void (*)(void))function);
248: return(0);
249: }
253: PetscErrorCode CharacteristicSetVelocityInterpolation(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal [], PetscInt, PetscInt [], PetscScalar [], void *), void *ctx)
254: {
256: c->velocityDA = da;
257: c->velocity = v;
258: c->velocityOld = vOld;
259: c->numVelocityComp = numComponents;
260: c->velocityComp = components;
261: c->velocityInterp = interp;
262: c->velocityCtx = ctx;
263: return(0);
264: }
268: PetscErrorCode CharacteristicSetVelocityInterpolationLocal(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void *, PetscReal [], PetscInt, PetscInt [], PetscScalar [], void *), void *ctx)
269: {
271: c->velocityDA = da;
272: c->velocity = v;
273: c->velocityOld = vOld;
274: c->numVelocityComp = numComponents;
275: c->velocityComp = components;
276: c->velocityInterpLocal = interp;
277: c->velocityCtx = ctx;
278: return(0);
279: }
283: PetscErrorCode CharacteristicSetFieldInterpolation(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal [], PetscInt, PetscInt [], PetscScalar [], void *), void *ctx)
284: {
286: #if 0
287: 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.");
288: #endif
289: c->fieldDA = da;
290: c->field = v;
291: c->numFieldComp = numComponents;
292: c->fieldComp = components;
293: c->fieldInterp = interp;
294: c->fieldCtx = ctx;
295: return(0);
296: }
300: PetscErrorCode CharacteristicSetFieldInterpolationLocal(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void *, PetscReal [], PetscInt, PetscInt [], PetscScalar [], void *), void *ctx)
301: {
303: #if 0
304: 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.");
305: #endif
306: c->fieldDA = da;
307: c->field = v;
308: c->numFieldComp = numComponents;
309: c->fieldComp = components;
310: c->fieldInterpLocal = interp;
311: c->fieldCtx = ctx;
312: return(0);
313: }
317: PetscErrorCode CharacteristicSolve(Characteristic c, PetscReal dt, Vec solution)
318: {
319: CharacteristicPointDA2D Qi;
320: DM da = c->velocityDA;
321: Vec velocityLocal, velocityLocalOld;
322: Vec fieldLocal;
323: DMDALocalInfo info;
324: PetscScalar **solArray;
325: void *velocityArray;
326: void *velocityArrayOld;
327: void *fieldArray;
328: PassiveScalar *interpIndices;
329: PassiveScalar *velocityValues, *velocityValuesOld;
330: PassiveScalar *fieldValues;
331: PetscMPIInt rank;
332: PetscInt dim;
333: PetscMPIInt neighbors[9];
334: PetscInt dof;
335: PetscInt gx, gy;
336: PetscInt n, is, ie, js, je, comp;
337: PetscErrorCode ierr;
340: c->queueSize = 0;
341: MPI_Comm_rank(((PetscObject)c)->comm, &rank);
342: DMDAGetNeighborsRank(da, neighbors);
343: CharacteristicSetNeighbors(c, 9, neighbors);
344: CharacteristicSetUp(c);
345: /* global and local grid info */
346: DMDAGetInfo(da, &dim, &gx, &gy, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
347: DMDAGetLocalInfo(da, &info);
348: is = info.xs; ie = info.xs+info.xm;
349: js = info.ys; je = info.ys+info.ym;
350: /* Allocation */
351: PetscMalloc(dim*sizeof(PetscScalar), &interpIndices);
352: PetscMalloc(c->numVelocityComp*sizeof(PetscScalar), &velocityValues);
353: PetscMalloc(c->numVelocityComp*sizeof(PetscScalar), &velocityValuesOld);
354: PetscMalloc(c->numFieldComp*sizeof(PetscScalar), &fieldValues);
355: PetscLogEventBegin(CHARACTERISTIC_Solve,0,0,0,0);
357: /* -----------------------------------------------------------------------
358: PART 1, AT t-dt/2
359: -----------------------------------------------------------------------*/
360: PetscLogEventBegin(CHARACTERISTIC_QueueSetup,0,0,0,0);
361: /* GET POSITION AT HALF TIME IN THE PAST */
362: if (c->velocityInterpLocal) {
363: DMGetLocalVector(c->velocityDA, &velocityLocal);
364: DMGetLocalVector(c->velocityDA, &velocityLocalOld);
365: DMGlobalToLocalBegin(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);
366: DMGlobalToLocalEnd(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);
367: DMGlobalToLocalBegin(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);
368: DMGlobalToLocalEnd(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);
369: DMDAVecGetArray(c->velocityDA, velocityLocal, &velocityArray);
370: DMDAVecGetArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);
371: }
372: PetscInfo(PETSC_NULL, "Calculating position at t_{n - 1/2}\n");
373: for(Qi.j = js; Qi.j < je; Qi.j++) {
374: for(Qi.i = is; Qi.i < ie; Qi.i++) {
375: interpIndices[0] = Qi.i;
376: interpIndices[1] = Qi.j;
377: if (c->velocityInterpLocal) {
378: c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);
379: } else {
380: c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);
381: }
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(PETSC_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(PETSC_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(PETSC_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(PETSC_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) {
496: c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);
497: } else {
498: c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);
499: }
500: for(comp = 0; comp < c->numFieldComp; comp++) {
501: c->queue[n].field[comp] = fieldValues[comp];
502: }
503: }
504: }
505: PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);
507: PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
508: CharacteristicSendCoordinatesEnd(c);
509: PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
511: /* GET VALUE AT FULL TIME IN THE PAST (REMOTE REQUESTS) */
512: PetscLogEventBegin(CHARACTERISTIC_FullTimeRemote,0,0,0,0);
513: PetscInfo1(PETSC_NULL, "Calculating %d remote field points at t_{n}\n", c->queueRemoteSize);
514: for(n = 0; n < c->queueRemoteSize; n++) {
515: interpIndices[0] = c->queueRemote[n].x;
516: interpIndices[1] = c->queueRemote[n].y;
518: /* for debugging purposes */
519: if (1) { /* hacked bounds test...let's do better */
520: PetscScalar im = interpIndices[0]; PetscScalar jm = interpIndices[1];
522: if (( im < (PetscScalar) is - 1.) || (im > (PetscScalar) ie) || (jm < (PetscScalar) js - 1.) || (jm > (PetscScalar) je)) {
523: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "Nonlocal point: (%g,%g)", im, jm);
524: }
525: }
527: if (c->fieldInterpLocal) {
528: c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);
529: } else {
530: c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);
531: }
532: for(comp = 0; comp < c->numFieldComp; comp++) {
533: c->queueRemote[n].field[comp] = fieldValues[comp];
534: }
535: }
536: PetscLogEventEnd(CHARACTERISTIC_FullTimeRemote,0,0,0,0);
538: PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
539: CharacteristicGetValuesBegin(c);
540: CharacteristicGetValuesEnd(c);
541: if (c->fieldInterpLocal) {
542: DMDAVecRestoreArray(c->fieldDA, fieldLocal, &fieldArray);
543: DMRestoreLocalVector(c->fieldDA, &fieldLocal);
544: }
545: PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);
547: /* Return field of characteristics at t_n-1 */
548: PetscLogEventBegin(CHARACTERISTIC_DAUpdate,0,0,0,0);
549: DMDAGetInfo(c->fieldDA,0,0,0,0,0,0,0,&dof,0,0,0,0,0);
550: DMDAVecGetArray(c->fieldDA, solution, &solArray);
551: for(n = 0; n < c->queueSize; n++) {
552: Qi = c->queue[n];
553: for(comp = 0; comp < c->numFieldComp; comp++) {
554: solArray[Qi.j][Qi.i*dof+c->fieldComp[comp]] = Qi.field[comp];
555: }
556: }
557: DMDAVecRestoreArray(c->fieldDA, solution, &solArray);
558: PetscLogEventEnd(CHARACTERISTIC_DAUpdate,0,0,0,0);
559: PetscLogEventEnd(CHARACTERISTIC_Solve,0,0,0,0);
561: /* Cleanup */
562: PetscFree(interpIndices);
563: PetscFree(velocityValues);
564: PetscFree(velocityValuesOld);
565: PetscFree(fieldValues);
566: return(0);
567: }
571: PetscErrorCode CharacteristicSetNeighbors(Characteristic c, PetscInt numNeighbors, PetscMPIInt neighbors[])
572: {
576: c->numNeighbors = numNeighbors;
577: PetscFree(c->neighbors);
578: PetscMalloc(numNeighbors * sizeof(PetscMPIInt), &c->neighbors);
579: PetscMemcpy(c->neighbors, neighbors, numNeighbors * sizeof(PetscMPIInt));
580: return(0);
581: }
585: PetscErrorCode CharacteristicAddPoint(Characteristic c, CharacteristicPointDA2D *point)
586: {
588: if (c->queueSize >= c->queueMax) {
589: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Exceeeded maximum queue size %d", c->queueMax);
590: }
591: c->queue[c->queueSize++] = *point;
592: return(0);
593: }
597: int CharacteristicSendCoordinatesBegin(Characteristic c)
598: {
599: PetscMPIInt rank, tag = 121;
600: PetscInt i, n;
604: MPI_Comm_rank(((PetscObject)c)->comm, &rank);
605: CharacteristicHeapSort(c, c->queue, c->queueSize);
606: PetscMemzero(c->needCount, c->numNeighbors * sizeof(PetscInt));
607: for(i = 0; i < c->queueSize; i++) {
608: c->needCount[c->queue[i].proc]++;
609: }
610: c->fillCount[0] = 0;
611: for(n = 1; n < c->numNeighbors; n++) {
612: MPI_Irecv(&(c->fillCount[n]), 1, MPIU_INT, c->neighbors[n], tag, ((PetscObject)c)->comm, &(c->request[n-1]));
613: }
614: for(n = 1; n < c->numNeighbors; n++) {
615: MPI_Send(&(c->needCount[n]), 1, MPIU_INT, c->neighbors[n], tag, ((PetscObject)c)->comm);
616: }
617: MPI_Waitall(c->numNeighbors-1, c->request, c->status);
618: /* Initialize the remote queue */
619: c->queueLocalMax = c->localOffsets[0] = 0;
620: c->queueRemoteMax = c->remoteOffsets[0] = 0;
621: for(n = 1; n < c->numNeighbors; n++) {
622: c->remoteOffsets[n] = c->queueRemoteMax;
623: c->queueRemoteMax += c->fillCount[n];
624: c->localOffsets[n] = c->queueLocalMax;
625: c->queueLocalMax += c->needCount[n];
626: }
627: /* HACK BEGIN */
628: for(n = 1; n < c->numNeighbors; n++) {
629: c->localOffsets[n] += c->needCount[0];
630: }
631: c->needCount[0] = 0;
632: /* HACK END */
633: if (c->queueRemoteMax) {
634: PetscMalloc(sizeof(CharacteristicPointDA2D) * c->queueRemoteMax, &c->queueRemote);
635: } else {
636: c->queueRemote = PETSC_NULL;
637: }
638: c->queueRemoteSize = c->queueRemoteMax;
640: /* Send and Receive requests for values at t_n+1/2, giving the coordinates for interpolation */
641: for(n = 1; n < c->numNeighbors; n++) {
642: PetscInfo2(PETSC_NULL, "Receiving %d requests for values from proc %d\n", c->fillCount[n], c->neighbors[n]);
643: MPI_Irecv(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm, &(c->request[n-1]));
644: }
645: for(n = 1; n < c->numNeighbors; n++) {
646: PetscInfo2(PETSC_NULL, "Sending %d requests for values from proc %d\n", c->needCount[n], c->neighbors[n]);
647: MPI_Send(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm);
648: }
649: return(0);
650: }
654: PetscErrorCode CharacteristicSendCoordinatesEnd(Characteristic c)
655: {
656: #if 0
657: PetscMPIInt rank;
658: PetscInt n;
659: #endif
663: MPI_Waitall(c->numNeighbors-1, c->request, c->status);
664: #if 0
665: MPI_Comm_rank(((PetscObject)c)->comm, &rank);
666: for(n = 0; n < c->queueRemoteSize; n++) {
667: if (c->neighbors[c->queueRemote[n].proc] == rank) {
668: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "This is fucked up, n = %d proc = %d", n, c->queueRemote[n].proc);
669: }
670: }
671: #endif
672: return(0);
673: }
677: PetscErrorCode CharacteristicGetValuesBegin(Characteristic c)
678: {
679: PetscMPIInt tag = 121;
680: PetscInt n;
684: /* SEND AND RECIEVE FILLED REQUESTS for velocities at t_n+1/2 */
685: for(n = 1; n < c->numNeighbors; n++) {
686: MPI_Irecv(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm, &(c->request[n-1]));
687: }
688: for(n = 1; n < c->numNeighbors; n++) {
689: MPI_Send(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm);
690: }
691: return(0);
692: }
696: PetscErrorCode CharacteristicGetValuesEnd(Characteristic c)
697: {
701: MPI_Waitall(c->numNeighbors-1, c->request, c->status);
702: /* Free queue of requests from other procs */
703: PetscFree(c->queueRemote);
704: return(0);
705: }
707: /*---------------------------------------------------------------------*/
710: /*
711: Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html
712: */
713: PetscErrorCode CharacteristicHeapSort(Characteristic c, Queue queue, PetscInt size)
714: /*---------------------------------------------------------------------*/
715: {
716: PetscErrorCode ierr;
717: CharacteristicPointDA2D temp;
718: PetscInt n;
719:
721: if (0) { /* Check the order of the queue before sorting */
722: PetscInfo(PETSC_NULL, "Before Heap sort\n");
723: for (n=0; n<size; n++) {
724: PetscInfo2(PETSC_NULL,"%d %d\n",n,queue[n].proc);
725: }
726: }
728: /* SORTING PHASE */
729: for (n = (size / 2)-1; n >= 0; n--) {
730: CharacteristicSiftDown(c, queue, n, size-1); /* Rich had size-1 here, Matt had size*/
731: }
732: for (n = size-1; n >= 1; n--) {
733: temp = queue[0];
734: queue[0] = queue[n];
735: queue[n] = temp;
736: CharacteristicSiftDown(c, queue, 0, n-1);
737: }
738: if (0) { /* Check the order of the queue after sorting */
739: PetscInfo(PETSC_NULL, "Avter Heap sort\n");
740: for (n=0; n<size; n++) {
741: PetscInfo2(PETSC_NULL,"%d %d\n",n,queue[n].proc);
742: }
743: }
744: return(0);
745: }
747: /*---------------------------------------------------------------------*/
750: /*
751: Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html
752: */
753: PetscErrorCode CharacteristicSiftDown(Characteristic c, Queue queue, PetscInt root, PetscInt bottom)
754: /*---------------------------------------------------------------------*/
755: {
756: PetscBool done = PETSC_FALSE;
757: PetscInt maxChild;
758: CharacteristicPointDA2D temp;
761: while ((root*2 <= bottom) && (!done)) {
762: if (root*2 == bottom) maxChild = root * 2;
763: else if (queue[root*2].proc > queue[root*2+1].proc) maxChild = root * 2;
764: else maxChild = root * 2 + 1;
766: if (queue[root].proc < queue[maxChild].proc) {
767: temp = queue[root];
768: queue[root] = queue[maxChild];
769: queue[maxChild] = temp;
770: root = maxChild;
771: } else
772: done = PETSC_TRUE;
773: }
774: return(0);
775: }
779: /* [center, left, top-left, top, top-right, right, bottom-right, bottom, bottom-left] */
780: PetscErrorCode DMDAGetNeighborsRank(DM da, PetscMPIInt neighbors[])
781: {
782: DMDABoundaryType bx, by;
783: PetscBool IPeriodic = PETSC_FALSE, JPeriodic = PETSC_FALSE;
784: MPI_Comm comm;
785: PetscMPIInt rank;
786: PetscInt **procs,pi,pj,pim,pip,pjm,pjp,PI,PJ;
790: PetscObjectGetComm((PetscObject) da, &comm);
791: MPI_Comm_rank(comm, &rank);
792: DMDAGetInfo(da, 0, 0, 0, 0, &PI,&PJ, 0, 0, 0, &bx, &by,0, 0);
794: if (bx == DMDA_BOUNDARY_PERIODIC) {
795: IPeriodic = PETSC_TRUE;
796: }
797: if (by == DMDA_BOUNDARY_PERIODIC) {
798: JPeriodic = PETSC_TRUE;
799: }
801: neighbors[0] = rank;
802: rank = 0;
803: PetscMalloc(sizeof(PetscInt*)*PJ,&procs);
804: for (pj=0;pj<PJ;pj++) {
805: PetscMalloc(sizeof(PetscInt)*PI,&(procs[pj]));
806: for (pi=0;pi<PI;pi++) {
807: procs[pj][pi] = rank;
808: rank++;
809: }
810: }
811:
812: pi = neighbors[0] % PI;
813: pj = neighbors[0] / PI;
814: pim = pi-1; if (pim<0) pim=PI-1;
815: pip = (pi+1)%PI;
816: pjm = pj-1; if (pjm<0) pjm=PJ-1;
817: pjp = (pj+1)%PJ;
819: neighbors[1] = procs[pj] [pim];
820: neighbors[2] = procs[pjp][pim];
821: neighbors[3] = procs[pjp][pi];
822: neighbors[4] = procs[pjp][pip];
823: neighbors[5] = procs[pj] [pip];
824: neighbors[6] = procs[pjm][pip];
825: neighbors[7] = procs[pjm][pi];
826: neighbors[8] = procs[pjm][pim];
828: if (!IPeriodic) {
829: if (pi==0) neighbors[1]=neighbors[2]=neighbors[8]=neighbors[0];
830: if (pi==PI-1) neighbors[4]=neighbors[5]=neighbors[6]=neighbors[0];
831: }
833: if (!JPeriodic) {
834: if (pj==0) neighbors[6]=neighbors[7]=neighbors[8]=neighbors[0];
835: if (pj==PJ-1) neighbors[2]=neighbors[3]=neighbors[4]=neighbors[0];
836: }
838: for(pj = 0; pj < PJ; pj++) {
839: PetscFree(procs[pj]);
840: }
841: PetscFree(procs);
842: return(0);
843: }
847: /*
848: SUBDOMAIN NEIGHBORHOOD PROCESS MAP:
849: 2 | 3 | 4
850: __|___|__
851: 1 | 0 | 5
852: __|___|__
853: 8 | 7 | 6
854: | |
855: */
856: PetscInt DMDAGetNeighborRelative(DM da, PassiveReal ir, PassiveReal jr)
857: {
858: DMDALocalInfo info;
859: PassiveReal is,ie,js,je;
861:
862: DMDAGetLocalInfo(da, &info);
863: is = (PassiveReal) info.xs - 0.5; ie = (PassiveReal) info.xs + info.xm - 0.5;
864: js = (PassiveReal) info.ys - 0.5; je = (PassiveReal) info.ys + info.ym - 0.5;
865:
866: if (ir >= is && ir <= ie) { /* center column */
867: if (jr >= js && jr <= je) {
868: return 0;
869: } else if (jr < js) {
870: return 7;
871: } else {
872: return 3;
873: }
874: } else if (ir < is) { /* left column */
875: if (jr >= js && jr <= je) {
876: return 1;
877: } else if (jr < js) {
878: return 8;
879: } else {
880: return 2;
881: }
882: } else { /* right column */
883: if (jr >= js && jr <= je) {
884: return 5;
885: } else if (jr < js) {
886: return 6;
887: } else {
888: return 4;
889: }
890: }
891: }