unitConverter.h

Go to the documentation of this file.

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2017 Max Gaedtke, Albert Mink
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/
 
#ifndef UNITCONVERTER_H
#define UNITCONVERTER_H
 
 
#include "utilities/omath.h"
#include "io/ostreamManager.h"
#include "core/util.h"
#include "io/xmlReader.h"
 
#include "descriptor/fields.h"
 
// known design issues
//    1. How can we prevent abuse of constructur by mixing up parameters?
//    2. physical problems may have different names for viscosity, e.g. diffusity,  temperature conductivity
//    4. Feedback about stability or comment the chosen discretization
//    5. Explain why Desctiptor as template
//    6. Is it worth to introduce invConversionDensity to avoid division
 
 
namespace olb {
 
namespace fields {
 
namespace converter {
 
struct PHYS_VELOCITY : public descriptors::FIELD_BASE<1> {
  template <typename T, typename DESCRIPTOR,typename FIELD>
  static constexpr auto isValid(FieldD<T,DESCRIPTOR,FIELD> value) {
    return value >= 0;
  }
};
struct PHYS_FORCE : public descriptors::FIELD_BASE<1> {
  template <typename T, typename DESCRIPTOR,typename FIELD>
  static constexpr auto isValid(FieldD<T,DESCRIPTOR,FIELD> value) {
    return value >= 0;
  }
};
struct PHYS_LENGTH : public descriptors::FIELD_BASE<1> {
  template <typename T, typename DESCRIPTOR,typename FIELD>
  static constexpr auto isValid(FieldD<T,DESCRIPTOR,FIELD> value) {
    return value > 0;
  }
 };
struct PHYS_DELTA_X : public descriptors::FIELD_BASE<1> {
  template <typename T, typename DESCRIPTOR,typename FIELD>
  static constexpr auto isValid(FieldD<T,DESCRIPTOR,FIELD> value) {
    return value > 0;
  }
};
struct PHYS_DELTA_T : public descriptors::FIELD_BASE<1> {
  template <typename T, typename DESCRIPTOR,typename FIELD>
  static constexpr auto isValid(FieldD<T,DESCRIPTOR,FIELD> value) {
    return value > 0;
  }
};
struct PHYS_PRESSURE : public descriptors::FIELD_BASE<1> {
  template <typename T, typename DESCRIPTOR,typename FIELD>
  static constexpr auto isValid(FieldD<T,DESCRIPTOR,FIELD> value) {
    return value >= 0;
  }
};
 
struct LATTICE_TIME : public descriptors::FIELD_BASE<1> {
  template <typename T, typename DESCRIPTOR,typename FIELD>
  static constexpr auto isValid(FieldD<T,DESCRIPTOR,FIELD> value) {
    return value >= 0;
  }
};
 
struct LATTICE_VISCOSITY : public descriptors::FIELD_BASE<1> { };
 
}
 
}
 
 
struct UnitConverterBase {
  virtual ~UnitConverterBase() = default;
 
  virtual void print() const = 0;
  virtual void print(std::ostream& fout) const = 0;
  virtual void write(std::string const& fileName = "unitConverter") const = 0;
};
 
 
template <typename T, typename DESCRIPTOR>
class UnitConverter : public UnitConverterBase {
public:
  constexpr UnitConverter( T physDeltaX, T physDeltaT, T charPhysLength, T charPhysVelocity,
                           T physViscosity, T physDensity, T charPhysPressure = 0 )
    : _conversionLength(physDeltaX),
      _conversionTime(physDeltaT),
      _conversionVelocity(_conversionLength / _conversionTime),
      _conversionDensity(physDensity),
      _conversionMass( _conversionDensity * util::pow(_conversionLength, 3) ),
      _conversionViscosity(_conversionLength * _conversionLength / _conversionTime),
      _conversionForce( _conversionMass * _conversionLength / (_conversionTime * _conversionTime) ),
      _conversionTorque( _conversionMass * _conversionLength * _conversionLength / (_conversionTime * _conversionTime) ),
      _conversionPressure( _conversionForce / util::pow(_conversionLength, 2) ),
      _charPhysLength(charPhysLength),
      _charPhysVelocity(charPhysVelocity),
      _physViscosity(physViscosity),
      _physDensity(physDensity),
      _charPhysPressure(charPhysPressure),
      _resolution((size_t)(_charPhysLength / _conversionLength + 0.5)),
      _latticeRelaxationTime( (_physViscosity / _conversionViscosity * descriptors::invCs2<T,DESCRIPTOR>()) + 0.5 ),
      _charLatticeVelocity( _charPhysVelocity / _conversionVelocity ),
      clout(std::cout,"UnitConverter")
  {
  }
 
