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 3480 by chuckv, Tue Nov 25 16:32:47 2008 UTC vs.
Revision 3481 by gezelter, Wed Nov 26 14:26:17 2008 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"
46 – #include "hydrodynamics/Sphere.hpp"
47 – #include "hydrodynamics/Ellipsoid.hpp"
48 – #include "utils/ElementsTable.hpp"
45   #include "math/ConvexHull.hpp"
46   #include "math/Triangle.hpp"
47  
52 –
48   namespace oopse {
49    
50 <  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
50 >  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
51 >
52      simParams = info->getSimParams();
57 –    veloMunge = new Velocitizer(info);
53      
54      // Create Hull, Convex Hull for now, other options later.
55 +    
56      surfaceMesh_ = new ConvexHull();
61 –
57      
63 –    
58      /* Check that the simulation has target pressure and target
59         temperature set*/
60      
61      if (!simParams->haveTargetTemp()) {
62 <      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
62 >      sprintf(painCave.errMsg,
63 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
64 >              "   without a targetTemp!\n");      
65        painCave.isFatal = 1;
66        painCave.severity = OOPSE_ERROR;
67        simError();
# Line 74 | Line 70 | namespace oopse {
70      }
71      
72      if (!simParams->haveTargetPressure()) {
73 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
74 <              "   without a targetPressure!\n");
75 <      
73 >      sprintf(painCave.errMsg,
74 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
75 >              "   without a targetPressure!\n");      
76        painCave.isFatal = 1;
77        simError();
78      } else {
79 <      /* Convert pressure from atm -> amu/(fs^2*Ang)*/
80 <      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
79 >      // Convert pressure from atm -> amu/(fs^2*Ang)
80 >      targetPressure_ = simParams->getTargetPressure() /
81 >        OOPSEConstant::pressureConvert;
82      }
86 –
83    
84      if (simParams->getUsePeriodicBoundaryConditions()) {
85 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
86 <              "   with periodic boundary conditions !\n");
87 <      
85 >      sprintf(painCave.errMsg,
86 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
87 >              "   with periodic boundary conditions!\n");    
88        painCave.isFatal = 1;
89        simError();
90      }
91 <
91 >    
92      if (!simParams->haveViscosity()) {
93 <      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
93 >      sprintf(painCave.errMsg,
94 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
95 >              "   without a viscosity!\n");
96        painCave.isFatal = 1;
97        painCave.severity = OOPSE_ERROR;
98        simError();
99      }else{
100        viscosity_ = simParams->getViscosity();
101      }
102 <
102 >    
103      dt_ = simParams->getDt();
104 +  
105 +    variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
106  
107 –    
108 –
109 –    //Compute initial hull
110 –    /*
111 –    surfaceMesh_->computeHull(localSites_);
112 –    Area0_ = surfaceMesh_->getArea();
113 –    int nTriangles = surfaceMesh_->getNMeshElements();
114 –    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
115 –    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
116 –      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
117 –    //RealType eta0 = gamma_0/
118 –    */
119 –
120 –
107      Molecule* mol;
108      StuntDouble* integrableObject;
109      SimInfo::MoleculeIterator i;
110 <    Molecule::IntegrableObjectIterator  j;              
110 >    Molecule::IntegrableObjectIterator  j;              
111  
112 <    
112 >    // Build a vector of integrable objects to determine if the are
113 >    // surface atoms
114  
115 <    /* Compute hull first time through to get the area of t=0*/
116 <
117 <    //Build a vector of integrable objects to determine if the are surface atoms
118 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
132 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
115 >    for (mol = info_->beginMolecule(i); mol != NULL;
116 >         mol = info_->nextMolecule(i)) {          
117 >      for (integrableObject = mol->beginIntegrableObject(j);
118 >           integrableObject != NULL;
119             integrableObject = mol->nextIntegrableObject(j)) {  
120          localSites_.push_back(integrableObject);
121        }
122 <    }
137 <
138 <
122 >    }  
123    }  
124 <
141 <  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
142 <    std::map<std::string, HydroProp*> props;
143 <    std::ifstream ifs(filename.c_str());
144 <    if (ifs.is_open()) {
145 <      
146 <    }
147 <    
148 <    const unsigned int BufferSize = 65535;
149 <    char buffer[BufferSize];  
150 <    while (ifs.getline(buffer, BufferSize)) {
151 <      HydroProp* currProp = new HydroProp(buffer);
152 <      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
153 <    }
154 <
155 <    return props;
156 <  }
157 <  
124 >  
125    void SMIPDForceManager::postCalculation(bool needStress){
126      SimInfo::MoleculeIterator i;
127      Molecule::IntegrableObjectIterator  j;
128      Molecule* mol;
129      StuntDouble* integrableObject;
163 –    RealType mass;
164 –    Vector3d pos;
165 –    Vector3d frc;
166 –    Mat3x3d A;
167 –    Mat3x3d Atrans;
168 –    Vector3d Tb;
169 –    Vector3d ji;
170 –    unsigned int index = 0;
171 –    int fdf;
130    
131 <    fdf = 0;
174 <  
175 <    /*Compute surface Mesh*/
131 >    // Compute surface Mesh
132      surfaceMesh_->computeHull(localSites_);
177 –
178 –    /* Get area and number of surface stunt doubles and compute new variance */
179 –     RealType area = surfaceMesh_->getArea();
180 –     int nSurfaceSDs = surfaceMesh_->getNs();
181 –
133      
134 +    // Get total area and number of surface stunt doubles
135 +    RealType area = surfaceMesh_->getArea();
136 +    int nSurfaceSDs = surfaceMesh_->getNs();        
137      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
138      int nTriangles = sMesh.size();
139 +        
140 +    // Generate all of the necessary random forces
141 +    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
142  
143 <
187 <
188 <     /* Compute variance for random forces */
189 <    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);
190 <
191 <    /* Loop over the mesh faces and apply random force to each of the faces*/    
143 >    // Loop over the mesh faces and apply random force to each of the faces
144      std::vector<Triangle>::iterator face;
145      std::vector<StuntDouble*>::iterator vertex;
146      int thisNumber = 0;
# Line 197 | Line 149 | namespace oopse {
149        Triangle thisTriangle = *face;
150        std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
151        RealType thisArea = thisTriangle.getArea();
200 –      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
201 –      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
202 –
203 –      /* Get Random Force */
152        Vector3d unitNormal = thisTriangle.getNormal();
153        unitNormal.normalize();
154 <      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
154 >
155        Vector3d centroid = thisTriangle.getCentroid();
156        Vector3d facetVel = thisTriangle.getFacetVelocity();
157 <      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/NumericConstant::PI;
157 >      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;
163  
164 <      RealType f_normal = viscosity_*hydroLength*OOPSEConstant::viscoConvert;
165 <      RealType extPressure = -(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;
213 <      RealType randomForce = randNums[thisNumber++] * sqrt(2.0 * f_normal * OOPSEConstant::kb*targetTemp_/dt_);
164 >      RealType extPressure = - (targetPressure_ * thisArea) /
165 >        OOPSEConstant::energyConvert;
166  
167 <      RealType dragForce = -f_normal * dot(facetVel, unitNormal);
167 >      RealType randomForce = randNums[thisNumber++] * sqrt(gamma);
168 >      RealType dragForce = -gamma * dot(facetVel, unitNormal);
169  
217 –
170        Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
171        
172 <      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
173 <      
174 <      // std::cout << " " << randomForce << " " << f_normal <<   std::endl;
175 <
176 <      /* Apply triangle force to stuntdouble vertices */
177 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
226 <        if ((*vertex) != NULL){
227 <          // mass = integrableObject->getMass();
228 <          Vector3d vertexForce = langevinForce/3.0;
229 <          (*vertex)->addFrc(vertexForce);
230 <
231 <          if ((*vertex)->isDirectional()){
232 <
172 >      // Apply triangle force to stuntdouble vertices
173 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
174 >        if ((*vertex) != NULL) {
175 >          Vector3d vertexForce = langevinForce / 3.0;
176 >          (*vertex)->addFrc(vertexForce);          
177 >          if ((*vertex)->isDirectional()) {
178              Vector3d vertexPos = (*vertex)->getPos();
179              Vector3d vertexCentroidVector = vertexPos - centroid;
180              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 237 | Line 182 | namespace oopse {
182          }  
183        }
184      }
185 <
185 >    
186      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
187      currSnapshot->setVolume(surfaceMesh_->getVolume());    
188      ForceManager::postCalculation(needStress);  
189    }
190  
191 <  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
191 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
192 >                                                             RealType variance) {
193  
248 –    
249 –    Vector<RealType, 6> Z;
250 –    Vector<RealType, 6> generalForce;
251 –    
252 –
253 –    Z[0] = randNumGen_.randNorm(0, variance);
254 –    Z[1] = randNumGen_.randNorm(0, variance);
255 –    Z[2] = randNumGen_.randNorm(0, variance);
256 –    Z[3] = randNumGen_.randNorm(0, variance);
257 –    Z[4] = randNumGen_.randNorm(0, variance);
258 –    Z[5] = randNumGen_.randNorm(0, variance);
259 –    
260 –    generalForce = hydroProps_[index]->getS()*Z;
261 –    
262 –    force[0] = generalForce[0];
263 –    force[1] = generalForce[1];
264 –    force[2] = generalForce[2];
265 –    torque[0] = generalForce[3];
266 –    torque[1] = generalForce[4];
267 –    torque[2] = generalForce[5];
268 –    
269 – }
270 –  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
271 –
194      // zero fill the random vector before starting:
195      std::vector<RealType> gaussRand;
196      gaussRand.resize(nTriangles);
197      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
198 <
277 <    
278 <
198 >  
199   #ifdef IS_MPI
200      if (worldRank == 0) {
201   #endif
202        for (int i = 0; i < nTriangles; i++) {
203 <        //gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
284 <        gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);    
203 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
204        }
205   #ifdef IS_MPI
206      }
# Line 291 | Line 210 | namespace oopse {
210  
211   #ifdef IS_MPI
212      if (worldRank == 0) {
213 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
213 >      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
214      } else {
215 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
215 >      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
216      }
217   #endif
218    
300 –    for (int i = 0; i < nTriangles; i++) {
301 –      gaussRand[i] = gaussRand[i] * variance;
302 –    }
303 –  
219      return gaussRand;
220    }
306 –
307 –
221   }

Diff Legend

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