WENO.C

Go to the documentation of this file.

// Functions for WENO schemes
//
// Source for basic scheme is:
// Jiang and Shu, J. Comput. Phys. 126 (1996) 202
// Deviations noted individually

#include "WENO.H"

namespace WENO {

  void ReconstPointVal(double vals[], CoeffsWENO& coeffs, int group, double& pVal) {
    // Computes point value reconstructions for an array of values on all
    //  stencils for a given group
    //
    // Inputs/Outputs
    //  vals    -- Field/flux values, index 0 corresponds to low offset for this group
    //  coeffs  -- WENO coefficients
    //  group   -- WENO stencil group to use
    //  pVal    -- Reconstructed value

    int nSten = coeffs.NSten(group);
    double pValSten[nSten];
    double siSten[nSten];
    for (int s=0; s<nSten; s++) {
      int i = coeffs.Lo(group, s) - coeffs.LoGroup(group);
      ReconstPointValSten(&vals[i], coeffs, group, s, pValSten[s]);
      SmoothInd(&vals[i], coeffs, group, s, siSten[s]);
    }

    double nlWeights[nSten];
    double nlWeightsSum = 0;
    for (int s=0; s<nSten; s++) {
      nlWeights[s] = coeffs.LinWeight(group, s)/std::pow(coeffs.epsW + siSten[s], coeffs.pW);
      nlWeightsSum += nlWeights[s];
    }

    pVal = 0;
    for (int s=0; s<nSten; s++) {
      pVal += nlWeights[s]*pValSten[s]/nlWeightsSum;
    }
  }

  void ReconstPointValSten(double vals[], CoeffsWENO& coeffs, int group, int sten, double& pVal) {
    // Computes point value reconstructions for an array of values on a single stencil
    //
    // Inputs/Outputs
    //  vals    -- Field/flux values, index 0 corresponds to low offset for this stencil
    //  coeffs  -- WENO coefficients
    //  group   -- WENO stencil group to use
    //  sten    -- Stencil to use within the group
    //  pVal    -- Reconstructed value

    StencilSet reconst = coeffs.Reconst(group);
    int lo = coeffs.Lo(group, sten);

    pVal = 0;
    for (int j=reconst.stencilStarts[sten]; j<reconst.stencilStarts[sten]+reconst.stencilSizes[sten]; j++) {
      pVal = pVal + reconst.stencilWeights[j]*vals[reconst.stencilOffsets[j]-lo];
    }
  }

  void SmoothInd(double vals[], CoeffsWENO& coeffs, int group, int sten, double& si) {
    // Computes smoothness indicator for an array of values on a single stencil
    //
    // Inputs/Outputs
    //  vals    -- Field/flux values, index 0 corresponds to low offset for this stencil
    //  coeffs  -- WENO coefficients
    //  group   -- WENO stencil group to use
    //  sten    -- Stencil to use within the group
    //  si    -- Smoothness indicator

    int lo = coeffs.Lo(group, sten);

    StencilSet si1 = coeffs.SmIndFirst(group);
    double comb = 0;
    for (int j=si1.stencilStarts[sten]; j<si1.stencilStarts[sten]+si1.stencilSizes[sten]; j++) {
      comb = comb + si1.stencilWeights[j]*vals[si1.stencilOffsets[j]-lo];
    }
    si = coeffs.siFirst*comb*comb;

    StencilSet si2 = coeffs.SmIndSecond(group);
    comb = 0;
    for (int j=si2.stencilStarts[sten]; j<si2.stencilStarts[sten]+si2.stencilSizes[sten]; j++) {
      comb = comb + si2.stencilWeights[j]*vals[si2.stencilOffsets[j]-lo];
    }
    si = si + coeffs.siSecond*comb*comb;
  }

  void Project(int n, int m, double T[], double vals[], double res[]) {
    // Computes characteristic projection for a given set of values
    //  with a given rotation matrix
    //
    // Inputs/Outputs
    //  n       -- size of rotation matrix
    //  m       -- number of columns in field/flux value matrix
    //  T       -- n x n rotation matrix (column-major ordering)
    //  vals    -- n x m matrix of field/flux values (row-major ordering)
    //  res     -- n x m result matrix (row-major ordering)

    for (int i=0; i<n; i++) {
      for (int j=0; j<m; j++) {
        int rIndex = i*m + j;   // (i,j) row-major
        res[rIndex] = 0;
        for (int k=0; k<n; k++) {
          int tIndex = i + k*n;   // (i,k) col-major
          int vIndex = k*m + j;   // (k,j) row-major
          res[rIndex] += T[tIndex]*vals[vIndex];
        }
      }
    }
  }

  double EntropyFixEta(int n, double lambdaL[], double lambdaR[]) {
    // Computes entropy fix threshold eta
    //
    // Inputs
    //  n       -- number of eigenvalues in problem
    //  lambdaL -- eigenvalues of left state (n x n diagonal matrix)
    //  lambdaR -- eigenvalues of right state (n x n diagonal matrix)

    double eta = 0;
    for (int i=0; i<n; i++) {
      int diag = i*(n+1);
      eta = std::max(eta, std::abs(lambdaR[diag] - lambdaL[diag]));
    }
    eta *= 0.5;

    return eta;
  }
}