  virtual ~UnitConverter() = default;
 
  constexpr int getResolution(  ) const
  {
    return _resolution;
  }
  constexpr T getLatticeRelaxationTime(  ) const
  {
    return _latticeRelaxationTime;
  }
  constexpr T getLatticeRelaxationFrequency(  ) const
  {
    return 1./_latticeRelaxationTime;
  }
  template <typename DESCRIPTOR_>
  constexpr T getLatticeRelaxationFrequencyFromDiffusivity( const T physDiffusivity ) const
  {
    return 1.0 / ( physDiffusivity / _conversionViscosity * descriptors::invCs2<T,DESCRIPTOR_>() + 0.5 );
  }
  constexpr T getCharPhysLength(  ) const
  {
    return _charPhysLength;
  }
  constexpr T getCharPhysVelocity(  ) const
  {
    return _charPhysVelocity;
  }
  constexpr T getCharLatticeVelocity(  ) const
  {
    return _charLatticeVelocity;
  }
  constexpr T getCharCFLnumber(  ) const
  {
    return _charLatticeVelocity;
  }
  constexpr T getPhysViscosity(  ) const
  {
    return _physViscosity;
  }
  constexpr T getPhysDensity(  ) const
  {
    return _physDensity;
  }
  constexpr T getCharPhysPressure(  ) const
  {
    return _charPhysPressure;
  }
  constexpr T getReynoldsNumber(  ) const
  {
    return _charPhysVelocity * _charPhysLength / _physViscosity;
  }
  constexpr T getMachNumber(  ) const
  {
    return getCharLatticeVelocity() * util::sqrt(descriptors::invCs2<T,DESCRIPTOR>());
  }
  constexpr T getKnudsenNumber(  ) const
  {
    // This calculates the lattice Knudsen number.
    // See e.g. (7.22) in "The Lattice Boltzmann Method: Principles and Practice" [kruger2017lattice].
    return getMachNumber() / getReynoldsNumber();
  }
  constexpr T getPhysLength( int latticeLength ) const
  {
    return _conversionLength * latticeLength;
  }
  constexpr int getLatticeLength( T physLength ) const
  {
    return int( physLength / _conversionLength + 0.5 );
  }
  constexpr T getConversionFactorLength() const
  {
    return _conversionLength;
  }
  constexpr T getPhysDeltaX() const
  {
    return _conversionLength;
  }
 
  constexpr T getPhysTime( size_t latticeTime ) const
  {
    return _conversionTime * latticeTime;
  }
  constexpr size_t getLatticeTime( T physTime ) const
  {
    return size_t(physTime / _conversionTime + 0.5);
  }
  constexpr T getConversionFactorTime() const
  {
    return _conversionTime;
  }
  constexpr T getPhysDeltaT() const
  {
    return _conversionTime;
  }
 
