step/stepcore/gas.cc

0001 /*
0002     SPDX-FileCopyrightText: 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
0003
0004     SPDX-License-Identifier: GPL-2.0-or-later
0005 */
0006
0007 #include "gas.h"
0008 #include "types.h"
0009 #include "constants.h"
0010 #include <cstdlib>
0011
0012 namespace StepCore
0013 {
0014
0015 STEPCORE_META_OBJECT(GasParticle, QT_TRANSLATE_NOOP("ObjectClass", "GasParticle"), QT_TRANSLATE_NOOP("ObjectDescription", "Gas particle"), 0, STEPCORE_SUPER_CLASS(Particle),)
0016
0017 STEPCORE_META_OBJECT(GasLJForce, QT_TRANSLATE_NOOP("ObjectClass", "GasLJForce"), QT_TRANSLATE_NOOP("ObjectDescription", "Lennard-Jones force"), 0,
0018     STEPCORE_SUPER_CLASS(Item) STEPCORE_SUPER_CLASS(Force),
0019     STEPCORE_PROPERTY_RW(double, depth, QT_TRANSLATE_NOOP("PropertyName", "depth"), QT_TRANSLATE_NOOP("Units", "J"), QT_TRANSLATE_NOOP("PropertyDescription", "Potential depth"), depth, setDepth)
0020     STEPCORE_PROPERTY_RW(double, rmin, QT_TRANSLATE_NOOP("PropertyName", "rmin"), QT_TRANSLATE_NOOP("Units", "m"), QT_TRANSLATE_NOOP("PropertyDescription", "Distance at which the force is zero"), rmin, setRmin)
0021     STEPCORE_PROPERTY_RW(double, cutoff, QT_TRANSLATE_NOOP("PropertyName", "cutoff"), QT_TRANSLATE_NOOP("Units", "m"), QT_TRANSLATE_NOOP("PropertyDescription", "Cut-off distance"), cutoff, setCutoff))
0022
0023 STEPCORE_META_OBJECT(GasLJForceErrors, QT_TRANSLATE_NOOP("ObjectClass", "GasLJForceErrors"), QT_TRANSLATE_NOOP("ObjectDescription", "Errors class for GasLJForce"), 0,
0024     STEPCORE_SUPER_CLASS(ObjectErrors),
0025     STEPCORE_PROPERTY_RW(double, depthVariance, QT_TRANSLATE_NOOP("PropertyName", "depthVariance"), QT_TRANSLATE_NOOP("Units", "J"),
0026             QT_TRANSLATE_NOOP("PropertyDescription", "Potential depth variance"), depthVariance, setDepthVariance)
0027     STEPCORE_PROPERTY_RW(double, rminVariance, QT_TRANSLATE_NOOP("PropertyName", "rminVariance"), QT_TRANSLATE_NOOP("Units", "m"),
0028             QT_TRANSLATE_NOOP("PropertyDescription", "Variance of the distance at which the force is zero"), rminVariance, setRminVariance))
0029
0030 // XXX: Check units for 2d
0031 // XXX: add cmPosition and cmVelocity
0032 STEPCORE_META_OBJECT(Gas, QT_TRANSLATE_NOOP("ObjectClass", "Gas"), QT_TRANSLATE_NOOP("ObjectDescription", "Particle gas"), 0, STEPCORE_SUPER_CLASS(ItemGroup),
0033     STEPCORE_PROPERTY_RW(StepCore::Vector2d, measureRectCenter, QT_TRANSLATE_NOOP("PropertyName", "measureRectCenter"), QT_TRANSLATE_NOOP("Units", "m"),
0034                 QT_TRANSLATE_NOOP("PropertyDescription", "Center of the rect for measurements"), measureRectCenter, setMeasureRectCenter)
0035     STEPCORE_PROPERTY_RW(StepCore::Vector2d, measureRectSize, QT_TRANSLATE_NOOP("PropertyName", "measureRectSize"), QT_TRANSLATE_NOOP("Units", "m"),
0036                 QT_TRANSLATE_NOOP("PropertyDescription", "Size of the rect for measurements"), measureRectSize, setMeasureRectSize)
0037     STEPCORE_PROPERTY_R_D(double, rectVolume, QT_TRANSLATE_NOOP("PropertyName", "rectVolume"), STEPCORE_FROM_UTF8(QT_TRANSLATE_NOOP("Units", "m²")),
0038                 QT_TRANSLATE_NOOP("PropertyDescription", "Volume of the measureRect"), rectVolume)
0039     STEPCORE_PROPERTY_R_D(double, rectParticleCount, QT_TRANSLATE_NOOP("PropertyName", "rectParticleCount"), STEPCORE_UNITS_1,
0040                 QT_TRANSLATE_NOOP("PropertyDescription", "Count of particles in the measureRect"), rectParticleCount)
0041     STEPCORE_PROPERTY_R_D(double, rectConcentration, QT_TRANSLATE_NOOP("PropertyName", "rectConcentration"), STEPCORE_FROM_UTF8(QT_TRANSLATE_NOOP("Units", "1/m²")),
0042                 QT_TRANSLATE_NOOP("PropertyDescription", "Concentration of particles in the measureRect"), rectConcentration)
0043     STEPCORE_PROPERTY_R_D(double, rectPressure, QT_TRANSLATE_NOOP("PropertyName", "rectPressure"), QT_TRANSLATE_NOOP("Units", "Pa"),
0044                 QT_TRANSLATE_NOOP("PropertyDescription", "Pressure of particles in the measureRect"), rectPressure)
0045     STEPCORE_PROPERTY_R_D(double, rectTemperature, QT_TRANSLATE_NOOP("PropertyName", "rectTemperature"), QT_TRANSLATE_NOOP("Units", "K"),
0046                 QT_TRANSLATE_NOOP("PropertyDescription", "Temperature of particles in the measureRect"), rectTemperature)
0047     STEPCORE_PROPERTY_R_D(double, rectMeanKineticEnergy, QT_TRANSLATE_NOOP("PropertyName", "rectMeanKineticEnergy"), QT_TRANSLATE_NOOP("Units", "J"),
0048                 QT_TRANSLATE_NOOP("PropertyDescription", "Mean kinetic energy of particles in the measureRect"), rectMeanKineticEnergy)
0049     STEPCORE_PROPERTY_R_D(StepCore::Vector2d, rectMeanVelocity, QT_TRANSLATE_NOOP("PropertyName", "rectMeanVelocity"), QT_TRANSLATE_NOOP("Units", "m/s"),
0050                 QT_TRANSLATE_NOOP("PropertyDescription", "Mean velocity of particles in the measureRect"), rectMeanVelocity)
0051     STEPCORE_PROPERTY_R_D(double, rectMeanParticleMass, QT_TRANSLATE_NOOP("PropertyName", "rectMeanParticleMass"), QT_TRANSLATE_NOOP("Units", "kg"),
0052                 QT_TRANSLATE_NOOP("PropertyDescription", "Mean mass of particles in the measureRect"), rectMeanParticleMass)
0053     STEPCORE_PROPERTY_R_D(double, rectMass, QT_TRANSLATE_NOOP("PropertyName", "rectMass"), QT_TRANSLATE_NOOP("Units", "kg"),
0054                 QT_TRANSLATE_NOOP("PropertyDescription", "Total mass of particles in the measureRect"), rectMass)
0055     )
0056
0057 STEPCORE_META_OBJECT(GasErrors, QT_TRANSLATE_NOOP("ObjectClass", "GasErrors"), QT_TRANSLATE_NOOP("ObjectDescription", "Errors class for Gas"), 0, STEPCORE_SUPER_CLASS(ObjectErrors),
0058     STEPCORE_PROPERTY_R_D(double, rectPressureVariance, QT_TRANSLATE_NOOP("PropertyName", "rectPressureVariance"), QT_TRANSLATE_NOOP("Units", "Pa"),
0059                 QT_TRANSLATE_NOOP("PropertyDescription", "Variance of pressure of particles in the measureRect"), rectPressureVariance)
0060     STEPCORE_PROPERTY_R_D(double, rectTemperatureVariance, QT_TRANSLATE_NOOP("PropertyName", "rectTemperatureVariance"), QT_TRANSLATE_NOOP("Units", "K"),
0061                 QT_TRANSLATE_NOOP("PropertyDescription", "Variance of temperature of particles in the measureRect"), rectTemperatureVariance)
0062     STEPCORE_PROPERTY_R_D(double, rectMeanKineticEnergyVariance, QT_TRANSLATE_NOOP("PropertyName", "rectMeanKineticEnergyVariance"), QT_TRANSLATE_NOOP("Units", "J"),
0063                 QT_TRANSLATE_NOOP("PropertyDescription", "Variance of mean kinetic energy of particles in the measureRect"), rectMeanKineticEnergyVariance)
0064     STEPCORE_PROPERTY_R_D(StepCore::Vector2d, rectMeanVelocityVariance, QT_TRANSLATE_NOOP("PropertyName", "rectMeanVelocityVariance"), QT_TRANSLATE_NOOP("Units", "m/s"),
0065                 QT_TRANSLATE_NOOP("PropertyDescription", "Variance of mean velocity of particles in the measureRect"), rectMeanVelocityVariance)
0066     STEPCORE_PROPERTY_R_D(double, rectMeanParticleMassVariance, QT_TRANSLATE_NOOP("PropertyName", "rectMeanParticleMassVariance"), QT_TRANSLATE_NOOP("Units", "kg"),
0067                 QT_TRANSLATE_NOOP("PropertyDescription", "Variance of mean mass of particles in the measureRect"), rectMeanParticleMassVariance)
0068     STEPCORE_PROPERTY_R_D(double, rectMassVariance, QT_TRANSLATE_NOOP("PropertyName", "rectMassVariance"), QT_TRANSLATE_NOOP("Units", "kg"),
0069                 QT_TRANSLATE_NOOP("PropertyDescription", "Variance of total mass of particles in the measureRect"), rectMassVariance)
0070     )
0071
0072 GasLJForce* GasLJForceErrors::gasLJForce() const
0073 {
0074     return static_cast<GasLJForce*>(owner());
0075 }
0076
0077 GasLJForce::GasLJForce(double depth, double rmin, double cutoff)
0078     : Force()
0079     , _depth(depth)
0080     , _rmin(rmin)
0081     , _cutoff(cutoff)
0082 {
0083     calcABC();
0084 }
0085
0086 void GasLJForce::calcABC()
0087 {
0088     double m = 12*_depth;
0089     _rmin6 = pow(_rmin, 6);
0090     _rmin12 = _rmin6*_rmin6;
0091     _a = m*_rmin12; _b = m*_rmin6;
0092     _c = _cutoff*_cutoff;
0093 }
0094
0095 void GasLJForce::calcForce(bool calcVariances)
0096 {
0097     if(!group()) return;
0098
0099     // NOTE: Currently we are handling only children of the same group
0100     const ItemList::const_iterator end = group()->items().end();
0101     for(ItemList::const_iterator i1 = group()->items().begin(); i1 != end; ++i1) {
0102         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0103         for(ItemList::const_iterator i2 = i1+1; i2 != end; ++i2) {
0104             if(!(*i2)->metaObject()->inherits<GasParticle>()) continue;
0105
0106             GasParticle* p1 = static_cast<GasParticle*>(*i1);
0107             GasParticle* p2 = static_cast<GasParticle*>(*i2);
0108             Vector2d r = p2->position() - p1->position();
0109             double rsquaredNorm = r.squaredNorm();
0110             if(rsquaredNorm < _c) {
0111                 double rnorm6 = rsquaredNorm*rsquaredNorm*rsquaredNorm;
0112                 double rnorm8 = rnorm6*rsquaredNorm;
0113                 Vector2d force = r * ((_a/rnorm6 - _b)/rnorm8);
0114                 p2->applyForce(force);
0115                 force = -force;
0116                 p1->applyForce(force);
0117
0118                 if(calcVariances) {
0119                     ParticleErrors* pe1 = p1->particleErrors();
0120                     ParticleErrors* pe2 = p2->particleErrors();
0121                     Vector2d rV = pe2->positionVariance() + pe1->positionVariance();
0122
0123                     GasLJForceErrors* ge = gasLJForceErrors();
0124                     Vector2d forceV = r.array().square().matrix() * (
0125                         ge->_rminVariance * square( (12*_a/_rmin/rnorm6 - 6*_b/_rmin)/rnorm8 ) +
0126                         ge->_depthVariance * square( 12*(_rmin12/rnorm6 - _rmin6)/rnorm8 ) );
0127
0128                     forceV[0] += rV[0] * square( (_a/rnorm6*( 1 - 14*r[0]*r[0]/rsquaredNorm ) -
0129                                                   _b*( 1 - 8*r[0]*r[0]/rsquaredNorm ))/rnorm8 ) +
0130                                  rV[1] * square( (_a/rnorm6*14 - _b*8)*r[0]*r[1]/(rnorm8*rsquaredNorm) );
0131                     forceV[1] += rV[1] * square( (_a/rnorm6*( 1 - 14*r[1]*r[1]/rsquaredNorm ) -
0132                                                   _b*( 1 - 8*r[1]*r[1]/rsquaredNorm ))/rnorm8 ) +
0133                                  rV[0] * square( (_a/rnorm6*14 - _b*8)*r[0]*r[1]/(rnorm8*rsquaredNorm) );
0134
0135                     pe1->applyForceVariance(forceV);
0136                     pe2->applyForceVariance(forceV);
0137                 }
0138             }
0139         }
0140     }
0141 }
0142
0143 Gas* GasErrors::gas() const
0144 {
0145     return static_cast<Gas*>(owner());
0146 }
0147
0148 double Gas::rectVolume() const
0149 {
0150     return _measureRectSize[0]*_measureRectSize[1];
0151 }
0152
0153 double Gas::randomUniform(double min, double max)
0154 {
0155     return double(std::rand())/RAND_MAX*(max-min) + min;
0156 }
0157
0158 double Gas::randomGauss(double mean, double deviation)
0159 {
0160     // Polar Box-Muller algorithm
0161     double x1, x2, w;
0162
0163     do {
0164         x1 = 2.0 * double(std::rand())/RAND_MAX - 1.0;
0165         x2 = 2.0 * double(std::rand())/RAND_MAX - 1.0;
0166         w = x1 * x1 + x2 * x2;
0167     } while( w >= 1.0 || w == 0 );
0168
0169     w = sqrt( (-2.0 * log( w ) ) / w );
0170     return x1 * w * deviation + mean;
0171 }
0172
0173 GasParticleList Gas::rectCreateParticles(int count,
0174                                 double mass, double temperature,
0175                                 const Vector2d& meanVelocity)
0176 {
0177     GasParticleList particles;
0178
0179     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0180     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0181     double deviation = std::sqrt( Constants::Boltzmann * temperature / mass );
0182
0183     for(int i=0; i<count; ++i) {
0184         Vector2d pos(randomUniform(r0[0], r1[0]), randomUniform(r0[1], r1[1]));
0185         Vector2d vel(randomGauss(meanVelocity[0], deviation), randomGauss(meanVelocity[1], deviation));
0186         GasParticle* particle = new GasParticle(pos, vel, mass);
0187         particles.push_back(particle);
0188     }
0189
0190     return particles;
0191 }
0192
0193 void Gas::addParticles(const GasParticleList& particles)
0194 {
0195     const GasParticleList::const_iterator end = particles.end();
0196     for(GasParticleList::const_iterator it = particles.begin(); it != end; ++it) {
0197         addItem(*it);
0198     }
0199 }
0200
0201 double Gas::rectParticleCount() const
0202 {
0203     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0204     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0205
0206     double count = 0;
0207     const ItemList::const_iterator end = items().end();
0208     for(ItemList::const_iterator i1 = items().begin(); i1 != end; ++i1) {
0209         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0210         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0211         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0212             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0213         ++count;
0214     }
0215
0216     return count;
0217 }
0218
0219 double Gas::rectMeanParticleMass() const
0220 {
0221     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0222     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0223
0224     double count = 0;
0225     double mass = 0;
0226
0227     const ItemList::const_iterator end = items().end();
0228     for(ItemList::const_iterator i1 = items().begin(); i1 != end; ++i1) {
0229         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0230         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0231         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0232             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0233         mass += p1->mass();
0234         ++count;
0235     }
0236
0237     mass /= count;
0238     return mass;
0239 }
0240
0241 double GasErrors::rectMeanParticleMassVariance() const
0242 {
0243     Vector2d r0 = gas()->_measureRectCenter-gas()->_measureRectSize/2.0;
0244     Vector2d r1 = gas()->_measureRectCenter+gas()->_measureRectSize/2.0;
0245
0246     double count = 0;
0247
0248     double mass = gas()->rectMeanParticleMass();
0249     double massVariance = 0;
0250
0251     const ItemList::const_iterator end = gas()->items().end();
0252     for(ItemList::const_iterator i1 = gas()->items().begin(); i1 != end; ++i1) {
0253         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0254         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0255         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0256             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0257
0258         massVariance += square(p1->mass() - mass);
0259
0260         ParticleErrors* pe1 = static_cast<ParticleErrors*>(p1->tryGetObjectErrors());
0261         if(pe1) massVariance += pe1->massVariance();
0262
0263         ++count;
0264     }
0265
0266     massVariance /= square(count);
0267     return massVariance;
0268 }
0269
0270 double Gas::rectMass() const
0271 {
0272     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0273     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0274
0275     double mass = 0;
0276
0277     const ItemList::const_iterator end = items().end();
0278     for(ItemList::const_iterator i1 = items().begin(); i1 != end; ++i1) {
0279         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0280         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0281         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0282             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0283         mass += p1->mass();
0284     }
0285
0286     return mass;
0287 }
0288
0289 double GasErrors::rectMassVariance() const
0290 {
0291     Vector2d r0 = gas()->_measureRectCenter-gas()->_measureRectSize/2.0;
0292     Vector2d r1 = gas()->_measureRectCenter+gas()->_measureRectSize/2.0;
0293
0294     double massVariance = 0;
0295
0296     const ItemList::const_iterator end = gas()->items().end();
0297     for(ItemList::const_iterator i1 = gas()->items().begin(); i1 != end; ++i1) {
0298         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0299         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0300         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0301             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0302
0303         ParticleErrors* pe1 = static_cast<ParticleErrors*>(p1->tryGetObjectErrors());
0304         if(pe1) massVariance += pe1->massVariance();
0305     }
0306
0307     return massVariance;
0308 }
0309
0310 double Gas::rectConcentration() const
0311 {
0312     return rectParticleCount() / rectVolume();
0313 }
0314
0315 Vector2d Gas::rectMeanVelocity() const
0316 {
0317     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0318     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0319
0320     double count = 0;
0321     Vector2d velocity(0.,0.);
0322
0323     const ItemList::const_iterator end = items().end();
0324     for(ItemList::const_iterator i1 = items().begin(); i1 != end; ++i1) {
0325         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0326         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0327         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0328             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0329         velocity += p1->velocity();
0330         ++count;
0331     }
0332
0333     velocity /= count;
0334     return velocity;
0335 }
0336
0337 Vector2d GasErrors::rectMeanVelocityVariance() const
0338 {
0339     Vector2d r0 = gas()->_measureRectCenter-gas()->_measureRectSize/2.0;
0340     Vector2d r1 = gas()->_measureRectCenter+gas()->_measureRectSize/2.0;
0341
0342     double count = 0;
0343
0344     Vector2d velocity = gas()->rectMeanVelocity();
0345     Vector2d velocityVariance(0.,0.);
0346
0347     const ItemList::const_iterator end = gas()->items().end();
0348     for(ItemList::const_iterator i1 = gas()->items().begin(); i1 != end; ++i1) {
0349         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0350         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0351         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0352             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0353
0354         velocityVariance += (p1->velocity() - velocity).array().square().matrix();
0355
0356         ParticleErrors* pe1 = static_cast<ParticleErrors*>(p1->tryGetObjectErrors());
0357         if(pe1) velocityVariance += pe1->velocityVariance();
0358
0359         ++count;
0360     }
0361
0362     velocityVariance /= square(count);
0363     return velocityVariance;
0364 }
0365
0366 double Gas::rectMeanKineticEnergy() const
0367 {
0368     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0369     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0370
0371     double count = 0;
0372     double energy = 0;
0373
0374     const ItemList::const_iterator end = items().end();
0375     for(ItemList::const_iterator i1 = items().begin(); i1 != end; ++i1) {
0376         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0377         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0378         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0379             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0380         energy += p1->mass() * p1->velocity().squaredNorm();
0381         ++count;
0382     }
0383
0384     energy /= (2.0*count);
0385     return energy;
0386 }
0387
0388 double GasErrors::rectMeanKineticEnergyVariance() const
0389 {
0390     Vector2d r0 = gas()->_measureRectCenter-gas()->_measureRectSize/2.0;
0391     Vector2d r1 = gas()->_measureRectCenter+gas()->_measureRectSize/2.0;
0392
0393     double count = 0;
0394
0395     double energy = 2*gas()->rectMeanKineticEnergy();
0396     double energyVariance = 0;
0397
0398     const ItemList::const_iterator end = gas()->items().end();
0399     for(ItemList::const_iterator i1 = gas()->items().begin(); i1 != end; ++i1) {
0400         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0401         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0402         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0403             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0404
0405         double pEnergy = p1->mass() * p1->velocity().squaredNorm();
0406         energyVariance += square(pEnergy - energy);
0407
0408         ParticleErrors* pe1 = static_cast<ParticleErrors*>(p1->tryGetObjectErrors());
0409         if(pe1) {
0410             energyVariance +=
0411                 pe1->massVariance() * square(p1->velocity().squaredNorm()) +
0412                 ((2*p1->mass()*p1->velocity()).array().square().matrix()).dot(pe1->velocityVariance());
0413         }
0414
0415         ++count;
0416     }
0417
0418     energyVariance /= square(2.0*count);
0419     return energyVariance;
0420 }
0421
0422 double Gas::rectTemperature() const
0423 {
0424     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0425     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0426
0427     double count = 0;
0428     double temperature = 0;
0429
0430     StepCore::Vector2d meanVelocity = rectMeanVelocity();
0431     const ItemList::const_iterator end = items().end();
0432     for(ItemList::const_iterator i1 = items().begin(); i1 != end; ++i1) {
0433         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0434         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0435         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0436             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0437         temperature += p1->mass() * (p1->velocity() - meanVelocity).squaredNorm();
0438         ++count;
0439     }
0440
0441     // no 3/2 factor since we live in 2d
0442     temperature /= (2.0*count*Constants::Boltzmann);
0443     return temperature;
0444 }
0445
0446 double GasErrors::rectTemperatureVariance() const
0447 {
0448     Vector2d r0 = gas()->_measureRectCenter-gas()->_measureRectSize/2.0;
0449     Vector2d r1 = gas()->_measureRectCenter+gas()->_measureRectSize/2.0;
0450
0451     double count = 0;
0452
0453     StepCore::Vector2d meanVelocity = gas()->rectMeanVelocity();
0454     double temperature = 2.0*Constants::Boltzmann*gas()->rectTemperature();
0455     double temperatureVariance = 0;
0456
0457     const ItemList::const_iterator end = gas()->items().end();
0458     for(ItemList::const_iterator i1 = gas()->items().begin(); i1 != end; ++i1) {
0459         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0460         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0461         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0462             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0463
0464         double pTemperature = p1->mass() * (p1->velocity() - meanVelocity).squaredNorm();
0465         temperatureVariance += square(pTemperature - temperature);
0466
0467         ParticleErrors* pe1 = static_cast<ParticleErrors*>(p1->tryGetObjectErrors());
0468         if(pe1) {
0469             temperatureVariance +=
0470                 pe1->massVariance() * square((p1->velocity() - meanVelocity).squaredNorm()) +
0471                 ((p1->mass()*(p1->velocity() - meanVelocity)).array().square().matrix()).dot(pe1->velocityVariance());
0472         }
0473
0474         ++count;
0475     }
0476
0477     temperatureVariance /= square(2.0*Constants::Boltzmann*count);
0478     // XXX: We could easily take into account BoltzmannError here
0479     // but this can confuse users so for now we don't do it
0480
0481     return temperatureVariance;
0482 }
0483
0484 // XXX: this formula is incorrect when forces are big
0485 // XXX: use better formula (for example from lammps)
0486 double Gas::rectPressure() const
0487 {
0488     Vector2d r0 = _measureRectCenter-_measureRectSize/2.0;
0489     Vector2d r1 = _measureRectCenter+_measureRectSize/2.0;
0490
0491     double pressure = 0;
0492
0493     StepCore::Vector2d meanVelocity = rectMeanVelocity();
0494     const ItemList::const_iterator end = items().end();
0495     for(ItemList::const_iterator i1 = items().begin(); i1 != end; ++i1) {
0496         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0497         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0498         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0499             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0500         pressure += p1->mass() * (p1->velocity() - meanVelocity).squaredNorm();
0501     }
0502
0503     pressure /= (2.0 * rectVolume());
0504     return pressure;
0505 }
0506
0507 double GasErrors::rectPressureVariance() const
0508 {
0509     Vector2d r0 = gas()->_measureRectCenter-gas()->_measureRectSize/2.0;
0510     Vector2d r1 = gas()->_measureRectCenter+gas()->_measureRectSize/2.0;
0511
0512     StepCore::Vector2d meanVelocity = gas()->rectMeanVelocity();
0513     double pressure = gas()->rectPressure();
0514     double pressureVariance = 0;
0515
0516     const ItemList::const_iterator end = gas()->items().end();
0517     for(ItemList::const_iterator i1 = gas()->items().begin(); i1 != end; ++i1) {
0518         if(!(*i1)->metaObject()->inherits<GasParticle>()) continue;
0519         GasParticle* p1 = static_cast<GasParticle*>(*i1);
0520         if(p1->position()[0] < r0[0] || p1->position()[0] > r1[0] ||
0521             p1->position()[1] < r0[1] || p1->position()[1] > r1[1]) continue;
0522
0523         double pPressure = p1->mass() * (p1->velocity() - meanVelocity).squaredNorm();
0524         pressureVariance += square(pPressure - pressure);
0525
0526         ParticleErrors* pe1 = static_cast<ParticleErrors*>(p1->tryGetObjectErrors());
0527         if(pe1) {
0528             pressureVariance +=
0529                 pe1->massVariance() * square((p1->velocity() - meanVelocity).squaredNorm()) +
0530                 ((p1->mass()*(p1->velocity() - meanVelocity)).array().square().matrix()).dot(pe1->velocityVariance());
0531         }
0532     }
0533
0534     pressureVariance /= square(2.0*gas()->rectVolume());
0535     return pressureVariance;
0536 }
0537
0538 }
0539