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 3458 by chuckv, Fri Oct 3 17:47:08 2008 UTC vs.
Revision 3461 by chuckv, Wed Oct 15 18:26:01 2008 UTC

# Line 121 | Line 121 | namespace oopse {
121          hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
122        } else {              
123          sprintf( painCave.errMsg,
124 <                 "HydroPropFile must be set to a file name if Langevin Dynamics\n"
124 >                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
125                   "\tis specified for rigidBodies which contain more than one atom\n"
126 <                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n");
126 >                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
127 >                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
128 >                 "\tset to 1.0 because the friction and random forces will be\n"
129 >                 "\tdynamically re-set assuming this is true.\n");
130          painCave.severity = OOPSE_ERROR;
131          painCave.isFatal = 1;
132          simError();  
# Line 215 | Line 218 | namespace oopse {
218                }
219              }
220            }
221 <          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
221 >          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
222            std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
223            if (iter != hydroPropMap.end())
224              hydroProps_.push_back(iter->second);
# Line 282 | Line 285 | namespace oopse {
285      surfaceMesh_->computeHull(localSites_);
286  
287      /* Get area and number of surface stunt doubles and compute new variance */
288 <     //RealType area = surfaceMesh_->getArea();
289 <     //RealType nSurfaceSDs = surfaceMesh_->getNs();
287 <
288 < //    std::cerr << "Surface Area is: " << area << " nSurfaceSDs is: " << nSurfaceSDs << std::endl;
288 >     RealType area = surfaceMesh_->getArea();
289 >     int nSurfaceSDs = surfaceMesh_->getNs();
290  
291 <    /* Compute variance for random forces */
291 >     /* Compute variance for random forces */
292      
293 < //    variance_ = sqrt(2.0*NumericConstant::PI)*(targetPressure_*area/nSurfaceSDs);
293 >     variance_ = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
294 >       /OOPSEConstant::energyConvert;
295      
296 < //    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
297 < //    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
298 <  
296 >    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
297 >    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
298 >    
299      /* Loop over the mesh faces and apply random force to each of the faces*/
300      
301      
302 < //    std::vector<Triangle*>::iterator face;
303 < //    std::vector<StuntDouble*>::iterator vertex;
304 < /*    
302 >    std::vector<Triangle*>::iterator face;
303 >    std::vector<StuntDouble*>::iterator vertex;
304 >    int thisNumber = 0;
305      for (face = sMesh.begin(); face != sMesh.end(); ++face){
306      
307        Triangle* thisTriangle = *face;
308        std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
309 +      
310 +      /* Get Random Force */
311 +      Vector3d unitNormal = thisTriangle->getNormal();
312 +      unitNormal.normalize();
313 +      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
314 +      Vector3d centroid = thisTriangle->getCentroid();
315  
316        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
309 –            std::cout << (*vertex)->getPos() << std::endl;
310 –         // mass = integrableObject->getMass();
317  
318 <         //      integrableObject->addFrc(randomForce);          
318 >         // mass = integrableObject->getMass();
319 >        Vector3d vertexForce = randomForce/3.0;
320 >        (*vertex)->addFrc(vertexForce);
321 >        if (integrableObject->isDirectional()){
322 >          Vector3d vertexPos = (*vertex)->getPos();
323 >          Vector3d vertexCentroidVector = vertexPos - centroid;
324 >          (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
325 >        }
326 >          
327        }
328 +    }
329  
330 +    /* Now loop over all surface particles and apply the drag*/
331  
332 <    }
333 <  */
334 <    /*
335 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
332 >    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
333 >    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
334 >      integrableObject = *vertex;
335 >      mass = integrableObject->getMass();
336 >      if (integrableObject->isDirectional()){
337 >        
338 >        // preliminaries for directional objects:
339 >        
340 >        A = integrableObject->getA();
341 >        Atrans = A.transpose();
342 >        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
343  
344 <      
345 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
346 <           integrableObject = mol->nextIntegrableObject(j)) {
347 <          
348 <          mass = integrableObject->getMass();
349 <          if (integrableObject->isDirectional()){
350 <
351 <            // preliminaries for directional objects:
352 <
353 <            A = integrableObject->getA();
354 <            Atrans = A.transpose();
355 <            Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
356 <
357 <            //apply random force and torque at center of resistance
358 <
359 <            Vector3d randomForceBody;
360 <            Vector3d randomTorqueBody;
361 <            genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
362 <            Vector3d randomForceLab = Atrans * randomForceBody;
363 <            Vector3d randomTorqueLab = Atrans * randomTorqueBody;
364 <            integrableObject->addFrc(randomForceLab);            
365 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
344 >        
345 >        Mat3x3d I = integrableObject->getI();
346 >        Vector3d omegaBody;
347 >        
348 >        // What remains contains velocity explicitly, but the velocity required
349 >        // is at the full step: v(t + h), while we have initially the velocity
350 >        // at the half step: v(t + h/2).  We need to iterate to converge the
351 >        // friction force and friction torque vectors.
352 >        
353 >        // this is the velocity at the half-step:
354 >        
355 >        Vector3d vel =integrableObject->getVel();
356 >        Vector3d angMom = integrableObject->getJ();
357 >        
358 >        //estimate velocity at full-step using everything but friction forces:          
359 >        
360 >        frc = integrableObject->getFrc();
361 >        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
362 >        
363 >        Tb = integrableObject->lab2Body(integrableObject->getTrq());
364 >        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
365 >        
366 >        Vector3d omegaLab;
367 >        Vector3d vcdLab;
368 >        Vector3d vcdBody;
369 >        Vector3d frictionForceBody;
370 >        Vector3d frictionForceLab(0.0);
371 >        Vector3d oldFFL;  // used to test for convergence
372 >        Vector3d frictionTorqueBody(0.0);
373 >        Vector3d oldFTB;  // used to test for convergence
374 >        Vector3d frictionTorqueLab;
375 >        RealType fdot;
376 >        RealType tdot;
377  
378 <            Mat3x3d I = integrableObject->getI();
379 <            Vector3d omegaBody;
380 <
381 <            // What remains contains velocity explicitly, but the velocity required
382 <            // is at the full step: v(t + h), while we have initially the velocity
383 <            // at the half step: v(t + h/2).  We need to iterate to converge the
384 <            // friction force and friction torque vectors.
385 <
386 <            // this is the velocity at the half-step:
378 >        //iteration starts here:
379 >        
380 >        for (int k = 0; k < maxIterNum_; k++) {
381 >          
382 >          if (integrableObject->isLinear()) {
383 >            int linearAxis = integrableObject->linearAxis();
384 >            int l = (linearAxis +1 )%3;
385 >            int m = (linearAxis +2 )%3;
386 >            omegaBody[l] = angMomStep[l] /I(l, l);
387 >            omegaBody[m] = angMomStep[m] /I(m, m);
388              
389 <            Vector3d vel =integrableObject->getVel();
390 <            Vector3d angMom = integrableObject->getJ();
391 <
392 <            //estimate velocity at full-step using everything but friction forces:          
393 <
394 <            frc = integrableObject->getFrc();
395 <            Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
396 <
397 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
398 <            Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
399 <
400 <            Vector3d omegaLab;
401 <            Vector3d vcdLab;
402 <            Vector3d vcdBody;
403 <            Vector3d frictionForceBody;
404 <            Vector3d frictionForceLab(0.0);
405 <            Vector3d oldFFL;  // used to test for convergence
406 <            Vector3d frictionTorqueBody(0.0);
407 <            Vector3d oldFTB;  // used to test for convergence
408 <            Vector3d frictionTorqueLab;
374 <            RealType fdot;
375 <            RealType tdot;
376 <
377 <            //iteration starts here:
378 <
379 <            for (int k = 0; k < maxIterNum_; k++) {
380 <                            
381 <              if (integrableObject->isLinear()) {
382 <                int linearAxis = integrableObject->linearAxis();
383 <                int l = (linearAxis +1 )%3;
384 <                int m = (linearAxis +2 )%3;
385 <                omegaBody[l] = angMomStep[l] /I(l, l);
386 <                omegaBody[m] = angMomStep[m] /I(m, m);
387 <                
388 <              } else {
389 <                omegaBody[0] = angMomStep[0] /I(0, 0);
390 <                omegaBody[1] = angMomStep[1] /I(1, 1);
391 <                omegaBody[2] = angMomStep[2] /I(2, 2);
392 <              }
389 >          } else {
390 >            omegaBody[0] = angMomStep[0] /I(0, 0);
391 >            omegaBody[1] = angMomStep[1] /I(1, 1);
392 >            omegaBody[2] = angMomStep[2] /I(2, 2);
393 >          }
394 >          
395 >          omegaLab = Atrans * omegaBody;
396 >          
397 >          // apply friction force and torque at center of resistance
398 >          
399 >          vcdLab = velStep + cross(omegaLab, rcrLab);      
400 >          vcdBody = A * vcdLab;
401 >          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
402 >          oldFFL = frictionForceLab;
403 >          frictionForceLab = Atrans * frictionForceBody;
404 >          oldFTB = frictionTorqueBody;
405 >          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
406 >          frictionTorqueLab = Atrans * frictionTorqueBody;
407 >          
408 >          // re-estimate velocities at full-step using friction forces:
409                
410 <              omegaLab = Atrans * omegaBody;
411 <              
396 <              // apply friction force and torque at center of resistance
397 <              
398 <              vcdLab = velStep + cross(omegaLab, rcrLab);      
399 <              vcdBody = A * vcdLab;
400 <              frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
401 <              oldFFL = frictionForceLab;
402 <              frictionForceLab = Atrans * frictionForceBody;
403 <              oldFTB = frictionTorqueBody;
404 <              frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
405 <              frictionTorqueLab = Atrans * frictionTorqueBody;
406 <              
407 <              // re-estimate velocities at full-step using friction forces:
408 <              
409 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
410 <              angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
410 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
411 >          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
412  
413 <              // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
413 >          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
414                
415 <              fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
416 <              tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
417 <              
418 <              if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
419 <                break; // iteration ends here
420 <            }
415 >          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
416 >          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
417 >          
418 >          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
419 >            break; // iteration ends here
420 >        }
421 >        
422 >        integrableObject->addFrc(frictionForceLab);
423 >        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
424  
421 –            integrableObject->addFrc(frictionForceLab);
422 –            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
423 –
425              
426 <          } else {
427 <            //spherical atom
426 >      } else {
427 >        //spherical atom
428  
429 <            Vector3d randomForce;
430 <            Vector3d randomTorque;
431 <            genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
432 <            integrableObject->addFrc(randomForce);            
433 <
434 <            // What remains contains velocity explicitly, but the velocity required
435 <            // is at the full step: v(t + h), while we have initially the velocity
436 <            // at the half step: v(t + h/2).  We need to iterate to converge the
437 <            // friction force vector.
438 <
439 <            // this is the velocity at the half-step:
440 <            
441 <            Vector3d vel =integrableObject->getVel();
442 <
443 <            //estimate velocity at full-step using everything but friction forces:          
444 <
445 <            frc = integrableObject->getFrc();
446 <            Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
447 <
448 <            Vector3d frictionForce(0.0);
449 <            Vector3d oldFF;  // used to test for convergence
450 <            RealType fdot;
451 <
452 <            //iteration starts here:
453 <
454 <            for (int k = 0; k < maxIterNum_; k++) {
455 <
455 <              oldFF = frictionForce;                            
456 <              frictionForce = -hydroProps_[index]->getXitt() * velStep;
457 <
458 <              // re-estimate velocities at full-step using friction forces:
459 <              
460 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
461 <
462 <              // check for convergence (if the vector has converged, fdot will be 1.0):
463 <              
464 <              fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
465 <              
466 <              if (fabs(1.0 - fdot) <= forceTolerance_)
467 <                break; // iteration ends here
468 <            }
469 <
470 <            integrableObject->addFrc(frictionForce);
471 <
429 >        
430 >        // What remains contains velocity explicitly, but the velocity required
431 >        // is at the full step: v(t + h), while we have initially the velocity
432 >        // at the half step: v(t + h/2).  We need to iterate to converge the
433 >        // friction force vector.
434 >        
435 >        // this is the velocity at the half-step:
436 >        
437 >        Vector3d vel =integrableObject->getVel();
438 >        
439 >        //estimate velocity at full-step using everything but friction forces:          
440 >        
441 >        frc = integrableObject->getFrc();
442 >        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
443 >        
444 >        Vector3d frictionForce(0.0);
445 >        Vector3d oldFF;  // used to test for convergence
446 >        RealType fdot;
447 >        
448 >        //iteration starts here:
449 >        
450 >        for (int k = 0; k < maxIterNum_; k++) {
451 >          
452 >          oldFF = frictionForce;                            
453 >          frictionForce = -hydroProps_[index]->getXitt() * velStep;
454 >          
455 >          // re-estimate velocities at full-step using friction forces:
456            
457 <          }
457 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
458 >          
459 >          // check for convergence (if the vector has converged, fdot will be 1.0):
460            
461 <        ++index;
462 <    
461 >          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
462 >          
463 >          if (fabs(1.0 - fdot) <= forceTolerance_)
464 >            break; // iteration ends here
465 >        }
466 >        
467 >        integrableObject->addFrc(frictionForce);
468 >        
469 >        
470        }
471 <    }    
472 <    */
473 <    // info_->setFdf(fdf);
474 <    // veloMunge->removeComDrift();
475 <    // Remove angular drift if we are not using periodic boundary conditions.
476 <    //if(!simParams->getUsePeriodicBoundaryConditions())
477 <    //  veloMunge->removeAngularDrift();
485 <
486 <    //ForceManager::postCalculation(needStress);  
471 >      
472 >      
473 >    }
474 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
475 >    currSnapshot->setVolume(surfaceMesh_->getVolume());
476 >    
477 >    ForceManager::postCalculation(needStress);  
478    }
479  
480 < void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
480 >  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
481  
482 <
482 >    
483      Vector<RealType, 6> Z;
484      Vector<RealType, 6> generalForce;
485 <        
485 >    
486 >
487      Z[0] = randNumGen_.randNorm(0, variance);
488      Z[1] = randNumGen_.randNorm(0, variance);
489      Z[2] = randNumGen_.randNorm(0, variance);

Diff Legend

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