  constexpr T getPhysVelocity( T latticeVelocity ) const
  {
    return _conversionVelocity * latticeVelocity;
  }
  constexpr T getLatticeVelocity( T physVelocity ) const
  {
    return physVelocity / _conversionVelocity;
  }
  template <unsigned D>
  constexpr Vector<T,D> getLatticeVelocity(Vector<T,D> physU) const
  {
    return Vector<T,D>([&](std::size_t iD) -> T {
      return this->getLatticeVelocity(physU[iD]);
    });
  }
  constexpr T getConversionFactorVelocity() const
  {
    return _conversionVelocity;
  }
 
  constexpr T getPhysDensity( T latticeDensity ) const
  {
    return _conversionDensity * latticeDensity;
  }
  constexpr T getLatticeDensity( T physDensity ) const
  {
    return physDensity / _conversionDensity;
  }
  constexpr T getLatticeDensityFromPhysPressure( T physPressure ) const
  {
    return util::densityFromPressure<T,DESCRIPTOR>( getLatticePressure( physPressure ) );
  }
  constexpr T getConversionFactorDensity() const
  {
    return _conversionDensity;
  }
 
  constexpr T getPhysMass( T latticeMass ) const
  {
    return _conversionMass * latticeMass;
  }
  constexpr T getLatticeMass( T physMass ) const
  {
    return physMass / _conversionMass;
  }
  constexpr T getConversionFactorMass() const
  {
    return _conversionMass;
  }
 
  constexpr T getPhysViscosity( T latticeViscosity ) const
  {
    return _conversionViscosity * latticeViscosity;
  }
  constexpr T getLatticeViscosity(  ) const
  {
    return _physViscosity / _conversionViscosity;
  }
  constexpr T getConversionFactorViscosity() const
  {
    return _conversionViscosity;
  }
 
  constexpr T getPhysForce( T latticeForce ) const
  {
    return _conversionForce * latticeForce;
  }
  template <unsigned D>
  constexpr Vector<T,D> getPhysForce(Vector<T,D> latticeForce) const
  {
    return Vector<T,D>([&](std::size_t iD) -> T {
      return this->getPhysForce(latticeForce[iD]);
    });
  }
  constexpr T getLatticeForce( T physForce ) const
  {
    return physForce / _conversionForce;
  }
  constexpr T getConversionFactorForce() const
  {
    return _conversionForce;
  }
 
  constexpr T getPhysTorque ( T latticeTorque ) const
  {
    return _conversionTorque * latticeTorque;
  }
  template <unsigned D>
  constexpr Vector<T,D> getPhysTorque(Vector<T,D> latticeTorque) const
  {
    return Vector<T,D>([&](std::size_t iD) -> T {
      return this->getPhysTorque(latticeTorque[iD]);
    });
  }
  constexpr T getLatticeTorque( T physTorque ) const
  {
    return physTorque / _conversionTorque;
  }
  constexpr T getConversionFactorTorque() const
  {
    return _conversionTorque;
  }
 
  constexpr T getPhysPressure( T latticePressure ) const
  {
    return _conversionPressure * latticePressure + _charPhysPressure;
  }
  constexpr T getLatticePressure( T physPressure ) const
  {
    return ( physPressure - _charPhysPressure ) / _conversionPressure;
  }
  constexpr T getConversionFactorPressure() const
  {
    return _conversionPressure;
  }
  virtual void print() const;
  void print(std::ostream& fout) const;
 
  void write(std::string const& fileName = "unitConverter") const;
 
protected:
  // conversion factors
  const T _conversionLength;      // m
  const T _conversionTime;        // s
  const T _conversionVelocity;    // m / s
  const T _conversionDensity;     // kg / m^3
  const T _conversionMass;        // kg
  const T _conversionViscosity;   // m^2 / s
  const T _conversionForce;       // kg m / s^2
  const T _conversionTorque;      // kg m^2 / s^2
  const T _conversionPressure;    // kg / m s^2
 
  // physical units, e.g characteristic or reference values
  const T _charPhysLength;        // m
  const T _charPhysVelocity;      // m / s
  const T _physViscosity;         // m^2 / s
  const T _physDensity;           // kg / m^3
  const T _charPhysPressure;      // kg / m s^2
 
