1 2 #include <petsc-private/characteristicimpl.h> /*I "petsccharacteristic.h" I*/ 3 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 = PETSC_NULL; 12 PetscBool CharacteristicRegisterAllCalled = PETSC_FALSE; 13 14 PetscErrorCode DMDAGetNeighborsRank(DM, PetscMPIInt []); 15 PetscInt DMDAGetNeighborRelative(DM, PassiveReal, PassiveReal); 16 PetscErrorCode DMDAMapToPeriodicDomain(DM, PetscScalar []); 17 18 PetscErrorCode CharacteristicHeapSort(Characteristic, Queue, PetscInt); 19 PetscErrorCode CharacteristicSiftDown(Characteristic, Queue, PetscInt, PetscInt); 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "CharacteristicView" 23 PetscErrorCode CharacteristicView(Characteristic c, PetscViewer viewer) 24 { 25 PetscBool iascii; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 PetscValidHeaderSpecific(c, CHARACTERISTIC_CLASSID, 1); 30 if (!viewer) { 31 ierr = PetscViewerASCIIGetStdout(((PetscObject)c)->comm,&viewer);CHKERRQ(ierr); 32 } 33 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 34 PetscCheckSameComm(c, 1, viewer, 2); 35 36 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 37 if (!iascii) { 38 if (c->ops->view) { 39 ierr = (*c->ops->view)(c, viewer);CHKERRQ(ierr); 40 } 41 } 42 PetscFunctionReturn(0); 43 } 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "CharacteristicDestroy" 47 PetscErrorCode CharacteristicDestroy(Characteristic *c) 48 { 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 if (!*c) PetscFunctionReturn(0); 53 PetscValidHeaderSpecific(*c, CHARACTERISTIC_CLASSID, 1); 54 if (--((PetscObject)(*c))->refct > 0) PetscFunctionReturn(0); 55 56 if ((*c)->ops->destroy) { 57 ierr = (*(*c)->ops->destroy)((*c));CHKERRQ(ierr); 58 } 59 ierr = MPI_Type_free(&(*c)->itemType);CHKERRQ(ierr); 60 ierr = PetscFree((*c)->queue);CHKERRQ(ierr); 61 ierr = PetscFree((*c)->queueLocal);CHKERRQ(ierr); 62 ierr = PetscFree((*c)->queueRemote);CHKERRQ(ierr); 63 ierr = PetscFree((*c)->neighbors);CHKERRQ(ierr); 64 ierr = PetscFree((*c)->needCount);CHKERRQ(ierr); 65 ierr = PetscFree((*c)->localOffsets);CHKERRQ(ierr); 66 ierr = PetscFree((*c)->fillCount);CHKERRQ(ierr); 67 ierr = PetscFree((*c)->remoteOffsets);CHKERRQ(ierr); 68 ierr = PetscFree((*c)->request);CHKERRQ(ierr); 69 ierr = PetscFree((*c)->status);CHKERRQ(ierr); 70 ierr = PetscHeaderDestroy(c);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "CharacteristicCreate" 76 PetscErrorCode CharacteristicCreate(MPI_Comm comm, Characteristic *c) 77 { 78 Characteristic newC; 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 PetscValidPointer(c, 2); 83 *c = PETSC_NULL; 84 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 85 ierr = CharacteristicInitializePackage(PETSC_NULL);CHKERRQ(ierr); 86 #endif 87 88 ierr = PetscHeaderCreate(newC, _p_Characteristic, struct _CharacteristicOps, CHARACTERISTIC_CLASSID, -1, "Characteristic", "Characteristic", "SemiLagrange", comm, CharacteristicDestroy, CharacteristicView);CHKERRQ(ierr); 89 ierr = PetscLogObjectCreate(newC);CHKERRQ(ierr); 90 *c = newC; 91 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 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "CharacteristicSetType" 132 /*@C 133 CharacteristicSetType - Builds Characteristic for a particular solver. 134 135 Logically Collective on Characteristic 136 137 Input Parameters: 138 + c - the method of characteristics context 139 - type - a known method 140 141 Options Database Key: 142 . -characteristic_type <method> - Sets the method; use -help for a list 143 of available methods 144 145 Notes: 146 See "include/petsccharacteristic.h" for available methods 147 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. 159 160 Level: intermediate 161 162 .keywords: Characteristic, set, method 163 164 .seealso: CharacteristicType 165 166 @*/ 167 PetscErrorCode CharacteristicSetType(Characteristic c, CharacteristicType type) 168 { 169 PetscErrorCode ierr, (*r)(Characteristic); 170 PetscBool match; 171 172 PetscFunctionBegin; 173 PetscValidHeaderSpecific(c, CHARACTERISTIC_CLASSID, 1); 174 PetscValidCharPointer(type, 2); 175 176 ierr = PetscObjectTypeCompare((PetscObject) c, type, &match);CHKERRQ(ierr); 177 if (match) PetscFunctionReturn(0); 178 179 if (c->data) { 180 /* destroy the old private Characteristic context */ 181 ierr = (*c->ops->destroy)(c);CHKERRQ(ierr); 182 c->ops->destroy = PETSC_NULL; 183 c->data = 0; 184 } 185 186 ierr = PetscFunctionListFind(((PetscObject)c)->comm,CharacteristicList, type,PETSC_TRUE, (void (**)(void)) &r);CHKERRQ(ierr); 187 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Characteristic type given: %s", type); 188 c->setupcalled = 0; 189 ierr = (*r)(c);CHKERRQ(ierr); 190 ierr = PetscObjectChangeTypeName((PetscObject) c, type);CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 194 #undef __FUNCT__ 195 #define __FUNCT__ "CharacteristicSetUp" 196 /*@ 197 CharacteristicSetUp - Sets up the internal data structures for the 198 later use of an iterative solver. 199 200 Collective on Characteristic 201 202 Input Parameter: 203 . ksp - iterative context obtained from CharacteristicCreate() 204 205 Level: developer 206 207 .keywords: Characteristic, setup 208 209 .seealso: CharacteristicCreate(), CharacteristicSolve(), CharacteristicDestroy() 210 @*/ 211 PetscErrorCode CharacteristicSetUp(Characteristic c) 212 { 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 PetscValidHeaderSpecific(c, CHARACTERISTIC_CLASSID, 1); 217 218 if (!((PetscObject)c)->type_name) { 219 ierr = CharacteristicSetType(c, CHARACTERISTICDA);CHKERRQ(ierr); 220 } 221 222 if (c->setupcalled == 2) PetscFunctionReturn(0); 223 224 ierr = PetscLogEventBegin(CHARACTERISTIC_SetUp,c,0,0,0);CHKERRQ(ierr); 225 if (!c->setupcalled) { 226 ierr = (*c->ops->setup)(c);CHKERRQ(ierr); 227 } 228 ierr = PetscLogEventEnd(CHARACTERISTIC_SetUp,c,0,0,0);CHKERRQ(ierr); 229 c->setupcalled = 2; 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "CharacteristicRegister" 235 /*@C 236 CharacteristicRegister - See CharacteristicRegisterDynamic() 237 238 Level: advanced 239 @*/ 240 PetscErrorCode CharacteristicRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Characteristic)) 241 { 242 PetscErrorCode ierr; 243 char fullname[PETSC_MAX_PATH_LEN]; 244 245 PetscFunctionBegin; 246 ierr = PetscFunctionListConcat(path,name,fullname);CHKERRQ(ierr); 247 ierr = PetscFunctionListAdd(PETSC_COMM_WORLD,&CharacteristicList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 #undef __FUNCT__ 252 #define __FUNCT__ "CharacteristicSetVelocityInterpolation" 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 { 255 PetscFunctionBegin; 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 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNCT__ 267 #define __FUNCT__ "CharacteristicSetVelocityInterpolationLocal" 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 { 270 PetscFunctionBegin; 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 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "CharacteristicSetFieldInterpolation" 283 PetscErrorCode CharacteristicSetFieldInterpolation(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal[], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx) 284 { 285 PetscFunctionBegin; 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 PetscFunctionReturn(0); 296 } 297 298 #undef __FUNCT__ 299 #define __FUNCT__ "CharacteristicSetFieldInterpolationLocal" 300 PetscErrorCode CharacteristicSetFieldInterpolationLocal(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void*, PetscReal[], PetscInt, PetscInt[], PetscScalar [], void*), void *ctx) 301 { 302 PetscFunctionBegin; 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 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "CharacteristicSolve" 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; 338 339 PetscFunctionBegin; 340 c->queueSize = 0; 341 ierr = MPI_Comm_rank(((PetscObject)c)->comm, &rank);CHKERRQ(ierr); 342 ierr = DMDAGetNeighborsRank(da, neighbors);CHKERRQ(ierr); 343 ierr = CharacteristicSetNeighbors(c, 9, neighbors);CHKERRQ(ierr); 344 ierr = CharacteristicSetUp(c);CHKERRQ(ierr); 345 /* global and local grid info */ 346 ierr = DMDAGetInfo(da, &dim, &gx, &gy, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);CHKERRQ(ierr); 347 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 348 is = info.xs; ie = info.xs+info.xm; 349 js = info.ys; je = info.ys+info.ym; 350 /* Allocation */ 351 ierr = PetscMalloc(dim*sizeof(PetscScalar), &interpIndices);CHKERRQ(ierr); 352 ierr = PetscMalloc(c->numVelocityComp*sizeof(PetscScalar), &velocityValues);CHKERRQ(ierr); 353 ierr = PetscMalloc(c->numVelocityComp*sizeof(PetscScalar), &velocityValuesOld);CHKERRQ(ierr); 354 ierr = PetscMalloc(c->numFieldComp*sizeof(PetscScalar), &fieldValues);CHKERRQ(ierr); 355 ierr = PetscLogEventBegin(CHARACTERISTIC_Solve,0,0,0,0);CHKERRQ(ierr); 356 357 /* ----------------------------------------------------------------------- 358 PART 1, AT t-dt/2 359 -----------------------------------------------------------------------*/ 360 ierr = PetscLogEventBegin(CHARACTERISTIC_QueueSetup,0,0,0,0);CHKERRQ(ierr); 361 /* GET POSITION AT HALF TIME IN THE PAST */ 362 if (c->velocityInterpLocal) { 363 ierr = DMGetLocalVector(c->velocityDA, &velocityLocal);CHKERRQ(ierr); 364 ierr = DMGetLocalVector(c->velocityDA, &velocityLocalOld);CHKERRQ(ierr); 365 ierr = DMGlobalToLocalBegin(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);CHKERRQ(ierr); 366 ierr = DMGlobalToLocalEnd(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);CHKERRQ(ierr); 367 ierr = DMGlobalToLocalBegin(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);CHKERRQ(ierr); 368 ierr = DMGlobalToLocalEnd(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);CHKERRQ(ierr); 369 ierr = DMDAVecGetArray(c->velocityDA, velocityLocal, &velocityArray);CHKERRQ(ierr); 370 ierr = DMDAVecGetArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);CHKERRQ(ierr); 371 } 372 ierr = PetscInfo(PETSC_NULL, "Calculating position at t_{n - 1/2}\n");CHKERRQ(ierr); 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) c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx); 378 else c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx); 379 Qi.x = Qi.i - velocityValues[0]*dt/2.0; 380 Qi.y = Qi.j - velocityValues[1]*dt/2.0; 381 382 /* Determine whether the position at t - dt/2 is local */ 383 Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y); 384 385 /* Check for Periodic boundaries and move all periodic points back onto the domain */ 386 ierr = DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));CHKERRQ(ierr); 387 ierr = CharacteristicAddPoint(c, &Qi);CHKERRQ(ierr); 388 } 389 } 390 ierr = PetscLogEventEnd(CHARACTERISTIC_QueueSetup,0,0,0,0);CHKERRQ(ierr); 391 392 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 393 ierr = CharacteristicSendCoordinatesBegin(c);CHKERRQ(ierr); 394 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 395 396 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);CHKERRQ(ierr); 397 /* Calculate velocity at t_n+1/2 (local values) */ 398 ierr = PetscInfo(PETSC_NULL, "Calculating local velocities at t_{n - 1/2}\n");CHKERRQ(ierr); 399 for (n = 0; n < c->queueSize; n++) { 400 Qi = c->queue[n]; 401 if (c->neighbors[Qi.proc] == rank) { 402 interpIndices[0] = Qi.x; 403 interpIndices[1] = Qi.y; 404 if (c->velocityInterpLocal) { 405 c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx); 406 c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx); 407 } else { 408 c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx); 409 c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx); 410 } 411 Qi.x = 0.5*(velocityValues[0] + velocityValuesOld[0]); 412 Qi.y = 0.5*(velocityValues[1] + velocityValuesOld[1]); 413 } 414 c->queue[n] = Qi; 415 } 416 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);CHKERRQ(ierr); 417 418 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 419 ierr = CharacteristicSendCoordinatesEnd(c);CHKERRQ(ierr); 420 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 421 422 423 /* Calculate velocity at t_n+1/2 (fill remote requests) */ 424 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);CHKERRQ(ierr); 425 ierr = PetscInfo1(PETSC_NULL, "Calculating %d remote velocities at t_{n - 1/2}\n", c->queueRemoteSize);CHKERRQ(ierr); 426 for (n = 0; n < c->queueRemoteSize; n++) { 427 Qi = c->queueRemote[n]; 428 interpIndices[0] = Qi.x; 429 interpIndices[1] = Qi.y; 430 if (c->velocityInterpLocal) { 431 c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx); 432 c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx); 433 } else { 434 c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx); 435 c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx); 436 } 437 Qi.x = 0.5*(velocityValues[0] + velocityValuesOld[0]); 438 Qi.y = 0.5*(velocityValues[1] + velocityValuesOld[1]); 439 c->queueRemote[n] = Qi; 440 } 441 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);CHKERRQ(ierr); 442 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 443 ierr = CharacteristicGetValuesBegin(c);CHKERRQ(ierr); 444 ierr = CharacteristicGetValuesEnd(c);CHKERRQ(ierr); 445 if (c->velocityInterpLocal) { 446 ierr = DMDAVecRestoreArray(c->velocityDA, velocityLocal, &velocityArray);CHKERRQ(ierr); 447 ierr = DMDAVecRestoreArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);CHKERRQ(ierr); 448 ierr = DMRestoreLocalVector(c->velocityDA, &velocityLocal);CHKERRQ(ierr); 449 ierr = DMRestoreLocalVector(c->velocityDA, &velocityLocalOld);CHKERRQ(ierr); 450 } 451 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 452 453 /* ----------------------------------------------------------------------- 454 PART 2, AT t-dt 455 -----------------------------------------------------------------------*/ 456 457 /* GET POSITION AT t_n (local values) */ 458 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 459 ierr = PetscInfo(PETSC_NULL, "Calculating position at t_{n}\n");CHKERRQ(ierr); 460 for (n = 0; n < c->queueSize; n++) { 461 Qi = c->queue[n]; 462 Qi.x = Qi.i - Qi.x*dt; 463 Qi.y = Qi.j - Qi.y*dt; 464 465 /* Determine whether the position at t-dt is local */ 466 Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y); 467 468 /* Check for Periodic boundaries and move all periodic points back onto the domain */ 469 ierr = DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));CHKERRQ(ierr); 470 471 c->queue[n] = Qi; 472 } 473 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 474 475 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 476 ierr = CharacteristicSendCoordinatesBegin(c);CHKERRQ(ierr); 477 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 478 479 /* GET VALUE AT FULL TIME IN THE PAST (LOCAL REQUESTS) */ 480 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 481 if (c->fieldInterpLocal) { 482 ierr = DMGetLocalVector(c->fieldDA, &fieldLocal);CHKERRQ(ierr); 483 ierr = DMGlobalToLocalBegin(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);CHKERRQ(ierr); 484 ierr = DMGlobalToLocalEnd(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);CHKERRQ(ierr); 485 ierr = DMDAVecGetArray(c->fieldDA, fieldLocal, &fieldArray);CHKERRQ(ierr); 486 } 487 ierr = PetscInfo(PETSC_NULL, "Calculating local field at t_{n}\n");CHKERRQ(ierr); 488 for (n = 0; n < c->queueSize; n++) { 489 if (c->neighbors[c->queue[n].proc] == rank) { 490 interpIndices[0] = c->queue[n].x; 491 interpIndices[1] = c->queue[n].y; 492 if (c->fieldInterpLocal) c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx); 493 else c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx); 494 for (comp = 0; comp < c->numFieldComp; comp++) c->queue[n].field[comp] = fieldValues[comp]; 495 } 496 } 497 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 498 499 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 500 ierr = CharacteristicSendCoordinatesEnd(c);CHKERRQ(ierr); 501 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 502 503 /* GET VALUE AT FULL TIME IN THE PAST (REMOTE REQUESTS) */ 504 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeRemote,0,0,0,0);CHKERRQ(ierr); 505 ierr = PetscInfo1(PETSC_NULL, "Calculating %d remote field points at t_{n}\n", c->queueRemoteSize);CHKERRQ(ierr); 506 for (n = 0; n < c->queueRemoteSize; n++) { 507 interpIndices[0] = c->queueRemote[n].x; 508 interpIndices[1] = c->queueRemote[n].y; 509 510 /* for debugging purposes */ 511 if (1) { /* hacked bounds test...let's do better */ 512 PetscScalar im = interpIndices[0]; PetscScalar jm = interpIndices[1]; 513 514 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); 515 } 516 517 if (c->fieldInterpLocal) c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx); 518 else c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx); 519 for (comp = 0; comp < c->numFieldComp; comp++) c->queueRemote[n].field[comp] = fieldValues[comp]; 520 } 521 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeRemote,0,0,0,0);CHKERRQ(ierr); 522 523 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 524 ierr = CharacteristicGetValuesBegin(c);CHKERRQ(ierr); 525 ierr = CharacteristicGetValuesEnd(c);CHKERRQ(ierr); 526 if (c->fieldInterpLocal) { 527 ierr = DMDAVecRestoreArray(c->fieldDA, fieldLocal, &fieldArray);CHKERRQ(ierr); 528 ierr = DMRestoreLocalVector(c->fieldDA, &fieldLocal);CHKERRQ(ierr); 529 } 530 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 531 532 /* Return field of characteristics at t_n-1 */ 533 ierr = PetscLogEventBegin(CHARACTERISTIC_DAUpdate,0,0,0,0);CHKERRQ(ierr); 534 ierr = DMDAGetInfo(c->fieldDA,0,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 535 ierr = DMDAVecGetArray(c->fieldDA, solution, &solArray);CHKERRQ(ierr); 536 for (n = 0; n < c->queueSize; n++) { 537 Qi = c->queue[n]; 538 for (comp = 0; comp < c->numFieldComp; comp++) solArray[Qi.j][Qi.i*dof+c->fieldComp[comp]] = Qi.field[comp]; 539 } 540 ierr = DMDAVecRestoreArray(c->fieldDA, solution, &solArray);CHKERRQ(ierr); 541 ierr = PetscLogEventEnd(CHARACTERISTIC_DAUpdate,0,0,0,0);CHKERRQ(ierr); 542 ierr = PetscLogEventEnd(CHARACTERISTIC_Solve,0,0,0,0);CHKERRQ(ierr); 543 544 /* Cleanup */ 545 ierr = PetscFree(interpIndices);CHKERRQ(ierr); 546 ierr = PetscFree(velocityValues);CHKERRQ(ierr); 547 ierr = PetscFree(velocityValuesOld);CHKERRQ(ierr); 548 ierr = PetscFree(fieldValues);CHKERRQ(ierr); 549 PetscFunctionReturn(0); 550 } 551 552 #undef __FUNCT__ 553 #define __FUNCT__ "CharacteristicSetNeighbors" 554 PetscErrorCode CharacteristicSetNeighbors(Characteristic c, PetscInt numNeighbors, PetscMPIInt neighbors[]) 555 { 556 PetscErrorCode ierr; 557 558 PetscFunctionBegin; 559 c->numNeighbors = numNeighbors; 560 ierr = PetscFree(c->neighbors);CHKERRQ(ierr); 561 ierr = PetscMalloc(numNeighbors * sizeof(PetscMPIInt), &c->neighbors);CHKERRQ(ierr); 562 ierr = PetscMemcpy(c->neighbors, neighbors, numNeighbors * sizeof(PetscMPIInt));CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 566 #undef __FUNCT__ 567 #define __FUNCT__ "CharacteristicAddPoint" 568 PetscErrorCode CharacteristicAddPoint(Characteristic c, CharacteristicPointDA2D *point) 569 { 570 PetscFunctionBegin; 571 if (c->queueSize >= c->queueMax) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Exceeeded maximum queue size %d", c->queueMax); 572 c->queue[c->queueSize++] = *point; 573 PetscFunctionReturn(0); 574 } 575 576 #undef __FUNCT__ 577 #define __FUNCT__ "CharacteristicSendCoordinatesBegin" 578 int CharacteristicSendCoordinatesBegin(Characteristic c) 579 { 580 PetscMPIInt rank, tag = 121; 581 PetscInt i, n; 582 PetscErrorCode ierr; 583 584 PetscFunctionBegin; 585 ierr = MPI_Comm_rank(((PetscObject)c)->comm, &rank);CHKERRQ(ierr); 586 ierr = CharacteristicHeapSort(c, c->queue, c->queueSize);CHKERRQ(ierr); 587 ierr = PetscMemzero(c->needCount, c->numNeighbors * sizeof(PetscInt));CHKERRQ(ierr); 588 for (i = 0; i < c->queueSize; i++) c->needCount[c->queue[i].proc]++; 589 c->fillCount[0] = 0; 590 for (n = 1; n < c->numNeighbors; n++) { 591 ierr = MPI_Irecv(&(c->fillCount[n]), 1, MPIU_INT, c->neighbors[n], tag, ((PetscObject)c)->comm, &(c->request[n-1]));CHKERRQ(ierr); 592 } 593 for (n = 1; n < c->numNeighbors; n++) { 594 ierr = MPI_Send(&(c->needCount[n]), 1, MPIU_INT, c->neighbors[n], tag, ((PetscObject)c)->comm);CHKERRQ(ierr); 595 } 596 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 597 /* Initialize the remote queue */ 598 c->queueLocalMax = c->localOffsets[0] = 0; 599 c->queueRemoteMax = c->remoteOffsets[0] = 0; 600 for (n = 1; n < c->numNeighbors; n++) { 601 c->remoteOffsets[n] = c->queueRemoteMax; 602 c->queueRemoteMax += c->fillCount[n]; 603 c->localOffsets[n] = c->queueLocalMax; 604 c->queueLocalMax += c->needCount[n]; 605 } 606 /* HACK BEGIN */ 607 for (n = 1; n < c->numNeighbors; n++) c->localOffsets[n] += c->needCount[0]; 608 c->needCount[0] = 0; 609 /* HACK END */ 610 if (c->queueRemoteMax) { 611 ierr = PetscMalloc(sizeof(CharacteristicPointDA2D) * c->queueRemoteMax, &c->queueRemote);CHKERRQ(ierr); 612 } else c->queueRemote = PETSC_NULL; 613 c->queueRemoteSize = c->queueRemoteMax; 614 615 /* Send and Receive requests for values at t_n+1/2, giving the coordinates for interpolation */ 616 for (n = 1; n < c->numNeighbors; n++) { 617 ierr = PetscInfo2(PETSC_NULL, "Receiving %d requests for values from proc %d\n", c->fillCount[n], c->neighbors[n]);CHKERRQ(ierr); 618 ierr = MPI_Irecv(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm, &(c->request[n-1]));CHKERRQ(ierr); 619 } 620 for (n = 1; n < c->numNeighbors; n++) { 621 ierr = PetscInfo2(PETSC_NULL, "Sending %d requests for values from proc %d\n", c->needCount[n], c->neighbors[n]);CHKERRQ(ierr); 622 ierr = MPI_Send(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm);CHKERRQ(ierr); 623 } 624 PetscFunctionReturn(0); 625 } 626 627 #undef __FUNCT__ 628 #define __FUNCT__ "CharacteristicSendCoordinatesEnd" 629 PetscErrorCode CharacteristicSendCoordinatesEnd(Characteristic c) 630 { 631 #if 0 632 PetscMPIInt rank; 633 PetscInt n; 634 #endif 635 PetscErrorCode ierr; 636 637 PetscFunctionBegin; 638 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 639 #if 0 640 ierr = MPI_Comm_rank(((PetscObject)c)->comm, &rank);CHKERRQ(ierr); 641 for (n = 0; n < c->queueRemoteSize; n++) { 642 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); 643 } 644 #endif 645 PetscFunctionReturn(0); 646 } 647 648 #undef __FUNCT__ 649 #define __FUNCT__ "CharacteristicGetValuesBegin" 650 PetscErrorCode CharacteristicGetValuesBegin(Characteristic c) 651 { 652 PetscMPIInt tag = 121; 653 PetscInt n; 654 PetscErrorCode ierr; 655 656 PetscFunctionBegin; 657 /* SEND AND RECIEVE FILLED REQUESTS for velocities at t_n+1/2 */ 658 for (n = 1; n < c->numNeighbors; n++) { 659 ierr = MPI_Irecv(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm, &(c->request[n-1]));CHKERRQ(ierr); 660 } 661 for (n = 1; n < c->numNeighbors; n++) { 662 ierr = MPI_Send(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, ((PetscObject)c)->comm);CHKERRQ(ierr); 663 } 664 PetscFunctionReturn(0); 665 } 666 667 #undef __FUNCT__ 668 #define __FUNCT__ "CharacteristicGetValuesEnd" 669 PetscErrorCode CharacteristicGetValuesEnd(Characteristic c) 670 { 671 PetscErrorCode ierr; 672 673 PetscFunctionBegin; 674 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 675 /* Free queue of requests from other procs */ 676 ierr = PetscFree(c->queueRemote);CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 /*---------------------------------------------------------------------*/ 681 #undef __FUNCT__ 682 #define __FUNCT__ "CharacteristicHeapSort" 683 /* 684 Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html 685 */ 686 PetscErrorCode CharacteristicHeapSort(Characteristic c, Queue queue, PetscInt size) 687 /*---------------------------------------------------------------------*/ 688 { 689 PetscErrorCode ierr; 690 CharacteristicPointDA2D temp; 691 PetscInt n; 692 693 PetscFunctionBegin; 694 if (0) { /* Check the order of the queue before sorting */ 695 PetscInfo(PETSC_NULL, "Before Heap sort\n"); 696 for (n=0; n<size; n++) { 697 ierr = PetscInfo2(PETSC_NULL,"%d %d\n",n,queue[n].proc);CHKERRQ(ierr); 698 } 699 } 700 701 /* SORTING PHASE */ 702 for (n = (size / 2)-1; n >= 0; n--) { 703 ierr = CharacteristicSiftDown(c, queue, n, size-1);CHKERRQ(ierr); /* Rich had size-1 here, Matt had size*/ 704 } 705 for (n = size-1; n >= 1; n--) { 706 temp = queue[0]; 707 queue[0] = queue[n]; 708 queue[n] = temp; 709 ierr = CharacteristicSiftDown(c, queue, 0, n-1);CHKERRQ(ierr); 710 } 711 if (0) { /* Check the order of the queue after sorting */ 712 ierr = PetscInfo(PETSC_NULL, "Avter Heap sort\n");CHKERRQ(ierr); 713 for (n=0; n<size; n++) { 714 ierr = PetscInfo2(PETSC_NULL,"%d %d\n",n,queue[n].proc);CHKERRQ(ierr); 715 } 716 } 717 PetscFunctionReturn(0); 718 } 719 720 /*---------------------------------------------------------------------*/ 721 #undef __FUNCT__ 722 #define __FUNCT__ "CharacteristicSiftDown" 723 /* 724 Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html 725 */ 726 PetscErrorCode CharacteristicSiftDown(Characteristic c, Queue queue, PetscInt root, PetscInt bottom) 727 /*---------------------------------------------------------------------*/ 728 { 729 PetscBool done = PETSC_FALSE; 730 PetscInt maxChild; 731 CharacteristicPointDA2D temp; 732 733 PetscFunctionBegin; 734 while ((root*2 <= bottom) && (!done)) { 735 if (root*2 == bottom) maxChild = root * 2; 736 else if (queue[root*2].proc > queue[root*2+1].proc) maxChild = root * 2; 737 else maxChild = root * 2 + 1; 738 739 if (queue[root].proc < queue[maxChild].proc) { 740 temp = queue[root]; 741 queue[root] = queue[maxChild]; 742 queue[maxChild] = temp; 743 root = maxChild; 744 } else done = PETSC_TRUE; 745 } 746 PetscFunctionReturn(0); 747 } 748 749 #undef __FUNCT__ 750 #define __FUNCT__ "DMDAGetNeighborsRank" 751 /* [center, left, top-left, top, top-right, right, bottom-right, bottom, bottom-left] */ 752 PetscErrorCode DMDAGetNeighborsRank(DM da, PetscMPIInt neighbors[]) 753 { 754 DMDABoundaryType bx, by; 755 PetscBool IPeriodic = PETSC_FALSE, JPeriodic = PETSC_FALSE; 756 MPI_Comm comm; 757 PetscMPIInt rank; 758 PetscInt **procs,pi,pj,pim,pip,pjm,pjp,PI,PJ; 759 PetscErrorCode ierr; 760 761 PetscFunctionBegin; 762 ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr); 763 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 764 ierr = DMDAGetInfo(da, 0, 0, 0, 0, &PI,&PJ, 0, 0, 0, &bx, &by,0, 0);CHKERRQ(ierr); 765 766 if (bx == DMDA_BOUNDARY_PERIODIC) IPeriodic = PETSC_TRUE; 767 if (by == DMDA_BOUNDARY_PERIODIC) JPeriodic = PETSC_TRUE; 768 769 neighbors[0] = rank; 770 rank = 0; 771 ierr = PetscMalloc(sizeof(PetscInt*)*PJ,&procs);CHKERRQ(ierr); 772 for (pj=0; pj<PJ; pj++) { 773 ierr = PetscMalloc(sizeof(PetscInt)*PI,&(procs[pj]));CHKERRQ(ierr); 774 for (pi=0; pi<PI; pi++) { 775 procs[pj][pi] = rank; 776 rank++; 777 } 778 } 779 780 pi = neighbors[0] % PI; 781 pj = neighbors[0] / PI; 782 pim = pi-1; if (pim<0) pim=PI-1; 783 pip = (pi+1)%PI; 784 pjm = pj-1; if (pjm<0) pjm=PJ-1; 785 pjp = (pj+1)%PJ; 786 787 neighbors[1] = procs[pj] [pim]; 788 neighbors[2] = procs[pjp][pim]; 789 neighbors[3] = procs[pjp][pi]; 790 neighbors[4] = procs[pjp][pip]; 791 neighbors[5] = procs[pj] [pip]; 792 neighbors[6] = procs[pjm][pip]; 793 neighbors[7] = procs[pjm][pi]; 794 neighbors[8] = procs[pjm][pim]; 795 796 if (!IPeriodic) { 797 if (pi==0) neighbors[1]=neighbors[2]=neighbors[8]=neighbors[0]; 798 if (pi==PI-1) neighbors[4]=neighbors[5]=neighbors[6]=neighbors[0]; 799 } 800 801 if (!JPeriodic) { 802 if (pj==0) neighbors[6]=neighbors[7]=neighbors[8]=neighbors[0]; 803 if (pj==PJ-1) neighbors[2]=neighbors[3]=neighbors[4]=neighbors[0]; 804 } 805 806 for (pj = 0; pj < PJ; pj++) { 807 ierr = PetscFree(procs[pj]);CHKERRQ(ierr); 808 } 809 ierr = PetscFree(procs);CHKERRQ(ierr); 810 PetscFunctionReturn(0); 811 } 812 813 #undef __FUNCT__ 814 #define __FUNCT__ "DMDAGetNeighborRelative" 815 /* 816 SUBDOMAIN NEIGHBORHOOD PROCESS MAP: 817 2 | 3 | 4 818 __|___|__ 819 1 | 0 | 5 820 __|___|__ 821 8 | 7 | 6 822 | | 823 */ 824 PetscInt DMDAGetNeighborRelative(DM da, PassiveReal ir, PassiveReal jr) 825 { 826 DMDALocalInfo info; 827 PassiveReal is,ie,js,je; 828 PetscErrorCode ierr; 829 830 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 831 is = (PassiveReal) info.xs - 0.5; ie = (PassiveReal) info.xs + info.xm - 0.5; 832 js = (PassiveReal) info.ys - 0.5; je = (PassiveReal) info.ys + info.ym - 0.5; 833 834 if (ir >= is && ir <= ie) { /* center column */ 835 if (jr >= js && jr <= je) return 0; 836 else if (jr < js) return 7; 837 else return 3; 838 } else if (ir < is) { /* left column */ 839 if (jr >= js && jr <= je) return 1; 840 else if (jr < js) return 8; 841 else return 2; 842 } else { /* right column */ 843 if (jr >= js && jr <= je) return 5; 844 else if (jr < js) return 6; 845 else return 4; 846 } 847 } 848