units.h

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library 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 Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef UNITS_H
#define UNITS_H

#include "core/globalDefs.h"
#include "core/globalDefs.h"
#include "parallelism/mpiManager.h"
#include "io/parallelIO.h"
#include <string>
#include <fstream>

namespace plb {

template<typename T>
class IncomprFlowParam {
public:

    IncomprFlowParam(T physicalU_, T latticeU_, T Re_, T physicalLength_, plint resolution_, T lx_, T ly_, T lz_=T() )
        : latticeU(latticeU_), physicalU(physicalU_), Re(Re_), physicalLength(physicalLength_), 
          resolution(resolution_), lx(lx_), ly(ly_), lz(lz_)
    { }
    
    IncomprFlowParam(T latticeU_, T Re_, plint resolution_, T lx_, T ly_, T lz_=T() )
    : latticeU(latticeU_), Re(Re_), resolution(resolution_), lx(lx_), ly(ly_), lz(lz_)
    { 
        physicalU      = (T)1;
        physicalLength = (T)1; 
    }
    T getLatticeU() const { return latticeU; }
    T getPhysicalU() const { return physicalU; }
    T getRe() const      { return Re; }
    T getPhysicalLength() const { return physicalLength; }
    plint getResolution() const { return resolution; }
    T getLx() const      { return getPhysicalLength()*lx; }
    T getLy() const      { return getPhysicalLength()*ly; }
    T getLz() const      { return getPhysicalLength()*lz; }
    T getDeltaX() const  { return (T)getPhysicalLength()/(T)getResolution(); }
    T getDeltaT() const  { return getDeltaX()*getLatticeU()/getPhysicalU(); }
    plint nCell(T l) const { return (int)(l/getDeltaX()+(T)0.5); }
    plint nStep(T t) const { return (int)(t/getDeltaT()+(T)0.5); }
    plint getNx(bool offLattice=false) const { return nCell(getLx())+1+(int)offLattice; }
    plint getNy(bool offLattice=false) const { return nCell(getLy())+1+(int)offLattice; }
    plint getNz(bool offLattice=false) const { return nCell(getLz())+1+(int)offLattice; }
    T getLatticeNu() const { return getLatticeU()*(T)getResolution()/Re; }
    T getTau() const       { return (T)3*getLatticeNu()+(T)0.5; }
    T getOmega() const     { return (T)1 / getTau(); }
private:
    T physicalU, latticeU, physicalLength, Re;
    plint resolution;
    T lx, ly, lz;
};

template<typename T>
void writeLogFile(IncomprFlowParam<T> const& parameters,
                  std::string const& title)
{
    std::string fullName = global::directories().getLogOutDir() + "plbLog.dat";
    plb_ofstream ofile(fullName.c_str());
    ofile << title << "\n\n";
    ofile << "Velocity in lattice units: u=" << parameters.getLatticeU() << "\n";
    ofile << "Reynolds number:           Re=" << parameters.getRe() << "\n";
    ofile << "Lattice resolution:        N=" << parameters.getResolution() << "\n";
    ofile << "Relaxation frequency:      omega=" << parameters.getOmega() << "\n";
    ofile << "Extent of the system:      lx=" << parameters.getLx() << "\n";
    ofile << "Extent of the system:      ly=" << parameters.getLy() << "\n";
    ofile << "Extent of the system:      lz=" << parameters.getLz() << "\n";
    ofile << "Grid spacing deltaX:       dx=" << parameters.getDeltaX() << "\n";
    ofile << "Time step deltaT:          dt=" << parameters.getDeltaT() << "\n";
}

template<typename T>
class ComprFlowParam {
public:

