| 11 |
|
#include "Integrator.hpp" |
| 12 |
|
#include "simError.h" |
| 13 |
|
#include "MatVec3.h" |
| 14 |
+ |
#include "ConstraintManager.hpp" |
| 15 |
+ |
#include "Mat3x3d.hpp" |
| 16 |
|
|
| 17 |
|
#ifdef IS_MPI |
| 18 |
|
#define __C |
| 19 |
|
#include "mpiSimulation.hpp" |
| 20 |
|
#endif // is_mpi |
| 21 |
|
|
| 22 |
+ |
inline double roundMe( double x ){ |
| 23 |
+ |
return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); |
| 24 |
+ |
} |
| 25 |
+ |
|
| 26 |
|
Thermo::Thermo( SimInfo* the_info ) { |
| 27 |
|
info = the_info; |
| 28 |
|
int baseSeed = the_info->getSeed(); |
| 29 |
|
|
| 30 |
|
gaussStream = new gaussianSPRNG( baseSeed ); |
| 31 |
+ |
|
| 32 |
+ |
cpIter = info->consMan->createPairIterator(); |
| 33 |
|
} |
| 34 |
|
|
| 35 |
|
Thermo::~Thermo(){ |
| 36 |
|
delete gaussStream; |
| 37 |
+ |
delete cpIter; |
| 38 |
|
} |
| 39 |
|
|
| 40 |
|
double Thermo::getKinetic(){ |
| 205 |
|
const double e_convert = 4.184e-4; |
| 206 |
|
|
| 207 |
|
double molmass, volume; |
| 208 |
< |
double vcom[3], pcom[3], fcom[3], scaled[3]; |
| 208 |
> |
double vcom[3]; |
| 209 |
|
double p_local[9], p_global[9]; |
| 210 |
< |
int i, j, k, nMols; |
| 202 |
< |
Molecule* molecules; |
| 210 |
> |
int i, j, k; |
| 211 |
|
|
| 204 |
– |
nMols = info->n_mol; |
| 205 |
– |
molecules = info->molecules; |
| 206 |
– |
//tau = info->tau; |
| 207 |
– |
|
| 208 |
– |
// use velocities of molecular centers of mass and molecular masses: |
| 212 |
|
for (i=0; i < 9; i++) { |
| 213 |
|
p_local[i] = 0.0; |
| 214 |
|
p_global[i] = 0.0; |
| 215 |
|
} |
| 216 |
|
|
| 217 |
+ |
// use velocities of integrableObjects and their masses: |
| 218 |
+ |
|
| 219 |
|
for (i=0; i < info->integrableObjects.size(); i++) { |
| 220 |
|
|
| 221 |
|
molmass = info->integrableObjects[i]->getMass(); |
| 222 |
|
|
| 223 |
|
info->integrableObjects[i]->getVel(vcom); |
| 219 |
– |
info->integrableObjects[i]->getPos(pcom); |
| 220 |
– |
info->integrableObjects[i]->getFrc(fcom); |
| 221 |
– |
|
| 222 |
– |
matVecMul3(info->HmatInv, pcom, scaled); |
| 223 |
– |
|
| 224 |
– |
for(j=0; j<3; j++) |
| 225 |
– |
scaled[j] -= roundMe(scaled[j]); |
| 226 |
– |
|
| 227 |
– |
// calc the wrapped real coordinates from the wrapped scaled coordinates |
| 228 |
– |
|
| 229 |
– |
matVecMul3(info->Hmat, scaled, pcom); |
| 224 |
|
|
| 225 |
< |
p_local[0] += molmass * (vcom[0] * vcom[0]) + fcom[0]*pcom[0]*eConvert; |
| 226 |
< |
p_local[1] += molmass * (vcom[0] * vcom[1]) + fcom[0]*pcom[1]*eConvert; |
| 227 |
< |
p_local[2] += molmass * (vcom[0] * vcom[2]) + fcom[0]*pcom[2]*eConvert; |
| 228 |
< |
p_local[3] += molmass * (vcom[1] * vcom[0]) + fcom[1]*pcom[0]*eConvert; |
| 229 |
< |
p_local[4] += molmass * (vcom[1] * vcom[1]) + fcom[1]*pcom[1]*eConvert; |
| 230 |
< |
p_local[5] += molmass * (vcom[1] * vcom[2]) + fcom[1]*pcom[2]*eConvert; |
| 231 |
< |
p_local[6] += molmass * (vcom[2] * vcom[0]) + fcom[2]*pcom[0]*eConvert; |
| 232 |
< |
p_local[7] += molmass * (vcom[2] * vcom[1]) + fcom[2]*pcom[1]*eConvert; |
| 233 |
< |
p_local[8] += molmass * (vcom[2] * vcom[2]) + fcom[2]*pcom[2]*eConvert; |
| 234 |
< |
|
| 225 |
> |
p_local[0] += molmass * (vcom[0] * vcom[0]); |
| 226 |
> |
p_local[1] += molmass * (vcom[0] * vcom[1]); |
| 227 |
> |
p_local[2] += molmass * (vcom[0] * vcom[2]); |
| 228 |
> |
p_local[3] += molmass * (vcom[1] * vcom[0]); |
| 229 |
> |
p_local[4] += molmass * (vcom[1] * vcom[1]); |
| 230 |
> |
p_local[5] += molmass * (vcom[1] * vcom[2]); |
| 231 |
> |
p_local[6] += molmass * (vcom[2] * vcom[0]); |
| 232 |
> |
p_local[7] += molmass * (vcom[2] * vcom[1]); |
| 233 |
> |
p_local[8] += molmass * (vcom[2] * vcom[2]); |
| 234 |
> |
|
| 235 |
|
} |
| 236 |
|
|
| 237 |
|
// Get total for entire system from MPI. |
| 238 |
< |
|
| 238 |
> |
|
| 239 |
|
#ifdef IS_MPI |
| 240 |
|
MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
| 241 |
|
#else |
| 246 |
|
|
| 247 |
|
volume = this->getVolume(); |
| 248 |
|
|
| 249 |
+ |
|
| 250 |
+ |
|
| 251 |
|
for(i = 0; i < 3; i++) { |
| 252 |
|
for (j = 0; j < 3; j++) { |
| 253 |
|
k = 3*i + j; |
| 254 |
< |
press[i][j] = p_global[k] / volume; |
| 259 |
< |
|
| 254 |
> |
press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume; |
| 255 |
|
} |
| 256 |
|
} |
| 257 |
|
} |
| 442 |
|
|
| 443 |
|
info->integrableObjects[vd]->setVel( aVel ); |
| 444 |
|
} |
| 445 |
< |
} |
| 445 |
> |
} |
| 446 |
> |
|
| 447 |
> |
void Thermo::removeAngularMomentum(){ |
| 448 |
> |
Vector3d vcom; |
| 449 |
> |
Vector3d qcom; |
| 450 |
> |
Vector3d pos; |
| 451 |
> |
Vector3d vel; |
| 452 |
> |
double mass; |
| 453 |
> |
double xx; |
| 454 |
> |
double yy; |
| 455 |
> |
double zz; |
| 456 |
> |
double xy; |
| 457 |
> |
double xz; |
| 458 |
> |
double yz; |
| 459 |
> |
Vector3d localAngMom; |
| 460 |
> |
Vector3d angMom; |
| 461 |
> |
Vector3d omega; |
| 462 |
> |
vector<StuntDouble *> integrableObjects; |
| 463 |
> |
double localInertiaVec[9]; |
| 464 |
> |
double inertiaVec[9]; |
| 465 |
> |
vector<Vector3d> qMinusQCom; |
| 466 |
> |
vector<Vector3d> vMinusVCom; |
| 467 |
> |
Mat3x3d inertiaMat; |
| 468 |
> |
Mat3x3d inverseInertiaMat; |
| 469 |
> |
|
| 470 |
> |
integrableObjects = info->integrableObjects; |
| 471 |
> |
qMinusQCom.resize(integrableObjects.size()); |
| 472 |
> |
vMinusVCom.resize(integrableObjects.size()); |
| 473 |
> |
|
| 474 |
> |
getCOM(qcom.vec); |
| 475 |
> |
getCOMVel(vcom.vec); |
| 476 |
> |
|
| 477 |
> |
//initialize components for inertia tensor |
| 478 |
> |
xx = 0.0; |
| 479 |
> |
yy = 0.0; |
| 480 |
> |
zz = 0.0; |
| 481 |
> |
xy = 0.0; |
| 482 |
> |
xz = 0.0; |
| 483 |
> |
yz = 0.0; |
| 484 |
> |
|
| 485 |
> |
//build components of Inertia tensor |
| 486 |
> |
// |
| 487 |
> |
// [ Ixx -Ixy -Ixz ] |
| 488 |
> |
// J = | -Iyx Iyy -Iyz | |
| 489 |
> |
// [ -Izx -Iyz Izz ] |
| 490 |
> |
//See Fowles and Cassidy Chapter 9 or Goldstein Chapter 5 |
| 491 |
> |
for(size_t i = 0; i < integrableObjects.size(); i++){ |
| 492 |
> |
integrableObjects[i]->getPos(pos.vec); |
| 493 |
> |
integrableObjects[i]->getVel(vel.vec); |
| 494 |
> |
mass = integrableObjects[i]->getMass(); |
| 495 |
> |
|
| 496 |
> |
qMinusQCom[i] = pos - qcom; |
| 497 |
> |
info->wrapVector(qMinusQCom[i].vec); |
| 498 |
> |
|
| 499 |
> |
vMinusVCom[i] = vel - vcom; |
| 500 |
> |
|
| 501 |
> |
//compute moment of inertia coefficents |
| 502 |
> |
xx += qMinusQCom[i].x * qMinusQCom[i].x * mass; |
| 503 |
> |
yy += qMinusQCom[i].y * qMinusQCom[i].y * mass; |
| 504 |
> |
zz += qMinusQCom[i].z * qMinusQCom[i].z * mass; |
| 505 |
> |
|
| 506 |
> |
// compute products of inertia |
| 507 |
> |
xy += qMinusQCom[i].x * qMinusQCom[i].y * mass; |
| 508 |
> |
xz += qMinusQCom[i].x * qMinusQCom[i].z * mass; |
| 509 |
> |
yz += qMinusQCom[i].y * qMinusQCom[i].z * mass; |
| 510 |
> |
|
| 511 |
> |
localAngMom += crossProduct(qMinusQCom[i] , vMinusVCom[i] ) * mass; |
| 512 |
> |
|
| 513 |
> |
} |
| 514 |
> |
|
| 515 |
> |
localInertiaVec[0] =yy+zz; |
| 516 |
> |
localInertiaVec[1] = -xy; |
| 517 |
> |
localInertiaVec[2] = -xz; |
| 518 |
> |
localInertiaVec[3] = -xy; |
| 519 |
> |
localInertiaVec[4] = xx+zz; |
| 520 |
> |
localInertiaVec[5] = -yz; |
| 521 |
> |
localInertiaVec[6] = -xz; |
| 522 |
> |
localInertiaVec[7] = -yz; |
| 523 |
> |
localInertiaVec[8] = xx+yy; |
| 524 |
> |
|
| 525 |
> |
//Sum and distribute inertia and angmom arrays |
| 526 |
> |
#ifdef MPI |
| 527 |
> |
|
| 528 |
> |
MPI_Allreduce(localInertiaVec, inertiaVec, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
| 529 |
> |
|
| 530 |
> |
MPI_Allreduce(localAngMom.vec, angMom.vec, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
| 531 |
> |
|
| 532 |
> |
inertiaMat.element[0][0] = inertiaVec[0]; |
| 533 |
> |
inertiaMat.element[0][1] = inertiaVec[1]; |
| 534 |
> |
inertiaMat.element[0][2] = inertiaVec[2]; |
| 535 |
> |
|
| 536 |
> |
inertiaMat.element[1][0] = inertiaVec[3]; |
| 537 |
> |
inertiaMat.element[1][1] = inertiaVec[4]; |
| 538 |
> |
inertiaMat.element[1][2] = inertiaVec[5]; |
| 539 |
> |
|
| 540 |
> |
inertiaMat.element[2][0] = inertiaVec[6]; |
| 541 |
> |
inertiaMat.element[2][1] = inertiaVec[7]; |
| 542 |
> |
inertiaMat.element[2][2] = inertiaVec[8]; |
| 543 |
> |
|
| 544 |
> |
#else |
| 545 |
> |
|
| 546 |
> |
inertiaMat.element[0][0] = localInertiaVec[0]; |
| 547 |
> |
inertiaMat.element[0][1] = localInertiaVec[1]; |
| 548 |
> |
inertiaMat.element[0][2] = localInertiaVec[2]; |
| 549 |
> |
|
| 550 |
> |
inertiaMat.element[1][0] = localInertiaVec[3]; |
| 551 |
> |
inertiaMat.element[1][1] = localInertiaVec[4]; |
| 552 |
> |
inertiaMat.element[1][2] = localInertiaVec[5]; |
| 553 |
> |
|
| 554 |
> |
inertiaMat.element[2][0] = localInertiaVec[6]; |
| 555 |
> |
inertiaMat.element[2][1] = localInertiaVec[7]; |
| 556 |
> |
inertiaMat.element[2][2] = localInertiaVec[8]; |
| 557 |
> |
|
| 558 |
> |
angMom = localAngMom; |
| 559 |
> |
#endif |
| 560 |
> |
|
| 561 |
> |
//invert the moment of inertia tensor by LU-decomposition / backsolving: |
| 562 |
> |
|
| 563 |
> |
inverseInertiaMat = inertiaMat.inverse(); |
| 564 |
> |
|
| 565 |
> |
//calculate the angular velocities: omega = I^-1 . L |
| 566 |
> |
|
| 567 |
> |
omega = inverseInertiaMat * angMom; |
| 568 |
> |
|
| 569 |
> |
//subtract out center of mass velocity and angular momentum from |
| 570 |
> |
//particle velocities |
| 571 |
> |
|
| 572 |
> |
for(size_t i = 0; i < integrableObjects.size(); i++){ |
| 573 |
> |
vel = vMinusVCom[i] - crossProduct(omega, qMinusQCom[i]); |
| 574 |
> |
integrableObjects[i]->setVel(vel.vec); |
| 575 |
> |
} |
| 576 |
> |
} |
| 577 |
> |
|
| 578 |
> |
double Thermo::getConsEnergy(){ |
| 579 |
> |
ConstraintPair* consPair; |
| 580 |
> |
double totConsEnergy; |
| 581 |
> |
double bondLen2; |
| 582 |
> |
double dist; |
| 583 |
> |
double lamda; |
| 584 |
> |
|
| 585 |
> |
totConsEnergy = 0; |
| 586 |
> |
|
| 587 |
> |
for(cpIter->first(); !cpIter->isEnd(); cpIter->next()){ |
| 588 |
> |
consPair = cpIter->currentItem(); |
| 589 |
> |
bondLen2 = consPair->getBondLength2(); |
| 590 |
> |
lamda = consPair->getLamda(); |
| 591 |
> |
//dist = consPair->getDistance(); |
| 592 |
> |
|
| 593 |
> |
//totConsEnergy += lamda * (dist*dist - bondLen2); |
| 594 |
> |
} |
| 595 |
> |
|
| 596 |
> |
return totConsEnergy; |
| 597 |
> |
} |
| 598 |
> |
|
| 599 |
> |
|
