ViscidUtil.H

Go to the documentation of this file.

#ifndef __VISCID_UTIL_H__
#define __VISCID_UTIL_H__

#include "PlasCom2.H"
#include "PCPPCommUtil.H"
#include "PCPPIntervalUtils.H"
#include "PCPPReport.H"

namespace viscid {
  namespace util {

    typedef simulation::grid::subregion GridRegionType;
    // Input:
    // inInterval:  (xStart xEnd yStart yEnd zStart zEnd)
    // bufferSize:  (numX numY numZ)
    // stateBuffer: [ *rho *rhoV *rhoE *scalarVars ]
    //
    // Output:
    // tvBuffer:    [mu lambda kappa]
    // 
    // simple implementation, ported from PlasComCM
    int ComputeTVBufferPower(const pcpp::IndexIntervalType &regionInterval,
                        const std::vector<size_t>     &bufferSize,
                        const double *                 temperatureBuffer,
                        std::vector<double *>         &tvBuffer,
                        const double                   beta,
                        const double                   power,
                        const double                   bulkViscFac,
                        const double                   specificHeat,
                        const double                   prandtlNumber);

    // Fortran-like interfaces 
    // Uses Fortran Kernels
    // MJA TBD, no yet implemented
    int ComputeTVBuffer(const int *numDimPtr,
                        const size_t *bufferSize,
                        const size_t *numPointsPtr,
                        const size_t *bufferInterval,
                        const double *rhoBuffer,
                        const double *rhoVBuffer,
                        const double *rhoEBuffer,
                        double *pressureBuffer,
                        double *tempBuffer,
                        double *rhom1Buffer,
                        double *velBuffer);

    // Fortran-like interfaces 
    int ComputeTauBuffer(const pcpp::IndexIntervalType &regionInterval,
                         const std::vector<size_t>     &bufferSize,
                         int                            gridType,
                         const std::vector<double>     &gridMetrics,
                         const std::vector<double>     &gridJacobian,
                         const double                  &Re,
                         const std::vector<double *>   &tvBuffers,
                         const std::vector<double *>   &velGradBuffers,
                         std::vector<double *>         &tauBuffers);

    // Fortran-like interfaces 
//     int ComputeHeatFluxBuffer(const pcpp::IndexIntervalType &regionInterval,
//                         const std::vector<size_t>     &bufferSize,
//                         const std::vector<double>     &gridMetrics,
//                         const double                  &gridJacobian,
//                         const double                  &Re,
//                         const std::vector<double *>   &tvBuffers,
//                         const std::vector<double *>   &temperatureGradBuffers,
//                         std::vector<double *>         &heatFluxBuffers);

    
    int ComputeHeatFluxBuffer(const pcpp::IndexIntervalType &regionInterval,
                              const std::vector<size_t>     &bufferSize,
                              int                            gridType,
                              const std::vector<double>     &gridMetrics,
                              const std::vector<double>     &gridJacobian,
                              const double                  &Re,
                              const std::vector<double *>   &tvBuffers,
                              const std::vector<double *>   &temperatureGradBuffers,
                              std::vector<double *>         &heatFluxBuffers);

