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 3481 by gezelter, Wed Nov 26 14:26:17 2008 UTC vs.
Revision 3482 by chuckv, Fri Dec 5 16:20:39 2008 UTC

# Line 41 | Line 41
41   #include <fstream>
42   #include <iostream>
43   #include "integrators/SMIPDForceManager.hpp"
44 + #include "math/CholeskyDecomposition.hpp"
45   #include "utils/OOPSEConstant.hpp"
46 + #include "hydrodynamics/Sphere.hpp"
47 + #include "hydrodynamics/Ellipsoid.hpp"
48 + #include "utils/ElementsTable.hpp"
49   #include "math/ConvexHull.hpp"
50   #include "math/Triangle.hpp"
51  
52   namespace oopse {
49 –  
50 –  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
53  
54 +  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info),
55 +                                                        forceTolerance_(1e-6),
56 +                                                        maxIterNum_(4) {
57      simParams = info->getSimParams();
58 +    veloMunge = new Velocitizer(info);
59      
60      // Create Hull, Convex Hull for now, other options later.
61      
62      surfaceMesh_ = new ConvexHull();
63      
64      /* Check that the simulation has target pressure and target
65 <       temperature set*/
65 >       temperature set */
66      
67      if (!simParams->haveTargetTemp()) {
68        sprintf(painCave.errMsg,
# Line 104 | Line 110 | namespace oopse {
110    
111      variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
112  
113 +    // Build the hydroProp map:
114 +    std::map<std::string, HydroProp*> hydroPropMap;
115 +
116      Molecule* mol;
117      StuntDouble* integrableObject;
118      SimInfo::MoleculeIterator i;
119 <    Molecule::IntegrableObjectIterator  j;              
119 >    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 >    }
135  
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 <
258 >    
259      for (mol = info_->beginMolecule(i); mol != NULL;
260           mol = info_->nextMolecule(i)) {          
261        for (integrableObject = mol->beginIntegrableObject(j);
# Line 121 | Line 265 | namespace oopse {
265        }
266      }  
267    }  
268 <  
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 >  }
285 >  
286    void SMIPDForceManager::postCalculation(bool needStress){
287      SimInfo::MoleculeIterator i;
288      Molecule::IntegrableObjectIterator  j;
289      Molecule* mol;
290      StuntDouble* integrableObject;
291 <  
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;
300 >  
301      // Compute surface Mesh
302      surfaceMesh_->computeHull(localSites_);
303 <    
303 >
304      // Get total area and number of surface stunt doubles
305      RealType area = surfaceMesh_->getArea();
306 <    int nSurfaceSDs = surfaceMesh_->getNs();        
306 >    int nSurfaceSDs = surfaceMesh_->getNs();
307      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
308      int nTriangles = sMesh.size();
139 –        
140 –    // Generate all of the necessary random forces
141 –    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
309  
310 <    // Loop over the mesh faces and apply random force to each of the faces
310 >    // Loop over the mesh faces and apply external pressure to each
311 >    // of the faces
312      std::vector<Triangle>::iterator face;
313      std::vector<StuntDouble*>::iterator vertex;
314      int thisNumber = 0;
# Line 151 | Line 319 | namespace oopse {
319        RealType thisArea = thisTriangle.getArea();
320        Vector3d unitNormal = thisTriangle.getNormal();
321        unitNormal.normalize();
154 –
322        Vector3d centroid = thisTriangle.getCentroid();
323 <      Vector3d facetVel = thisTriangle.getFacetVelocity();
324 <      RealType hydroLength = thisTriangle.getIncircleRadius() * 2.0 /
158 <        NumericConstant::PI;
323 >      Vector3d extPressure = - unitNormal * (targetPressure_ * thisArea)
324 >        / OOPSEConstant::energyConvert;
325        
326 <      // gamma is the drag coefficient normal to the face of the triangle      
327 <      RealType gamma = viscosity_ * hydroLength *
328 <        OOPSEConstant::viscoConvert;
329 <
330 <      RealType extPressure = - (targetPressure_ * thisArea) /
331 <        OOPSEConstant::energyConvert;
166 <
167 <      RealType randomForce = randNums[thisNumber++] * sqrt(gamma);
168 <      RealType dragForce = -gamma * dot(facetVel, unitNormal);
169 <
170 <      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
171 <      
172 <      // Apply triangle force to stuntdouble vertices
173 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
174 <        if ((*vertex) != NULL) {
175 <          Vector3d vertexForce = langevinForce / 3.0;
176 <          (*vertex)->addFrc(vertexForce);          
177 <          if ((*vertex)->isDirectional()) {
326 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
327 >        if ((*vertex) != NULL){
328 >          Vector3d vertexForce = extPressure/3.0;
329 >          (*vertex)->addFrc(vertexForce);
330 >          
331 >          if ((*vertex)->isDirectional()){          
332              Vector3d vertexPos = (*vertex)->getPos();
333              Vector3d vertexCentroidVector = vertexPos - centroid;
334              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 182 | Line 336 | namespace oopse {
336          }  
337        }
338      }
339 +
340 +    // Now loop over all surface particles and apply the drag
341 +    // and random forces
342      
343 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
344 <    currSnapshot->setVolume(surfaceMesh_->getVolume());    
345 <    ForceManager::postCalculation(needStress);  
346 <  }
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 <  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
192 <                                                             RealType variance) {
356 >        // apply random force and torque at center of resistance
357  
358 <    // zero fill the random vector before starting:
359 <    std::vector<RealType> gaussRand;
360 <    gaussRand.resize(nTriangles);
361 <    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
362 <  
363 < #ifdef IS_MPI
364 <    if (worldRank == 0) {
365 < #endif
366 <      for (int i = 0; i < nTriangles; i++) {
203 <        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
204 <      }
205 < #ifdef IS_MPI
206 <    }
207 < #endif
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 <    // push these out to the other processors
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 < #ifdef IS_MPI
402 <    if (worldRank == 0) {
403 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
404 <    } else {
405 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
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 < #endif
499 <  
500 <    return gaussRand;
498 >        
499 >    veloMunge->removeComDrift();
500 >    veloMunge->removeAngularDrift();
501 >
502 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
503 >    currSnapshot->setVolume(surfaceMesh_->getVolume());    
504 >    ForceManager::postCalculation(needStress);  
505    }
506 +
507 +  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force,
508 +                                                  Vector3d& torque,
509 +                                                  unsigned int index,
510 +                                                  RealType variance) {
511 +        
512 +    Vector<RealType, 6> Z;
513 +    Vector<RealType, 6> generalForce;
514 +    
515 +    Z[0] = randNumGen_.randNorm(0, variance);
516 +    Z[1] = randNumGen_.randNorm(0, variance);
517 +    Z[2] = randNumGen_.randNorm(0, variance);
518 +    Z[3] = randNumGen_.randNorm(0, variance);
519 +    Z[4] = randNumGen_.randNorm(0, variance);
520 +    Z[5] = randNumGen_.randNorm(0, variance);
521 +    
522 +
523 +    generalForce = hydroProps_[index]->getS()*Z;
524 +    
525 +    force[0] = generalForce[0];
526 +    force[1] = generalForce[1];
527 +    force[2] = generalForce[2];
528 +    torque[0] = generalForce[3];
529 +    torque[1] = generalForce[4];
530 +    torque[2] = generalForce[5];    
531 +  }
532   }

Diff Legend

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