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