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 3463 by chuckv, Thu Oct 16 18:25:36 2008 UTC vs.
Revision 3473 by chuckv, Fri Nov 14 15:44:34 2008 UTC

# Line 58 | Line 58 | namespace oopse {
58      
59      // Create Hull, Convex Hull for now, other options later.
60      surfaceMesh_ = new ConvexHull();
61 +
62      
63      
64      /* Check that the simulation has target pressure and target
# Line 79 | Line 80 | namespace oopse {
80        painCave.isFatal = 1;
81        simError();
82      } else {
83 <      targetPressure_ = simParams->getTargetPressure();
83 >      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
84      }
85  
86    
# Line 91 | 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 +    
104 +
105 +    //Compute initial hull
106 +    /*
107 +    surfaceMesh_->computeHull(localSites_);
108 +    Area0_ = surfaceMesh_->getArea();
109 +    int nTriangles = surfaceMesh_->getNMeshElements();
110 +    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
111 +    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
112 +      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
113 +    //RealType eta0 = gamma_0/
114 +    */
115  
116      // Build the hydroProp map:
117      std::map<std::string, HydroProp*> hydroPropMap;
# Line 195 | Line 216 | namespace oopse {
216                    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
217                    if (ljData != NULL) {
218                      LJParam ljParam = ljData->getData();
219 <                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
219 >                    currShape = new Sphere(atom->getPos(), 2.0);
220                    } else {
221                      sprintf( painCave.errMsg,
222                               "Can not cast GenericData to LJParam\n");
# Line 207 | Line 228 | namespace oopse {
228                } else {
229                  int aNum = etab.GetAtomicNum((atom->getType()).c_str());
230                  if (aNum != 0) {
231 <                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
231 >                  currShape = new Sphere(atom->getPos(), 2.0);
232                  } else {
233                    sprintf( painCave.errMsg,
234                             "Could not find atom type in default element.txt\n");
# Line 218 | Line 239 | namespace oopse {
239                }
240              }
241            }
242 <          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
242 >          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
243            std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
244            if (iter != hydroPropMap.end())
245              hydroProps_.push_back(iter->second);
# Line 233 | Line 254 | namespace oopse {
254  
255      /* Compute hull first time through to get the area of t=0*/
256  
257 <    /* Build a vector of integrable objects to determine if the are surface atoms */
257 >    //Build a vector of integrable objects to determine if the are surface atoms
258      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
259        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
260             integrableObject = mol->nextIntegrableObject(j)) {  
# Line 241 | Line 262 | namespace oopse {
262        }
263      }
264  
265 <    surfaceMesh_->computeHull(localSites_);
245 <    Area0_ = surfaceMesh_->getArea();
246 <    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
247 <    
265 >
266    }  
267  
268    std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
# Line 288 | Line 306 | namespace oopse {
306       RealType area = surfaceMesh_->getArea();
307       int nSurfaceSDs = surfaceMesh_->getNs();
308  
291 –     /* Compute variance for random forces */
309      
310 <     RealType TD_variance = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
311 <       /OOPSEConstant::energyConvert;
312 <    
313 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
314 <    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),TD_variance);
315 <    
310 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
311 >    int nTriangles = sMesh.size();
312 >
313 >
314 >
315 >     /* Compute variance for random forces */
316 >  
317 >
318 >
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      
325 <    std::vector<Triangle*>::iterator face;
325 >    std::vector<Triangle>::iterator face;
326      std::vector<StuntDouble*>::iterator vertex;
327      int thisNumber = 0;
328      for (face = sMesh.begin(); face != sMesh.end(); ++face){
329      
330 <      Triangle* thisTriangle = *face;
331 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
332 <      
330 >      Triangle thisTriangle = *face;
331 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
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();
337 >      Vector3d unitNormal = thisTriangle.getNormal();
338        unitNormal.normalize();
339 <      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
340 <      Vector3d centroid = thisTriangle->getCentroid();
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 +      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  
361 <         // mass = integrableObject->getMass();
362 <        Vector3d vertexForce = randomForce/3.0;
363 <        (*vertex)->addFrc(vertexForce);
364 <        if (integrableObject->isDirectional()){
365 <          Vector3d vertexPos = (*vertex)->getPos();
366 <          Vector3d vertexCentroidVector = vertexPos - centroid;
367 <          (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
325 <        }
326 <          
361 >          if ((*vertex)->isDirectional()){
362 >
363 >            Vector3d vertexPos = (*vertex)->getPos();
364 >            Vector3d vertexCentroidVector = vertexPos - centroid;
365 >            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
366 >          }
367 >        }  
368        }
369      }
370  
371      /* Now loop over all surface particles and apply the drag*/
372 <
372 >    /*
373      std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
374      for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
375        integrableObject = *vertex;
# Line 341 | Line 382 | namespace oopse {
382          Atrans = A.transpose();
383          Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
384          //apply random force and torque at center of resistance
344 –
345 –        Vector3d randomForceBody;
346 –        Vector3d randomTorqueBody;
347 –        genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
348 –        Vector3d randomForceLab = Atrans * randomForceBody;
349 –        Vector3d randomTorqueLab = Atrans * randomTorqueBody;
350 –        integrableObject->addFrc(randomForceLab);            
351 –        integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
352 –        
353 –
354 –        
385          Mat3x3d I = integrableObject->getI();
386          Vector3d omegaBody;
387          
# Line 434 | Line 464 | namespace oopse {
464  
465              
466        } else {
467 <        //spherical atom
467 >        //spherical atom
468  
439 –        Vector3d randomForce;
440 –        Vector3d randomTorque;
441 –        genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
442 –        integrableObject->addFrc(randomForce);  
443 –        
469          // What remains contains velocity explicitly, but the velocity required
470          // is at the full step: v(t + h), while we have initially the velocity
471          // at the half step: v(t + h/2).  We need to iterate to converge the
# Line 465 | Line 490 | namespace oopse {
490            
491            oldFF = frictionForce;                            
492            frictionForce = -hydroProps_[index]->getXitt() * velStep;
493 <          
493 >          //frictionForce = -gamma_t*velStep;
494            // re-estimate velocities at full-step using friction forces:
495            
496            velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
# Line 482 | Line 507 | namespace oopse {
507          
508          
509        }
510 +  
511        
512 <      
513 <    }
512 >  }
513 >    */
514      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
515 <    currSnapshot->setVolume(surfaceMesh_->getVolume());
490 <    
515 >    currSnapshot->setVolume(surfaceMesh_->getVolume());    
516      ForceManager::postCalculation(needStress);  
517    }
518  
# Line 527 | 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 542 | 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