NavierStokesBC.C

Go to the documentation of this file.

#include "NavierStokesBC.H"
#include <cmath>

namespace simulation {
  namespace domain {
    namespace boundary {
      namespace bc {
        namespace navierstokes {
          
          int ResolveBCName(const std::string &inName)
          {
            for(int iType = 0;iType < (int)NUMBCTYPE;iType++){
              const std::string &bcName(bcNames[iType]);
              if(inName == bcName) return(iType);
            }
            return((int)NUMBCTYPE);
          }
          
          std::string ResolveBCType(int inType)
          {
            if(inType >= 0 && inType < (int)NUMBCTYPE)
              return(std::string(bcNames[inType]));
            return(std::string(bcNames[(int)NUMBCTYPE]));

          }
          

          void LinearSponge(const std::vector<size_t>     &bufferSizes,int normalDir,
                            const pcpp::IndexIntervalType &spongeRegion,
                            const pcpp::IndexIntervalType &partitionOpInterval,
                            const pcpp::IndexIntervalType &opBufferInterval,
                            const double *bcParams,const int *bcFlags,
                            const int numEquations,
                            const std::vector<double *>   &stateBuffers,
                            const std::vector<double *>   &targetBuffers,
                            std::vector<double *>         &rhsBuffers){
            
            int numDim = bufferSizes.size();
            double spongeAmp   = bcParams[0];
            double spongePower = bcParams[1];
            
            int opDir = normalDir;
            if(opDir < 0) opDir *= -1;
            opDir -= 1;
            
            size_t spongeStart = 0;
            size_t spongeEnd   = 0;
            size_t spongeSize  = 0;
            size_t sponge_end1 = spongeRegion[opDir].first;
            size_t sponge_end2 = spongeRegion[opDir].second;
            
            if(normalDir > 0){
              spongeStart = sponge_end2;
              spongeEnd   = sponge_end1;
              spongeSize = spongeStart - spongeEnd;
            } else {
              spongeStart = sponge_end1;
              spongeEnd   = sponge_end2;
              spongeSize = spongeEnd - spongeStart;
            }
            assert(spongeSize > 0);
            double spongeScale = 1.0/((double)spongeSize);            

            if(numDim == 3){
              
              size_t kPartStart = partitionOpInterval[2].first;
              size_t kPartEnd   = partitionOpInterval[2].second;
              size_t jPartStart = partitionOpInterval[1].first;
              size_t jPartEnd   = partitionOpInterval[1].second;
              size_t iPartStart = partitionOpInterval[0].first;
              size_t iPartEnd   = partitionOpInterval[0].second;
              
              size_t kBufferStart = opBufferInterval[2].first;
              size_t kBufferEnd   = opBufferInterval[2].second;
              size_t jBufferStart = opBufferInterval[1].first;
              size_t jBufferEnd   = opBufferInterval[1].second;
              size_t iBufferStart = opBufferInterval[0].first;
              size_t iBufferEnd   = opBufferInterval[0].second;
              
              size_t nPlane = bufferSizes[0]*bufferSizes[1];
              size_t kBufferIndex,jBufferIndex,iBufferIndex;
              
              size_t partIndex[] = {0,0,0};
              for(partIndex[2] = kPartStart,kBufferIndex = kBufferStart;
                  partIndex[2] <= kPartEnd;
                  partIndex[2]++,kBufferIndex++){
                size_t zIndex = kBufferIndex*nPlane;
                for(partIndex[1] = jPartStart,jBufferIndex = jBufferStart;
                    partIndex[1] <= jPartEnd;
                    partIndex[1]++,jBufferIndex++){
                  size_t yzIndex = jBufferIndex*bufferSizes[0] + zIndex;
                  for(partIndex[0] = iPartStart,iBufferIndex = iBufferStart;
                      partIndex[0] <= iPartEnd;
                      partIndex[0]++,iBufferIndex++){
                    size_t bufferIndex = yzIndex + iBufferIndex;
                    double spongeDistance = std::abs(double(partIndex[opDir]) 
                                                     - double(spongeStart));
                    double eta = spongeDistance*spongeScale;
                    double spongeFac = spongeAmp*std::pow(eta,spongePower);
                      
                    // rhs = rhs - spongeFac*(state - target)
                    for(int iEquation = 0;iEquation < numEquations;iEquation++){
                      rhsBuffers[iEquation][bufferIndex] = rhsBuffers[iEquation][bufferIndex] -
                        spongeFac*(stateBuffers[iEquation][bufferIndex] - 
                                   targetBuffers[iEquation][bufferIndex]);
                    }
                  }
                }
              }
            } else if(numDim == 2){
              
              size_t jPartStart = partitionOpInterval[1].first;
              size_t jPartEnd   = partitionOpInterval[1].second;
              size_t iPartStart = partitionOpInterval[0].first;
              size_t iPartEnd   = partitionOpInterval[0].second;
              
              size_t jBufferStart = opBufferInterval[1].first;
              size_t jBufferEnd   = opBufferInterval[1].second;
              size_t iBufferStart = opBufferInterval[0].first;
              size_t iBufferEnd   = opBufferInterval[0].second;
              
              size_t jBufferIndex,iBufferIndex;
              
              size_t partIndex[] = {0,0};
              for(partIndex[1] = jPartStart,jBufferIndex = jBufferStart;
                  partIndex[1] <= jPartEnd;
                  partIndex[1]++,jBufferIndex++){
                size_t yIndex = jBufferIndex*bufferSizes[0];
                for(partIndex[0] = iPartStart,iBufferIndex = iBufferStart;
                    partIndex[0] <= iPartEnd;
                    partIndex[0]++,iBufferIndex++){
                  size_t bufferIndex = yIndex + iBufferIndex;
                  double spongeDistance = std::abs(double(partIndex[opDir]) - 
                                                   double(spongeStart));
                  double eta = spongeDistance*spongeScale;
                  double spongeFac = spongeAmp*std::pow(eta,spongePower);

                    // rhs = rhs - spongeFac*(state - target)
                  for(int iEquation = 0;iEquation < numEquations;iEquation++){
                    rhsBuffers[iEquation][bufferIndex] = rhsBuffers[iEquation][bufferIndex] -
                      spongeFac*(stateBuffers[iEquation][bufferIndex] - 
                                 targetBuffers[iEquation][bufferIndex]);
                  }
                }
              }
            } else if(numDim == 1){
              
              size_t iPartStart = partitionOpInterval[0].first;
              size_t iPartEnd   = partitionOpInterval[0].second;
              
              size_t iBufferStart = opBufferInterval[0].first;
              size_t iBufferEnd   = opBufferInterval[0].second;
              
              size_t iBufferIndex;
              
              for(size_t partIndex = iPartStart,iBufferIndex = iBufferStart;
                  partIndex <= iPartEnd;
                  partIndex++,iBufferIndex++){
                
                double spongeDistance = std::abs(double(partIndex) - double(spongeStart));
                double eta = spongeDistance*spongeScale;
                double spongeFac = spongeAmp*std::pow(eta,spongePower);
                
                // rhs = rhs - spongeFac*(state - target)
                for(int iEquation = 0;iEquation < numEquations;iEquation++){
                  
                  rhsBuffers[iEquation][iBufferIndex] = rhsBuffers[iEquation][iBufferIndex] -
                    spongeFac*(stateBuffers[iEquation][iBufferIndex] - 
                               targetBuffers[iEquation][iBufferIndex]);
                }
              }
            }
            return;
          } // LinearSponge BC