  // lattice units, discretization parameters
  const size_t _resolution;
  const T _latticeRelaxationTime;
  const T _charLatticeVelocity;   //
 
private:
  mutable OstreamManager clout;
};
 
template <typename T, typename DESCRIPTOR>
UnitConverter<T,DESCRIPTOR> convectivelyRefineUnitConverter(
  const UnitConverter<T,DESCRIPTOR>& converter,
  unsigned scale = 2)
{
  const T refinementFactor = T{1} / scale;
  return UnitConverter<T,DESCRIPTOR>(
    converter.getPhysDeltaX() * refinementFactor,
    converter.getPhysDeltaT() * refinementFactor,
    converter.getCharPhysLength(),
    converter.getCharPhysVelocity(),
    converter.getPhysViscosity(),
    converter.getPhysDensity(),
    converter.getCharPhysPressure()
  );
}
 
template <typename T, typename DESCRIPTOR>
UnitConverter<T, DESCRIPTOR>* createUnitConverter(XMLreader const& params);
 
template <typename T, typename DESCRIPTOR>
class UnitConverterFromResolutionAndRelaxationTime : public UnitConverter<T, DESCRIPTOR> {
public:
  constexpr UnitConverterFromResolutionAndRelaxationTime(
    size_t resolution,
    T latticeRelaxationTime,
    T charPhysLength,
    T charPhysVelocity,
    T physViscosity,
    T physDensity,
    T charPhysPressure = 0) : UnitConverter<T, DESCRIPTOR>(
        (charPhysLength/resolution),
        (latticeRelaxationTime - 0.5) / descriptors::invCs2<T,DESCRIPTOR>() * util::pow((charPhysLength/resolution),2) / physViscosity,
        charPhysLength,
        charPhysVelocity,
        physViscosity,
        physDensity,
        charPhysPressure)
  { }
};
 
template <typename T, typename DESCRIPTOR>
class UnitConverterFromResolutionAndLatticeVelocity : public UnitConverter<T, DESCRIPTOR> {
public:
  constexpr UnitConverterFromResolutionAndLatticeVelocity(
    size_t resolution,
    T charLatticeVelocity,
    T charPhysLength,
    T charPhysVelocity,
    T physViscosity,
    T physDensity,
    T charPhysPressure = 0) : UnitConverter<T, DESCRIPTOR>(
        (charPhysLength/resolution),
        (charLatticeVelocity / charPhysVelocity * charPhysLength / resolution),
        charPhysLength,
        charPhysVelocity,
        physViscosity,
        physDensity,
        charPhysPressure)
  { }
};
 
template <typename T, typename DESCRIPTOR>
class UnitConverterFromRelaxationTimeAndLatticeVelocity : public UnitConverter<T, DESCRIPTOR> {
public:
  constexpr UnitConverterFromRelaxationTimeAndLatticeVelocity(
    T latticeRelaxationTime,
    T charLatticeVelocity,
    T charPhysLength,
    T charPhysVelocity,
    T physViscosity,
    T physDensity,
    T charPhysPressure = 0) : UnitConverter<T, DESCRIPTOR>(
        (physViscosity * charLatticeVelocity / charPhysVelocity * descriptors::invCs2<T,DESCRIPTOR>() / (latticeRelaxationTime - 0.5)),
        (physViscosity * charLatticeVelocity * charLatticeVelocity / charPhysVelocity / charPhysVelocity * descriptors::invCs2<T,DESCRIPTOR>() / (latticeRelaxationTime - 0.5)),
        charPhysLength,
        charPhysVelocity,
        physViscosity,
        physDensity,
        charPhysPressure)
  {
  }
};
 
}  // namespace olb
 
#include "thermalUnitConverter.h"
#include "multiPhaseUnitConverter.h"
#endif