WENO.H

Go to the documentation of this file.

#ifndef __WENO_H__
#define __WENO_H__

#include "Stencil.H"

#include <cmath>
#include <algorithm>
#include <vector>

namespace WENO {

  typedef plascom2::operators::stencilset StencilSet;

  class CoeffsWENO {

    // Private data
    std::vector<StencilSet> reconstCoeffs;        // StencilSets with reconstruction coefficients
    std::vector<StencilSet> siCoeffsFirst;        // StencilSets with first set of smoothness indicator coefficients
    std::vector<StencilSet> siCoeffsSecond;       // StencilSets with second set of smoothness indicator coefficients
    std::vector<std::vector<double> > linWeights; // linear weights for non-linear combination step

    std::vector<int> nSten;                       // number of stencils in each group
    std::vector<int> loGroup;                     // smallest offset from center point in group
    std::vector<std::vector<int> > loOffset;      // smallest offset from center point in each stencil

    public:

    // Public data
    int loAll;        // smallest offset from center point in group of stencils
    int hiAll;        // largest offset from center point in group of stencils
    int nGroup;       // number of groups

    // coefficients for WENO combination
    double siFirst;
    double siSecond;
    double epsW;
    double pW;

    // Constructor
    CoeffsWENO() {

      nGroup = 2;
      // Group 0 = negative upwinding
      // Group 1 = positive upwinding

      loOffset.resize(nGroup);

      // Smoothness indicator constant coefficients
      siFirst = 13/12.0;
      siSecond = 1/4.0;

      // Non-linear combination coefficients
      epsW = 1e-6;
      pW = 2;
      linWeights.resize(nGroup);

      // fill vectors with blank stencilsets
      for (int g=0; g<nGroup; g++) {
        StencilSet ss;
        reconstCoeffs.push_back(ss);
      }
      for (int g=0; g<nGroup; g++) {
        StencilSet ss;
        siCoeffsFirst.push_back(ss);
      }
      for (int g=0; g<nGroup; g++) {
        StencilSet ss;
        siCoeffsSecond.push_back(ss);
      }

      StencilSet temp;  // temporary StencilSet for adding values
      temp.ownData = true;

      // Reconstruction coefficients -- group 0
      temp.numStencils = 3;
      temp.numValues = 9;
      temp.stencilSizes = new int[temp.numStencils];
      temp.stencilStarts = new int[temp.numStencils];
      temp.stencilOffsets = new int[temp.numValues];
      temp.stencilWeights = new double[temp.numValues];

      temp.stencilSizes[0] = 3;
      temp.stencilSizes[1] = 3;
      temp.stencilSizes[2] = 3;

      temp.stencilStarts[0] = 0;
      for (int i=1; i<temp.numStencils; i++)
        temp.stencilStarts[i] = temp.stencilStarts[i-1] + temp.stencilSizes[i-1];

      temp.stencilOffsets = new int[temp.numValues];
      temp.stencilOffsets[0] = 3;
      temp.stencilOffsets[1] = 2;
      temp.stencilOffsets[2] = 1;
      temp.stencilOffsets[3] = 2;
      temp.stencilOffsets[4] = 1;
      temp.stencilOffsets[5] = 0;
      temp.stencilOffsets[6] = 1;
      temp.stencilOffsets[7] = 0;
      temp.stencilOffsets[8] = -1;

      loAll = temp.stencilOffsets[0];
      hiAll = temp.stencilOffsets[0];
      loGroup.push_back(temp.stencilOffsets[0]);
      for (int i=0; i<temp.numValues; i++) {
        if (loGroup[0] > temp.stencilOffsets[i]) loGroup[0] = temp.stencilOffsets[i];
        if (loAll > temp.stencilOffsets[i]) loAll = temp.stencilOffsets[i];
        if (hiAll < temp.stencilOffsets[i]) hiAll = temp.stencilOffsets[i];
      }
      for (int s=0; s<temp.numStencils; s++) {
        loOffset[0].push_back(temp.stencilOffsets[temp.stencilStarts[s]]);
        for (int i=temp.stencilStarts[s]+1; i<temp.stencilStarts[s]+temp.stencilSizes[s]; i++) {
          if (loOffset[0][s] > temp.stencilOffsets[i]) loOffset[0][s] = temp.stencilOffsets[i];
        }
      }

      temp.stencilWeights[0] = 1/3.0;
      temp.stencilWeights[1] = -7/6.0;
      temp.stencilWeights[2] = 11/6.0;
      temp.stencilWeights[3] = -1/6.0;
      temp.stencilWeights[4] = 5/6.0;
      temp.stencilWeights[5] = 1/3.0;
      temp.stencilWeights[6] = 1/3.0;
      temp.stencilWeights[7] = 5/6.0;
      temp.stencilWeights[8] = -1/6.0;

      reconstCoeffs[0].Copy(temp);
      nSten.push_back(temp.numStencils);
      linWeights[0].push_back(0.1);
      linWeights[0].push_back(0.6);
      linWeights[0].push_back(0.3);

      // Reconstruction coefficients -- group 1
      // Note that group 1 is symmetric to group 0 about i+1/2, so
      //  sizes/starts/weights do not change
      temp.stencilOffsets[0] = -2;
      temp.stencilOffsets[1] = -1;
      temp.stencilOffsets[2] = 0;
      temp.stencilOffsets[3] = -1;
      temp.stencilOffsets[4] = 0;
      temp.stencilOffsets[5] = 1;
      temp.stencilOffsets[6] = 0;
      temp.stencilOffsets[7] = 1;
      temp.stencilOffsets[8] = 2;

      loGroup.push_back(temp.stencilOffsets[0]);
      for (int i=0; i<temp.numValues; i++) {
        if (loGroup[1] > temp.stencilOffsets[i]) loGroup[1] = temp.stencilOffsets[i];
        if (loAll > temp.stencilOffsets[i]) loAll = temp.stencilOffsets[i];
        if (hiAll < temp.stencilOffsets[i]) hiAll = temp.stencilOffsets[i];
      }
      for (int s=0; s<temp.numStencils; s++) {
        loOffset[1].push_back(temp.stencilOffsets[temp.stencilStarts[s]]);
        for (int i=temp.stencilStarts[s]+1; i<temp.stencilStarts[s]+temp.stencilSizes[s]; i++) {
          if (loOffset[1][s] > temp.stencilOffsets[i]) loOffset[1][s] = temp.stencilOffsets[i];
        }
      }

      reconstCoeffs[1].Copy(temp);
      nSten.push_back(temp.numStencils);
      linWeights[1].push_back(0.1);
      linWeights[1].push_back(0.6);
      linWeights[1].push_back(0.3);

      // First set of smoothness indicator coefficients -- group 0
      // Same sizes as reconstruction coefficients, but since we just did
      //  group 1, need to enter the offsets again
      temp.stencilOffsets[0] = 3;
      temp.stencilOffsets[1] = 2;
      temp.stencilOffsets[2] = 1;
      temp.stencilOffsets[3] = 2;
      temp.stencilOffsets[4] = 1;
      temp.stencilOffsets[5] = 0;
      temp.stencilOffsets[6] = 1;
      temp.stencilOffsets[7] = 0;
      temp.stencilOffsets[8] = -1;

      temp.stencilWeights[0] = 1;
      temp.stencilWeights[1] = -2;
      temp.stencilWeights[2] = 1;
      temp.stencilWeights[3] = 1;
      temp.stencilWeights[4] = -2;
      temp.stencilWeights[5] = 1;
      temp.stencilWeights[6] = 1;
      temp.stencilWeights[7] = -2;
      temp.stencilWeights[8] = 1;

      siCoeffsFirst[0].Copy(temp);

      // First set of smoothness indicator coefficients -- group 1
      // Note that group 1 is symmetric to group 0 about i+1/2, so
      //  sizes/starts/weights do not change
      temp.stencilOffsets[0] = -2;
      temp.stencilOffsets[1] = -1;
      temp.stencilOffsets[2] = 0;
      temp.stencilOffsets[3] = -1;
      temp.stencilOffsets[4] = 0;
      temp.stencilOffsets[5] = 1;
      temp.stencilOffsets[6] = 0;
      temp.stencilOffsets[7] = 1;
      temp.stencilOffsets[8] = 2;

      siCoeffsFirst[1].Copy(temp);

      // Second set of smoothness indicator coefficients -- group 0
      temp.numStencils = 3;
      temp.numValues = 8;
      temp.stencilSizes = new int[temp.numStencils];
      temp.stencilStarts = new int[temp.numStencils];
      temp.stencilOffsets = new int[temp.numValues];
      temp.stencilWeights = new double[temp.numValues];

      temp.stencilSizes[0] = 3;
      temp.stencilSizes[1] = 2;
      temp.stencilSizes[2] = 3;

      temp.stencilStarts[0] = 0;
      for (int i=1; i<temp.numStencils; i++)
        temp.stencilStarts[i] = temp.stencilStarts[i-1] + temp.stencilSizes[i-1];

      temp.stencilOffsets = new int[temp.numValues];
      temp.stencilOffsets[0] = 3;
      temp.stencilOffsets[1] = 2;
      temp.stencilOffsets[2] = 1;
      temp.stencilOffsets[3] = 2;
      temp.stencilOffsets[4] = 0;
      temp.stencilOffsets[5] = 1;
      temp.stencilOffsets[6] = 0;
      temp.stencilOffsets[7] = -1;

      temp.stencilWeights[0] = 1;
      temp.stencilWeights[1] = -4;
      temp.stencilWeights[2] = 3;
      temp.stencilWeights[3] = 1;
      temp.stencilWeights[4] = -1;
      temp.stencilWeights[5] = 3;
      temp.stencilWeights[6] = -4;
      temp.stencilWeights[7] = 1;

      siCoeffsSecond[0].Copy(temp);

      // Second set of smoothness indicator coefficients -- group 1
      // Note that group 1 is symmetric to group 0 about i+1/2, so
      //  sizes/starts/weights do not change
      temp.stencilOffsets[0] = -2;
      temp.stencilOffsets[1] = -1;
      temp.stencilOffsets[2] = 0;
      temp.stencilOffsets[3] = -1;
      temp.stencilOffsets[4] = 1;
      temp.stencilOffsets[5] = 0;
      temp.stencilOffsets[6] = 1;
      temp.stencilOffsets[7] = 2;

      siCoeffsSecond[1].Copy(temp);
    }

