Actual source code: characteristic.c

petsc-3.8.2 2017-11-09
Report Typos and Errors

  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: }