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 3504 by gezelter, Wed May 13 22:27:29 2009 UTC vs.
Revision 3516 by gezelter, Wed Jul 22 15:00:21 2009 UTC

# Line 41 | Line 41
41   #include <fstream>
42   #include <iostream>
43   #include "integrators/SMIPDForceManager.hpp"
44 – #include "math/CholeskyDecomposition.hpp"
44   #include "utils/OOPSEConstant.hpp"
45   #include "math/ConvexHull.hpp"
46   #include "math/Triangle.hpp"
# Line 51 | Line 50 | namespace oopse {
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.
# Line 63 | Line 63 | namespace oopse {
63      if (!simParams->haveTargetTemp()) {
64        sprintf(painCave.errMsg,
65                "SMIPDynamics error: You can't use the SMIPD integrator\n"
66 <              "   without a targetTemp!\n");      
66 >              "   without a targetTemp (K)!\n");      
67        painCave.isFatal = 1;
68        painCave.severity = OOPSE_ERROR;
69        simError();
# Line 74 | Line 74 | namespace oopse {
74      if (!simParams->haveTargetPressure()) {
75        sprintf(painCave.errMsg,
76                "SMIPDynamics error: You can't use the SMIPD integrator\n"
77 <              "   without a targetPressure!\n");      
77 >              "   without a targetPressure (atm)!\n");      
78        painCave.isFatal = 1;
79        simError();
80      } else {
# Line 86 | 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 (W m^-1 K^-1)!\n");
98        painCave.isFatal = 1;
99        painCave.severity = OOPSE_ERROR;
100        simError();
101      }else{
102 <      viscosity_ = simParams->getViscosity();
102 >      thermalConductivity_ = simParams->getThermalConductivity() *
103 >        OOPSEConstant::thermalConductivityConvert;
104      }
105 +
106 +    if (!simParams->haveThermalLength()) {
107 +      sprintf(painCave.errMsg,
108 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
109 +              "   without a thermalLength (Angstroms)!\n");
110 +      painCave.isFatal = 1;
111 +      painCave.severity = OOPSE_ERROR;
112 +      simError();
113 +    }else{
114 +      thermalLength_ = simParams->getThermalLength();
115 +    }
116      
117      dt_ = simParams->getDt();
118    
# Line 139 | Line 151 | namespace oopse {
151      int nTriangles = sMesh.size();
152  
153      // Generate all of the necessary random forces
154 <    std::vector<Vector3d>  randNums = genTriangleForces(nTriangles, variance_);
154 >    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
155  
156 +    RealType instaTemp = thermo->getTemperature();
157 +
158      // Loop over the mesh faces and apply external pressure to each
159      // of the faces
160      std::vector<Triangle>::iterator face;
# Line 155 | Line 169 | namespace oopse {
169        unitNormal.normalize();
170        Vector3d centroid = thisTriangle.getCentroid();
171        Vector3d facetVel = thisTriangle.getFacetVelocity();
172 +      RealType thisMass = thisTriangle.getFacetMass();
173  
174 <      Mat3x3d hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);
174 >      // gamma is the drag coefficient normal to the face of the triangle      
175 >      RealType gamma = thermalConductivity_ * thisMass * thisArea  
176 >        / (2.0 * thermalLength_ * OOPSEConstant::kB);
177 >      
178 >      gamma *= fabs(1.0 - targetTemp_/instaTemp);      
179 >      
180 >      RealType extPressure = - (targetPressure_ * thisArea) /
181 >        OOPSEConstant::energyConvert;
182  
183 <      hydroTensor *= OOPSEConstant::viscoConvert;
184 <      Mat3x3d S;
163 <      CholeskyDecomposition(hydroTensor, S);
183 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
184 >      RealType dragForce = -gamma * dot(facetVel, unitNormal);
185  
186 <      Vector3d extPressure = -unitNormal*(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;      
187 <      Vector3d randomForce = S * randNums[thisFacet++];
188 <      Vector3d dragForce = -hydroTensor * facetVel;
168 <
169 <      Vector3d langevinForce = (extPressure + randomForce + dragForce);
170 <
186 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
187 >        unitNormal;
188 >      
189        // Apply triangle force to stuntdouble vertices
190        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
191          if ((*vertex) != NULL){
# Line 181 | Line 199 | namespace oopse {
199          }  
200        }
201      }
202 <        
202 >    
203      veloMunge->removeComDrift();
204      veloMunge->removeAngularDrift();
205 <
205 >    
206      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
207      currSnapshot->setVolume(surfaceMesh_->getVolume());    
208      ForceManager::postCalculation(needStress);  
209    }
192 –
210    
211 <  std::vector<Vector3d> SMIPDForceManager::genTriangleForces(int nTriangles,
211 >  
212 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
213                                                               RealType variance)
214    {
215      
216      // zero fill the random vector before starting:
217 <    std::vector<Vector3d> gaussRand;
217 >    std::vector<RealType> gaussRand;
218      gaussRand.resize(nTriangles);
219 <    std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
219 >    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
220      
221   #ifdef IS_MPI
222      if (worldRank == 0) {
223   #endif
206 –      RealType rx, ry, rz;
224        for (int i = 0; i < nTriangles; i++) {
225 <        rx = randNumGen_.randNorm(0.0, variance);
209 <        ry = randNumGen_.randNorm(0.0, variance);
210 <        rz = randNumGen_.randNorm(0.0, variance);
211 <        gaussRand[i][0] = rx;
212 <        gaussRand[i][1] = ry;
213 <        gaussRand[i][2] = rz;
225 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
226        }
227   #ifdef IS_MPI
228      }
# Line 220 | Line 232 | namespace oopse {
232      
233   #ifdef IS_MPI
234      if (worldRank == 0) {
235 <      MPI_Bcast(&gaussRand[0], nTriangles*3 , MPI_REALTYPE, 0, MPI_COMM_WORLD);
235 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
236      } else {
237 <      MPI_Bcast(&gaussRand[0], nTriangles*3, MPI_REALTYPE, 0, MPI_COMM_WORLD);
237 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
238      }
239   #endif
240      

Diff Legend

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