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 3515 by gezelter, Mon Jul 20 17:50:40 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 +    thermo = new Thermo(info);
54 +    veloMunge = new Velocitizer(info);
55      
56      // Create Hull, Convex Hull for now, other options later.
57      
58      surfaceMesh_ = new ConvexHull();
59      
60      /* Check that the simulation has target pressure and target
61 <       temperature set*/
61 >       temperature set */
62      
63      if (!simParams->haveTargetTemp()) {
64        sprintf(painCave.errMsg,
# Line 84 | 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 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      surfaceMesh_->computeHull(localSites_);
145 <    
145 >
146      // Get total area and number of surface stunt doubles
147      RealType area = surfaceMesh_->getArea();
148 <    int nSurfaceSDs = surfaceMesh_->getNs();        
148 >    int nSurfaceSDs = surfaceMesh_->getNs();
149      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
150      int nTriangles = sMesh.size();
151 <        
151 >
152      // Generate all of the necessary random forces
153      std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
154  
155 <    // Loop over the mesh faces and apply random force to each of the faces
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 151 | Line 166 | namespace oopse {
166        RealType thisArea = thisTriangle.getArea();
167        Vector3d unitNormal = thisTriangle.getNormal();
168        unitNormal.normalize();
154 –
169        Vector3d centroid = thisTriangle.getCentroid();
170        Vector3d facetVel = thisTriangle.getFacetVelocity();
171 <      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;
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[thisNumber++] * sqrt(gamma);
182 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
183        RealType dragForce = -gamma * dot(facetVel, unitNormal);
184  
185 <      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
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) {
189 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
190 >        if ((*vertex) != NULL){
191            Vector3d vertexForce = langevinForce / 3.0;
192 <          (*vertex)->addFrc(vertexForce);          
193 <          if ((*vertex)->isDirectional()) {
192 >          (*vertex)->addFrc(vertexForce);          
193 >          if ((*vertex)->isDirectional()){          
194              Vector3d vertexPos = (*vertex)->getPos();
195              Vector3d vertexCentroidVector = vertexPos - centroid;
196              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 183 | Line 199 | namespace oopse {
199        }
200      }
201      
202 +    veloMunge->removeComDrift();
203 +    veloMunge->removeAngularDrift();
204 +    
205      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
206      currSnapshot->setVolume(surfaceMesh_->getVolume());    
207      ForceManager::postCalculation(needStress);  
208    }
209 <
209 >  
210 >  
211    std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
212 <                                                             RealType variance) {
213 <
212 >                                                             RealType variance)
213 >  {
214 >    
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 <  
219 >    
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);
224 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
225        }
226   #ifdef IS_MPI
227      }
228   #endif
229 <
229 >    
230      // push these out to the other processors
231 <
231 >    
232   #ifdef IS_MPI
233      if (worldRank == 0) {
234 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
234 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
235      } else {
236 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
236 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
237      }
238   #endif
239 <  
239 >    
240      return gaussRand;
241    }
242   }

Diff Legend

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