ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/hydrodynamics/ApproximationModel.cpp
(Generate patch)

Comparing trunk/src/applications/hydrodynamics/ApproximationModel.cpp (file contents):
Revision 972 by gezelter, Wed May 24 18:31:12 2006 UTC vs.
Revision 1177 by xsun, Tue Aug 14 17:40:33 2007 UTC

# Line 69 | Line 69 | namespace oopse {
69      
70    }
71    
72 <  bool ApproximationModel::calcHydroProps(Shape* shape, double viscosity, double temperature) {
72 >  bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73      
74      bool ret = true;
75 <    HydroProps cr;
76 <    HydroProps cd;
75 >    HydroProp* cr = new HydroProp();
76 >    HydroProp* cd = new HydroProp();
77      calcHydroPropsAtCR(beads_, viscosity, temperature, cr);
78 <    //calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
78 >    calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
79      setCR(cr);
80      setCD(cd);
81 –    
81      return true;    
82    }
83    
84 <  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
84 >  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cr) {
85      
86      int nbeads = beads.size();
87 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
88 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
87 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
88 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
89      Mat3x3d I;
90      I(0, 0) = 1.0;
91      I(1, 1) = 1.0;
# Line 97 | Line 96 | namespace oopse {
96          Mat3x3d Tij;
97              if (i != j ) {
98                Vector3d Rij = beads[i].pos - beads[j].pos;
99 <              double rij = Rij.length();
100 <              double rij2 = rij * rij;
101 <              double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
99 >              RealType rij = Rij.length();
100 >              RealType rij2 = rij * rij;
101 >              RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
102                Mat3x3d tmpMat;
103                tmpMat = outProduct(Rij, Rij) / rij2;
104 <              double constant = 8.0 * NumericConstant::PI * viscosity * rij;
105 <              Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
104 >              RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
105 >              RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
106 >              RealType tmp2 = 1.0 - sumSigma2OverRij2;
107 >              Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
108              }else {
109 <              double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
109 >              RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
110                Tij(0, 0) = constant;
111                Tij(1, 1) = constant;
112                Tij(2, 2) = constant;
# Line 132 | Line 133 | namespace oopse {
133      
134      //calculate the total volume
135      
136 <    double volume = 0.0;
136 >    RealType volume = 0.0;
137      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
138        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
139      }
# Line 149 | Line 150 | namespace oopse {
150        }
151      }
152      
153 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
153 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
154      Xiott *= convertConstant;
155      Xiotr *= convertConstant;
156      Xiorr *= convertConstant;
# Line 185 | Line 186 | namespace oopse {
186      Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose();
187      
188  
189 <    SquareMatrix<double,6> Xir6x6;
190 <    SquareMatrix<double,6> Dr6x6;
189 >    SquareMatrix<RealType,6> Xir6x6;
190 >    SquareMatrix<RealType,6> Dr6x6;
191  
192      Xir6x6.setSubMatrix(0, 0, Xirtt);
193      Xir6x6.setSubMatrix(0, 3, Xirtr.transpose());
# Line 202 | Line 203 | namespace oopse {
203      Dr6x6.getSubMatrix(0, 3, Drrt);
204      Dr6x6.getSubMatrix(3, 0, Drtr);
205      Dr6x6.getSubMatrix(3, 3, Drrr);
206 <    double kt = OOPSEConstant::kB * temperature ;
206 >    RealType kt = OOPSEConstant::kB * temperature ;
207      Drtt *= kt;
208      Drrt *= kt;
209      Drtr *= kt;
# Line 211 | Line 212 | namespace oopse {
212      Xirtr *= OOPSEConstant::kb * temperature;
213      Xirrr *= OOPSEConstant::kb * temperature;
214      
215 +    Mat6x6d Xi, D;
216  
217 <    cr.center = ror;
218 <    cr.Xi.setSubMatrix(0, 0, Xirtt);
219 <    cr.Xi.setSubMatrix(0, 3, Xirtr);
220 <    cr.Xi.setSubMatrix(3, 0, Xirtr);
221 <    cr.Xi.setSubMatrix(3, 3, Xirrr);
222 <    cr.D.setSubMatrix(0, 0, Drtt);
223 <    cr.D.setSubMatrix(0, 3, Drrt);
224 <    cr.D.setSubMatrix(3, 0, Drtr);
225 <    cr.D.setSubMatrix(3, 3, Drrr);    
217 >    cr->setCOR(ror);
218 >
219 >    Xi.setSubMatrix(0, 0, Xirtt);
220 >    Xi.setSubMatrix(0, 3, Xirtr);
221 >    Xi.setSubMatrix(3, 0, Xirtr);
222 >    Xi.setSubMatrix(3, 3, Xirrr);
223 >
224 >    cr->setXi(Xi);
225 >
226 >    D.setSubMatrix(0, 0, Drtt);
227 >    D.setSubMatrix(0, 3, Drrt);
228 >    D.setSubMatrix(3, 0, Drtr);
229 >    D.setSubMatrix(3, 3, Drrr);    
230 >
231 >    cr->setD(D);
232      
233      std::cout << "-----------------------------------------\n";
234      std::cout << "center of resistance :" << std::endl;
# Line 246 | Line 254 | namespace oopse {
254      return true;
255   }
256    
257 <  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
257 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cd) {
258      
259      int nbeads = beads.size();
260 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
261 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
260 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
261 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
262      Mat3x3d I;
263      I(0, 0) = 1.0;
264      I(1, 1) = 1.0;
# Line 261 | Line 269 | namespace oopse {
269          Mat3x3d Tij;
270          if (i != j ) {
271            Vector3d Rij = beads[i].pos - beads[j].pos;
272 <          double rij = Rij.length();
273 <          double rij2 = rij * rij;
274 <          double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
272 >          RealType rij = Rij.length();
273 >          RealType rij2 = rij * rij;
274 >          RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
275            Mat3x3d tmpMat;
276            tmpMat = outProduct(Rij, Rij) / rij2;
277 <          double constant = 8.0 * NumericConstant::PI * viscosity * rij;
278 <          Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
277 >          RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
278 >          RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
279 >          RealType tmp2 = 1.0 - sumSigma2OverRij2;
280 >          Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
281          }else {
282 <          double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
282 >          RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
283            Tij(0, 0) = constant;
284            Tij(1, 1) = constant;
285            Tij(2, 2) = constant;
# Line 296 | Line 306 | namespace oopse {
306  
307      //calculate the total volume
308  
309 <    double volume = 0.0;
309 >    RealType volume = 0.0;
310      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
311        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
312      }
# Line 313 | Line 323 | namespace oopse {
323        }
324      }
325      
326 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
326 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
327      Xitt *= convertConstant;
328      Xitr *= convertConstant;
329      Xirr *= convertConstant;
330      
331 <    double kt = OOPSEConstant::kB * temperature;
331 >    RealType kt = OOPSEConstant::kB * temperature;
332      
333      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
334      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
# Line 377 | Line 387 | namespace oopse {
387      Ddrr = Dorr;
388      Ddtr = Dotr + Dorr * Uod;
389  
390 <    SquareMatrix<double, 6> Dd;
390 >    SquareMatrix<RealType, 6> Dd;
391      Dd.setSubMatrix(0, 0, Ddtt);
392      Dd.setSubMatrix(0, 3, Ddtr.transpose());
393      Dd.setSubMatrix(3, 0, Ddtr);
394      Dd.setSubMatrix(3, 3, Ddrr);    
395 <    SquareMatrix<double, 6> Xid;
395 >    SquareMatrix<RealType, 6> Xid;
396      Ddtt *= kt;
397      Ddtr *=kt;
398      Ddrr *= kt;
# Line 394 | Line 404 | namespace oopse {
404      //Xid /= OOPSEConstant::energyConvert;
405      Xid *= OOPSEConstant::kb * temperature;
406  
407 <    cr.center = rod;
398 <    cr.D.setSubMatrix(0, 0, Ddtt);
399 <    cr.D.setSubMatrix(0, 3, Ddtr);
400 <    cr.D.setSubMatrix(3, 0, Ddtr);
401 <    cr.D.setSubMatrix(3, 3, Ddrr);
402 <    cr.Xi = Xid;
407 >    Mat6x6d Xi, D;
408  
409 +    cd->setCOR(rod);
410 +
411 +    cd->setXi(Xid);
412 +
413 +    D.setSubMatrix(0, 0, Ddtt);
414 +    D.setSubMatrix(0, 3, Ddtr);
415 +    D.setSubMatrix(3, 0, Ddtr);
416 +    D.setSubMatrix(3, 3, Ddrr);
417 +
418 +    cd->setD(D);
419 +
420      std::cout << "viscosity = " << viscosity << std::endl;
421      std::cout << "temperature = " << temperature << std::endl;
422      std::cout << "center of diffusion :" << std::endl;

Diff Legend

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