    // initialization for Poiseulle problems
    //template<typename GridType,typename StateType,typename BoundaryType>
    template<typename GridType,typename StateType>
    int InitializePoiseuille(const GridType &inGrid,StateType &inState, StateType &inParamState,
                             const std::vector<double> &inParams,int threadId,
                             std::ostream *messageStream = NULL)
    {
      
      if(!messageStream)
        messageStream = &std::cout;

      double inletX = 0.0; // inlet x coordinate
      double inletY = 0.0; // inlet y coordinate
      double inletZ = 0.0; // inlet z coordinate
      double deltaP = .01; // pressure differential across pipe length
      
      int numParams = inParams.size();
      if(numParams > 0)
        deltaP = inParams[0];
      if(numParams > 1)
        inletX = inParams[1];
      if(numParams > 2)
        inletY = inParams[2];
      

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      int numDim = gridSizes.size();
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &readGridExtent(inGrid.PhysicalExtent());

      if(numDim > 2 && numParams>3)
        inletZ = inParams[3];

      // check on dimensionality
      int nonDimensionalMode=0;
      {
        int paramHandle = inParamState.GetDataIndex("nonDimensional");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(inParamState.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          nonDimensionalMode = *paramPtr;
        } else {
          nonDimensionalMode = 0;
        }
      }

      // get reference conditions
      double *gammaPtr  = inParamState.template GetFieldBuffer<double>("gamma");
      double *rhoRefPtr  = inParamState.template GetFieldBuffer<double>("refRho");
      double *refLengthPtr  = inParamState.template GetFieldBuffer<double>("refLength");
      double *presRefPtr  = inParamState.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr  = inParamState.template GetFieldBuffer<double>("refTemperature");

      double gamma = *gammaPtr;
      double refLength = *refLengthPtr;
      double refRho = *rhoRefPtr;
      double refPressure = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef = sqrt(gamma*refPressure/refRho);

      // adjust input parameters
      if(nonDimensionalMode == 1) {
        inletX /= refLength;
        inletY /= refLength;
        inletZ /= refLength;
      } else if (nonDimensionalMode ==2) {
        inletX /= refLength;
        inletY /= refLength;
        inletZ /= refLength;
      }

      // determine the flow direction 
      int direction = 0;
      double xMin = readGridExtent[0];
      double xMax = readGridExtent[1];
      double yMin = readGridExtent[2];
      double yMax = readGridExtent[3];
      double zMin = 0.;
      double zMax = 0.;
      if(numDim > 2 ) {
        zMin = readGridExtent[4];
        zMax = readGridExtent[5];
      }
      double height = 0.;
      double length = 0.;
      double width = 0.;

      //*messageStream << "xMin " << xMin << std::endl;
      //*messageStream << "xMax " << xMax << std::endl;
      //*messageStream << "yMin " << yMin << std::endl;
      //*messageStream << "yMax " << yMax << std::endl;
      //*messageStream << "inletX " << inletX << std::endl;
      //*messageStream << "inletY " << inletY << std::endl;
      if(inletX == xMin) {
        // positive X
        if((inletY == yMin || inletY == yMax) ||
          ((numDim > 2) && (inletZ == zMin || inletZ == zMax))) {
          *messageStream << "Invalid inlet specification, point must live on one boundary only! " << std::endl;
          return(1);
        }
        *messageStream << "Initializing Poiseulle flow in the positive X-direction" << std::endl;
        height = yMax - yMin;
        length = xMax - xMin;
        if(numDim > 2)
          width = zMax - zMin;
        direction = 1;
      } else if (inletX == xMax) {
        // negative X
        if((inletY == yMin || inletY == yMax) || 
          ((numDim > 2) && (inletZ == zMin || inletZ == zMax))) {
          *messageStream << "Invalid inlet specification, point must live on one boundary only! " << std::endl;
          return(1);
        }
        *messageStream << "Initializing Poiseulle flow in the negative X-direction" << std::endl;
        height = yMax - yMin;
        length = xMax - xMin;
        if(numDim > 2)
          width = zMax - zMin;
        direction = -1;
      } else if (inletY == yMin) {
        // positive Y
        if((inletX == xMin || inletX == xMax) ||
          ((numDim > 2) && (inletZ == zMin || inletZ == zMax))) {
          *messageStream << "Invalid inlet specification, point must live on one boundary only! " << std::endl;
          return(1);
        }
        *messageStream << "Initializing Poiseulle flow in the positive Y-direction" << std::endl;
        direction = 2;
        height = xMax - xMin;
        length = yMax - yMin;
        if(numDim > 2)
          width = zMax - zMin;
      } else if (inletY == yMax) {
        // negative Y
        if((inletX == xMin || inletX == xMax) ||
          ((numDim > 2) && (inletZ == zMin || inletZ == zMax))) {
          *messageStream << "Invalid inlet specification, point must live on one boundary only! " << std::endl;
          return(1);
        }
        *messageStream << "Initializing Poiseulle flow in the negative Y-direction" << std::endl;
        height = xMax - xMin;
        length = yMax - yMin;
        if(numDim > 2)
          width = zMax - zMin;
        direction = -2;
      } else if (inletZ == zMin) {
        // positive Z, 3D
        if(inletX == xMin || inletX == xMax || inletY == yMin || inletY == yMax) {
          *messageStream << "Invalid inlet specification, point must live on one boundary only! " << std::endl;
          return(1);
        }
        *messageStream << "Initializing 3D Poiseulle flow in the positive Z-direction" << std::endl;
        height = xMax - xMin;
        width = yMax - yMin;
        length = zMax - zMin;
        direction = 3;
      } else if (inletZ == zMax) {
        // negative Z, 3D
        if(inletX == xMin || inletX == xMax || inletY == yMin || inletY == yMax) {
          *messageStream << "Invalid inlet specification, point must live on one boundary only! " << std::endl;
          return(1);
        }
        *messageStream << "Initializing 3D Poiseulle flow in the positive Z-direction" << std::endl;
        height = xMax - xMin;
        width = yMax - yMin;
        length = zMax - zMin;
        direction = -3;
      } else {
        *messageStream << "Could not dermine flow direction from inlet specification" << std::endl;
        return(1);
      }
        
      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);

