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 3459 by chuckv, Tue Oct 7 17:12:48 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();
288 >     RealType area = surfaceMesh_->getArea();
289 >     int nSurfaceSDs = surfaceMesh_->getNs();
290  
291 < //    std::cerr << "Surface Area is: " << area << " nSurfaceSDs is: " << nSurfaceSDs << std::endl;
289 <
290 <    /* 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      
295 < //    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
296 < //    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
297 <  
295 >    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
296 >    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
297 >    
298      /* Loop over the mesh faces and apply random force to each of the faces*/
299      
300      
301 < //    std::vector<Triangle*>::iterator face;
302 < //    std::vector<StuntDouble*>::iterator vertex;
303 < /*    
301 >    std::vector<Triangle*>::iterator face;
302 >    std::vector<StuntDouble*>::iterator vertex;
303 >    int thisNumber = 0;
304      for (face = sMesh.begin(); face != sMesh.end(); ++face){
305      
306        Triangle* thisTriangle = *face;
307        std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
308 <
308 >      
309 >      /* Get Random Force */
310 >      Vector3d unitNormal = thisTriangle->getNormal();
311 >      unitNormal.normalize();
312 >      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
313 >      
314        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
315 <            std::cout << (*vertex)->getPos() << std::endl;
315 >
316           // mass = integrableObject->getMass();
317  
318 <         //      integrableObject->addFrc(randomForce);          
318 >                    (*vertex)->addFrc(randomForce/3.0);          
319        }
320 +    }
321  
322 +    /* Now loop over all surface particles and apply the drag*/
323  
324 <    }
325 <  */
326 <    /*
327 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
324 >    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
325 >    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
326 >      integrableObject = *vertex;
327 >      mass = integrableObject->getMass();
328 >      if (integrableObject->isDirectional()){
329 >        
330 >        // preliminaries for directional objects:
331 >        
332 >        A = integrableObject->getA();
333 >        Atrans = A.transpose();
334 >        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
335  
336 <      
337 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
338 <           integrableObject = mol->nextIntegrableObject(j)) {
339 <          
340 <          mass = integrableObject->getMass();
341 <          if (integrableObject->isDirectional()){
342 <
343 <            // preliminaries for directional objects:
344 <
345 <            A = integrableObject->getA();
346 <            Atrans = A.transpose();
347 <            Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
348 <
349 <            //apply random force and torque at center of resistance
350 <
351 <            Vector3d randomForceBody;
352 <            Vector3d randomTorqueBody;
353 <            genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
354 <            Vector3d randomForceLab = Atrans * randomForceBody;
355 <            Vector3d randomTorqueLab = Atrans * randomTorqueBody;
356 <            integrableObject->addFrc(randomForceLab);            
357 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
358 <
359 <            Mat3x3d I = integrableObject->getI();
360 <            Vector3d omegaBody;
336 >        
337 >        Mat3x3d I = integrableObject->getI();
338 >        Vector3d omegaBody;
339 >        
340 >        // What remains contains velocity explicitly, but the velocity required
341 >        // is at the full step: v(t + h), while we have initially the velocity
342 >        // at the half step: v(t + h/2).  We need to iterate to converge the
343 >        // friction force and friction torque vectors.
344 >        
345 >        // this is the velocity at the half-step:
346 >        
347 >        Vector3d vel =integrableObject->getVel();
348 >        Vector3d angMom = integrableObject->getJ();
349 >        
350 >        //estimate velocity at full-step using everything but friction forces:          
351 >        
352 >        frc = integrableObject->getFrc();
353 >        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
354 >        
355 >        Tb = integrableObject->lab2Body(integrableObject->getTrq());
356 >        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
357 >        
358 >        Vector3d omegaLab;
359 >        Vector3d vcdLab;
360 >        Vector3d vcdBody;
361 >        Vector3d frictionForceBody;
362 >        Vector3d frictionForceLab(0.0);
363 >        Vector3d oldFFL;  // used to test for convergence
364 >        Vector3d frictionTorqueBody(0.0);
365 >        Vector3d oldFTB;  // used to test for convergence
366 >        Vector3d frictionTorqueLab;
367 >        RealType fdot;
368 >        RealType tdot;
369  
370 <            // What remains contains velocity explicitly, but the velocity required
371 <            // is at the full step: v(t + h), while we have initially the velocity
372 <            // at the half step: v(t + h/2).  We need to iterate to converge the
373 <            // friction force and friction torque vectors.
374 <
375 <            // this is the velocity at the half-step:
370 >        //iteration starts here:
371 >        
372 >        for (int k = 0; k < maxIterNum_; k++) {
373 >          
374 >          if (integrableObject->isLinear()) {
375 >            int linearAxis = integrableObject->linearAxis();
376 >            int l = (linearAxis +1 )%3;
377 >            int m = (linearAxis +2 )%3;
378 >            omegaBody[l] = angMomStep[l] /I(l, l);
379 >            omegaBody[m] = angMomStep[m] /I(m, m);
380              
381 <            Vector3d vel =integrableObject->getVel();
382 <            Vector3d angMom = integrableObject->getJ();
383 <
384 <            //estimate velocity at full-step using everything but friction forces:          
385 <
386 <            frc = integrableObject->getFrc();
387 <            Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
388 <
389 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
390 <            Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
391 <
392 <            Vector3d omegaLab;
393 <            Vector3d vcdLab;
394 <            Vector3d vcdBody;
395 <            Vector3d frictionForceBody;
396 <            Vector3d frictionForceLab(0.0);
397 <            Vector3d oldFFL;  // used to test for convergence
398 <            Vector3d frictionTorqueBody(0.0);
399 <            Vector3d oldFTB;  // used to test for convergence
400 <            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 <              }
393 <              
394 <              omegaLab = Atrans * omegaBody;
395 <              
396 <              // apply friction force and torque at center of resistance
381 >          } else {
382 >            omegaBody[0] = angMomStep[0] /I(0, 0);
383 >            omegaBody[1] = angMomStep[1] /I(1, 1);
384 >            omegaBody[2] = angMomStep[2] /I(2, 2);
385 >          }
386 >          
387 >          omegaLab = Atrans * omegaBody;
388 >          
389 >          // apply friction force and torque at center of resistance
390 >          
391 >          vcdLab = velStep + cross(omegaLab, rcrLab);      
392 >          vcdBody = A * vcdLab;
393 >          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
394 >          oldFFL = frictionForceLab;
395 >          frictionForceLab = Atrans * frictionForceBody;
396 >          oldFTB = frictionTorqueBody;
397 >          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
398 >          frictionTorqueLab = Atrans * frictionTorqueBody;
399 >          
400 >          // re-estimate velocities at full-step using friction forces:
401                
402 <              vcdLab = velStep + cross(omegaLab, rcrLab);      
403 <              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);
402 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
403 >          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
404  
405 <              // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
405 >          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
406                
407 <              fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
408 <              tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
409 <              
410 <              if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
411 <                break; // iteration ends here
412 <            }
407 >          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
408 >          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
409 >          
410 >          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
411 >            break; // iteration ends here
412 >        }
413 >        
414 >        integrableObject->addFrc(frictionForceLab);
415 >        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
416  
421 –            integrableObject->addFrc(frictionForceLab);
422 –            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
423 –
417              
418 <          } else {
419 <            //spherical atom
418 >      } else {
419 >        //spherical atom
420  
421 <            Vector3d randomForce;
422 <            Vector3d randomTorque;
423 <            genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
424 <            integrableObject->addFrc(randomForce);            
425 <
426 <            // What remains contains velocity explicitly, but the velocity required
427 <            // is at the full step: v(t + h), while we have initially the velocity
428 <            // at the half step: v(t + h/2).  We need to iterate to converge the
429 <            // friction force vector.
430 <
431 <            // this is the velocity at the half-step:
432 <            
433 <            Vector3d vel =integrableObject->getVel();
434 <
435 <            //estimate velocity at full-step using everything but friction forces:          
436 <
437 <            frc = integrableObject->getFrc();
438 <            Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
439 <
440 <            Vector3d frictionForce(0.0);
441 <            Vector3d oldFF;  // used to test for convergence
442 <            RealType fdot;
443 <
444 <            //iteration starts here:
445 <
446 <            for (int k = 0; k < maxIterNum_; k++) {
447 <
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 <
421 >        
422 >        // What remains contains velocity explicitly, but the velocity required
423 >        // is at the full step: v(t + h), while we have initially the velocity
424 >        // at the half step: v(t + h/2).  We need to iterate to converge the
425 >        // friction force vector.
426 >        
427 >        // this is the velocity at the half-step:
428 >        
429 >        Vector3d vel =integrableObject->getVel();
430 >        
431 >        //estimate velocity at full-step using everything but friction forces:          
432 >        
433 >        frc = integrableObject->getFrc();
434 >        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
435 >        
436 >        Vector3d frictionForce(0.0);
437 >        Vector3d oldFF;  // used to test for convergence
438 >        RealType fdot;
439 >        
440 >        //iteration starts here:
441 >        
442 >        for (int k = 0; k < maxIterNum_; k++) {
443 >          
444 >          oldFF = frictionForce;                            
445 >          frictionForce = -hydroProps_[index]->getXitt() * velStep;
446 >          
447 >          // re-estimate velocities at full-step using friction forces:
448            
449 <          }
449 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
450 >          
451 >          // check for convergence (if the vector has converged, fdot will be 1.0):
452            
453 <        ++index;
454 <    
453 >          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
454 >          
455 >          if (fabs(1.0 - fdot) <= forceTolerance_)
456 >            break; // iteration ends here
457 >        }
458 >        
459 >        integrableObject->addFrc(frictionForce);
460 >        
461 >        
462        }
463 <    }    
464 <    */
465 <    // info_->setFdf(fdf);
466 <    // veloMunge->removeComDrift();
467 <    // Remove angular drift if we are not using periodic boundary conditions.
468 <    //if(!simParams->getUsePeriodicBoundaryConditions())
469 <    //  veloMunge->removeAngularDrift();
485 <
486 <    //ForceManager::postCalculation(needStress);  
463 >      
464 >      
465 >    }
466 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
467 >    currSnapshot->setVolume(surfaceMesh_->getVolume());
468 >    
469 >    ForceManager::postCalculation(needStress);  
470    }
471  
472 < void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
472 >  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
473  
474 <
474 >    
475      Vector<RealType, 6> Z;
476      Vector<RealType, 6> generalForce;
477 <        
477 >    
478 >
479      Z[0] = randNumGen_.randNorm(0, variance);
480      Z[1] = randNumGen_.randNorm(0, variance);
481      Z[2] = randNumGen_.randNorm(0, variance);

Diff Legend

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