    // Get necessary width of halo region
    int HaloWidth() {
      // WENO computes results at face i+1/2 for stencils centered on cell i,
      //  so the +1 below handles the case where WENO needs to compute the
      //  left face of the leftmost cell
      return std::max(std::abs(loAll)+1, std::abs(hiAll));
    }

    int NSten(int group) {
      return nSten[group];
    }

    int LoGroup(int group) {
      return loGroup[group];
    }

    int Lo(int group, int sten) {
      return loOffset[group][sten];
    }

    double LinWeight(int group, int sten) {
      return linWeights[group][sten];
    }

    StencilSet Reconst(int group) {
      StencilSet res;
      res.Copy(reconstCoeffs[group]);
      return res;
    }

    StencilSet SmIndFirst(int group) {
      StencilSet res;
      res.Copy(siCoeffsFirst[group]);
      return res;
    }

    StencilSet SmIndSecond(int group) {
      StencilSet res;
      res.Copy(siCoeffsSecond[group]);
      return res;
    }
  };

  void ReconstPointVal(double vals[], CoeffsWENO& coeffs, int group, double& pVal);
  void ReconstPointValSten(double vals[], CoeffsWENO& coeffs, int group, int sten, double& pVal);
  void SmoothInd(double vals[], CoeffsWENO& coeffs, int group, int sten, double& si);
  void Project(int n, int m, double T[], double vals[], double res[]);
  double EntropyFixEta(int n, double lambdaL[], double lambdaR[]);
}

#endif