Initialize at openlb-1-3

1 files changed, 173 insertions, 0 deletions

diff --git a/src/core/powerLawUnitConverter.h b/src/core/powerLawUnitConverter.h
new file mode 100644
index 0000000..241a1e5
--- /dev/null
+++ b/src/core/powerLawUnitConverter.h
@@ -0,0 +1,173 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2017-2018 Max Gaedtke, Albert Mink, Davide Dapelo
+ *  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.
+*/
+
+/** \file
+ * Unit conversion handling -- header file.
+ */
+
+#ifndef PL_UNITCONVERTER_H
+#define PL_UNITCONVERTER_H
+
+
+#include <math.h>
+#include "io/ostreamManager.h"
+#include "core/util.h"
+#include "io/xmlReader.h"
+#include "unitConverter.h"
+
+
+/// All OpenLB code is contained in this namespace.
+namespace olb {
+
+
+
+/** Conversion between physical and lattice units, as well as discretization specialized for power-law rheology: \f$\nu=m^{n-1}\f$, with \f$m$ being the consistency coefficient and \f$n\f$ the power-law index.
+ *  The Newtonian case is recovered for \f$n=1\f$.
+ *  A characteristic (apparent) kinematic viscosity for problem similarity (e.g., to compute the Reynolds number) is obtained as __physViscosity = poweLawConsistencyCoeff * [physVelocity / (2*physLength)]^(n-1)__.
+* Be aware of the nomenclature:
+* We distingish between physical (dimensioned) and lattice (dimensionless) values.
+* A specific conversion factor maps the two different scopes,
+* e.g. __physLength = conversionLength * latticeLength__
+*
+* For pressure and temperature we first shift the physical values by a characteristic value to asure a lattice pressure and lattice temperature between 0 and 1, e.g. __physPressure - charPhysPressure = conversionPressure * latticePressure__
+*
+*  \param latticeRelaxationTime   relaxation time, have to be greater than 0.5!
+*  - - -
+*  \param physViscosity         physical kinematic viscosity in __m^2 / s__
+*  \param poweLawIndex          power-law index
+*  \param physConsistencyCoeff  physicsl consistency coefficient __m^2 s^(n-2)__
+*  \param physDensity           physical density in __kg / m^3__
+*  - - -
+*  \param conversionLength      conversion factor for length __m__
+*  \param conversionTime        conversion factor for time __s__
+*  \param conversionMass        conversion factor for mass __kg__
+*  - - -
+*  \param conversionVelocity    conversion velocity __m / s__
+*  \param conversionViscosity   conversion kinematic viscosity __m^2 / s__
+*  \param conversionConsistencyCoeff conversion power-law consistency coefficient __m^2 s^(n-2)__
+*  \param conversionDensity     conversion density __kg / m^3__
+*  \param conversionForce       conversion force __kg m / s^2__
+*  \param conversionPressure    conversion pressure __kg / m s^2__
+*  - - -
+*  \param resolution            number of grid points per charPhysLength
+*  - - -
+*  \param charLatticeVelocity
+*/
+template <typename T, typename DESCRIPTOR>
+class PowerLawUnitConverter : public UnitConverter<T, DESCRIPTOR> {
+public:
+  /** Documentation of constructor:
+    *  \param physDeltaX              spacing between two lattice cells in __m__
+    *  \param physDeltaT              time step in __s__
+    *  \param charPhysLength          reference/characteristic length of simulation geometry in __m__
+    *  \param charPhysVelocity        maximal or highest expected velocity during simulation in __m / s__
+    *  \param physConsistencyCoeff    physical power-law consistency coefficient in __m^2 s^(n-2)__
+    *  \param powerLawIndex           Power-law index
+    *  \param physDensity             physical density in __kg / m^3__
+    *  \param charPhysPressure        reference/characteristic physical pressure in Pa = kg / m s^2
+    */
+  constexpr PowerLawUnitConverter( T physDeltaX, T physDeltaT, T charPhysLength,
+                                   T charPhysVelocity, T physConsistencyCoeff, T powerLawIndex,
+                                   T physDensity, T charPhysPressure = 0 )
+    : UnitConverter<T, DESCRIPTOR>( physDeltaX, physDeltaT, charPhysLength, charPhysVelocity,
+                                 physConsistencyCoeff * pow(charPhysVelocity / (2*charPhysLength), powerLawIndex-1),
+                                 physDensity, charPhysPressure ),
+      _conversionConsistencyCoeff(pow( physDeltaT,powerLawIndex-2 ) * pow( physDeltaX,2 ) ),
+      _powerLawIndex(powerLawIndex),
+      _physConsistencyCoeff(physConsistencyCoeff),
+      clout(std::cout,"PowerLawUnitConverter")
+  {
+  }
+
+  /// return consistency coefficient in physical units
+  constexpr T getPhysConsistencyCoeff( ) const
+  {
+    return _physConsistencyCoeff;
+  }
+  /// conversion from lattice to  physical consistency coefficient
+  constexpr T getPhysConsistencyCoeff( T latticeConsistencyCoeff ) const
+  {
+    return _conversionConsistencyCoeff * latticeConsistencyCoeff;
+  }
+  /// conversion from physical to lattice consistency coefficient
+  constexpr T getLatticeConsistencyCoeff(  ) const
+  {
+    return _physConsistencyCoeff / _conversionConsistencyCoeff;
+  }
+  /// access (read-only) to private member variable
+  constexpr T getConversionFactorConsistencyCoeff() const
+  {
+    return _conversionConsistencyCoeff;
+  }
+
+  /// access (read-only) to private member variable
+  constexpr T getPowerLawIndex() const
+  {
+    return _powerLawIndex;
+  }
+
+  /// nice terminal output for conversion factors, characteristical and physical data
+  virtual void print() const;
+
+  void write(std::string const& fileName = "unitConverter") const;
+
+protected:
+  // conversion factors
+  const T _conversionConsistencyCoeff;         // m^2 s^(n-2)
+
+  // physical units, e.g characteristic or reference values
+  const T _powerLawIndex;
+  const T _physConsistencyCoeff;         // m^2 s^(n-2)
+
+private:
+  mutable OstreamManager clout;
+};
+
+/// creator function with data given by a XML file
+template <typename T, typename DESCRIPTOR>
+PowerLawUnitConverter<T, DESCRIPTOR>* createPowerLawUnitConverter(XMLreader const& params);
+
+template <typename T, typename DESCRIPTOR>
+class PowerLawUnitConverterFrom_Resolution_RelaxationTime_Reynolds_PLindex : public PowerLawUnitConverter<T, DESCRIPTOR> {
+public:
+  constexpr PowerLawUnitConverterFrom_Resolution_RelaxationTime_Reynolds_PLindex(
+    int resolution,
+    T latticeRelaxationTime,
+    T charPhysLength,
+    T charPhysVelocity,
+    T Re,
+    T powerLawIndex,
+    T physDensity,
+    T charPhysPressure = 0)
+    : PowerLawUnitConverter<T, DESCRIPTOR>( (charPhysLength/resolution),
+                                         (latticeRelaxationTime - 0.5) / descriptors::invCs2<T,DESCRIPTOR>() * pow((charPhysLength/resolution),2) / ( ( charPhysLength * charPhysVelocity * pow( charPhysVelocity / ( 2 * charPhysLength ), 1 - powerLawIndex ) / Re ) * pow( charPhysVelocity / (2 * charPhysLength ), powerLawIndex - 1 ) ),
+                                         charPhysLength, charPhysVelocity,
+                                         charPhysLength * charPhysVelocity * pow( charPhysVelocity / ( 2 * charPhysLength ), 1 - powerLawIndex ) / Re, powerLawIndex, physDensity, charPhysPressure )
+  {
+  }
+
+};
+
+}  // namespace olb
+
+#endif