ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/SMIPDForceManager.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/integrators/SMIPDForceManager.cpp (file contents):
Revision 3465 by chuckv, Tue Oct 21 16:44:00 2008 UTC vs.
Revision 3473 by chuckv, Fri Nov 14 15:44:34 2008 UTC

# Line 92 | Line 92 | namespace oopse {
92        simError();
93      }
94  
95 +    if (!simParams->haveViscosity()) {
96 +      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
97 +      painCave.isFatal = 1;
98 +      painCave.severity = OOPSE_ERROR;
99 +      simError();
100 +    }
101  
102  
103      
# Line 308 | Line 314 | namespace oopse {
314  
315       /* Compute variance for random forces */
316    
311 –    RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*area/nTriangles)
312 –       /OOPSEConstant::energyConvert;
317  
314 –    gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * area * area * simParams->getDt()) /
315 –      (4.0 * nTriangles * nTriangles* OOPSEConstant::kB*simParams->getTargetTemp()) /OOPSEConstant::energyConvert;
318  
317 –    std::vector<RealType>  randNums = genTriangleForces(nTriangles, sigma_t);
319  
320 +    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);
321 +
322      /* Loop over the mesh faces and apply random force to each of the faces*/
323      
324      
# Line 326 | Line 329 | namespace oopse {
329      
330        Triangle thisTriangle = *face;
331        std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
332 <      
332 >      RealType thisArea = thisTriangle.getArea();
333 >      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
334 >      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
335 >
336        /* Get Random Force */
337        Vector3d unitNormal = thisTriangle.getNormal();
338        unitNormal.normalize();
339 <      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
339 >      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
340        Vector3d centroid = thisTriangle.getCentroid();
341 +      Vector3d facetVel = thisTriangle.getFacetVelocity();
342 +      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/3.14;
343  
344 <      Vector3d langevinForce = randomForce - gamma_t_*thisTriangle.getFacetVelocity();
344 >      RealType f_normal = simParams->getViscosity()*hydroLength*1.439326479e4*OOPSEConstant::energyConvert;
345 >      RealType extPressure = -(targetPressure_ * thisArea);
346 >      RealType randomForce = randNums[thisNumber] * f_normal * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
347 >      RealType dragForce = -f_normal * dot(facetVel, unitNormal);      
348 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
349        
350 +      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
351 +      
352 +      //std::cout << "randomForce " << randomForce << " dragForce " << dragForce <<  " hydro  " << hydroLength << std::endl;
353 +
354 +
355        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
356          if ((*vertex) != NULL){
357            // mass = integrableObject->getMass();
358            Vector3d vertexForce = langevinForce/3.0;
359            (*vertex)->addFrc(vertexForce);
360 <          if (integrableObject->isDirectional()){
360 >
361 >          if ((*vertex)->isDirectional()){
362 >
363              Vector3d vertexPos = (*vertex)->getPos();
364              Vector3d vertexCentroidVector = vertexPos - centroid;
365              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 493 | Line 512 | namespace oopse {
512    }
513      */
514      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
515 <    currSnapshot->setVolume(surfaceMesh_->getVolume());
497 <    
515 >    currSnapshot->setVolume(surfaceMesh_->getVolume());    
516      ForceManager::postCalculation(needStress);  
517    }
518  
# Line 534 | Line 552 | namespace oopse {
552      if (worldRank == 0) {
553   #endif
554        for (int i = 0; i < nTriangles; i++) {
555 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
555 >        //gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
556 >        gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);    
557        }
558   #ifdef IS_MPI
559      }
# Line 549 | Line 568 | namespace oopse {
568        MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
569      }
570   #endif
571 <
571 >  
572      for (int i = 0; i < nTriangles; i++) {
573        gaussRand[i] = gaussRand[i] * variance;
574      }
575 <
575 >  
576      return gaussRand;
577    }
578  

Diff Legend

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