ErrorEstimation.h

#ifndef _Error_ESTIMATION_H_
#define _Error_ESTIMATION_H_
 
namespace coupling {
namespace error {
 
class ErrorEstimation;
 
}
} // namespace coupling

class coupling::error::ErrorEstimation {
public:
  ErrorEstimation(double velocity, double temperature, double numberOfParticle, double particleMass, double soundSpeed, double numberOfSamples,
                  double cellVolume)
      : _velocity(velocity), _temperature(temperature), _numberOfParticle(numberOfParticle), _particleMass(particleMass), _soundSpeed(soundSpeed),
        _numberOfSamples(numberOfSamples), _cellVolume(cellVolume), _desiredAbsErrorVelocity(0.08), _desiredRelErrorVelocity(0.1),
        _desiredAbsErrorDensity(0.05), _desiredRelErrorDensity(0.05), _desiredAbsErrorTemperature(0.05), _desiredRelErrorTemperature(0.05),
        _desiredAbsErrorPressure(0.05), _desiredRelErrorPressure(0.05), _gamma(1.667), _C_v(1.4973), _k(1) {}

  ~ErrorEstimation() {
#ifdef DEBUG_Error
    std::cout << " Error estimation deconstructed" << std::endl;
#endif
  }

  enum errorBaseQuantity {
    Velocity = 0 ,
    Density = 1 ,
    Temperature = 2 ,
    Pressure = 3 
  };

  enum errorType {
    Relative = 0 ,
    Absolute = 1 
  };

  double getError(errorBaseQuantity BaseQuantity, errorType Type) {
 
    double error;
 
    if (BaseQuantity == Velocity) {
 
      error = getErrorVelocity(_numberOfSamples, Type == Relative ? _velocity : 1, _temperature, _numberOfParticle, _particleMass);
 
    } else if (BaseQuantity == Density) {
 
      error = getErrorDensity(_numberOfSamples, _soundSpeed, _temperature, Type == Relative ? _numberOfParticle : 1 / _numberOfParticle, _particleMass);
 
    } else if (BaseQuantity == Temperature) {
 
      error = getErrorTemperature(_numberOfSamples, _numberOfParticle, Type == Relative ? 1 : _temperature);
 
    } else if (BaseQuantity == Pressure) {
 
      error = getErrorPressure(_numberOfSamples, _numberOfParticle, _temperature, _soundSpeed, _particleMass, _cellVolume, 1);
 
    } else {
 
      std::cout << "ERROR coupling::ErrorEstimation(): Base quantity invalid! " << std::endl;
      exit(EXIT_FAILURE);
    }
 
    return error;
  }

  double getCorrectorNumberOfSamples(errorBaseQuantity BaseQuantity, errorType Type) {
 
    if (BaseQuantity == Velocity) {
 
      return requiredSamplesV(Type == Relative ? _desiredRelErrorVelocity : _desiredAbsErrorVelocity, _temperature, _soundSpeed,
                              Type == Relative ? _velocity : 1, _numberOfParticle, _particleMass);
 
    } else if (BaseQuantity == Density) {
 
      return requiredSamplesD(_desiredRelErrorDensity, _temperature, _soundSpeed, _velocity, _particleMass, _numberOfParticle);
 
    } else if (BaseQuantity == Temperature) {
 
      return 3.0;
 
    } else if (BaseQuantity == Pressure) {
 
      return 4.0;
    }
 
    std::cout << "Error Estimation::getCorrectorNumberOfSamples not a Valid quantity" << std::endl;
    std::exit(EXIT_FAILURE);
    return 0.0;
  }

  double getErrorVelocity(double numberOfSamples, double velocity, double temperature, double numberOfParticle, double particleMass) {
    double error = 1 / (velocitySNR(velocity, temperature, numberOfParticle, particleMass) * std::sqrt(numberOfSamples));
    return error;
  }

  double getErrorDensity(double numberOfSamples, double soundSpeed, double temperature, double numberOfParticle, double particleMass) {
    double refSP = referenceSoundSpeed(temperature, particleMass);
    double Ac = acousticNumber(soundSpeed, refSP);
    double error = 1 / std::sqrt(numberOfParticle * numberOfSamples) / Ac;
    //  std::cout << "numberOfParticle "  << numberOfParticle <<
    //"numberOfSamples  " << numberOfSamples << "Ac  "  << Ac <<  std::endl;
    return error;
  }

  double getErrorTemperature(double numberOfSamples, double numberOfParticle, double temperature) {
 
    //      double deltaT = k*temperature*temperature/_C_v/numberOfParticle;
 
    double error = std::sqrt(_k / (_C_v * numberOfParticle * numberOfSamples)) * temperature;
    return error;
  }

