1 2 #include <petsc/private/characteristicimpl.h> /*I "petsccharacteristic.h" I*/ 3 #include <petscdmda.h> 4 #include <petscviewer.h> 5 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; 15 16 PetscErrorCode DMDAGetNeighborsRank(DM, PetscMPIInt []); 17 PetscInt DMDAGetNeighborRelative(DM, PetscReal, PetscReal); 18 PetscErrorCode DMDAMapToPeriodicDomain(DM, PetscScalar []); 19 20 PetscErrorCode CharacteristicHeapSort(Characteristic, Queue, PetscInt); 21 PetscErrorCode CharacteristicSiftDown(Characteristic, Queue, PetscInt, PetscInt); 22 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(PetscObjectComm((PetscObject)c),&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 PetscErrorCode CharacteristicDestroy(Characteristic *c) 46 { 47 PetscErrorCode ierr; 48 49 PetscFunctionBegin; 50 if (!*c) PetscFunctionReturn(0); 51 PetscValidHeaderSpecific(*c, CHARACTERISTIC_CLASSID, 1); 52 if (--((PetscObject)(*c))->refct > 0) PetscFunctionReturn(0); 53 54 if ((*c)->ops->destroy) { 55 ierr = (*(*c)->ops->destroy)((*c));CHKERRQ(ierr); 56 } 57 ierr = MPI_Type_free(&(*c)->itemType);CHKERRQ(ierr); 58 ierr = PetscFree((*c)->queue);CHKERRQ(ierr); 59 ierr = PetscFree((*c)->queueLocal);CHKERRQ(ierr); 60 ierr = PetscFree((*c)->queueRemote);CHKERRQ(ierr); 61 ierr = PetscFree((*c)->neighbors);CHKERRQ(ierr); 62 ierr = PetscFree((*c)->needCount);CHKERRQ(ierr); 63 ierr = PetscFree((*c)->localOffsets);CHKERRQ(ierr); 64 ierr = PetscFree((*c)->fillCount);CHKERRQ(ierr); 65 ierr = PetscFree((*c)->remoteOffsets);CHKERRQ(ierr); 66 ierr = PetscFree((*c)->request);CHKERRQ(ierr); 67 ierr = PetscFree((*c)->status);CHKERRQ(ierr); 68 ierr = PetscHeaderDestroy(c);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 PetscErrorCode CharacteristicCreate(MPI_Comm comm, Characteristic *c) 73 { 74 Characteristic newC; 75 PetscErrorCode ierr; 76 77 PetscFunctionBegin; 78 PetscValidPointer(c, 2); 79 *c = NULL; 80 ierr = CharacteristicInitializePackage();CHKERRQ(ierr); 81 82 ierr = PetscHeaderCreate(newC, CHARACTERISTIC_CLASSID, "Characteristic", "Characteristic", "Characteristic", comm, CharacteristicDestroy, CharacteristicView);CHKERRQ(ierr); 83 *c = newC; 84 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 PetscFunctionReturn(0); 121 } 122 123 /*@C 124 CharacteristicSetType - Builds Characteristic for a particular solver. 125 126 Logically Collective on Characteristic 127 128 Input Parameters: 129 + c - the method of characteristics context 130 - type - a known method 131 132 Options Database Key: 133 . -characteristic_type <method> - Sets the method; use -help for a list 134 of available methods 135 136 Notes: 137 See "include/petsccharacteristic.h" for available methods 138 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. 150 151 Level: intermediate 152 153 .keywords: Characteristic, set, method 154 155 .seealso: CharacteristicType 156 157 @*/ 158 PetscErrorCode CharacteristicSetType(Characteristic c, CharacteristicType type) 159 { 160 PetscErrorCode ierr, (*r)(Characteristic); 161 PetscBool match; 162 163 PetscFunctionBegin; 164 PetscValidHeaderSpecific(c, CHARACTERISTIC_CLASSID, 1); 165 PetscValidCharPointer(type, 2); 166 167 ierr = PetscObjectTypeCompare((PetscObject) c, type, &match);CHKERRQ(ierr); 168 if (match) PetscFunctionReturn(0); 169 170 if (c->data) { 171 /* destroy the old private Characteristic context */ 172 ierr = (*c->ops->destroy)(c);CHKERRQ(ierr); 173 c->ops->destroy = NULL; 174 c->data = 0; 175 } 176 177 ierr = PetscFunctionListFind(CharacteristicList,type,&r);CHKERRQ(ierr); 178 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Characteristic type given: %s", type); 179 c->setupcalled = 0; 180 ierr = (*r)(c);CHKERRQ(ierr); 181 ierr = PetscObjectChangeTypeName((PetscObject) c, type);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 /*@ 186 CharacteristicSetUp - Sets up the internal data structures for the 187 later use of an iterative solver. 188 189 Collective on Characteristic 190 191 Input Parameter: 192 . ksp - iterative context obtained from CharacteristicCreate() 193 194 Level: developer 195 196 .keywords: Characteristic, setup 197 198 .seealso: CharacteristicCreate(), CharacteristicSolve(), CharacteristicDestroy() 199 @*/ 200 PetscErrorCode CharacteristicSetUp(Characteristic c) 201 { 202 PetscErrorCode ierr; 203 204 PetscFunctionBegin; 205 PetscValidHeaderSpecific(c, CHARACTERISTIC_CLASSID, 1); 206 207 if (!((PetscObject)c)->type_name) { 208 ierr = CharacteristicSetType(c, CHARACTERISTICDA);CHKERRQ(ierr); 209 } 210 211 if (c->setupcalled == 2) PetscFunctionReturn(0); 212 213 ierr = PetscLogEventBegin(CHARACTERISTIC_SetUp,c,0,0,0);CHKERRQ(ierr); 214 if (!c->setupcalled) { 215 ierr = (*c->ops->setup)(c);CHKERRQ(ierr); 216 } 217 ierr = PetscLogEventEnd(CHARACTERISTIC_SetUp,c,0,0,0);CHKERRQ(ierr); 218 c->setupcalled = 2; 219 PetscFunctionReturn(0); 220 } 221 222 /*@C 223 CharacteristicRegister - Adds a solver to the method of characteristics package. 224 225 Not Collective 226 227 Input Parameters: 228 + name_solver - name of a new user-defined solver 229 - routine_create - routine to create method context 230 231 Sample usage: 232 .vb 233 CharacteristicRegister("my_char", MyCharCreate); 234 .ve 235 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 245 246 Notes: 247 CharacteristicRegister() may be called multiple times to add several user-defined solvers. 248 249 .keywords: Characteristic, register 250 251 .seealso: CharacteristicRegisterAll(), CharacteristicRegisterDestroy() 252 253 Level: advanced 254 @*/ 255 PetscErrorCode CharacteristicRegister(const char sname[],PetscErrorCode (*function)(Characteristic)) 256 { 257 PetscErrorCode ierr; 258 259 PetscFunctionBegin; 260 ierr = PetscFunctionListAdd(&CharacteristicList,sname,function);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 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 { 266 PetscFunctionBegin; 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 PetscFunctionReturn(0); 275 } 276 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 { 279 PetscFunctionBegin; 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 PetscFunctionReturn(0); 288 } 289 290 PetscErrorCode CharacteristicSetFieldInterpolation(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal[], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx) 291 { 292 PetscFunctionBegin; 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 PetscFunctionReturn(0); 303 } 304 305 PetscErrorCode CharacteristicSetFieldInterpolationLocal(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void*, PetscReal[], PetscInt, PetscInt[], PetscScalar [], void*), void *ctx) 306 { 307 PetscFunctionBegin; 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 PetscFunctionReturn(0); 318 } 319 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; 341 342 PetscFunctionBegin; 343 c->queueSize = 0; 344 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);CHKERRQ(ierr); 345 ierr = DMDAGetNeighborsRank(da, neighbors);CHKERRQ(ierr); 346 ierr = CharacteristicSetNeighbors(c, 9, neighbors);CHKERRQ(ierr); 347 ierr = CharacteristicSetUp(c);CHKERRQ(ierr); 348 /* global and local grid info */ 349 ierr = DMDAGetInfo(da, &dim, &gx, &gy, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);CHKERRQ(ierr); 350 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 351 is = info.xs; ie = info.xs+info.xm; 352 js = info.ys; je = info.ys+info.ym; 353 /* Allocation */ 354 ierr = PetscMalloc1(dim, &interpIndices);CHKERRQ(ierr); 355 ierr = PetscMalloc1(c->numVelocityComp, &velocityValues);CHKERRQ(ierr); 356 ierr = PetscMalloc1(c->numVelocityComp, &velocityValuesOld);CHKERRQ(ierr); 357 ierr = PetscMalloc1(c->numFieldComp, &fieldValues);CHKERRQ(ierr); 358 ierr = PetscLogEventBegin(CHARACTERISTIC_Solve,0,0,0,0);CHKERRQ(ierr); 359 360 /* ----------------------------------------------------------------------- 361 PART 1, AT t-dt/2 362 -----------------------------------------------------------------------*/ 363 ierr = PetscLogEventBegin(CHARACTERISTIC_QueueSetup,0,0,0,0);CHKERRQ(ierr); 364 /* GET POSITION AT HALF TIME IN THE PAST */ 365 if (c->velocityInterpLocal) { 366 ierr = DMGetLocalVector(c->velocityDA, &velocityLocal);CHKERRQ(ierr); 367 ierr = DMGetLocalVector(c->velocityDA, &velocityLocalOld);CHKERRQ(ierr); 368 ierr = DMGlobalToLocalBegin(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);CHKERRQ(ierr); 369 ierr = DMGlobalToLocalEnd(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);CHKERRQ(ierr); 370 ierr = DMGlobalToLocalBegin(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);CHKERRQ(ierr); 371 ierr = DMGlobalToLocalEnd(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);CHKERRQ(ierr); 372 ierr = DMDAVecGetArray(c->velocityDA, velocityLocal, &velocityArray);CHKERRQ(ierr); 373 ierr = DMDAVecGetArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);CHKERRQ(ierr); 374 } 375 ierr = PetscInfo(NULL, "Calculating position at t_{n - 1/2}\n");CHKERRQ(ierr); 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) {ierr = c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr);} 381 else {ierr = c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr);} 382 Qi.x = Qi.i - velocityValues[0]*dt/2.0; 383 Qi.y = Qi.j - velocityValues[1]*dt/2.0; 384 385 /* Determine whether the position at t - dt/2 is local */ 386 Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y); 387 388 /* Check for Periodic boundaries and move all periodic points back onto the domain */ 389 ierr = DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));CHKERRQ(ierr); 390 ierr = CharacteristicAddPoint(c, &Qi);CHKERRQ(ierr); 391 } 392 } 393 ierr = PetscLogEventEnd(CHARACTERISTIC_QueueSetup,0,0,0,0);CHKERRQ(ierr); 394 395 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 396 ierr = CharacteristicSendCoordinatesBegin(c);CHKERRQ(ierr); 397 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 398 399 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);CHKERRQ(ierr); 400 /* Calculate velocity at t_n+1/2 (local values) */ 401 ierr = PetscInfo(NULL, "Calculating local velocities at t_{n - 1/2}\n");CHKERRQ(ierr); 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 ierr = c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 409 ierr = c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 410 } else { 411 ierr = c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 412 ierr = c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 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 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);CHKERRQ(ierr); 420 421 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 422 ierr = CharacteristicSendCoordinatesEnd(c);CHKERRQ(ierr); 423 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 424 425 426 /* Calculate velocity at t_n+1/2 (fill remote requests) */ 427 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);CHKERRQ(ierr); 428 ierr = PetscInfo1(NULL, "Calculating %d remote velocities at t_{n - 1/2}\n", c->queueRemoteSize);CHKERRQ(ierr); 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 ierr = c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 435 ierr = c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 436 } else { 437 ierr = c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 438 ierr = c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 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 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);CHKERRQ(ierr); 445 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 446 ierr = CharacteristicGetValuesBegin(c);CHKERRQ(ierr); 447 ierr = CharacteristicGetValuesEnd(c);CHKERRQ(ierr); 448 if (c->velocityInterpLocal) { 449 ierr = DMDAVecRestoreArray(c->velocityDA, velocityLocal, &velocityArray);CHKERRQ(ierr); 450 ierr = DMDAVecRestoreArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);CHKERRQ(ierr); 451 ierr = DMRestoreLocalVector(c->velocityDA, &velocityLocal);CHKERRQ(ierr); 452 ierr = DMRestoreLocalVector(c->velocityDA, &velocityLocalOld);CHKERRQ(ierr); 453 } 454 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 455 456 /* ----------------------------------------------------------------------- 457 PART 2, AT t-dt 458 -----------------------------------------------------------------------*/ 459 460 /* GET POSITION AT t_n (local values) */ 461 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 462 ierr = PetscInfo(NULL, "Calculating position at t_{n}\n");CHKERRQ(ierr); 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; 467 468 /* Determine whether the position at t-dt is local */ 469 Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y); 470 471 /* Check for Periodic boundaries and move all periodic points back onto the domain */ 472 ierr = DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));CHKERRQ(ierr); 473 474 c->queue[n] = Qi; 475 } 476 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 477 478 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 479 ierr = CharacteristicSendCoordinatesBegin(c);CHKERRQ(ierr); 480 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 481 482 /* GET VALUE AT FULL TIME IN THE PAST (LOCAL REQUESTS) */ 483 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 484 if (c->fieldInterpLocal) { 485 ierr = DMGetLocalVector(c->fieldDA, &fieldLocal);CHKERRQ(ierr); 486 ierr = DMGlobalToLocalBegin(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);CHKERRQ(ierr); 487 ierr = DMGlobalToLocalEnd(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);CHKERRQ(ierr); 488 ierr = DMDAVecGetArray(c->fieldDA, fieldLocal, &fieldArray);CHKERRQ(ierr); 489 } 490 ierr = PetscInfo(NULL, "Calculating local field at t_{n}\n");CHKERRQ(ierr); 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) {ierr = c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 496 else {ierr = c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 497 for (comp = 0; comp < c->numFieldComp; comp++) c->queue[n].field[comp] = fieldValues[comp]; 498 } 499 } 500 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 501 502 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 503 ierr = CharacteristicSendCoordinatesEnd(c);CHKERRQ(ierr); 504 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 505 506 /* GET VALUE AT FULL TIME IN THE PAST (REMOTE REQUESTS) */ 507 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeRemote,0,0,0,0);CHKERRQ(ierr); 508 ierr = PetscInfo1(NULL, "Calculating %d remote field points at t_{n}\n", c->queueRemoteSize);CHKERRQ(ierr); 509 for (n = 0; n < c->queueRemoteSize; n++) { 510 interpIndices[0] = c->queueRemote[n].x; 511 interpIndices[1] = c->queueRemote[n].y; 512 513 /* for debugging purposes */ 514 if (1) { /* hacked bounds test...let's do better */ 515 PetscScalar im = interpIndices[0]; PetscScalar jm = interpIndices[1]; 516 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 } 519 520 if (c->fieldInterpLocal) {ierr = c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 521 else {ierr = c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 522 for (comp = 0; comp < c->numFieldComp; comp++) c->queueRemote[n].field[comp] = fieldValues[comp]; 523 } 524 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeRemote,0,0,0,0);CHKERRQ(ierr); 525 526 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 527 ierr = CharacteristicGetValuesBegin(c);CHKERRQ(ierr); 528 ierr = CharacteristicGetValuesEnd(c);CHKERRQ(ierr); 529 if (c->fieldInterpLocal) { 530 ierr = DMDAVecRestoreArray(c->fieldDA, fieldLocal, &fieldArray);CHKERRQ(ierr); 531 ierr = DMRestoreLocalVector(c->fieldDA, &fieldLocal);CHKERRQ(ierr); 532 } 533 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 534 535 /* Return field of characteristics at t_n-1 */ 536 ierr = PetscLogEventBegin(CHARACTERISTIC_DAUpdate,0,0,0,0);CHKERRQ(ierr); 537 ierr = DMDAGetInfo(c->fieldDA,0,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 538 ierr = DMDAVecGetArray(c->fieldDA, solution, &solArray);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(c->fieldDA, solution, &solArray);CHKERRQ(ierr); 544 ierr = PetscLogEventEnd(CHARACTERISTIC_DAUpdate,0,0,0,0);CHKERRQ(ierr); 545 ierr = PetscLogEventEnd(CHARACTERISTIC_Solve,0,0,0,0);CHKERRQ(ierr); 546 547 /* Cleanup */ 548 ierr = PetscFree(interpIndices);CHKERRQ(ierr); 549 ierr = PetscFree(velocityValues);CHKERRQ(ierr); 550 ierr = PetscFree(velocityValuesOld);CHKERRQ(ierr); 551 ierr = PetscFree(fieldValues);CHKERRQ(ierr); 552 PetscFunctionReturn(0); 553 } 554 555 PetscErrorCode CharacteristicSetNeighbors(Characteristic c, PetscInt numNeighbors, PetscMPIInt neighbors[]) 556 { 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 c->numNeighbors = numNeighbors; 561 ierr = PetscFree(c->neighbors);CHKERRQ(ierr); 562 ierr = PetscMalloc1(numNeighbors, &c->neighbors);CHKERRQ(ierr); 563 ierr = PetscMemcpy(c->neighbors, neighbors, numNeighbors * sizeof(PetscMPIInt));CHKERRQ(ierr); 564 PetscFunctionReturn(0); 565 } 566 567 PetscErrorCode CharacteristicAddPoint(Characteristic c, CharacteristicPointDA2D *point) 568 { 569 PetscFunctionBegin; 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 PetscFunctionReturn(0); 573 } 574 575 int CharacteristicSendCoordinatesBegin(Characteristic c) 576 { 577 PetscMPIInt rank, tag = 121; 578 PetscInt i, n; 579 PetscErrorCode ierr; 580 581 PetscFunctionBegin; 582 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);CHKERRQ(ierr); 583 ierr = CharacteristicHeapSort(c, c->queue, c->queueSize);CHKERRQ(ierr); 584 ierr = PetscMemzero(c->needCount, c->numNeighbors * sizeof(PetscInt));CHKERRQ(ierr); 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 ierr = MPI_Irecv(&(c->fillCount[n]), 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &(c->request[n-1]));CHKERRQ(ierr); 589 } 590 for (n = 1; n < c->numNeighbors; n++) { 591 ierr = MPI_Send(&(c->needCount[n]), 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));CHKERRQ(ierr); 592 } 593 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 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 ierr = PetscMalloc1(c->queueRemoteMax, &c->queueRemote);CHKERRQ(ierr); 609 } else c->queueRemote = NULL; 610 c->queueRemoteSize = c->queueRemoteMax; 611 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 ierr = PetscInfo2(NULL, "Receiving %d requests for values from proc %d\n", c->fillCount[n], c->neighbors[n]);CHKERRQ(ierr); 615 ierr = MPI_Irecv(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &(c->request[n-1]));CHKERRQ(ierr); 616 } 617 for (n = 1; n < c->numNeighbors; n++) { 618 ierr = PetscInfo2(NULL, "Sending %d requests for values from proc %d\n", c->needCount[n], c->neighbors[n]);CHKERRQ(ierr); 619 ierr = MPI_Send(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));CHKERRQ(ierr); 620 } 621 PetscFunctionReturn(0); 622 } 623 624 PetscErrorCode CharacteristicSendCoordinatesEnd(Characteristic c) 625 { 626 #if 0 627 PetscMPIInt rank; 628 PetscInt n; 629 #endif 630 PetscErrorCode ierr; 631 632 PetscFunctionBegin; 633 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 634 #if 0 635 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 641 } 642 643 PetscErrorCode CharacteristicGetValuesBegin(Characteristic c) 644 { 645 PetscMPIInt tag = 121; 646 PetscInt n; 647 PetscErrorCode ierr; 648 649 PetscFunctionBegin; 650 /* SEND AND RECIEVE FILLED REQUESTS for velocities at t_n+1/2 */ 651 for (n = 1; n < c->numNeighbors; n++) { 652 ierr = MPI_Irecv(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &(c->request[n-1]));CHKERRQ(ierr); 653 } 654 for (n = 1; n < c->numNeighbors; n++) { 655 ierr = MPI_Send(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));CHKERRQ(ierr); 656 } 657 PetscFunctionReturn(0); 658 } 659 660 PetscErrorCode CharacteristicGetValuesEnd(Characteristic c) 661 { 662 PetscErrorCode ierr; 663 664 PetscFunctionBegin; 665 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 666 /* Free queue of requests from other procs */ 667 ierr = PetscFree(c->queueRemote);CHKERRQ(ierr); 668 PetscFunctionReturn(0); 669 } 670 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; 681 682 PetscFunctionBegin; 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 ierr = PetscInfo2(NULL,"%d %d\n",n,queue[n].proc);CHKERRQ(ierr); 687 } 688 } 689 690 /* SORTING PHASE */ 691 for (n = (size / 2)-1; n >= 0; n--) { 692 ierr = CharacteristicSiftDown(c, queue, n, size-1);CHKERRQ(ierr); /* 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 ierr = CharacteristicSiftDown(c, queue, 0, n-1);CHKERRQ(ierr); 699 } 700 if (0) { /* Check the order of the queue after sorting */ 701 ierr = PetscInfo(NULL, "Avter Heap sort\n");CHKERRQ(ierr); 702 for (n=0; n<size; n++) { 703 ierr = PetscInfo2(NULL,"%d %d\n",n,queue[n].proc);CHKERRQ(ierr); 704 } 705 } 706 PetscFunctionReturn(0); 707 } 708 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; 719 720 PetscFunctionBegin; 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; 725 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 PetscFunctionReturn(0); 734 } 735 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; 745 746 PetscFunctionBegin; 747 ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr); 748 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 749 ierr = DMDAGetInfo(da, 0, 0, 0, 0, &PI,&PJ, 0, 0, 0, &bx, &by,0, 0);CHKERRQ(ierr); 750 751 if (bx == DM_BOUNDARY_PERIODIC) IPeriodic = PETSC_TRUE; 752 if (by == DM_BOUNDARY_PERIODIC) JPeriodic = PETSC_TRUE; 753 754 neighbors[0] = rank; 755 rank = 0; 756 ierr = PetscMalloc1(PJ,&procs);CHKERRQ(ierr); 757 for (pj=0; pj<PJ; pj++) { 758 ierr = PetscMalloc1(PI,&(procs[pj]));CHKERRQ(ierr); 759 for (pi=0; pi<PI; pi++) { 760 procs[pj][pi] = rank; 761 rank++; 762 } 763 } 764 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; 771 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]; 780 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 } 785 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 } 790 791 for (pj = 0; pj < PJ; pj++) { 792 ierr = PetscFree(procs[pj]);CHKERRQ(ierr); 793 } 794 ierr = PetscFree(procs);CHKERRQ(ierr); 795 PetscFunctionReturn(0); 796 } 797 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; 811 PetscErrorCode ierr; 812 813 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 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; 816 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 } 831