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 3517 by gezelter, Thu Jul 23 18:49:45 2009 UTC

# Line 46 | Line 46 | namespace oopse {
46   #include "math/Triangle.hpp"
47  
48   namespace oopse {
49 <  
49 >
50    SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
51  
52      simParams = info->getSimParams();
53 +    veloMunge = new Velocitizer(info);
54      
55      // Create Hull, Convex Hull for now, other options later.
56      
57      surfaceMesh_ = new ConvexHull();
58      
59      /* Check that the simulation has target pressure and target
60 <       temperature set*/
60 >       temperature set */
61      
62      if (!simParams->haveTargetTemp()) {
63        sprintf(painCave.errMsg,
64                "SMIPDynamics error: You can't use the SMIPD integrator\n"
65 <              "   without a targetTemp!\n");      
65 >              "\twithout a targetTemp (K)!\n");      
66        painCave.isFatal = 1;
67        painCave.severity = OOPSE_ERROR;
68        simError();
# Line 72 | Line 73 | namespace oopse {
73      if (!simParams->haveTargetPressure()) {
74        sprintf(painCave.errMsg,
75                "SMIPDynamics error: You can't use the SMIPD integrator\n"
76 <              "   without a targetPressure!\n");      
76 >              "\twithout a targetPressure (atm)!\n");      
77        painCave.isFatal = 1;
78        simError();
79      } else {
# Line 84 | Line 85 | namespace oopse {
85      if (simParams->getUsePeriodicBoundaryConditions()) {
86        sprintf(painCave.errMsg,
87                "SMIPDynamics error: You can't use the SMIPD integrator\n"
88 <              "   with periodic boundary conditions!\n");    
88 >              "\twith periodic boundary conditions!\n");    
89        painCave.isFatal = 1;
90        simError();
91      }
92      
93 <    if (!simParams->haveViscosity()) {
93 >    if (!simParams->haveThermalConductivity()) {
94        sprintf(painCave.errMsg,
95                "SMIPDynamics error: You can't use the SMIPD integrator\n"
96 <              "   without a viscosity!\n");
96 >              "\twithout a thermalConductivity (W m^-1 K^-1)!\n");
97        painCave.isFatal = 1;
98        painCave.severity = OOPSE_ERROR;
99        simError();
100      }else{
101 <      viscosity_ = simParams->getViscosity();
101 >      thermalConductivity_ = simParams->getThermalConductivity() *
102 >        OOPSEConstant::thermalConductivityConvert;
103      }
104 +
105 +    if (!simParams->haveThermalLength()) {
106 +      sprintf(painCave.errMsg,
107 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
108 +              "\twithout a thermalLength (Angstroms)!\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 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;              
125 >    Molecule::IntegrableObjectIterator  j;
126  
112 –    // Build a vector of integrable objects to determine if the are
113 –    // surface atoms
114 –
127      for (mol = info_->beginMolecule(i); mol != NULL;
128           mol = info_->nextMolecule(i)) {          
129        for (integrableObject = mol->beginIntegrableObject(j);
# Line 121 | Line 133 | namespace oopse {
133        }
134      }  
135    }  
136 <  
136 >  
137    void SMIPDForceManager::postCalculation(bool needStress){
138      SimInfo::MoleculeIterator i;
139      Molecule::IntegrableObjectIterator  j;
140      Molecule* mol;
141      StuntDouble* integrableObject;
142 <  
142 >  
143      // Compute surface Mesh
144 +
145      surfaceMesh_->computeHull(localSites_);
146 <    
146 >
147      // Get total area and number of surface stunt doubles
148      RealType area = surfaceMesh_->getArea();
149 <    int nSurfaceSDs = surfaceMesh_->getNs();        
149 >    int nSurfaceSDs = surfaceMesh_->getNs();
150      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
151      int nTriangles = sMesh.size();
152 <        
152 >
153      // Generate all of the necessary random forces
154      std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
155  
156 <    // Loop over the mesh faces and apply random force to each of the faces
156 >    // Loop over the mesh faces and apply external pressure to each
157 >    // of the faces
158      std::vector<Triangle>::iterator face;
159      std::vector<StuntDouble*>::iterator vertex;
160 <    int thisNumber = 0;
160 >    int thisFacet = 0;
161      for (face = sMesh.begin(); face != sMesh.end(); ++face){
162      
163        Triangle thisTriangle = *face;
# Line 151 | Line 165 | namespace oopse {
165        RealType thisArea = thisTriangle.getArea();
166        Vector3d unitNormal = thisTriangle.getNormal();
167        unitNormal.normalize();
154 –
168        Vector3d centroid = thisTriangle.getCentroid();
169        Vector3d facetVel = thisTriangle.getFacetVelocity();
170 <      RealType hydroLength = thisTriangle.getIncircleRadius() * 2.0 /
158 <        NumericConstant::PI;
159 <      
160 <      // gamma is the drag coefficient normal to the face of the triangle      
161 <      RealType gamma = viscosity_ * hydroLength *
162 <        OOPSEConstant::viscoConvert;
170 >      RealType thisMass = thisTriangle.getFacetMass();
171  
172 +      // gamma is the drag coefficient normal to the face of the triangle      
173 +      RealType gamma = thermalConductivity_ * thisMass * thisArea
174 +        / (2.0 * thermalLength_ * OOPSEConstant::kB);
175 +      
176        RealType extPressure = - (targetPressure_ * thisArea) /
177          OOPSEConstant::energyConvert;
178  
179 <      RealType randomForce = randNums[thisNumber++] * sqrt(gamma);
179 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
180        RealType dragForce = -gamma * dot(facetVel, unitNormal);
181  
182 <      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
182 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
183 >        unitNormal;
184        
185        // Apply triangle force to stuntdouble vertices
186 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
187 <        if ((*vertex) != NULL) {
186 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
187 >        if ((*vertex) != NULL){
188            Vector3d vertexForce = langevinForce / 3.0;
189 <          (*vertex)->addFrc(vertexForce);          
190 <          if ((*vertex)->isDirectional()) {
189 >          (*vertex)->addFrc(vertexForce);          
190 >          if ((*vertex)->isDirectional()){          
191              Vector3d vertexPos = (*vertex)->getPos();
192              Vector3d vertexCentroidVector = vertexPos - centroid;
193              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 183 | Line 196 | namespace oopse {
196        }
197      }
198      
199 +    veloMunge->removeComDrift();
200 +    veloMunge->removeAngularDrift();
201 +    
202      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
203      currSnapshot->setVolume(surfaceMesh_->getVolume());    
204      ForceManager::postCalculation(needStress);  
205    }
206 <
206 >  
207 >  
208    std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
209 <                                                             RealType variance) {
210 <
209 >                                                             RealType variance)
210 >  {
211 >    
212      // zero fill the random vector before starting:
213      std::vector<RealType> gaussRand;
214      gaussRand.resize(nTriangles);
215      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
216 <  
216 >    
217   #ifdef IS_MPI
218      if (worldRank == 0) {
219   #endif
220        for (int i = 0; i < nTriangles; i++) {
221 <        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
221 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
222        }
223   #ifdef IS_MPI
224      }
225   #endif
226 <
226 >    
227      // push these out to the other processors
228 <
228 >    
229   #ifdef IS_MPI
230      if (worldRank == 0) {
231 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
231 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
232      } else {
233 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
233 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
234      }
235   #endif
236 <  
236 >    
237      return gaussRand;
238    }
239   }

Diff Legend

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