  double getErrorPressure(double numberOfSamples, double numberOfParticle, double temperature, double soundSpeed, double particleMass, double cellVolume,
                          double pressure) {
 
    double Ac = acousticNumber(soundSpeed, referenceSoundSpeed(temperature, particleMass));
    double referenceP = referencePressure(temperature, numberOfParticle, cellVolume);
 
    double error = referenceP / pressure * Ac * std::sqrt(_gamma / (numberOfParticle * numberOfSamples));
    return error;
  }

  double requiredSamplesV(double desiredError, double temperature, double soundSpeed, double velocity, double numberOfParticle, double particleMass) {
    //  double refSP= referenceSoundSpeed( temperature, particleMass);
    //  double Ac = acousticNumber(soundSpeed, refSP);
    double SNR = velocitySNR(velocity, temperature, numberOfParticle, particleMass);
    double desiredNumber = 1 / ((SNR * desiredError) * (SNR * desiredError));
 
    return desiredNumber;
  }

  double requiredSamplesD(double desiredError, double temperature, double soundSpeed, double velocity, double particleMass, double numberOfParticle) {
    double refSP = referenceSoundSpeed(temperature, particleMass);
    double Ac = acousticNumber(soundSpeed, refSP);
    double desiredNumber = 1 / numberOfParticle / (desiredError * desiredError) / (Ac * Ac);
 
    return desiredNumber;
  }

  double reqiredSamplesT(double desiredError) {
 
    double desiredNumber = (_k / (_C_v * _numberOfParticle * desiredError)) * (_k / (_C_v * _numberOfParticle * desiredError));
 
    return desiredNumber;
  }

  double velocitySNR(double velocity, double temperature, double numberOfParticle, double particleMass) {
 
    return velocity / veloyityFluctuation(temperature, numberOfParticle, particleMass);
  }

  double veloyityFluctuation(double temperature, double numberOfParticle, double particleMass) {
 
    return std::sqrt(_k * temperature / particleMass / numberOfParticle);
  }

  void setAbsVelocityError(double error) { _desiredAbsErrorVelocity = error; }

  void setRelVelocityError(double error) { _desiredRelErrorVelocity = error; }

  void setAbsDensityError(double error) { _desiredAbsErrorDensity = error; }

  void setRelDensityError(double error) { _desiredRelErrorDensity = error; }

  void setAbsTemperatureError(double error) { _desiredAbsErrorTemperature = error; }

  void setRelTemperatureError(double error) { _desiredRelErrorTemperature = error; }

  void setAbsPressureError(double error) { _desiredAbsErrorPressure = error; }

  void setRelPressureError(double error) { _desiredRelErrorPressure = error; }
 
private:
  double getErrorVelocity() {
    double er = getErrorVelocity(_numberOfSamples, _velocity, _temperature, _numberOfParticle, _particleMass);
    return er;
  }

  double getErrorDensity() {
    double er = getErrorDensity(_numberOfSamples, _soundSpeed, _temperature, _numberOfParticle, _particleMass);
    return er;
  }

  double getErrorTemperature() {
    double er = getErrorTemperature(_numberOfSamples, _numberOfParticle, _temperature);
    return er;
  }
 
  double acousticNumber(double soundSpeed, double soundSpeed_ref) { return (soundSpeed / soundSpeed_ref); }
 
  //------------------------- Reference Properties
  //---------------------------------------------------
  double referenceSoundSpeed(double temperature, double particleMass) { return std::sqrt(_gamma * _k * temperature / particleMass); }

  double referencePressure(double temperature, double numberOfParticle, double cellVolume) {
    return std::sqrt(numberOfParticle * _k * temperature / cellVolume);
  }
 
  double _velocity;
  double _temperature;
  double _numberOfParticle;
  double _particleMass;
  double _soundSpeed;
  double _numberOfSamples;
  double _cellVolume;
 
  double _desiredAbsErrorVelocity; //=0.08;
  double _desiredRelErrorVelocity; // =0.1;
 
  double _desiredAbsErrorDensity; //=0.05;
  double _desiredRelErrorDensity; //=0.05;
 
  double _desiredAbsErrorTemperature; //=0.05;
  double _desiredRelErrorTemperature; //=0.05;
 
  double _desiredAbsErrorPressure; //=0.05;
  double _desiredRelErrorPressure; //=0.05;
 
  double _gamma;
  double _C_v;
  double _k;
};
 
#endif