      *messageStream << "Poiseulle: Grid Sizes: (";
      pcpp::io::DumpContents(*messageStream,gridSizes,",");
      *messageStream << ")" << std::endl
                << "Poiseulle: Partition Interval: ";
      partitionInterval.PrettyPrint(*messageStream);
      *messageStream << std::endl << "Poiseulle: Partition Buffer Interval: ";
      partitionBufferInterval.PrettyPrint(*messageStream);
      *messageStream << std::endl;
      *messageStream << "Poiseulle: pipe height " << height << std::endl;
      if(numDim > 2)
        *messageStream << "Poiseulle: pipe width " << width << std::endl;
      *messageStream << "Poiseulle: pipe length " << length << std::endl;

      const std::vector<double> &coordinateData(inGrid.CoordinateData());

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }

      // nonDimensional use Re, otherwise calculate mu from the transport model
      // Re will be set to 1 for dimensionalMode
      double *RePtr  = inParamState.template GetFieldBuffer<double>("refRe");
      double *betaPtr  = inParamState.template GetFieldBuffer<double>("beta");
      double *powerPtr  = inParamState.template GetFieldBuffer<double>("power");
      // assume Power law, dimensional
      double mu = (*betaPtr)*pow(refTemperature,(*powerPtr));
      double Re = 1;

      if(nonDimensionalMode == 1){
        mu = 1.0;
        Re = *RePtr;
      } if(nonDimensionalMode == 2){
        mu = 1.0;
        Re = *RePtr;
      }

      // flow conditions
      // working in dimensional units
      double rho = refRho;
      std::vector<double> velocity(numDim,0);
      double pLo = refPressure; 
      double pHi = pLo + deltaP;

      double tempLo = refTemperature;

      if(nonDimensionalMode == 1){
        rho /= refRho;
        pLo /= refPressure;
        pHi /= refPressure;
        tempLo /= refTemperature;
      } else if (nonDimensionalMode == 2) {
        rho /= refRho;
        pLo /= gamma*refPressure;
        pHi /= gamma*refPressure;
        tempLo /= (gamma-1)*refTemperature;
      }

      double rhoE = pLo/(gamma-1)/rho;

      double dPdX = (pHi - pLo)/length;
      //*messageStream << "dPdX " << dPdX << std::endl;