    ComprFlowParam(T latticeU_, T latticeRho_, T latticeTemp_, T physU_, T physRho_, T physTemp_, 
                   T Re_, T Pe_, plint resolution_, T lx_, T ly_, T lz_=T() )
        : latticeU(latticeU_), latticeRho(latticeRho_), latticeTemp(latticeTemp_), 
          physRho(physRho_), physU(physU_), physTemp(physTemp_), 
          Re(Re_), Pe(Pe_), resolution(resolution_), lx(lx_), ly(ly_), lz(lz_)
    { }
    T getLatticeU() const { return latticeU; }
    T getLatticeRho() const { return latticeRho; }
    T getLatticeTemp() const { return latticeTemp; }
    T getPhysicalU() const { return physU; }
    T getPhysicalRho() const { return physRho; }
    T getPhysicalTemp() const { return physTemp; }
    T getRe() const      { return Re; }
    T getPe() const      { return Pe; }
    plint getResolution() const { return resolution; }
    T getLx() const      { return lx; }
    T getLy() const      { return ly; }
    T getLz() const      { return lz; }
    T getDeltaX() const  { return (T)1/(T)getResolution(); }
    T getDeltaT() const  { return getDeltaX() * getLatticeU() / getPhysicalU(); }
    T getDeltaRho() const  { return getPhysicalRho() / getLatticeRho(); }
    T getDeltaTemp() const  { return getPhysicalTemp() / getLatticeTemp(); }
    plint nCell(T l) const { return (int)(l/getDeltaX()+(T)0.5); }
    plint nStep(T t) const { return (int)(t/getDeltaT()+(T)0.5); }
    plint getNx(bool offLattice=false) const { return nCell(lx)+1+(int)offLattice; }
    plint getNy(bool offLattice=false) const { return nCell(ly)+1+(int)offLattice; }
    plint getNz(bool offLattice=false) const { return nCell(lz)+1+(int)offLattice; }
    T getLatticeMu() const { return getLatticeU()*getResolution()*getLatticeRho() / Re; }
    T getLatticeKappa() const { return getLatticeU()*getResolution()*getLatticeRho() / Pe; }
    T getTau() const       { return getLatticeMu() / getLatticeRho() / getLatticeTemp()+(T)0.5; }
    // TODO for the moment only Pr = 1 (Pe = Re) fluids are simulable....
    T getOmega() const     { return (T)1 / getTau(); }
private:
    T latticeU, latticeRho, latticeTemp, physRho, physU, physTemp, Re, Pe;
    plint resolution;
    T lx, ly, lz;
};

template<typename T>
void writeLogFile(ComprFlowParam<T> const& parameters,
                  std::string const& title)
{
    std::string fullName = global::directories().getLogOutDir() + "plbLog.dat";
    plb_ofstream ofile(fullName.c_str());
    ofile << title << "\n\n";
    ofile << "Velocity in lattice units:    u=" << parameters.getLatticeU() << "\n";
    ofile << "Density in lattice units:     rho=" << parameters.getLatticeRho() << "\n";
    ofile << "Temperature in lattice units: T=" << parameters.getLatticeTemp() << "\n";
    ofile << "Reynolds number:              Re=" << parameters.getRe() << "\n";
    ofile << "Peclet number:                Pe=" << parameters.getPe() << "\n";
    ofile << "Lattice resolution:           N=" << parameters.getResolution() << "\n";
    ofile << "Extent of the system:         lx=" << parameters.getLx() << "\n";
    ofile << "Extent of the system:         ly=" << parameters.getLy() << "\n";
    ofile << "Extent of the system:         lz=" << parameters.getLz() << "\n";
    ofile << "Grid spacing deltaX:          dx=" << parameters.getDeltaX() << "\n";
    ofile << "Time step deltaT:             dt=" << parameters.getDeltaT() << "\n";
    ofile << "Density deltaRho:             dRho=" << parameters.getDeltaRho() << "\n";
    ofile << "Temp deltaTemp:               dTemp=" << parameters.getDeltaTemp() << "\n";
    ofile << "Relaxation time:              Tau=" << parameters.getTau() << "\n";
}

}  // namespace plb

#endif