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 = CharacteristicInitializePackage();CHKERRQ(ierr); 261 ierr = PetscFunctionListAdd(&CharacteristicList,sname,function);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 PetscErrorCode CharacteristicSetVelocityInterpolation(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal[], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx) 266 { 267 PetscFunctionBegin; 268 c->velocityDA = da; 269 c->velocity = v; 270 c->velocityOld = vOld; 271 c->numVelocityComp = numComponents; 272 c->velocityComp = components; 273 c->velocityInterp = interp; 274 c->velocityCtx = ctx; 275 PetscFunctionReturn(0); 276 } 277 278 PetscErrorCode CharacteristicSetVelocityInterpolationLocal(Characteristic c, DM da, Vec v, Vec vOld, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void*, PetscReal [], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx) 279 { 280 PetscFunctionBegin; 281 c->velocityDA = da; 282 c->velocity = v; 283 c->velocityOld = vOld; 284 c->numVelocityComp = numComponents; 285 c->velocityComp = components; 286 c->velocityInterpLocal = interp; 287 c->velocityCtx = ctx; 288 PetscFunctionReturn(0); 289 } 290 291 PetscErrorCode CharacteristicSetFieldInterpolation(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(Vec, PetscReal[], PetscInt, PetscInt[], PetscScalar[], void*), void *ctx) 292 { 293 PetscFunctionBegin; 294 #if 0 295 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."); 296 #endif 297 c->fieldDA = da; 298 c->field = v; 299 c->numFieldComp = numComponents; 300 c->fieldComp = components; 301 c->fieldInterp = interp; 302 c->fieldCtx = ctx; 303 PetscFunctionReturn(0); 304 } 305 306 PetscErrorCode CharacteristicSetFieldInterpolationLocal(Characteristic c, DM da, Vec v, PetscInt numComponents, PetscInt components[], PetscErrorCode (*interp)(void*, PetscReal[], PetscInt, PetscInt[], PetscScalar [], void*), void *ctx) 307 { 308 PetscFunctionBegin; 309 #if 0 310 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."); 311 #endif 312 c->fieldDA = da; 313 c->field = v; 314 c->numFieldComp = numComponents; 315 c->fieldComp = components; 316 c->fieldInterpLocal = interp; 317 c->fieldCtx = ctx; 318 PetscFunctionReturn(0); 319 } 320 321 PetscErrorCode CharacteristicSolve(Characteristic c, PetscReal dt, Vec solution) 322 { 323 CharacteristicPointDA2D Qi; 324 DM da = c->velocityDA; 325 Vec velocityLocal, velocityLocalOld; 326 Vec fieldLocal; 327 DMDALocalInfo info; 328 PetscScalar **solArray; 329 void *velocityArray; 330 void *velocityArrayOld; 331 void *fieldArray; 332 PetscScalar *interpIndices; 333 PetscScalar *velocityValues, *velocityValuesOld; 334 PetscScalar *fieldValues; 335 PetscMPIInt rank; 336 PetscInt dim; 337 PetscMPIInt neighbors[9]; 338 PetscInt dof; 339 PetscInt gx, gy; 340 PetscInt n, is, ie, js, je, comp; 341 PetscErrorCode ierr; 342 343 PetscFunctionBegin; 344 c->queueSize = 0; 345 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);CHKERRQ(ierr); 346 ierr = DMDAGetNeighborsRank(da, neighbors);CHKERRQ(ierr); 347 ierr = CharacteristicSetNeighbors(c, 9, neighbors);CHKERRQ(ierr); 348 ierr = CharacteristicSetUp(c);CHKERRQ(ierr); 349 /* global and local grid info */ 350 ierr = DMDAGetInfo(da, &dim, &gx, &gy, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);CHKERRQ(ierr); 351 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 352 is = info.xs; ie = info.xs+info.xm; 353 js = info.ys; je = info.ys+info.ym; 354 /* Allocation */ 355 ierr = PetscMalloc1(dim, &interpIndices);CHKERRQ(ierr); 356 ierr = PetscMalloc1(c->numVelocityComp, &velocityValues);CHKERRQ(ierr); 357 ierr = PetscMalloc1(c->numVelocityComp, &velocityValuesOld);CHKERRQ(ierr); 358 ierr = PetscMalloc1(c->numFieldComp, &fieldValues);CHKERRQ(ierr); 359 ierr = PetscLogEventBegin(CHARACTERISTIC_Solve,0,0,0,0);CHKERRQ(ierr); 360 361 /* ----------------------------------------------------------------------- 362 PART 1, AT t-dt/2 363 -----------------------------------------------------------------------*/ 364 ierr = PetscLogEventBegin(CHARACTERISTIC_QueueSetup,0,0,0,0);CHKERRQ(ierr); 365 /* GET POSITION AT HALF TIME IN THE PAST */ 366 if (c->velocityInterpLocal) { 367 ierr = DMGetLocalVector(c->velocityDA, &velocityLocal);CHKERRQ(ierr); 368 ierr = DMGetLocalVector(c->velocityDA, &velocityLocalOld);CHKERRQ(ierr); 369 ierr = DMGlobalToLocalBegin(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);CHKERRQ(ierr); 370 ierr = DMGlobalToLocalEnd(c->velocityDA, c->velocity, INSERT_VALUES, velocityLocal);CHKERRQ(ierr); 371 ierr = DMGlobalToLocalBegin(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);CHKERRQ(ierr); 372 ierr = DMGlobalToLocalEnd(c->velocityDA, c->velocityOld, INSERT_VALUES, velocityLocalOld);CHKERRQ(ierr); 373 ierr = DMDAVecGetArray(c->velocityDA, velocityLocal, &velocityArray);CHKERRQ(ierr); 374 ierr = DMDAVecGetArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);CHKERRQ(ierr); 375 } 376 ierr = PetscInfo(NULL, "Calculating position at t_{n - 1/2}\n");CHKERRQ(ierr); 377 for (Qi.j = js; Qi.j < je; Qi.j++) { 378 for (Qi.i = is; Qi.i < ie; Qi.i++) { 379 interpIndices[0] = Qi.i; 380 interpIndices[1] = Qi.j; 381 if (c->velocityInterpLocal) {ierr = c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr);} 382 else {ierr = c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr);} 383 Qi.x = Qi.i - velocityValues[0]*dt/2.0; 384 Qi.y = Qi.j - velocityValues[1]*dt/2.0; 385 386 /* Determine whether the position at t - dt/2 is local */ 387 Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y); 388 389 /* Check for Periodic boundaries and move all periodic points back onto the domain */ 390 ierr = DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));CHKERRQ(ierr); 391 ierr = CharacteristicAddPoint(c, &Qi);CHKERRQ(ierr); 392 } 393 } 394 ierr = PetscLogEventEnd(CHARACTERISTIC_QueueSetup,0,0,0,0);CHKERRQ(ierr); 395 396 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 397 ierr = CharacteristicSendCoordinatesBegin(c);CHKERRQ(ierr); 398 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 399 400 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);CHKERRQ(ierr); 401 /* Calculate velocity at t_n+1/2 (local values) */ 402 ierr = PetscInfo(NULL, "Calculating local velocities at t_{n - 1/2}\n");CHKERRQ(ierr); 403 for (n = 0; n < c->queueSize; n++) { 404 Qi = c->queue[n]; 405 if (c->neighbors[Qi.proc] == rank) { 406 interpIndices[0] = Qi.x; 407 interpIndices[1] = Qi.y; 408 if (c->velocityInterpLocal) { 409 ierr = c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 410 ierr = c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 411 } else { 412 ierr = c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 413 ierr = c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 414 } 415 Qi.x = 0.5*(velocityValues[0] + velocityValuesOld[0]); 416 Qi.y = 0.5*(velocityValues[1] + velocityValuesOld[1]); 417 } 418 c->queue[n] = Qi; 419 } 420 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeLocal,0,0,0,0);CHKERRQ(ierr); 421 422 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 423 ierr = CharacteristicSendCoordinatesEnd(c);CHKERRQ(ierr); 424 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 425 426 427 /* Calculate velocity at t_n+1/2 (fill remote requests) */ 428 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);CHKERRQ(ierr); 429 ierr = PetscInfo1(NULL, "Calculating %d remote velocities at t_{n - 1/2}\n", c->queueRemoteSize);CHKERRQ(ierr); 430 for (n = 0; n < c->queueRemoteSize; n++) { 431 Qi = c->queueRemote[n]; 432 interpIndices[0] = Qi.x; 433 interpIndices[1] = Qi.y; 434 if (c->velocityInterpLocal) { 435 ierr = c->velocityInterpLocal(velocityArray, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 436 ierr = c->velocityInterpLocal(velocityArrayOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 437 } else { 438 ierr = c->velocityInterp(c->velocity, interpIndices, c->numVelocityComp, c->velocityComp, velocityValues, c->velocityCtx);CHKERRQ(ierr); 439 ierr = c->velocityInterp(c->velocityOld, interpIndices, c->numVelocityComp, c->velocityComp, velocityValuesOld, c->velocityCtx);CHKERRQ(ierr); 440 } 441 Qi.x = 0.5*(velocityValues[0] + velocityValuesOld[0]); 442 Qi.y = 0.5*(velocityValues[1] + velocityValuesOld[1]); 443 c->queueRemote[n] = Qi; 444 } 445 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeRemote,0,0,0,0);CHKERRQ(ierr); 446 ierr = PetscLogEventBegin(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 447 ierr = CharacteristicGetValuesBegin(c);CHKERRQ(ierr); 448 ierr = CharacteristicGetValuesEnd(c);CHKERRQ(ierr); 449 if (c->velocityInterpLocal) { 450 ierr = DMDAVecRestoreArray(c->velocityDA, velocityLocal, &velocityArray);CHKERRQ(ierr); 451 ierr = DMDAVecRestoreArray(c->velocityDA, velocityLocalOld, &velocityArrayOld);CHKERRQ(ierr); 452 ierr = DMRestoreLocalVector(c->velocityDA, &velocityLocal);CHKERRQ(ierr); 453 ierr = DMRestoreLocalVector(c->velocityDA, &velocityLocalOld);CHKERRQ(ierr); 454 } 455 ierr = PetscLogEventEnd(CHARACTERISTIC_HalfTimeExchange,0,0,0,0);CHKERRQ(ierr); 456 457 /* ----------------------------------------------------------------------- 458 PART 2, AT t-dt 459 -----------------------------------------------------------------------*/ 460 461 /* GET POSITION AT t_n (local values) */ 462 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 463 ierr = PetscInfo(NULL, "Calculating position at t_{n}\n");CHKERRQ(ierr); 464 for (n = 0; n < c->queueSize; n++) { 465 Qi = c->queue[n]; 466 Qi.x = Qi.i - Qi.x*dt; 467 Qi.y = Qi.j - Qi.y*dt; 468 469 /* Determine whether the position at t-dt is local */ 470 Qi.proc = DMDAGetNeighborRelative(da, Qi.x, Qi.y); 471 472 /* Check for Periodic boundaries and move all periodic points back onto the domain */ 473 ierr = DMDAMapCoordsToPeriodicDomain(da,&(Qi.x),&(Qi.y));CHKERRQ(ierr); 474 475 c->queue[n] = Qi; 476 } 477 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 478 479 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 480 ierr = CharacteristicSendCoordinatesBegin(c);CHKERRQ(ierr); 481 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 482 483 /* GET VALUE AT FULL TIME IN THE PAST (LOCAL REQUESTS) */ 484 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 485 if (c->fieldInterpLocal) { 486 ierr = DMGetLocalVector(c->fieldDA, &fieldLocal);CHKERRQ(ierr); 487 ierr = DMGlobalToLocalBegin(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);CHKERRQ(ierr); 488 ierr = DMGlobalToLocalEnd(c->fieldDA, c->field, INSERT_VALUES, fieldLocal);CHKERRQ(ierr); 489 ierr = DMDAVecGetArray(c->fieldDA, fieldLocal, &fieldArray);CHKERRQ(ierr); 490 } 491 ierr = PetscInfo(NULL, "Calculating local field at t_{n}\n");CHKERRQ(ierr); 492 for (n = 0; n < c->queueSize; n++) { 493 if (c->neighbors[c->queue[n].proc] == rank) { 494 interpIndices[0] = c->queue[n].x; 495 interpIndices[1] = c->queue[n].y; 496 if (c->fieldInterpLocal) {ierr = c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 497 else {ierr = c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 498 for (comp = 0; comp < c->numFieldComp; comp++) c->queue[n].field[comp] = fieldValues[comp]; 499 } 500 } 501 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeLocal,0,0,0,0);CHKERRQ(ierr); 502 503 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 504 ierr = CharacteristicSendCoordinatesEnd(c);CHKERRQ(ierr); 505 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 506 507 /* GET VALUE AT FULL TIME IN THE PAST (REMOTE REQUESTS) */ 508 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeRemote,0,0,0,0);CHKERRQ(ierr); 509 ierr = PetscInfo1(NULL, "Calculating %d remote field points at t_{n}\n", c->queueRemoteSize);CHKERRQ(ierr); 510 for (n = 0; n < c->queueRemoteSize; n++) { 511 interpIndices[0] = c->queueRemote[n].x; 512 interpIndices[1] = c->queueRemote[n].y; 513 514 /* for debugging purposes */ 515 if (1) { /* hacked bounds test...let's do better */ 516 PetscScalar im = interpIndices[0]; PetscScalar jm = interpIndices[1]; 517 518 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); 519 } 520 521 if (c->fieldInterpLocal) {ierr = c->fieldInterpLocal(fieldArray, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 522 else {ierr = c->fieldInterp(c->field, interpIndices, c->numFieldComp, c->fieldComp, fieldValues, c->fieldCtx);CHKERRQ(ierr);} 523 for (comp = 0; comp < c->numFieldComp; comp++) c->queueRemote[n].field[comp] = fieldValues[comp]; 524 } 525 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeRemote,0,0,0,0);CHKERRQ(ierr); 526 527 ierr = PetscLogEventBegin(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 528 ierr = CharacteristicGetValuesBegin(c);CHKERRQ(ierr); 529 ierr = CharacteristicGetValuesEnd(c);CHKERRQ(ierr); 530 if (c->fieldInterpLocal) { 531 ierr = DMDAVecRestoreArray(c->fieldDA, fieldLocal, &fieldArray);CHKERRQ(ierr); 532 ierr = DMRestoreLocalVector(c->fieldDA, &fieldLocal);CHKERRQ(ierr); 533 } 534 ierr = PetscLogEventEnd(CHARACTERISTIC_FullTimeExchange,0,0,0,0);CHKERRQ(ierr); 535 536 /* Return field of characteristics at t_n-1 */ 537 ierr = PetscLogEventBegin(CHARACTERISTIC_DAUpdate,0,0,0,0);CHKERRQ(ierr); 538 ierr = DMDAGetInfo(c->fieldDA,0,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 539 ierr = DMDAVecGetArray(c->fieldDA, solution, &solArray);CHKERRQ(ierr); 540 for (n = 0; n < c->queueSize; n++) { 541 Qi = c->queue[n]; 542 for (comp = 0; comp < c->numFieldComp; comp++) solArray[Qi.j][Qi.i*dof+c->fieldComp[comp]] = Qi.field[comp]; 543 } 544 ierr = DMDAVecRestoreArray(c->fieldDA, solution, &solArray);CHKERRQ(ierr); 545 ierr = PetscLogEventEnd(CHARACTERISTIC_DAUpdate,0,0,0,0);CHKERRQ(ierr); 546 ierr = PetscLogEventEnd(CHARACTERISTIC_Solve,0,0,0,0);CHKERRQ(ierr); 547 548 /* Cleanup */ 549 ierr = PetscFree(interpIndices);CHKERRQ(ierr); 550 ierr = PetscFree(velocityValues);CHKERRQ(ierr); 551 ierr = PetscFree(velocityValuesOld);CHKERRQ(ierr); 552 ierr = PetscFree(fieldValues);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 PetscErrorCode CharacteristicSetNeighbors(Characteristic c, PetscInt numNeighbors, PetscMPIInt neighbors[]) 557 { 558 PetscErrorCode ierr; 559 560 PetscFunctionBegin; 561 c->numNeighbors = numNeighbors; 562 ierr = PetscFree(c->neighbors);CHKERRQ(ierr); 563 ierr = PetscMalloc1(numNeighbors, &c->neighbors);CHKERRQ(ierr); 564 ierr = PetscMemcpy(c->neighbors, neighbors, numNeighbors * sizeof(PetscMPIInt));CHKERRQ(ierr); 565 PetscFunctionReturn(0); 566 } 567 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 int CharacteristicSendCoordinatesBegin(Characteristic c) 577 { 578 PetscMPIInt rank, tag = 121; 579 PetscInt i, n; 580 PetscErrorCode ierr; 581 582 PetscFunctionBegin; 583 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);CHKERRQ(ierr); 584 ierr = CharacteristicHeapSort(c, c->queue, c->queueSize);CHKERRQ(ierr); 585 ierr = PetscMemzero(c->needCount, c->numNeighbors * sizeof(PetscInt));CHKERRQ(ierr); 586 for (i = 0; i < c->queueSize; i++) c->needCount[c->queue[i].proc]++; 587 c->fillCount[0] = 0; 588 for (n = 1; n < c->numNeighbors; n++) { 589 ierr = MPI_Irecv(&(c->fillCount[n]), 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c), &(c->request[n-1]));CHKERRQ(ierr); 590 } 591 for (n = 1; n < c->numNeighbors; n++) { 592 ierr = MPI_Send(&(c->needCount[n]), 1, MPIU_INT, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));CHKERRQ(ierr); 593 } 594 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 595 /* Initialize the remote queue */ 596 c->queueLocalMax = c->localOffsets[0] = 0; 597 c->queueRemoteMax = c->remoteOffsets[0] = 0; 598 for (n = 1; n < c->numNeighbors; n++) { 599 c->remoteOffsets[n] = c->queueRemoteMax; 600 c->queueRemoteMax += c->fillCount[n]; 601 c->localOffsets[n] = c->queueLocalMax; 602 c->queueLocalMax += c->needCount[n]; 603 } 604 /* HACK BEGIN */ 605 for (n = 1; n < c->numNeighbors; n++) c->localOffsets[n] += c->needCount[0]; 606 c->needCount[0] = 0; 607 /* HACK END */ 608 if (c->queueRemoteMax) { 609 ierr = PetscMalloc1(c->queueRemoteMax, &c->queueRemote);CHKERRQ(ierr); 610 } else c->queueRemote = NULL; 611 c->queueRemoteSize = c->queueRemoteMax; 612 613 /* Send and Receive requests for values at t_n+1/2, giving the coordinates for interpolation */ 614 for (n = 1; n < c->numNeighbors; n++) { 615 ierr = PetscInfo2(NULL, "Receiving %d requests for values from proc %d\n", c->fillCount[n], c->neighbors[n]);CHKERRQ(ierr); 616 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); 617 } 618 for (n = 1; n < c->numNeighbors; n++) { 619 ierr = PetscInfo2(NULL, "Sending %d requests for values from proc %d\n", c->needCount[n], c->neighbors[n]);CHKERRQ(ierr); 620 ierr = MPI_Send(&(c->queue[c->localOffsets[n]]), c->needCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));CHKERRQ(ierr); 621 } 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode CharacteristicSendCoordinatesEnd(Characteristic c) 626 { 627 #if 0 628 PetscMPIInt rank; 629 PetscInt n; 630 #endif 631 PetscErrorCode ierr; 632 633 PetscFunctionBegin; 634 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 635 #if 0 636 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)c), &rank);CHKERRQ(ierr); 637 for (n = 0; n < c->queueRemoteSize; n++) { 638 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); 639 } 640 #endif 641 PetscFunctionReturn(0); 642 } 643 644 PetscErrorCode CharacteristicGetValuesBegin(Characteristic c) 645 { 646 PetscMPIInt tag = 121; 647 PetscInt n; 648 PetscErrorCode ierr; 649 650 PetscFunctionBegin; 651 /* SEND AND RECIEVE FILLED REQUESTS for velocities at t_n+1/2 */ 652 for (n = 1; n < c->numNeighbors; n++) { 653 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); 654 } 655 for (n = 1; n < c->numNeighbors; n++) { 656 ierr = MPI_Send(&(c->queueRemote[c->remoteOffsets[n]]), c->fillCount[n], c->itemType, c->neighbors[n], tag, PetscObjectComm((PetscObject)c));CHKERRQ(ierr); 657 } 658 PetscFunctionReturn(0); 659 } 660 661 PetscErrorCode CharacteristicGetValuesEnd(Characteristic c) 662 { 663 PetscErrorCode ierr; 664 665 PetscFunctionBegin; 666 ierr = MPI_Waitall(c->numNeighbors-1, c->request, c->status);CHKERRQ(ierr); 667 /* Free queue of requests from other procs */ 668 ierr = PetscFree(c->queueRemote);CHKERRQ(ierr); 669 PetscFunctionReturn(0); 670 } 671 672 /*---------------------------------------------------------------------*/ 673 /* 674 Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html 675 */ 676 PetscErrorCode CharacteristicHeapSort(Characteristic c, Queue queue, PetscInt size) 677 /*---------------------------------------------------------------------*/ 678 { 679 PetscErrorCode ierr; 680 CharacteristicPointDA2D temp; 681 PetscInt n; 682 683 PetscFunctionBegin; 684 if (0) { /* Check the order of the queue before sorting */ 685 ierr = PetscInfo(NULL, "Before Heap sort\n");CHKERRQ(ierr); 686 for (n=0; n<size; n++) { 687 ierr = PetscInfo2(NULL,"%d %d\n",n,queue[n].proc);CHKERRQ(ierr); 688 } 689 } 690 691 /* SORTING PHASE */ 692 for (n = (size / 2)-1; n >= 0; n--) { 693 ierr = CharacteristicSiftDown(c, queue, n, size-1);CHKERRQ(ierr); /* Rich had size-1 here, Matt had size*/ 694 } 695 for (n = size-1; n >= 1; n--) { 696 temp = queue[0]; 697 queue[0] = queue[n]; 698 queue[n] = temp; 699 ierr = CharacteristicSiftDown(c, queue, 0, n-1);CHKERRQ(ierr); 700 } 701 if (0) { /* Check the order of the queue after sorting */ 702 ierr = PetscInfo(NULL, "Avter Heap sort\n");CHKERRQ(ierr); 703 for (n=0; n<size; n++) { 704 ierr = PetscInfo2(NULL,"%d %d\n",n,queue[n].proc);CHKERRQ(ierr); 705 } 706 } 707 PetscFunctionReturn(0); 708 } 709 710 /*---------------------------------------------------------------------*/ 711 /* 712 Based on code from http://linux.wku.edu/~lamonml/algor/sort/heap.html 713 */ 714 PetscErrorCode CharacteristicSiftDown(Characteristic c, Queue queue, PetscInt root, PetscInt bottom) 715 /*---------------------------------------------------------------------*/ 716 { 717 PetscBool done = PETSC_FALSE; 718 PetscInt maxChild; 719 CharacteristicPointDA2D temp; 720 721 PetscFunctionBegin; 722 while ((root*2 <= bottom) && (!done)) { 723 if (root*2 == bottom) maxChild = root * 2; 724 else if (queue[root*2].proc > queue[root*2+1].proc) maxChild = root * 2; 725 else maxChild = root * 2 + 1; 726 727 if (queue[root].proc < queue[maxChild].proc) { 728 temp = queue[root]; 729 queue[root] = queue[maxChild]; 730 queue[maxChild] = temp; 731 root = maxChild; 732 } else done = PETSC_TRUE; 733 } 734 PetscFunctionReturn(0); 735 } 736 737 /* [center, left, top-left, top, top-right, right, bottom-right, bottom, bottom-left] */ 738 PetscErrorCode DMDAGetNeighborsRank(DM da, PetscMPIInt neighbors[]) 739 { 740 DMBoundaryType bx, by; 741 PetscBool IPeriodic = PETSC_FALSE, JPeriodic = PETSC_FALSE; 742 MPI_Comm comm; 743 PetscMPIInt rank; 744 PetscInt **procs,pi,pj,pim,pip,pjm,pjp,PI,PJ; 745 PetscErrorCode ierr; 746 747 PetscFunctionBegin; 748 ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr); 749 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 750 ierr = DMDAGetInfo(da, 0, 0, 0, 0, &PI,&PJ, 0, 0, 0, &bx, &by,0, 0);CHKERRQ(ierr); 751 752 if (bx == DM_BOUNDARY_PERIODIC) IPeriodic = PETSC_TRUE; 753 if (by == DM_BOUNDARY_PERIODIC) JPeriodic = PETSC_TRUE; 754 755 neighbors[0] = rank; 756 rank = 0; 757 ierr = PetscMalloc1(PJ,&procs);CHKERRQ(ierr); 758 for (pj=0; pj<PJ; pj++) { 759 ierr = PetscMalloc1(PI,&(procs[pj]));CHKERRQ(ierr); 760 for (pi=0; pi<PI; pi++) { 761 procs[pj][pi] = rank; 762 rank++; 763 } 764 } 765 766 pi = neighbors[0] % PI; 767 pj = neighbors[0] / PI; 768 pim = pi-1; if (pim<0) pim=PI-1; 769 pip = (pi+1)%PI; 770 pjm = pj-1; if (pjm<0) pjm=PJ-1; 771 pjp = (pj+1)%PJ; 772 773 neighbors[1] = procs[pj] [pim]; 774 neighbors[2] = procs[pjp][pim]; 775 neighbors[3] = procs[pjp][pi]; 776 neighbors[4] = procs[pjp][pip]; 777 neighbors[5] = procs[pj] [pip]; 778 neighbors[6] = procs[pjm][pip]; 779 neighbors[7] = procs[pjm][pi]; 780 neighbors[8] = procs[pjm][pim]; 781 782 if (!IPeriodic) { 783 if (pi==0) neighbors[1]=neighbors[2]=neighbors[8]=neighbors[0]; 784 if (pi==PI-1) neighbors[4]=neighbors[5]=neighbors[6]=neighbors[0]; 785 } 786 787 if (!JPeriodic) { 788 if (pj==0) neighbors[6]=neighbors[7]=neighbors[8]=neighbors[0]; 789 if (pj==PJ-1) neighbors[2]=neighbors[3]=neighbors[4]=neighbors[0]; 790 } 791 792 for (pj = 0; pj < PJ; pj++) { 793 ierr = PetscFree(procs[pj]);CHKERRQ(ierr); 794 } 795 ierr = PetscFree(procs);CHKERRQ(ierr); 796 PetscFunctionReturn(0); 797 } 798 799 /* 800 SUBDOMAIN NEIGHBORHOOD PROCESS MAP: 801 2 | 3 | 4 802 __|___|__ 803 1 | 0 | 5 804 __|___|__ 805 8 | 7 | 6 806 | | 807 */ 808 PetscInt DMDAGetNeighborRelative(DM da, PetscReal ir, PetscReal jr) 809 { 810 DMDALocalInfo info; 811 PetscReal is,ie,js,je; 812 PetscErrorCode ierr; 813 814 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 815 is = (PetscReal) info.xs - 0.5; ie = (PetscReal) info.xs + info.xm - 0.5; 816 js = (PetscReal) info.ys - 0.5; je = (PetscReal) info.ys + info.ym - 0.5; 817 818 if (ir >= is && ir <= ie) { /* center column */ 819 if (jr >= js && jr <= je) return 0; 820 else if (jr < js) return 7; 821 else return 3; 822 } else if (ir < is) { /* left column */ 823 if (jr >= js && jr <= je) return 1; 824 else if (jr < js) return 8; 825 else return 2; 826 } else { /* right column */ 827 if (jr >= js && jr <= je) return 5; 828 else if (jr < js) return 6; 829 else return 4; 830 } 831 } 832