          void SpongeZone(const std::vector<size_t>     &bufferSizes,int normalDir,
                          const pcpp::IndexIntervalType &spongeRegion,
                          const pcpp::IndexIntervalType &partitionOpInterval,
                          const pcpp::IndexIntervalType &opBufferInterval,
                          const double *bcParams,const int *bcFlags,
                          const int *iMask,const int disabledMask,
                          const int numEquations,
                          const std::vector<double *>   &stateBuffers,
                          const std::vector<double *>   &targetBuffers,
                          std::vector<double *>         &rhsBuffers){
            
            int numDim = bufferSizes.size();
            double spongeAmp   = bcParams[0];
            double spongeWidth = bcParams[1];
            double spongeScale = 1.0/spongeWidth;
            
//             std::cout << "SpongeAmp   = " << spongeAmp   << std::endl
//                       << "SpongeWidth = " << spongeWidth << std::endl;

//            std::cout << "In sponge zone " << std::endl;
            if(numDim == 3){
              
              double spongeSize[] = {
                (double(spongeRegion[0].second)-double(spongeRegion[0].first)+1),
                (double(spongeRegion[1].second)-double(spongeRegion[1].first)+1),
                (double(spongeRegion[2].second)-double(spongeRegion[2].first)+1)
              };

              double spongeCenter[] = {
                double(spongeRegion[0].first)+(spongeSize[0]-1)/2.0,
                double(spongeRegion[1].first)+(spongeSize[1]-1)/2.0,
                double(spongeRegion[2].first)+(spongeSize[2]-1)/2.0
              };
              
              size_t kPartStart = partitionOpInterval[2].first;
              size_t kPartEnd   = partitionOpInterval[2].second;
              size_t jPartStart = partitionOpInterval[1].first;
              size_t jPartEnd   = partitionOpInterval[1].second;
              size_t iPartStart = partitionOpInterval[0].first;
              size_t iPartEnd   = partitionOpInterval[0].second;
              
              size_t kBufferStart = opBufferInterval[2].first;
              size_t kBufferEnd   = opBufferInterval[2].second;
              size_t jBufferStart = opBufferInterval[1].first;
              size_t jBufferEnd   = opBufferInterval[1].second;
              size_t iBufferStart = opBufferInterval[0].first;
              size_t iBufferEnd   = opBufferInterval[0].second;
              
              size_t nPlane = bufferSizes[0]*bufferSizes[1];
              size_t kBufferIndex,jBufferIndex,iBufferIndex;
              
              size_t partIndex[] = {0,0,0};
              for(partIndex[2] = kPartStart,kBufferIndex = kBufferStart;
                  partIndex[2] <= kPartEnd;
                  partIndex[2]++,kBufferIndex++){
                size_t zIndex = kBufferIndex*nPlane;
                double spongeDistanceZ  = double(partIndex[2])-spongeCenter[2];
                double spongeDistanceZ2 = spongeDistanceZ*spongeDistanceZ;
                spongeDistanceZ  = std::sqrt(spongeDistanceZ2);
                for(partIndex[1] = jPartStart,jBufferIndex = jBufferStart;
                    partIndex[1] <= jPartEnd;
                    partIndex[1]++,jBufferIndex++){
                  size_t yzIndex = jBufferIndex*bufferSizes[0] + zIndex;
                  double spongeDistanceY  = double(partIndex[1])-spongeCenter[1];
                  double spongeDistanceY2 = spongeDistanceY*spongeDistanceY;
                  spongeDistanceY  = std::sqrt(spongeDistanceY2);
                  for(partIndex[0] = iPartStart,iBufferIndex = iBufferStart;
                      partIndex[0] <= iPartEnd;
                      partIndex[0]++,iBufferIndex++){
                    size_t bufferIndex = yzIndex + iBufferIndex;
                    
                    if(iMask[bufferIndex]&disabledMask)
                      continue;
                    
                    double spongeDistanceX  = double(partIndex[0])-spongeCenter[0];
                    double spongeDistanceX2 = spongeDistanceX*spongeDistanceX;
                    
                    double spongeDistance   = spongeDistanceX2+
                      spongeDistanceY2+ 
                      spongeDistanceZ2;
                    
                    double spongeDistance2  = -1.0*spongeDistance*spongeScale;
                    double spongeFac = spongeAmp*std::exp(spongeDistance2);
                    
                    // rhs = rhs - spongeFac*(state - target)
                    for(int iEquation = 0;iEquation < numEquations;iEquation++){
                      rhsBuffers[iEquation][bufferIndex] = rhsBuffers[iEquation][bufferIndex] -
                        spongeFac*(stateBuffers[iEquation][bufferIndex] - 
                                   targetBuffers[iEquation][bufferIndex]);
                    }
                  }
                }
              }
            } else if(numDim == 2){
              double spongeSize[] = {
                (double(spongeRegion[0].second)-double(spongeRegion[0].first)+1),
                (double(spongeRegion[1].second)-double(spongeRegion[1].first)+1)
              };
              if(bcFlags[1] > 0){
                spongeSize[0] = bcFlags[1];
              }
              if(bcFlags[3] > 0){
                spongeSize[1] = bcFlags[3];
              }
              double spongeCenter[] = {
                double(spongeRegion[0].first)+(spongeSize[0]-1)/2.0,
                double(spongeRegion[1].first)+(spongeSize[1]-1)/2.0
              };
              if(bcFlags[0] > 0)
                spongeCenter[0] = bcFlags[0] - 1;
              if(bcFlags[2] > 0)
                spongeCenter[1] = bcFlags[2] - 1;

//               std::cout << "Sponge: (" << spongeCenter[0] << "," 
//                         << spongeCenter[1] << ") (" << spongeSize[0] << "," 
//                         << spongeSize[1] << ")" << std::endl;
              
              size_t jPartStart = partitionOpInterval[1].first;
              size_t jPartEnd   = partitionOpInterval[1].second;
              size_t iPartStart = partitionOpInterval[0].first;
              size_t iPartEnd   = partitionOpInterval[0].second;
              
              size_t jBufferStart = opBufferInterval[1].first;
              size_t jBufferEnd   = opBufferInterval[1].second;
              size_t iBufferStart = opBufferInterval[0].first;
              size_t iBufferEnd   = opBufferInterval[0].second;
              
              size_t jBufferIndex,iBufferIndex;
              
              size_t partIndex[] = {0,0};
              for(partIndex[1] = jPartStart, jBufferIndex = jBufferStart;
                  partIndex[1] <= jPartEnd;
                  partIndex[1]++,jBufferIndex++){
                size_t yIndex = jBufferIndex*bufferSizes[0];
                double spongeDistanceY  = double(partIndex[1])-spongeCenter[1];
                double spongeDistanceY2 = spongeDistanceY*spongeDistanceY;
                spongeDistanceY  = std::sqrt(spongeDistanceY2);
                for(partIndex[0] = iPartStart,iBufferIndex = iBufferStart;
                    partIndex[0] <= iPartEnd;
                    partIndex[0]++,iBufferIndex++){
                  size_t bufferIndex = yIndex + iBufferIndex;
                  
                  if(iMask[bufferIndex]&disabledMask)
                    continue;
                  
                  double spongeDistanceX  = double(partIndex[0])-spongeCenter[0];
                  double spongeDistanceX2 = spongeDistanceX*spongeDistanceX;
                  double spongeDistance   = spongeDistanceX2 + 
                    spongeDistanceY2;
                  double spongeDistance2  = -1.0*spongeDistance*spongeScale;
                  double spongeFac = spongeAmp*std::exp(spongeDistance2);
                  spongeDistance = std::sqrt(spongeDistance);
                  
                  for(int iEquation = 0;iEquation < numEquations;iEquation++){
                    rhsBuffers[iEquation][bufferIndex] = rhsBuffers[iEquation][bufferIndex] -
                      spongeFac*(stateBuffers[iEquation][bufferIndex] - 
                                 targetBuffers[iEquation][bufferIndex]);
                    
                  }
                }
              }
            } else if(numDim == 1){
              
              double spongeSize[] = {
                (double(spongeRegion[0].second)-double(spongeRegion[0].first)+1),
              };

              double spongeCenter[] = {
                (double(spongeRegion[0].first)+(spongeSize[0]-1))/2.0,
              };
              
              size_t iPartStart = partitionOpInterval[0].first;
              size_t iPartEnd   = partitionOpInterval[0].second;
              
              size_t iBufferStart = opBufferInterval[0].first;
              size_t iBufferEnd   = opBufferInterval[0].second;
              
              size_t iBufferIndex;
              
              for(size_t partIndex = iPartStart,iBufferIndex = iBufferStart;
                  partIndex <= iPartEnd;
                  partIndex++,iBufferIndex++){
                
                
                if(iMask[iBufferIndex]&disabledMask)
                  continue;
                
                double spongeDistanceX  = double(partIndex)-spongeCenter[0];
                double spongeDistanceX2 = spongeDistanceX*spongeDistanceX;
                double spongeDistance   = std::sqrt(spongeDistanceX2);
                double spongeDistance2  = -1.0*spongeDistance*spongeDistance*spongeScale;
                double spongeFac = spongeAmp*std::exp(spongeDistance2);
                
                // rhs = rhs - spongeFac*(state - target)
                for(int iEquation = 0;iEquation < numEquations;iEquation++){
                  rhsBuffers[iEquation][iBufferIndex] = rhsBuffers[iEquation][iBufferIndex] -
                    spongeFac*(stateBuffers[iEquation][iBufferIndex] - 
                               targetBuffers[iEquation][iBufferIndex]);
                }
              }
            }
            return;
          } // SpongeZone BC
          
          // Namespaces close
        }
      }
    }
  }
}