Initialize at openlb-1-3

1 files changed, 209 insertions, 0 deletions

diff --git a/src/core/radiativeUnitConverter.h b/src/core/radiativeUnitConverter.h
new file mode 100644
index 0000000..08551ad
--- /dev/null
+++ b/src/core/radiativeUnitConverter.h
@@ -0,0 +1,209 @@
+/*  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.
+*/
+
+/** \file
+ * Unit conversion handling -- header file.
+ */
+
+#ifndef RADIATIVEUNITCONVERTER_H
+#define RADIATIVEUNITCONVERTER_H
+
+
+#include "io/ostreamManager.h"
+#include "core/unitConverter.h"
+
+
+/// All OpenLB code is contained in this namespace.
+namespace olb {
+
+double getThetaRefracted(double const& thetaIncident, double const& refractiveRelative);
+double getFresnelFunction(double const& theta, double const& refractiveRelative);
+double R_phi_diff (double const& theta, double const& refractiveRelative);
+double R_j_diff (double const& theta, double const& refractiveRelative);
+double getRefractionFunction(const double& refractiveRelative);
+double getRefractionFunction(const double& refractiveRelative);
+double getPartialBBCoefficient(double const& latticeDiffusionCoefficient, double const& relativeRefractiveIndex );
+
+// forward declaration
+template<typename T, typename DESCRIPTOR> class RadiativeUnitConverter;
+
+// wrapper for above function
+template <typename T, typename DESCRIPTOR>
+double getPartialBBCoefficient(RadiativeUnitConverter<T,DESCRIPTOR> const& converter)
+{
+  return getPartialBBCoefficient( converter.getLatticeDiffusion(), converter.getRefractiveRelative() );
+};
+
+// wrapper for above function
+template <typename T, typename DESCRIPTOR>
+double getRefractionFunction(RadiativeUnitConverter<T,DESCRIPTOR> const& converter)
+{
+  return getRefractionFunction(converter.getRefractiveRelative());
+};
+
+/** Conversion between physical and lattice units, as well as discretization.
+* 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__
+*
+*/
+template <typename T, typename DESCRIPTOR>
+class RadiativeUnitConverter : public UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR> {
+public:
+  /** Documentation of constructor:
+    *   \param resolution   is number of voxel per 1 meter
+    *   \param latticeRelaxationTime    see class UnitConverterFromResolutionAndRelaxationTime
+    *   \param physAbsorption   physical absorption in 1/meter
+    *   \param physScattering   physical scattering in 1/meter
+    */
+  constexpr RadiativeUnitConverter( int resolution, T latticeRelaxationTime, T physAbsorption, T physScattering, T anisotropyFactor=0, T charPhysLength=1, T refractiveMedia=1, T refractiveAmbient=1 )
+    : UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR>( resolution, latticeRelaxationTime, charPhysLength, T(1), T(1), T(1) ),
+      clout(std::cout, "RadiativeUnitConverter"),
+      _physAbsorption(physAbsorption),
+      _physScattering(physScattering),
+      _anisotropyFactor(anisotropyFactor),
+      _extinction( physAbsorption+physScattering ),
+      _scatteringAlbedo( physScattering/(physAbsorption+physScattering) ),
+      _physDiffusion( 1.0 / (3.0*(physAbsorption+physScattering)) ),
+      _effectiveAttenuation( std::sqrt(3*physAbsorption*(physAbsorption+physScattering)) ),
+      _refractiveRelative(refractiveMedia/refractiveAmbient),
+      _latticeAbsorption( physAbsorption*this->getConversionFactorLength() ),
+      _latticeScattering( physScattering*this->getConversionFactorLength() ),
+      _latticeDiffusion(_physDiffusion/this->getConversionFactorLength())
+  { };
+
+  constexpr T getPhysAbsorption() const
+  {
+    return _physAbsorption;
+  };
+
+  constexpr T getPhysScattering() const
+  {
+    return _physScattering;
+  };
+
+  constexpr T getAnisotropyFactor() const
+  {
+    return _anisotropyFactor;
+  };
+
+  constexpr T getExtinction() const
+  {
+    return _extinction;
+  };
+
+  constexpr T getScatteringAlbedo() const
+  {
+    return _scatteringAlbedo;
+  };
+
+  constexpr T getPhysDiffusion() const
+  {
+    return _physDiffusion;
+  };
+
+  constexpr T getEffectiveAttenuation() const
+  {
+    return _effectiveAttenuation;
+  };
+
+  constexpr T getLatticeAbsorption() const
+  {
+    return _latticeAbsorption;
+  };
+
+  constexpr T getLatticeScattering() const
+  {
+    return _latticeScattering;
+  };
+
+  constexpr T getLatticeDiffusion() const
+  {
+    return _latticeDiffusion;
+  };
+
+  constexpr T getRefractiveRelative() const
+  {
+    return _refractiveRelative;
+  };
+
+  void print() const override;
+
+private:
+  mutable OstreamManager clout;
+
+  double _physAbsorption;
+  double _physScattering;
+  double _anisotropyFactor;
+  double _extinction;
+  double _scatteringAlbedo;
+  double _physDiffusion;
+  double _effectiveAttenuation;
+
+  double _refractiveRelative;
+
+  double _latticeAbsorption;
+  double _latticeScattering;
+  double _latticeDiffusion;
+};
+
+template <typename T, typename DESCRIPTOR>
+void RadiativeUnitConverter<T, DESCRIPTOR>::print() const
+{
+  clout << "----------------- UnitConverter information -----------------" << std::endl;
+  clout << "-- Parameters:" << std::endl;
+  clout << "Resolution:                       N=              " << this->getResolution() << std::endl;
+  clout << "Lattice relaxation frequency:     omega=          " << this->getLatticeRelaxationFrequency() << std::endl;
+  clout << "Lattice relaxation time:          tau=            " << this->getLatticeRelaxationTime() << std::endl;
+  clout << "Characteristical length(m):       charL=          " << this->getCharPhysLength() << std::endl;
+  clout << "Phys. density(kg/m^d):            charRho=        " << this->getPhysDensity() << std::endl;
+  clout << "Phys. absorption(1/m):            mu_a=           " << getPhysAbsorption() << std::endl;
+  clout << "Phys. scattering(1/m):            mu_s=           " << getPhysScattering() << std::endl;
+  clout << "Extinction(1/m):                  mu_t=           " << getExtinction() << std::endl;
+  clout << "Effective attenuation(1/m):       mu_eff=         " << getEffectiveAttenuation() << std::endl;
+  clout << "Phys. diffusion(m):               D=              " << getPhysDiffusion() << std::endl;
+  clout << "Single scattering albedo:         albedo=         " << getScatteringAlbedo() << std::endl;
+  clout << "Anisotropy factor:                g=              " << getAnisotropyFactor() << std::endl;
+
+  clout << std::endl;
+  clout << "Lattice diffusion:                D^*=            " << getLatticeDiffusion() << std::endl;
+  clout << "Lattice absorption:               absorption=     " << getLatticeAbsorption() << std::endl;
+  clout << "Lattice scattering:               scattering=     " << getLatticeScattering() << std::endl;
+  clout << "Lattice sink:                     sink=           " << 3./8.*getLatticeAbsorption()*(getLatticeScattering()+getLatticeAbsorption()) << std::endl;
+  clout << "C_R: " << getRefractionFunction(getRefractiveRelative()) << std::endl;
+  clout << "r_F: " << getPartialBBCoefficient(getLatticeDiffusion(),getRefractiveRelative()) << std::endl;
+
+  clout << std::endl;
+  clout << "-- Conversion factors:" << std::endl;
+  clout << "Voxel length(m):                  physDeltaX=     " << this->getConversionFactorLength() << std::endl;
+  clout << "Time step(s):                     physDeltaT=     " << this->getConversionFactorTime() << std::endl;
+  clout << "Density factor(kg/m^3):           physDensity=    " << this->getConversionFactorDensity() <<  std::endl;
+  clout << "-------------------------------------------------------------" << std::endl;
+
+}
+
+
+}  // namespace olb
+
+#endif