      *messageStream << "Poiseulle: Initial state" << std::endl;
      *messageStream << "\t nonDimensional: " << nonDimensionalMode << std::endl;
      //*messageStream << "\t Reynolds Number: " << Re << std::endl;
      *messageStream << "\t rho: " << rho << std::endl;
      *messageStream << "\t Inlet pressure: " << pLo << std::endl;
      *messageStream << "\t Outlet pressure: " << pHi << std::endl;
      *messageStream << "\t temperature: " << tempLo << std::endl;
      *messageStream << "\t rhoE: " << rhoE << std::endl;
      *messageStream << "\t dPdX: " << dPdX << std::endl;
      //*messageStream << "\t u1: " << velocity[0] << std::endl;
      //*messageStream << "\t u2: " << velocity[1] << std::endl;

      // initialize the interior points
      // apply the linear pressure gradient in the flow direction
      double pi = 3.14159265359;
      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;
      size_t kStart = 0;
      size_t kEnd   = 0;
      size_t nPlane = 0;
      if(numDim > 2 ){
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
        nPlane = bufferSizes[0]*bufferSizes[1];
      }

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        size_t kPartIndex = 0;
        double z = 0;
        if (numDim > 2){
          kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          z = kPartIndex*dX[2]; 
        }

        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
          double y = jPartIndex*dX[1];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
            double x = iPartIndex*dX[0];

            double pressure = pLo;
            if (direction == 1) {
              pressure += dPdX*(xMax-x);
              if(x == xMin || x == xMax)
                velocity[0] = Re*dPdX*y*(height-y)/(2*mu);
              else
                velocity[0]=0.;
            } else if (direction == -1) {
              pressure += dPdX*(x-xMin);
              if(x == xMin || x == xMax)
                velocity[0] = Re*dPdX*y*(y-height)/(2*mu);
              else
                velocity[0]=0.;
            } else if (direction == 2) {
              pressure += dPdX*(yMax-y);
              if(y == yMin || y == yMax)
                velocity[1] = Re*dPdX*x*(height-x)/(2*mu);
              else
                velocity[1]=0.;
            } else if (direction == -2) {
              pressure += dPdX*(y-yMin);
              if(y == yMin || y == yMax)
                velocity[1] = Re*dPdX*x*(x-height)/(2*mu);
              else
                velocity[1]=0.;
            } else if (direction == 3) {
              pressure += dPdX*(zMax-z);
              if(z == zMin || z == zMax) {

                velocity[2] = Re*dPdX*x*(height-x)/(2*mu);

                // height = x
                // width = y
                // length = z
                double sumTerm = 0;
                int nMax=10;
                for(int n=1;n<nMax;n++){
                  double beta = (2*n-1)*pi/height;
                  sumTerm += sin(beta*x)/pow((2*n-1),3)*(sinh(beta*y)+sinh(beta*(width-y)))/sinh(beta*width);
                }
                sumTerm *= 4*dPdX*height*height/mu/pow(pi,3);
                velocity[2] -= sumTerm;
              } else {
                velocity[2]=0;
              }
             
            } else if (direction == -3) {
              pressure += dPdX*(z-zMin);
              if(z == zMin || z == zMax) {

                velocity[2] = Re*dPdX*x*(x-height)/(2*mu);

                // height = x
                // width = y
                // length = z
                double sumTerm = 0;
                int nMax=10;
                for(int n=1;n<nMax;n++){
                  double beta = (2*n-1)*pi/height;
                  sumTerm += sin(beta*x)/pow((2*n-1),3)*(sinh(beta*y)+sinh(beta*(width-y)))/sinh(beta*width);
                }
                sumTerm *= 4*dPdX*height*height/mu/pow(pi,3);
                velocity[2] += sumTerm;
              } else {
                velocity[2]=0;
              }
            }

            rhoPtr[bufferIndex]  = rho;
            rhoEPtr[bufferIndex] = pressure/(gamma-1);
            for (int iDim=0;iDim<numDim;iDim++){
              rhoVPtr[bufferIndex+numPointsBuffer*iDim] = rhoPtr[bufferIndex]*velocity[iDim];
              rhoEPtr[bufferIndex] += 0.5*rho*velocity[iDim]*velocity[iDim];
            }
          }
        }
      }

      return(0);
    };

  }
}
#endif