ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NPTi.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NPTi.cpp (file contents):
Revision 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC vs.
Revision 763 by tim, Mon Sep 15 16:52:02 2003 UTC

# Line 9 | Line 9
9   #include "Integrator.hpp"
10   #include "simError.h"
11  
12 + #ifdef IS_MPI
13 + #include "mpiSimulation.hpp"
14 + #endif
15  
16 +
17   // Basic isotropic thermostating and barostating via the Melchionna
18   // modification of the Hoover algorithm:
19   //
# Line 20 | Line 24 | NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
24   //
25   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
26  
27 < NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
28 <  Integrator( theInfo, the_ff )
27 > template<typename T> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
28 >  T( theInfo, the_ff )
29   {
30    chi = 0.0;
31    eta = 0.0;
32 +  integralOfChidt = 0.0;
33    have_tau_thermostat = 0;
34    have_tau_barostat = 0;
35    have_target_temp = 0;
36    have_target_pressure = 0;
37 +  have_chi_tolerance = 0;
38 +  have_eta_tolerance = 0;
39 +  have_pos_iter_tolerance = 0;
40 +
41 +  oldPos = new double[3*nAtoms];
42 +  oldVel = new double[3*nAtoms];
43 +  oldJi = new double[3*nAtoms];
44 + #ifdef IS_MPI
45 +  Nparticles = mpiSim->getTotAtoms();
46 + #else
47 +  Nparticles = theInfo->n_atoms;
48 + #endif
49 +
50   }
51  
52 < void NPTi::moveA() {
52 > template<typename T> NPTi<T>::~NPTi() {
53 >  delete[] oldPos;
54 >  delete[] oldVel;
55 >  delete[] oldJi;
56 > }
57 >
58 > template<typename T> void NPTi<T>::moveA() {
59 >
60 >
61 > //   int i, j;
62 > //   DirectionalAtom* dAtom;
63 > //   double Tb[3], ji[3];
64 > //   double A[3][3], I[3][3];
65 > //   double angle, mass;
66 > //   double vel[3], pos[3], frc[3];
67 >
68 > //   double rj[3];
69 > //   double instaTemp, instaPress, instaVol;
70 > //   double tt2, tb2, scaleFactor;
71 >
72 > //   tt2 = tauThermostat * tauThermostat;
73 > //   tb2 = tauBarostat * tauBarostat;
74 >
75 > //   instaTemp = tStats->getTemperature();
76 > //   instaPress = tStats->getPressure();
77 > //   instaVol = tStats->getVolume();
78 >  
79 > //    // first evolve chi a half step
80    
81 <  int i,j,k;
82 <  int atomIndex, aMatIndex;
81 > //   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
82 > //   eta += dt2 * ( instaVol * (instaPress - targetPressure) /
83 > //               (p_convert*NkBT*tb2));
84 >
85 > //   integralOfChidt += dt2* chi;
86 >
87 > //   for( i=0; i<nAtoms; i++ ){
88 > //     atoms[i]->getVel( vel );
89 > //     atoms[i]->getPos( pos );
90 > //     atoms[i]->getFrc( frc );
91 >
92 > //     mass = atoms[i]->getMass();
93 >
94 > //     for (j=0; j < 3; j++) {
95 > //       vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
96 > //       rj[j] = pos[j];
97 > //     }
98 >
99 > //     atoms[i]->setVel( vel );
100 >
101 > //     info->wrapVector(rj);
102 >
103 > //     for (j = 0; j < 3; j++)
104 > //       pos[j] += dt * (vel[j] + eta*rj[j]);
105 >
106 > //     atoms[i]->setPos( pos );
107 >
108 > //     if( atoms[i]->isDirectional() ){
109 >
110 > //       dAtom = (DirectionalAtom *)atoms[i];
111 >          
112 > //       // get and convert the torque to body frame
113 >      
114 > //       dAtom->getTrq( Tb );
115 > //       dAtom->lab2Body( Tb );
116 >      
117 > //       // get the angular momentum, and propagate a half step
118 >
119 > //       dAtom->getJ( ji );
120 >
121 > //       for (j=0; j < 3; j++)
122 > //         ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
123 >      
124 > //       // use the angular velocities to propagate the rotation matrix a
125 > //       // full time step
126 >
127 > //       dAtom->getA(A);
128 > //       dAtom->getI(I);
129 >    
130 > //       // rotate about the x-axis      
131 > //       angle = dt2 * ji[0] / I[0][0];
132 > //       this->rotate( 1, 2, angle, ji, A );
133 >
134 > //       // rotate about the y-axis
135 > //       angle = dt2 * ji[1] / I[1][1];
136 > //       this->rotate( 2, 0, angle, ji, A );
137 >      
138 > //       // rotate about the z-axis
139 > //       angle = dt * ji[2] / I[2][2];
140 > //       this->rotate( 0, 1, angle, ji, A);
141 >      
142 > //       // rotate about the y-axis
143 > //       angle = dt2 * ji[1] / I[1][1];
144 > //       this->rotate( 2, 0, angle, ji, A );
145 >      
146 > //        // rotate about the x-axis
147 > //       angle = dt2 * ji[0] / I[0][0];
148 > //       this->rotate( 1, 2, angle, ji, A );
149 >      
150 > //       dAtom->setJ( ji );
151 > //       dAtom->setA( A  );    
152 > //     }                
153 >
154 > //   }
155 >
156 > //   // Scale the box after all the positions have been moved:
157 >  
158 > //   scaleFactor = exp(dt*eta);
159 >
160 > //   if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
161 > //     sprintf( painCave.errMsg,
162 > //              "NPTi error: Attempting a Box scaling of more than 10 percent"
163 > //              " check your tauBarostat, as it is probably too small!\n"
164 > //              " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
165 > //              );
166 > //     painCave.isFatal = 1;
167 > //     simError();
168 > //   } else {        
169 > //     info->scaleBox(exp(dt*eta));      
170 > //   }
171 >  
172 >
173 >  //new version of NPTi
174 >  int i, j, k;
175    DirectionalAtom* dAtom;
176 <  double Tb[3];
177 <  double ji[3];
176 >  double Tb[3], ji[3];
177 >  double A[3][3], I[3][3];
178 >  double angle, mass;
179 >  double vel[3], pos[3], frc[3];
180 >
181    double rj[3];
182    double instaTemp, instaPress, instaVol;
183 <  double tt2, tb2;
184 <  double angle;
183 >  double tt2, tb2, scaleFactor;
184 >  double COM[3];
185  
46 –
186    tt2 = tauThermostat * tauThermostat;
187    tb2 = tauBarostat * tauBarostat;
188  
189    instaTemp = tStats->getTemperature();
190    instaPress = tStats->getPressure();
191    instaVol = tStats->getVolume();
53 –  
54 –   // first evolve chi a half step
192    
193 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
194 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
195 <                 (p_convert*NkBT*tb2));
59 <
193 >  tStats->getCOM(COM);
194 >  
195 >  //evolve velocity half step
196    for( i=0; i<nAtoms; i++ ){
61 –    atomIndex = i * 3;
62 –    aMatIndex = i * 9;
63 –    
64 –    // velocity half step
65 –    for( j=atomIndex; j<(atomIndex+3); j++ )
66 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
67 –                       - vel[j]*(chi+eta));
197  
198 <    // position whole step    
198 >    atoms[i]->getVel( vel );
199 >    atoms[i]->getFrc( frc );
200  
201 <    rj[0] = pos[atomIndex];
72 <    rj[1] = pos[atomIndex+1];
73 <    rj[2] = pos[atomIndex+2];
74 <    
75 <    info->wrapVector(rj);
201 >    mass = atoms[i]->getMass();
202  
203 <    pos[atomIndex] += dt * (vel[atomIndex] + eta*rj[0]);
204 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + eta*rj[1]);
205 <    pos[atomIndex+2] += dt * (vel[atomIndex+2] + eta*rj[2]);
203 >    for (j=0; j < 3; j++) {
204 >      // velocity half step  (use chi from previous step here):
205 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
206 >  
207 >    }
208 >
209 >    atoms[i]->setVel( vel );
210    
211      if( atoms[i]->isDirectional() ){
212  
213        dAtom = (DirectionalAtom *)atoms[i];
214 <          
214 >
215        // get and convert the torque to body frame
216        
217 <      Tb[0] = dAtom->getTx();
88 <      Tb[1] = dAtom->getTy();
89 <      Tb[2] = dAtom->getTz();
90 <      
217 >      dAtom->getTrq( Tb );
218        dAtom->lab2Body( Tb );
219        
220        // get the angular momentum, and propagate a half step
221  
222 <      ji[0] = dAtom->getJx();
223 <      ji[1] = dAtom->getJy();
224 <      ji[2] = dAtom->getJz();
222 >      dAtom->getJ( ji );
223 >
224 >      for (j=0; j < 3; j++)
225 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
226        
99 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
100 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
101 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
102 –      
227        // use the angular velocities to propagate the rotation matrix a
228        // full time step
229 <      
229 >
230 >      dAtom->getA(A);
231 >      dAtom->getI(I);
232 >    
233        // rotate about the x-axis      
234 <      angle = dt2 * ji[0] / dAtom->getIxx();
235 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
236 <      
234 >      angle = dt2 * ji[0] / I[0][0];
235 >      this->rotate( 1, 2, angle, ji, A );
236 >
237        // rotate about the y-axis
238 <      angle = dt2 * ji[1] / dAtom->getIyy();
239 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
238 >      angle = dt2 * ji[1] / I[1][1];
239 >      this->rotate( 2, 0, angle, ji, A );
240        
241        // rotate about the z-axis
242 <      angle = dt * ji[2] / dAtom->getIzz();
243 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
242 >      angle = dt * ji[2] / I[2][2];
243 >      this->rotate( 0, 1, angle, ji, A);
244        
245        // rotate about the y-axis
246 <      angle = dt2 * ji[1] / dAtom->getIyy();
247 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
246 >      angle = dt2 * ji[1] / I[1][1];
247 >      this->rotate( 2, 0, angle, ji, A );
248        
249         // rotate about the x-axis
250 <      angle = dt2 * ji[0] / dAtom->getIxx();
251 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
250 >      angle = dt2 * ji[0] / I[0][0];
251 >      this->rotate( 1, 2, angle, ji, A );
252        
253 <      dAtom->setJx( ji[0] );
254 <      dAtom->setJy( ji[1] );
255 <      dAtom->setJz( ji[2] );
256 <    }
253 >      dAtom->setJ( ji );
254 >      dAtom->setA( A  );    
255 >    }    
256 >  }
257 >
258 >  // evolve chi and eta  half step
259 >  
260 >  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
261 >  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2));
262 >
263 >  //calculate the integral of chidt
264 >  integralOfChidt += dt2*chi;
265 >
266 >  //save the old positions
267 >  for(i = 0; i < nAtoms; i++){
268 >    atoms[i]->getPos(pos);
269 >    for(j = 0; j < 3; j++)
270 >      oldPos[i*3 + j] = pos[j];
271 >  }
272 >  
273 >  //the first estimation of r(t+dt) is equal to  r(t)
274      
275 +  for(k = 0; k < 4; k ++){
276 +
277 +    for(i =0 ; i < nAtoms; i++){
278 +
279 +      atoms[i]->getVel(vel);
280 +      atoms[i]->getPos(pos);
281 +
282 +      for(j = 0; j < 3; j++)
283 +        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
284 +
285 +      
286 +      //wrapVector(r(t)) = r(t)-R0
287 +      //info->wrapVector(rj);
288 +      
289 +      for(j = 0; j < 3; j++)
290 +        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]);
291 +
292 +      atoms[i]->setPos( pos );
293 +
294 +    }
295 +
296    }
297 +    
298 +
299    // Scale the box after all the positions have been moved:
300 +  
301 +  scaleFactor = exp(dt*eta);
302  
303 <  cerr << "eta = " << eta
304 <       << "; exp(dt*eta) = " << exp(eta*dt) << "\n";
305 <
306 <  info->scaleBox(exp(dt*eta));
303 >  if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
304 >    sprintf( painCave.errMsg,
305 >             "NPTi error: Attempting a Box scaling of more than 10 percent"
306 >             " check your tauBarostat, as it is probably too small!\n"
307 >             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
308 >             );
309 >    painCave.isFatal = 1;
310 >    simError();
311 >  } else {        
312 >    info->scaleBox(scaleFactor);      
313 >  }  
314  
315 +  //advance volume;
316 +  volume = volume * exp(dt*eta);
317   }
318  
319 < void NPTi::moveB( void ){
320 <  int i,j,k;
321 <  int atomIndex;
319 > template<typename T> void NPTi<T>::moveB( void ){
320 >
321 > /*
322 >  int i, j;
323    DirectionalAtom* dAtom;
324 <  double Tb[3];
325 <  double ji[3];
324 >  double Tb[3], ji[3];
325 >  double vel[3], frc[3];
326 >  double mass;
327 >
328    double instaTemp, instaPress, instaVol;
329    double tt2, tb2;
330    
# Line 157 | Line 338 | void NPTi::moveB( void ){
338    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
339    eta += dt2 * ( instaVol * (instaPress - targetPressure) /
340                   (p_convert*NkBT*tb2));
341 +  integralOfChidt += dt2*chi;
342    
343    for( i=0; i<nAtoms; i++ ){
344 <    atomIndex = i * 3;
345 <    
344 >
345 >    atoms[i]->getVel( vel );
346 >    atoms[i]->getFrc( frc );
347 >
348 >    mass = atoms[i]->getMass();
349 >
350      // velocity half step
351 <    for( j=atomIndex; j<(atomIndex+3); j++ )
352 <    for( j=atomIndex; j<(atomIndex+3); j++ )
167 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
168 <                       - vel[j]*(chi+eta));
351 >    for (j=0; j < 3; j++)
352 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
353      
354 +    atoms[i]->setVel( vel );
355 +
356      if( atoms[i]->isDirectional() ){
357 <      
357 >
358        dAtom = (DirectionalAtom *)atoms[i];
359 <      
360 <      // get and convert the torque to body frame
361 <      
362 <      Tb[0] = dAtom->getTx();
177 <      Tb[1] = dAtom->getTy();
178 <      Tb[2] = dAtom->getTz();
179 <      
359 >
360 >      // get and convert the torque to body frame      
361 >
362 >      dAtom->getTrq( Tb );
363        dAtom->lab2Body( Tb );
364 +
365 +      // get the angular momentum, and propagate a half step
366 +
367 +      dAtom->getJ( ji );
368 +
369 +      for (j=0; j < 3; j++)
370 +        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);    
371 +
372 +      dAtom->setJ( ji );
373 +    }
374 +  }
375 +
376 + */
377 +  
378 +  //new version of NPTi
379 +  int i, j, k;
380 +  DirectionalAtom* dAtom;
381 +  double Tb[3], ji[3];
382 +  double vel[3], frc[3];
383 +  double mass;
384 +
385 +  double instTemp, instPress, instVol;
386 +  double tt2, tb2;
387 +  double oldChi, prevChi;
388 +  double oldEta, preEta;
389 +  
390 +  tt2 = tauThermostat * tauThermostat;
391 +  tb2 = tauBarostat * tauBarostat;
392 +
393 +
394 +  // Set things up for the iteration:
395 +
396 +  oldChi = chi;
397 +  oldEta = eta;
398 +
399 +  for( i=0; i<nAtoms; i++ ){
400 +
401 +    atoms[i]->getVel( vel );
402 +
403 +    for (j=0; j < 3; j++)
404 +      oldVel[3*i + j]  = vel[j];
405 +
406 +    if( atoms[i]->isDirectional() ){
407 +
408 +      dAtom = (DirectionalAtom *)atoms[i];
409 +
410 +      dAtom->getJ( ji );
411 +
412 +      for (j=0; j < 3; j++)
413 +        oldJi[3*i + j] = ji[j];
414 +
415 +    }
416 +  }
417 +
418 +  // do the iteration:
419 +
420 +  instVol = tStats->getVolume();
421 +  
422 +  for (k=0; k < 4; k++) {
423 +    
424 +    instTemp = tStats->getTemperature();
425 +    instPress = tStats->getPressure();
426 +
427 +    // evolve chi another half step using the temperature at t + dt/2
428 +
429 +    prevChi = chi;
430 +    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
431 +      (tauThermostat*tauThermostat);
432 +
433 +    preEta = eta;
434 +    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
435 +       (p_convert*NkBT*tb2));
436 +
437 +  
438 +    for( i=0; i<nAtoms; i++ ){
439 +
440 +      atoms[i]->getFrc( frc );
441 +      atoms[i]->getVel(vel);
442        
443 <      // get the angular momentum, and complete the angular momentum
183 <      // half step
443 >      mass = atoms[i]->getMass();
444        
445 <      ji[0] = dAtom->getJx();
446 <      ji[1] = dAtom->getJy();
447 <      ji[2] = dAtom->getJz();
445 >      // velocity half step
446 >      for (j=0; j < 3; j++)
447 >        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
448        
449 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
190 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
191 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
449 >      atoms[i]->setVel( vel );
450        
451 <      dAtom->setJx( ji[0] );
452 <      dAtom->setJy( ji[1] );
453 <      dAtom->setJz( ji[2] );
451 >      if( atoms[i]->isDirectional() ){
452 >
453 >        dAtom = (DirectionalAtom *)atoms[i];
454 >  
455 >        // get and convert the torque to body frame      
456 >  
457 >        dAtom->getTrq( Tb );
458 >        dAtom->lab2Body( Tb );      
459 >            
460 >        for (j=0; j < 3; j++)
461 >          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
462 >      
463 >          dAtom->setJ( ji );
464 >      }
465      }
466 +
467 +    if (fabs(prevChi - chi) <= chiTolerance && fabs(preEta -eta) <= etaTolerance)
468 +      break;
469    }
470 +
471 +  //calculate integral of chida
472 +  integralOfChidt += dt2*chi;
473 +
474 +
475   }
476  
477 < int NPTi::readyCheck() {
477 > template<typename T> void NPTi<T>::resetIntegrator() {
478 >  chi = 0.0;
479 >  eta = 0.0;
480 > }
481 >
482 > template<typename T> int NPTi<T>::readyCheck() {
483 >
484 >  //check parent's readyCheck() first
485 >  if (T::readyCheck() == -1)
486 >    return -1;
487  
488    // First check to see if we have a target temperature.
489    // Not having one is fatal.
# Line 244 | Line 530 | int NPTi::readyCheck() {
530      return -1;
531    }    
532  
533 +  if (!have_chi_tolerance) {
534 +    sprintf( painCave.errMsg,
535 +             "NPTi warning: setting chi tolerance to 1e-6\n");
536 +    chiTolerance = 1e-6;
537 +    have_chi_tolerance = 1;
538 +    painCave.isFatal = 0;
539 +    simError();
540 +  }
541 +
542 +    if (!have_eta_tolerance) {
543 +    sprintf( painCave.errMsg,
544 +             "NPTi warning: setting eta tolerance to 1e-6\n");
545 +    etaTolerance = 1e-6;
546 +    have_eta_tolerance = 1;
547 +    painCave.isFatal = 0;
548 +    simError();
549 +  }
550    // We need NkBT a lot, so just set it here:
551  
552 <  NkBT = (double)info->ndf * kB * targetTemp;
552 >  NkBT = (double)Nparticles * kB * targetTemp;
553 >  fkBT = (double)info->ndf * kB * targetTemp;
554  
555    return 1;
556   }
557 +
558 + template<typename T> double NPTi<T>::getConservedQuantity(void){
559 +
560 +  double conservedQuantity;
561 +  double tb2;
562 +  double eta2;  
563 +  double E_NPT;
564 +  double U;
565 +  double TS;
566 +  double PV;
567 +  double extra;
568 +
569 +  static double pre_U;
570 +  static double pre_TS;
571 +  static double pre_PV;
572 +  static double pre_extra;
573 +  static int hackCount = 0;
574 +
575 +  double delta_U;
576 +  double delta_TS;
577 +  double delta_PV;
578 +  double delta_extra;
579 +
580 +  U = tStats->getTotalE();
581 +
582 +  TS = fkBT *
583 +    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
584 +
585 +  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
586 +
587 +  tb2 = tauBarostat * tauBarostat;
588 +  eta2 = eta * eta;
589 +
590 +  extra = (fkBT * tb2 * eta2 / 2.0 ) / eConvert;
591 +  /*
592 +  if(hackCount == 0){
593 +    pre_U = U;
594 +    pre_TS =TS;
595 +    pre_PV = PV;
596 +    pre_extra =extra;
597 +    hackCount ++;
598 +  }
599 +
600 +  delta_U = U - pre_U;
601 +  delta_TS = TS - pre_TS;
602 +  delta_PV = PV - pre_PV;
603 +  delta_extra = extra - pre_extra;
604 + */
605 +  cout.width(8);
606 +  cout.precision(8);
607 +
608 +  
609 +  cout << info->getTime() << "\t"
610 +       << chi << "\t"
611 +       << eta << "\t"
612 +       << U << "\t"
613 +       << TS << "\t"
614 +       << PV << "\t"
615 +       << extra << "\t"
616 +       << U+TS+PV+extra << endl;
617 +
618 + /*
619 +    pre_U = U;
620 +    pre_TS =TS;
621 +    pre_PV = PV;
622 +    pre_extra =extra;
623 +
624 +
625 +  cout << info->getTime() << "\t"
626 +       << U << "\t"
627 +       << U+TS << "\t"
628 +       << U+TS+PV << "\t"
629 +       << U+TS+PV+extra << endl;
630 + */
631 +  conservedQuantity = U+TS+PV+extra;
632 +  return conservedQuantity;
633 + }

Diff Legend

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