ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/hydrodynamics/HydrodynamicsModel.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 2596 by tim, Wed Feb 22 20:35:16 2006 UTC vs.
Revision 2598 by tim, Fri Feb 24 21:17:05 2006 UTC

# Line 43 | Line 43 | namespace oopse {
43   #include "math/LU.hpp"
44   #include "math/DynamicRectMatrix.hpp"
45   #include "math/SquareMatrix3.hpp"
46 + #include "utils/OOPSEConstant.hpp"
47   namespace oopse {
48   /**
49   * Reference:
# Line 50 | Line 51 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
51   * Comparison of Different Modeling and Computational Procedures.
52   * Biophysical Journal, 75(6), 3044, 1999
53   */
54 < bool HydrodynamicsModel::calcHydrodyanmicsProps(double eta) {
54 >
55 > HydrodynamicsModel::HydrodynamicsModel(StuntDouble* sd, const DynamicProperty& extraParams) : sd_(sd){
56 >    DynamicProperty::const_iterator iter;
57 >
58 >    iter = extraParams.find("Viscosity");
59 >    if (iter != extraParams.end()) {
60 >        boost::any param = iter->second;
61 >        viscosity_ = boost::any_cast<double>(param);
62 >    }else {
63 >        std::cout << "HydrodynamicsModel Error\n" ;
64 >    }
65 >
66 >    iter = extraParams.find("Temperature");
67 >    if (iter != extraParams.end()) {
68 >        boost::any param = iter->second;
69 >        temperature_ = boost::any_cast<double>(param);
70 >    }else {
71 >        std::cout << "HydrodynamicsModel Error\n" ;
72 >    }    
73 > }
74 >
75 > bool HydrodynamicsModel::calcHydrodyanmicsProps() {
76      if (!createBeads(beads_)) {
77          std::cout << "can not create beads" << std::endl;
78          return false;
# Line 60 | Line 82 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
82      DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
83      DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
84      Mat3x3d I;
85 +    I(0, 0) = 1.0;
86 +    I(1, 1) = 1.0;
87 +    I(2, 2) = 1.0;
88 +    
89      for (std::size_t i = 0; i < nbeads; ++i) {
90          for (std::size_t j = 0; j < nbeads; ++j) {
91              Mat3x3d Tij;
# Line 69 | Line 95 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
95                  double rij2 = rij * rij;
96                  double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
97                  Mat3x3d tmpMat;
98 <                tmpMat = outProduct(beads_[i].pos, beads_[j].pos) / rij2;
99 <                double constant = 8.0 * NumericConstant::PI * eta * rij;
98 >                tmpMat = outProduct(Rij, Rij) / rij2;
99 >                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
100                  Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
101              }else {
102 <                double constant = 1.0 / (6.0 * NumericConstant::PI * eta * beads_[i].radius);
102 >                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
103                  Tij(0, 0) = constant;
104                  Tij(1, 1) = constant;
105                  Tij(2, 2) = constant;
106              }
107              B.setSubMatrix(i*3, j*3, Tij);
108 +            std::cout << Tij << std::endl;
109          }
110      }
111  
112 +    std::cout << "B=\n"
113 +                  << B << std::endl;
114      //invert B Matrix
115      invertMatrix(B, C);
116 <    
116 >
117 >    std::cout << "C=\n"
118 >                  << C << std::endl;
119 >
120      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
121      std::vector<Mat3x3d> U;
122      for (int i = 0; i < nbeads; ++i) {
# Line 97 | Line 129 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
129      Mat3x3d Xitt;
130      Mat3x3d Xirr;
131      Mat3x3d Xitr;
132 +
133 +    //calculate the total volume
134 +
135 +    double volume = 0.0;
136 +    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
137 +        volume = 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
138 +    }
139          
140      for (std::size_t i = 0; i < nbeads; ++i) {
141          for (std::size_t j = 0; j < nbeads; ++j) {
# Line 104 | Line 143 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
143              C.getSubMatrix(i*3, j*3, Cij);
144              
145              Xitt += Cij;
146 <            Xirr += U[i] * Cij;
147 <            Xitr += U[i] * Cij * U[j];            
146 >            Xitr += U[i] * Cij;
147 >            Xirr += -U[i] * Cij * U[j];            
148 >            //Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;            
149          }
150      }
151  
# Line 115 | Line 155 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
155      Xi.setSubMatrix(0, 0, Xitt);
156      Xi.setSubMatrix(0, 3, Xitr.transpose());
157      Xi.setSubMatrix(3, 0, Xitr);
158 <    Xi.setSubMatrix(3, 3, Xitt);
159 <    invertMatrix(Xi, Do);
158 >    Xi.setSubMatrix(3, 3, Xirr);
159 >    //invertMatrix(Xi, Do);
160 >    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
161 >    //Do *= kt;    
162  
163 +
164      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
165      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
166      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
124 –    Do.getSubMatrix(0, 0 , Dott);
125 –    Do.getSubMatrix(3, 0, Dotr);
126 –    Do.getSubMatrix(3, 3, Dorr);
167  
168 +    const static Mat3x3d zeroMat(0.0);
169 +    
170 +    Mat3x3d XittInv(0.0);
171 +    XittInv = Xitt.inverse();
172 +    
173 +    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
174 +    //( a11 a12 0)
175 +    //( a21 a22 0)
176 +    //( 0    0    0)
177 +    Mat3x3d XirrInv;
178 +    XirrInv = Xirr.inverse();
179 +
180 +    Mat3x3d tmp;
181 +    Mat3x3d tmpInv;
182 +    tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
183 +    tmpInv = tmp.inverse();
184 +
185 +    Dott = kt * tmpInv;
186 +    Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8;
187 +
188 +    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
189 +    tmpInv = tmp.inverse();
190 +    
191 +    Dorr = kt * tmpInv*1.0E16;
192 +
193 +    //Do.getSubMatrix(0, 0 , Dott);
194 +    //Do.getSubMatrix(3, 0, Dotr);
195 +    //Do.getSubMatrix(3, 3, Dorr);
196 +
197      //calculate center of diffusion
198 <    Mat3x3d tmpMat;
199 <    tmpMat(0, 0) = Dorr(1, 1) + Dorr(2, 2);
200 <    tmpMat(0, 1) = - Dorr(0, 1);
201 <    tmpMat(0, 2) = -Dorr(0, 2);
202 <    tmpMat(1, 0) = -Dorr(0, 1);
203 <    tmpMat(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
204 <    tmpMat(1, 2) = -Dorr(1, 2);
205 <    tmpMat(2, 0) = -Dorr(0, 2);
206 <    tmpMat(2, 1) = -Dorr(1, 2);
138 <    tmpMat(2, 2) = Dorr(1, 1) + Dorr(0, 0);
198 >    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
199 >    tmp(0, 1) = - Dorr(0, 1);
200 >    tmp(0, 2) = -Dorr(0, 2);
201 >    tmp(1, 0) = -Dorr(0, 1);
202 >    tmp(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
203 >    tmp(1, 2) = -Dorr(1, 2);
204 >    tmp(2, 0) = -Dorr(0, 2);
205 >    tmp(2, 1) = -Dorr(1, 2);
206 >    tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
207  
208      Vector3d tmpVec;
209      tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
210      tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
211      tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
144 –        
145 –    Vector3d rod = tmpMat.inverse() * tmpVec;
212  
213 +    tmpInv = tmp.inverse();
214 +    
215 +    Vector3d rod = tmpInv * tmpVec;
216 +
217      //calculate Diffusion Tensor at center of diffusion
218      Mat3x3d Uod;
219      Uod.setupSkewMat(rod);
# Line 165 | Line 235 | void HydrodynamicsModel::writeBeads(std::ostream& os)
235   }
236  
237   void HydrodynamicsModel::writeBeads(std::ostream& os) {
238 +    std::vector<BeadParam>::iterator iter;
239 +    os << beads_.size() << std::endl;
240 +    os << "Generated by Hydro" << std::endl;
241 +    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
242 +        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
243 +    }
244  
245   }
246  
247   void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
248 +    os << "//viscosity = " << viscosity_ << std::endl;
249 +    os << "//temperature = " << temperature_<< std::endl;
250 +    std::vector<BeadParam>::iterator iter;
251 +    os << sd_->getType() << "\n";
252  
253 +    os << "//diffusion center" << std::endl;
254 +    os << props_.diffCenter << std::endl;
255 +
256 +    os << "//translational diffusion tensor" << std::endl;
257 +    os << props_.transDiff << std::endl;
258 +
259 +    os << "//translation-rotation coupling diffusion tensor" << std::endl;
260 +    os << props_.transRotDiff << std::endl;
261 +
262 +    os << "//rotational diffusion tensor" << std::endl;
263 +    os << props_.rotDiff << std::endl;
264 +    
265 +    /*
266 +    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
267 +
268 +    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
269 +        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
270 +        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
271 +    
272 +    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
273 +        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
274 +        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
275 +
276 +    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
277 +        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
278 +        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
279 +        << std::endl;
280 +    */
281   }
282  
283   }

Diff Legend

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