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 3482 by chuckv, Fri Dec 5 16:20:39 2008 UTC vs.
Revision 3515 by gezelter, Mon Jul 20 17:50:40 2009 UTC

# Line 41 | Line 41
41   #include <fstream>
42   #include <iostream>
43   #include "integrators/SMIPDForceManager.hpp"
44 – #include "math/CholeskyDecomposition.hpp"
44   #include "utils/OOPSEConstant.hpp"
46 – #include "hydrodynamics/Sphere.hpp"
47 – #include "hydrodynamics/Ellipsoid.hpp"
48 – #include "utils/ElementsTable.hpp"
45   #include "math/ConvexHull.hpp"
46   #include "math/Triangle.hpp"
47  
48   namespace oopse {
49  
50 <  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info),
51 <                                                        forceTolerance_(1e-6),
56 <                                                        maxIterNum_(4) {
50 >  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
51 >
52      simParams = info->getSimParams();
53 +    thermo = new Thermo(info);
54      veloMunge = new Velocitizer(info);
55      
56      // Create Hull, Convex Hull for now, other options later.
# Line 90 | Line 86 | namespace oopse {
86      if (simParams->getUsePeriodicBoundaryConditions()) {
87        sprintf(painCave.errMsg,
88                "SMIPDynamics error: You can't use the SMIPD integrator\n"
89 <              "   with periodic boundary conditions!\n");    
89 >              "   with periodic boundary conditions!\n");    
90        painCave.isFatal = 1;
91        simError();
92      }
93      
94 <    if (!simParams->haveViscosity()) {
94 >    if (!simParams->haveThermalConductivity()) {
95        sprintf(painCave.errMsg,
96                "SMIPDynamics error: You can't use the SMIPD integrator\n"
97 <              "   without a viscosity!\n");
97 >              "   without a thermalConductivity!\n");
98        painCave.isFatal = 1;
99        painCave.severity = OOPSE_ERROR;
100        simError();
101      }else{
102 <      viscosity_ = simParams->getViscosity();
102 >      thermalConductivity_ = simParams->getThermalConductivity();
103      }
104 +
105 +    if (!simParams->haveThermalLength()) {
106 +      sprintf(painCave.errMsg,
107 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
108 +              "   without a thermalLength!\n");
109 +      painCave.isFatal = 1;
110 +      painCave.severity = OOPSE_ERROR;
111 +      simError();
112 +    }else{
113 +      thermalLength_ = simParams->getThermalLength();
114 +    }
115      
116      dt_ = simParams->getDt();
117    
118      variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
119  
120 <    // Build the hydroProp map:
121 <    std::map<std::string, HydroProp*> hydroPropMap;
115 <
120 >    // Build a vector of integrable objects to determine if the are
121 >    // surface atoms
122      Molecule* mol;
123      StuntDouble* integrableObject;
124      SimInfo::MoleculeIterator i;
125 <    Molecule::IntegrableObjectIterator  j;              
120 <    bool needHydroPropFile = false;
121 <    
122 <    for (mol = info->beginMolecule(i); mol != NULL;
123 <         mol = info->nextMolecule(i)) {
124 <      for (integrableObject = mol->beginIntegrableObject(j);
125 <           integrableObject != NULL;
126 <           integrableObject = mol->nextIntegrableObject(j)) {
127 <        
128 <        if (integrableObject->isRigidBody()) {
129 <          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
130 <          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
131 <        }
132 <        
133 <      }
134 <    }
125 >    Molecule::IntegrableObjectIterator  j;
126  
136 –    hydroProps_.resize(info->getNIntegrableObjects());
137 –    
138 –    if (needHydroPropFile) {              
139 –      if (simParams->haveHydroPropFile()) {
140 –        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
141 –      } else {              
142 –        sprintf( painCave.errMsg,
143 –                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
144 –                 "\tis specified for rigidBodies which contain more than one atom\n"
145 –                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n");
146 –        painCave.severity = OOPSE_ERROR;
147 –        painCave.isFatal = 1;
148 –        simError();  
149 –      }      
150 –      
151 –
152 –
153 –
154 –
155 –
156 –      for (mol = info->beginMolecule(i); mol != NULL;
157 –           mol = info->nextMolecule(i)) {
158 –        for (integrableObject = mol->beginIntegrableObject(j);
159 –             integrableObject != NULL;
160 –             integrableObject = mol->nextIntegrableObject(j)) {
161 –
162 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
163 –          if (iter != hydroPropMap.end()) {
164 –            hydroProps_[integrableObject->getLocalIndex()] = iter->second;
165 –          } else {
166 –            sprintf( painCave.errMsg,
167 –                     "Can not find resistance tensor for atom [%s]\n",
168 –                     integrableObject->getType().c_str());
169 –            painCave.severity = OOPSE_ERROR;
170 –            painCave.isFatal = 1;
171 –            simError();  
172 –          }        
173 –        }
174 –      }
175 –    } else {
176 –      
177 –      std::map<std::string, HydroProp*> hydroPropMap;
178 –      for (mol = info->beginMolecule(i); mol != NULL;
179 –           mol = info->nextMolecule(i)) {
180 –        for (integrableObject = mol->beginIntegrableObject(j);
181 –             integrableObject != NULL;
182 –             integrableObject = mol->nextIntegrableObject(j)) {
183 –          Shape* currShape = NULL;
184 –
185 –          if (integrableObject->isAtom()){
186 –            Atom* atom = static_cast<Atom*>(integrableObject);
187 –            AtomType* atomType = atom->getAtomType();
188 –            if (atomType->isGayBerne()) {
189 –              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
190 –              GenericData* data = dAtomType->getPropertyByName("GayBerne");
191 –              if (data != NULL) {
192 –                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
193 –                
194 –                if (gayBerneData != NULL) {  
195 –                  GayBerneParam gayBerneParam = gayBerneData->getData();
196 –                  currShape = new Ellipsoid(V3Zero,
197 –                                            gayBerneParam.GB_l / 2.0,
198 –                                            gayBerneParam.GB_d / 2.0,
199 –                                            Mat3x3d::identity());
200 –                } else {
201 –                  sprintf( painCave.errMsg,
202 –                           "Can not cast GenericData to GayBerneParam\n");
203 –                  painCave.severity = OOPSE_ERROR;
204 –                  painCave.isFatal = 1;
205 –                  simError();  
206 –                }
207 –              } else {
208 –                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
209 –                painCave.severity = OOPSE_ERROR;
210 –                painCave.isFatal = 1;
211 –                simError();    
212 –              }
213 –            } else {
214 –              if (atomType->isLennardJones()){
215 –                GenericData* data = atomType->getPropertyByName("LennardJones");
216 –                if (data != NULL) {
217 –                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
218 –                  if (ljData != NULL) {
219 –                    LJParam ljParam = ljData->getData();
220 –                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
221 –                  } else {
222 –                    sprintf( painCave.errMsg,
223 –                             "Can not cast GenericData to LJParam\n");
224 –                    painCave.severity = OOPSE_ERROR;
225 –                    painCave.isFatal = 1;
226 –                    simError();          
227 –                  }      
228 –                }
229 –              } else {
230 –                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
231 –                if (aNum != 0) {
232 –                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
233 –                } else {
234 –                  sprintf( painCave.errMsg,
235 –                           "Could not find atom type in default element.txt\n");
236 –                  painCave.severity = OOPSE_ERROR;
237 –                  painCave.isFatal = 1;
238 –                  simError();          
239 –                }
240 –              }
241 –            }
242 –          }
243 –          HydroProp* currHydroProp = currShape->getHydroProp(viscosity_, targetTemp_);
244 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
245 –          if (iter != hydroPropMap.end()) {
246 –            hydroProps_[integrableObject->getLocalIndex()] = iter->second;
247 –          } else {
248 –            currHydroProp->complete();
249 –            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
250 –            hydroProps_[integrableObject->getLocalIndex()] = currHydroProp;
251 –          }
252 –        }
253 –      }
254 –    }
255 –
256 –    // Build a vector of integrable objects to determine if the are
257 –    // surface atoms
258 –    
127      for (mol = info_->beginMolecule(i); mol != NULL;
128           mol = info_->nextMolecule(i)) {          
129        for (integrableObject = mol->beginIntegrableObject(j);
# Line 265 | Line 133 | namespace oopse {
133        }
134      }  
135    }  
268 –
269 –  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
270 –    std::map<std::string, HydroProp*> props;
271 –    std::ifstream ifs(filename.c_str());
272 –    if (ifs.is_open()) {
273 –      
274 –    }
275 –    
276 –    const unsigned int BufferSize = 65535;
277 –    char buffer[BufferSize];  
278 –    while (ifs.getline(buffer, BufferSize)) {
279 –      HydroProp* currProp = new HydroProp(buffer);
280 –      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
281 –    }
282 –    
283 –    return props;
284 –  }
136    
137    void SMIPDForceManager::postCalculation(bool needStress){
138      SimInfo::MoleculeIterator i;
139      Molecule::IntegrableObjectIterator  j;
140      Molecule* mol;
141      StuntDouble* integrableObject;
291 –
292 –    RealType mass;
293 –    Vector3d pos;
294 –    Vector3d frc;
295 –    Mat3x3d A;
296 –    Mat3x3d Atrans;
297 –    Vector3d Tb;
298 –    Vector3d ji;
299 –    int index = 0;
142    
143      // Compute surface Mesh
144      surfaceMesh_->computeHull(localSites_);
# Line 307 | Line 149 | namespace oopse {
149      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
150      int nTriangles = sMesh.size();
151  
152 +    // Generate all of the necessary random forces
153 +    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
154 +
155 +    RealType instaTemp = thermo->getTemperature();
156 +
157      // Loop over the mesh faces and apply external pressure to each
158      // of the faces
159      std::vector<Triangle>::iterator face;
160      std::vector<StuntDouble*>::iterator vertex;
161 <    int thisNumber = 0;
161 >    int thisFacet = 0;
162      for (face = sMesh.begin(); face != sMesh.end(); ++face){
163      
164        Triangle thisTriangle = *face;
# Line 320 | Line 167 | namespace oopse {
167        Vector3d unitNormal = thisTriangle.getNormal();
168        unitNormal.normalize();
169        Vector3d centroid = thisTriangle.getCentroid();
170 <      Vector3d extPressure = - unitNormal * (targetPressure_ * thisArea)
171 <        / OOPSEConstant::energyConvert;
170 >      Vector3d facetVel = thisTriangle.getFacetVelocity();
171 >      RealType thisMass = thisTriangle.getFacetMass();
172 >
173 >      // gamma is the drag coefficient normal to the face of the triangle      
174 >      RealType gamma = thermalConductivity_ * thisMass * thisArea  
175 >        / (2.0 * thermalLength_ * OOPSEConstant::kb);
176 >      
177 >      gamma *= fabs(1.0 - targetTemp_/instaTemp);      
178 >      
179 >      RealType extPressure = - (targetPressure_ * thisArea) /
180 >        OOPSEConstant::energyConvert;
181 >
182 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
183 >      RealType dragForce = -gamma * dot(facetVel, unitNormal);
184 >
185 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
186 >        unitNormal;
187        
188 +      // Apply triangle force to stuntdouble vertices
189        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
190          if ((*vertex) != NULL){
191 <          Vector3d vertexForce = extPressure/3.0;
192 <          (*vertex)->addFrc(vertexForce);
330 <          
191 >          Vector3d vertexForce = langevinForce / 3.0;
192 >          (*vertex)->addFrc(vertexForce);          
193            if ((*vertex)->isDirectional()){          
194              Vector3d vertexPos = (*vertex)->getPos();
195              Vector3d vertexCentroidVector = vertexPos - centroid;
# Line 336 | Line 198 | namespace oopse {
198          }  
199        }
200      }
339 –
340 –    // Now loop over all surface particles and apply the drag
341 –    // and random forces
201      
343 –    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
344 –    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
345 –      integrableObject = *vertex;      
346 –      index = integrableObject->getLocalIndex();
347 –      mass = integrableObject->getMass();
348 –      if (integrableObject->isDirectional()){
349 –        
350 –        // preliminaries for directional objects:
351 –        
352 –        A = integrableObject->getA();
353 –        Atrans = A.transpose();
354 –        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
355 –
356 –        // apply random force and torque at center of resistance
357 –
358 –        Vector3d randomForceBody;
359 –        Vector3d randomTorqueBody;
360 –        genRandomForceAndTorque(randomForceBody, randomTorqueBody,
361 –                                index, variance_);
362 –        Vector3d randomForceLab = Atrans * randomForceBody;
363 –        Vector3d randomTorqueLab = Atrans * randomTorqueBody;
364 –        integrableObject->addFrc(randomForceLab);            
365 –        integrableObject->addTrq(randomTorqueLab + cross(rcrLab,
366 –                                                         randomForceLab ));
367 –
368 –        Mat3x3d I = integrableObject->getI();
369 –        Vector3d omegaBody;
370 –        
371 –        // What remains contains velocity explicitly, but the velocity required
372 –        // is at the full step: v(t + h), while we have initially the velocity
373 –        // at the half step: v(t + h/2).  We need to iterate to converge the
374 –        // friction force and friction torque vectors.
375 –        
376 –        // this is the velocity at the half-step:
377 –        
378 –        Vector3d vel =integrableObject->getVel();
379 –        Vector3d angMom = integrableObject->getJ();
380 –        
381 –        // estimate velocity at full-step using everything but friction forces:
382 –        
383 –        frc = integrableObject->getFrc();
384 –        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
385 –        
386 –        Tb = integrableObject->lab2Body(integrableObject->getTrq());
387 –        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
388 –        
389 –        Vector3d omegaLab;
390 –        Vector3d vcdLab;
391 –        Vector3d vcdBody;
392 –        Vector3d frictionForceBody;
393 –        Vector3d frictionForceLab(0.0);
394 –        Vector3d oldFFL;  // used to test for convergence
395 –        Vector3d frictionTorqueBody(0.0);
396 –        Vector3d oldFTB;  // used to test for convergence
397 –        Vector3d frictionTorqueLab;
398 –        RealType fdot;
399 –        RealType tdot;
400 –
401 –        //iteration starts here:
402 –        
403 –        for (int k = 0; k < maxIterNum_; k++) {
404 –          
405 –          if (integrableObject->isLinear()) {
406 –            int linearAxis = integrableObject->linearAxis();
407 –            int l = (linearAxis +1 )%3;
408 –            int m = (linearAxis +2 )%3;
409 –            omegaBody[l] = angMomStep[l] /I(l, l);
410 –            omegaBody[m] = angMomStep[m] /I(m, m);
411 –            
412 –          } else {
413 –            omegaBody[0] = angMomStep[0] /I(0, 0);
414 –            omegaBody[1] = angMomStep[1] /I(1, 1);
415 –            omegaBody[2] = angMomStep[2] /I(2, 2);
416 –          }
417 –          
418 –          omegaLab = Atrans * omegaBody;
419 –          
420 –          // apply friction force and torque at center of resistance
421 –          
422 –          vcdLab = velStep + cross(omegaLab, rcrLab);      
423 –          vcdBody = A * vcdLab;
424 –          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
425 –          oldFFL = frictionForceLab;
426 –          frictionForceLab = Atrans * frictionForceBody;
427 –          oldFTB = frictionTorqueBody;
428 –          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
429 –          frictionTorqueLab = Atrans * frictionTorqueBody;
430 –          
431 –          // re-estimate velocities at full-step using friction forces:
432 –              
433 –          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
434 –          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
435 –
436 –          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
437 –              
438 –          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
439 –          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
440 –          
441 –          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
442 –            break; // iteration ends here
443 –        }
444 –        
445 –        integrableObject->addFrc(frictionForceLab);
446 –        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
447 –
448 –            
449 –      } else {
450 –        //spherical atom
451 –
452 –        Vector3d randomForce;
453 –        Vector3d randomTorque;
454 –        genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
455 –        integrableObject->addFrc(randomForce);            
456 –        
457 –        // What remains contains velocity explicitly, but the velocity required
458 –        // is at the full step: v(t + h), while we have initially the velocity
459 –        // at the half step: v(t + h/2).  We need to iterate to converge the
460 –        // friction force vector.
461 –        
462 –        // this is the velocity at the half-step:
463 –        
464 –        Vector3d vel =integrableObject->getVel();
465 –        
466 –        //estimate velocity at full-step using everything but friction forces:          
467 –        
468 –        frc = integrableObject->getFrc();
469 –        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
470 –        
471 –        Vector3d frictionForce(0.0);
472 –        Vector3d oldFF;  // used to test for convergence
473 –        RealType fdot;
474 –        
475 –        //iteration starts here:
476 –        
477 –        for (int k = 0; k < maxIterNum_; k++) {
478 –          
479 –          oldFF = frictionForce;                            
480 –          frictionForce = -hydroProps_[index]->getXitt() * velStep;
481 –          //frictionForce = -gamma_t*velStep;
482 –          // re-estimate velocities at full-step using friction forces:
483 –          
484 –          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
485 –          
486 –          // check for convergence (if the vector has converged, fdot will be 1.0):
487 –          
488 –          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
489 –          
490 –          if (fabs(1.0 - fdot) <= forceTolerance_)
491 –            break; // iteration ends here
492 –        }
493 –        
494 –        integrableObject->addFrc(frictionForce);
495 –                
496 –      }      
497 –    }
498 –        
202      veloMunge->removeComDrift();
203      veloMunge->removeAngularDrift();
204 <
204 >    
205      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
206      currSnapshot->setVolume(surfaceMesh_->getVolume());    
207      ForceManager::postCalculation(needStress);  
208    }
209 <
210 <  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force,
211 <                                                  Vector3d& torque,
212 <                                                  unsigned int index,
213 <                                                  RealType variance) {
511 <        
512 <    Vector<RealType, 6> Z;
513 <    Vector<RealType, 6> generalForce;
209 >  
210 >  
211 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
212 >                                                             RealType variance)
213 >  {
214      
215 <    Z[0] = randNumGen_.randNorm(0, variance);
216 <    Z[1] = randNumGen_.randNorm(0, variance);
217 <    Z[2] = randNumGen_.randNorm(0, variance);
218 <    Z[3] = randNumGen_.randNorm(0, variance);
519 <    Z[4] = randNumGen_.randNorm(0, variance);
520 <    Z[5] = randNumGen_.randNorm(0, variance);
215 >    // zero fill the random vector before starting:
216 >    std::vector<RealType> gaussRand;
217 >    gaussRand.resize(nTriangles);
218 >    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
219      
220 <
221 <    generalForce = hydroProps_[index]->getS()*Z;
220 > #ifdef IS_MPI
221 >    if (worldRank == 0) {
222 > #endif
223 >      for (int i = 0; i < nTriangles; i++) {
224 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
225 >      }
226 > #ifdef IS_MPI
227 >    }
228 > #endif
229      
230 <    force[0] = generalForce[0];
231 <    force[1] = generalForce[1];
232 <    force[2] = generalForce[2];
233 <    torque[0] = generalForce[3];
234 <    torque[1] = generalForce[4];
235 <    torque[2] = generalForce[5];    
236 <  }
230 >    // push these out to the other processors
231 >    
232 > #ifdef IS_MPI
233 >    if (worldRank == 0) {
234 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
235 >    } else {
236 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
237 >    }
238 > #endif
239 >    
240 >    return gaussRand;
241 >  }
242   }

Diff Legend

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