ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Thermo.cpp
(Generate patch)

Comparing branches/development/src/brains/Thermo.cpp (file contents):
Revision 1715 by gezelter, Tue May 22 21:55:31 2012 UTC vs.
Revision 1723 by gezelter, Thu May 24 20:59:54 2012 UTC

# Line 246 | Line 246 | namespace OpenMD {
246  
247      RealType volume = this->getVolume();
248      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
249 <    Mat3x3d tau = curSnapshot->getTau();
249 >    Mat3x3d stressTensor = curSnapshot->getStressTensor();
250  
251 <    pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
251 >    pressureTensor =  (p_global +
252 >                       PhysicalConstants::energyConvert * stressTensor)/volume;
253      
254      return pressureTensor;
255    }
# Line 285 | Line 286 | namespace OpenMD {
286      }
287  
288      Globals* simParams = info_->getSimParams();
289 +    // grab the heat flux if desired
290 +    if (simParams->havePrintHeatFlux()) {
291 +      if (simParams->getPrintHeatFlux()){
292 +        Vector3d heatFlux = getHeatFlux();
293 +        stat[Stats::HEATFLUX_X] = heatFlux(0);
294 +        stat[Stats::HEATFLUX_Y] = heatFlux(1);
295 +        stat[Stats::HEATFLUX_Z] = heatFlux(2);
296 +      }
297 +    }
298  
299      if (simParams->haveTaggedAtomPair() &&
300          simParams->havePrintTaggedPairDistance()) {
# Line 457 | Line 467 | namespace OpenMD {
467      boxDipole += (pPos - nPos) * chg_value;
468  
469      return boxDipole;
470 +  }
471 +
472 +  // Returns the Heat Flux Vector for the system
473 +  Vector3d Thermo::getHeatFlux(){
474 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
475 +    SimInfo::MoleculeIterator miter;
476 +    std::vector<StuntDouble*>::iterator iiter;
477 +    Molecule* mol;
478 +    StuntDouble* integrableObject;    
479 +    RigidBody::AtomIterator ai;
480 +    Atom* atom;      
481 +    Vector3d vel;
482 +    Vector3d angMom;
483 +    Mat3x3d I;
484 +    int i;
485 +    int j;
486 +    int k;
487 +    RealType mass;
488 +
489 +    Vector3d x_a;
490 +    RealType kinetic;
491 +    RealType potential;
492 +    RealType eatom;
493 +    RealType AvgE_a_ = 0;
494 +    // Convective portion of the heat flux
495 +    Vector3d heatFluxJc = V3Zero;
496 +
497 +    /* Calculate convective portion of the heat flux */
498 +    for (mol = info_->beginMolecule(miter); mol != NULL;
499 +         mol = info_->nextMolecule(miter)) {
500 +      
501 +      for (integrableObject = mol->beginIntegrableObject(iiter);
502 +           integrableObject != NULL;
503 +           integrableObject = mol->nextIntegrableObject(iiter)) {
504 +        
505 +        mass = integrableObject->getMass();
506 +        vel = integrableObject->getVel();
507 +
508 +        kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
509 +        
510 +        if (integrableObject->isDirectional()) {
511 +          angMom = integrableObject->getJ();
512 +          I = integrableObject->getI();
513 +
514 +          if (integrableObject->isLinear()) {
515 +            i = integrableObject->linearAxis();
516 +            j = (i + 1) % 3;
517 +            k = (i + 2) % 3;
518 +            kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
519 +          } else {                        
520 +            kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
521 +              + angMom[2]*angMom[2]/I(2, 2);
522 +          }
523 +        }
524 +
525 +        potential = 0.0;
526 +
527 +        if (integrableObject->isRigidBody()) {
528 +          RigidBody* rb = dynamic_cast<RigidBody*>(integrableObject);
529 +          for (atom = rb->beginAtom(ai); atom != NULL;
530 +               atom = rb->nextAtom(ai)) {
531 +            potential +=  atom->getParticlePot();
532 +          }          
533 +        } else {
534 +          potential = integrableObject->getParticlePot();
535 +          cerr << "ppot = "  << potential << "\n";
536 +        }
537 +
538 +        potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2
539 +        // The potential may not be a 1/2 factor
540 +        eatom = (kinetic + potential)/2.0;  // amu A^2/fs^2
541 +        heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3
542 +        heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3
543 +        heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3
544 +      }
545 +    }
546 +
547 +    std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl;
548 +
549 +    /* The J_v vector is reduced in fortan so everyone has the global
550 +     *  Jv. Jc is computed over the local atoms and must be reduced
551 +     *  among all processors.
552 +     */
553 + #ifdef IS_MPI
554 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
555 +                              MPI::SUM);
556 + #endif
557 +    
558 +    // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
559 +
560 +    Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() *
561 +      PhysicalConstants::energyConvert;
562 +    
563 +    std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl;
564 +    std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl;
565 +    
566 +    // Correct for the fact the flux is 1/V (Jc + Jv)
567 +    return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
568    }
569   } //end namespace OpenMD

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines