TestWENOSerial.C

Go to the documentation of this file.

// Serial unit tests for WENO routines and classes

#include "Testing.H"
#include "Simulation.H"
#include "RKTestFixtures.H"
#include "EulerRHS.H"
#include "Stencil.H"
#include "PCPPCommUtil.H"
#include "PCPPReport.H"
#include "PCPPIntervalUtils.H"
#include "TestFixtures.H"
#include "EulerTestFixtures.H"
#include "SATKernels.H"
#include "WENO.H"


void TestWENO_Stencils(ix::test::results &serialUnitResults) {
  // Tests WENO stencil collection classes/objects

  bool result;
  double eps = 1e-8;

  WENO::CoeffsWENO coeffs;

  // --- nGroup ---
  int nGroup = 2;

  result = (coeffs.nGroup == nGroup);
  serialUnitResults.UpdateResult("WENO:Stencils:nGroup", result);

  // --- NSten ---
  int nSten = 3;   // both groups have 3 stencils

  result = true;
  for (int g=0; g<nGroup; g++) {
    result = result && (coeffs.NSten(g) == nSten);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:nGroup", result);

  // --- siFirst ---
  double siFirst = 13/12.0;

  result = (std::abs(coeffs.siFirst - siFirst) < eps);
  serialUnitResults.UpdateResult("WENO:Stencils:siFirst", result);

  // --- siSecond ---
  double siSecond = 1/4.0;

  result = (std::abs(coeffs.siSecond - siSecond) < eps);
  serialUnitResults.UpdateResult("WENO:Stencils:siSecond", result);

  // --- epsW ---
  double epsW = 1e-6;

  result = (std::abs(coeffs.epsW - epsW) < eps);
  serialUnitResults.UpdateResult("WENO:Stencils:epsW", result);

  // --- pW ---
  double pW = 2;

  result = (std::abs(coeffs.pW - pW) < eps);
  serialUnitResults.UpdateResult("WENO:Stencils:pW", result);

  // --- linWeights ---
  double linWeights[2][3] = {{0.1, 0.6, 0.3}, {0.1, 0.6, 0.3}};

  result = true;
  for (int g=0; g<nGroup; g++) {
    for (int i=0; i<nSten; i++) {
      result = (std::abs(coeffs.LinWeight(g, i) - linWeights[g][i]) < eps);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:linWeights", result);

  // --- HaloWidth ---
  int loAll = -2;
  int hiAll = 3;
  int expHaloWidth = 3;  // max(abs(loAll) + 1, abs(hiAll)), +1 is because we may need the flux on the left face of the leftmost cell

  result = (coeffs.HaloWidth() == expHaloWidth);
  serialUnitResults.UpdateResult("WENO:Stencils:HaloWidth", result);

  // --- LoGroup ---
  // Test that coeffs.LoGroup returns correct lo offset for each stencil group
  // For now, 2 groups
  // Group 0 = negative upwinding
  // Group 1 = positive upwinding
  int expLoGroup[2] = {-1, -2};

  result = true;
  for (int g=0; g<nGroup; g++) {
    result = result && (coeffs.LoGroup(g) == expLoGroup[g]);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:LoGroup", result);

  // --- Lo ---
  // Test that coeffs.Lo returns correct lo offset for each stencil
  // Assumes stencils are entered in a prespecified order
  // For now, 2 groups, 3 stencils in each
  // Group 0 = negative upwinding
  // Group 1 = positive upwinding
  int expLo[2][3] = {{1, 0, -1}, {-2, -1, 0}};

  result = true;
  for (int g=0; g<nGroup; g++) {
    for (int s=0; s<nSten; s++) {
      result = result && (coeffs.Lo(g, s) == expLo[g][s]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:Lo", result);

  // --- Reconst ---
  // Test that coeffs.Reconst returns stencilset with reconstruction coefficients
  // For now, 2 groups

  // Stencils (expected values)
  // Group 0 = negative upwinding
  // Group 1 = positive upwinding
  int nVal = 9;
  int rSizes[2][3] = {{3, 3, 3}, {3, 3, 3}};
  int rStarts[2][3] = {{0, 3, 6}, {0, 3, 6}};
  int rOffsets[2][9] = {{3, 2, 1, 2, 1, 0, 1, 0, -1}, {-2, -1, 0, -1, 0, 1, 0, 1, 2}};
  double rWeights[2][9] = {{1/3.0, -7/6.0, 11/6.0, -1/6.0, 5/6.0, 1/3.0, 1/3.0, 5/6.0, -1/6.0},
                           {1/3.0, -7/6.0, 11/6.0, -1/6.0, 5/6.0, 1/3.0, 1/3.0, 5/6.0, -1/6.0}};

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.Reconst(g);

    result = result && (test.numStencils == nSten);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:Reconst:nSten", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.Reconst(g);

    result = result && (test.numValues == nVal);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:Reconst:nVal", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.Reconst(g);

    for (int i=0; i<nSten; i++) {
      result = result && (test.stencilSizes[i] == rSizes[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:Reconst:Sizes", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.Reconst(g);

    for (int i=0; i<nSten; i++) {
      result = result && (test.stencilStarts[i] == rStarts[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:Reconst:Starts", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.Reconst(g);

    for (int i=0; i<nVal; i++) {
      //std::cout << "Stencils:Reconst:Offsets -- g = " << g << ", i = " << i << ", Actual   = " << test.stencilOffsets[i] << std::endl;
      //std::cout << "Stencils:Reconst:Offsets -- g = " << g << ", i = " << i << ", Expected = " << rOffsets[g][i] << std::endl;
      result = result && (test.stencilOffsets[i] == rOffsets[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:Reconst:Offsets", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.Reconst(g);

    for (int i=0; i<nVal; i++) {
      result = result && (std::abs(test.stencilWeights[i] - rWeights[g][i]) < eps);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:Reconst:Weights", result);

  // --- SmIndFirst ---
  // Test that coeffs.SmIndFirst returns stencilset with first set of smoothness indicator coefficients
  // For now, 2 groups

  // Stencils (expected values)
  // Group 0 = negative upwinding
  // Group 1 = positive upwinding
  int si1Sizes[2][3] = {{3, 3, 3}, {3, 3, 3}};
  int si1Starts[2][3] = {{0, 3, 6}, {0, 3, 6}};
  int si1Offsets[2][9] = {{3, 2, 1, 2, 1, 0, 1, 0, -1}, {-2, -1, 0, -1, 0, 1, 0, 1, 2}};
  double si1Weights[2][9] = {{1, -2, 1, 1, -2, 1, 1, -2, 1}, {1, -2, 1, 1, -2, 1, 1, -2, 1}};

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndFirst(g);

    result = result && (test.numStencils == nSten);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndFirst:nSten", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndFirst(g);

    result = result && (test.numValues == nVal);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndFirst:nVal", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndFirst(g);

    for (int i=0; i<nSten; i++) {
      result = result && (test.stencilSizes[i] == si1Sizes[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndFirst:Sizes", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndFirst(g);

    for (int i=0; i<nSten; i++) {
      result = result && (test.stencilStarts[i] == si1Starts[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndFirst:Starts", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndFirst(g);

    for (int i=0; i<nVal; i++) {
      //std::cout << "Stencils:SmIndFirst:Offsets -- g = " << g << ", i = " << i << ", Actual   = " << test.stencilOffsets[i] << std::endl;
      //std::cout << "Stencils:SmIndFirst:Offsets -- g = " << g << ", i = " << i << ", Expected = " << si1Offsets[g][i] << std::endl;
      result = result && (test.stencilOffsets[i] == si1Offsets[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndFirst:Offsets", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndFirst(g);

    for (int i=0; i<nVal; i++) {
      result = result && (std::abs(test.stencilWeights[i] - si1Weights[g][i]) < eps);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndFirst:Weights", result);

  // --- SmIndSecond ---
  // Test that coeffs.SmIndSecond returns stencilset with second set of smoothness indicator coefficients
  // For now, 2 groups

  // Stencils (expected values)
  // Group 0 = negative upwinding
  // Group 1 = positive upwinding
  nVal = 8;
  int si2Sizes[2][3] = {{3, 2, 3}, {3, 2, 3}};
  int si2Starts[2][3] = {{0, 3, 5}, {0, 3, 5}};
  int si2Offsets[2][8] = {{3, 2, 1, 2, 0, 1, 0, -1}, {-2, -1, 0, -1, 1, 0, 1, 2}};
  double si2Weights[2][8] = {{1, -4, 3, 1, -1, 3, -4, 1}, {1, -4, 3, 1, -1, 3, -4, 1}};

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndSecond(g);

    result = result && (test.numStencils == nSten);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndSecond:nSten", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndSecond(g);

    result = result && (test.numValues == nVal);
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndSecond:nVal", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndSecond(g);

    for (int i=0; i<nSten; i++) {
      result = result && (test.stencilSizes[i] == si2Sizes[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndSecond:Sizes", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndSecond(g);

    for (int i=0; i<nSten; i++) {
      result = result && (test.stencilStarts[i] == si2Starts[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndSecond:Starts", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndSecond(g);

    for (int i=0; i<nVal; i++) {
      //std::cout << "Stencils:SmIndSecond:Offsets -- g = " << g << ", i = " << i << ", Actual   = " << test.stencilOffsets[i] << std::endl;
      //std::cout << "Stencils:SmIndSecond:Offsets -- g = " << g << ", i = " << i << ", Expected = " << si2Offsets[g][i] << std::endl;
      result = result && (test.stencilOffsets[i] == si2Offsets[g][i]);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndSecond:Offsets", result);

  result = true;
  for (int g=0; g<nGroup; g++) {
    WENO::StencilSet test = coeffs.SmIndSecond(g);

    for (int i=0; i<nVal; i++) {
      result = result && (std::abs(test.stencilWeights[i] - si2Weights[g][i]) < eps);
    }
  }
  serialUnitResults.UpdateResult("WENO:Stencils:SmIndSecond:Weights", result);

}

void TestWENO_ReconstPointVal(ix::test::results &serialUnitResults) {
  // Tests for WENO point value reconstruction over a group of stencils
  // Assumes formulation given in
  //  Jiang and Shu, J. Comput. Phys. 126 (1996) 202

  bool result;
  double eps = 1e-8;
  //std::cout << "WENO:ReconstPointVal" << std::endl;

  WENO::CoeffsWENO cw;

  // --- Equal ---
  // Test for correct result if all values are equal
  // Apply once for each group, should be same on both

  int nGroup = 2;
  double vals[5] = {1, 1, 1, 1, 1};
  double pvalExpRes = 1.0;

  result = true;
  for (int g=0; g<nGroup; g++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    WENO::ReconstPointVal(vals, cw, g, pval);

    result = (std::abs(pval - pvalExpRes) < eps);
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:Equal", result);

  // --- Discontinuous ---
  // Test for correct result if there is a discontinuity
  // In all cases point values are reconstructed at face i+1/2, where central
  //  point is i
  // Test for a discontinuity at each face
  // Using group 1 -- positive upwinding

  // > Face-3/2 <
  // Discontinuity occurs at face i-3/2, where central point is i

  double valsD[5] = {2, 1, 1, 1, 1};
  double pvalExpResD = 1;

  double pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 1, pval);

  //std::cout << "Discontinuous:Face-3/2 -- Actual:   " << pval << std::endl;
  //std::cout << "Discontinuous:Face-3/2 -- Expected: " << pvalExpResD << std::endl;
  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:Discontinuous:Face-3/2", result);

  // > Face-1/2 <
  // Discontinuity occurs at face i-1/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 1;
  valsD[3] = 1;
  valsD[4] = 1;
  pvalExpResD = 1;

  pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 1, pval);

  //std::cout << "Discontinuous:Face-1/2 -- Actual:   " << pval << std::endl;
  //std::cout << "Discontinuous:Face-1/2 -- Expected: " << pvalExpResD << std::endl;
  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:Discontinuous:Face-1/2", result);

  // > Face+1/2 <
  // Discontinuity occurs at face i+1/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 1;
  valsD[4] = 1;
  pvalExpResD = 2;

  pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 1, pval);

  //std::cout << "Discontinuous:Face+1/2 -- Actual:   " << pval << std::endl;
  //std::cout << "Discontinuous:Face+1/2 -- Expected: " << pvalExpResD << std::endl;
  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:Discontinuous:Face+1/2", result);

  // > Face+3/2 <
  // Discontinuity occurs at face i+3/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 1;
  pvalExpResD = 2;

  pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 1, pval);

  //std::cout << "Discontinuous:Face+3/2 -- Actual:   " << pval << std::endl;
  //std::cout << "Discontinuous:Face+3/2 -- Expected: " << pvalExpResD << std::endl;
  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:Discontinuous:Face+3/2", result);

  // --- ReverseDiscont ---
  // Test for correct result if there is a discontinuity
  // In all cases point values are reconstructed at face i+1/2, where central
  //  point is i+1
  // Test for a discontinuity at each face
  // Using group 0 -- negative upwinding

  // > Face+5/2 <
  // Discontinuity occurs at face i+5/2, where central point is i

  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 1;
  valsD[3] = 1;
  valsD[4] = 2;
  pvalExpResD = 1;

  pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 0, pval);

  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:ReverseDiscont:Face+5/2", result);

  // > Face+3/2 <
  // Discontinuity occurs at face i+3/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 1;
  valsD[3] = 2;
  valsD[4] = 2;
  pvalExpResD = 1;

  pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 0, pval);

  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:ReverseDiscont:Face+3/2", result);

  // > Face+1/2 <
  // Discontinuity occurs at face i+1/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 2;
  pvalExpResD = 2;

  pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 0, pval);

  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:ReverseDiscont:Face+1/2", result);

  // > Face-1/2 <
  // Discontinuity occurs at face i-1/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 2;
  pvalExpResD = 2;

  pval = -1000; // we don't want test to pass because this is initialized to 0
  WENO::ReconstPointVal(valsD, cw, 0, pval);

  result = (std::abs(pval - pvalExpResD) < eps);
  serialUnitResults.UpdateResult("WENO:ReconstPointVal:ReverseDiscont:Face-1/2", result);
}

void TestWENO_ReconstPointValSten(ix::test::results &serialUnitResults) {
  // Tests for WENO point value reconstruction on a stencil
  // Assumes formulation given in
  //  Jiang and Shu, J. Comput. Phys. 126 (1996) 202

  bool result;
  double eps = 1e-8;
  //std::cout << "WENO:ReconstPointValSten" << std::endl;

  WENO::CoeffsWENO cw;

  // --- Equal ---
  // Test for correct result on each stencil if all values are equal
  // Test all stencil groups
  // Since stencils are all of length 3 and all values are equal, we can pass
  //  the same array of values for each stencil
  int nGroups = 2;
  int nSten = 3;
  double vals[3] = {1, 1, 1};
  double pValExpRes = 1;

  result = true;
  for (int g=0; g<nGroups; g++) {
    for (int i=0; i<nSten; i++) {
      double pVal = -1000; // we don't want test to pass because this is initialized to 0
      WENO::ReconstPointValSten(vals, cw, g, i, pVal);

      //std::cout << "ReconstPointValSten:Equal -- g = " << g << ", i = " << i << ", Actual   = " << pVal << std::endl;
      //std::cout << "ReconstPointValSten:Equal -- g = " << g << ", i = " << i << ", Expected = " << pValExpRes << std::endl;
      result = result && (std::abs(pVal - pValExpRes) < eps);
    }
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:Equal", result);

  // --- Discontinuous ---
  // Test for correct result on each stencil if there is a discontinuity
  // In all cases point values are reconstructed at face i+1/2, where central
  //  point is i
  // Test for a discontinuity at each face
  // For this test vals will be long enough to cover all stencils
  // Using group 1 -- positive upwinding

  // > Face-3/2 <
  // Discontinuity occurs at face i-3/2, where central point is i
  double valsD[5] = {2, 1, 1, 1, 1};
  double pvalsExpResD[3] = {4/3.0, 1, 1};

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[i], cw, 1, i, pval);

    //std::cout << "Discontinuous:Face-3/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "Discontinuous:Face-3/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:Discontinuous:Face-3/2", result);

  // > Face-1/2 <
  // Discontinuity occurs at face i-1/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 1;
  valsD[3] = 1;
  valsD[4] = 1;
  pvalsExpResD[0] = 1/6.0;
  pvalsExpResD[1] = 5/6.0;
  pvalsExpResD[2] = 1;

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[i], cw, 1, i, pval);

    //std::cout << "Discontinuous:Face-1/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "Discontinuous:Face-1/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:Discontinuous:Face-1/2", result);

  // > Face+1/2 <
  // Discontinuity occurs at face i+1/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 1;
  valsD[4] = 1;
  pvalsExpResD[0] = 2;
  pvalsExpResD[1] = 5/3.0;
  pvalsExpResD[2] = 4/3.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[i], cw, 1, i, pval);

    //std::cout << "Discontinuous:Face+1/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "Discontinuous:Face+1/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:Discontinuous:Face+1/2", result);

  // > Face+3/2 <
  // Discontinuity occurs at face i+3/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 1;
  pvalsExpResD[0] = 2;
  pvalsExpResD[1] = 2;
  pvalsExpResD[2] = 13/6.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[i], cw, 1, i, pval);

    //std::cout << "Discontinuous:Face+3/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "Discontinuous:Face+3/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:Discontinuous:Face+3/2", result);

  // --- ReverseDiscont ---
  // Same as Discontinuous above, but uses negative upwinding
  // Test for correct result on each stencil if there is a discontinuity
  // In all cases point values are reconstructed at face i+1/2, where central
  //  point is i+1
  // Test for a discontinuity at each face
  // For this test vals will be long enough to cover all stencils
  // Using group 0 -- negative upwinding

  // > Face+5/2 <
  // Discontinuity occurs at face i+5/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 1;
  valsD[3] = 1;
  valsD[4] = 2;
  pvalsExpResD[0] = 4/3.0;
  pvalsExpResD[1] = 1;
  pvalsExpResD[2] = 1;

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[2-i], cw, 0, i, pval);

    //std::cout << "ReverseDiscont:Face+5/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "ReverseDiscont:Face+5/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
    //std::cout << "ReverseDiscont:Face+5/2 (" << i << ") result:   " << result << std::endl;
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:ReverseDiscont:Face+5/2", result);

  // > Face+3/2 <
  // Discontinuity occurs at face i+3/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 1;
  valsD[3] = 2;
  valsD[4] = 2;
  pvalsExpResD[0] = 1/6.0;
  pvalsExpResD[1] = 5/6.0;
  pvalsExpResD[2] = 1;

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[2-i], cw, 0, i, pval);

    //std::cout << "ReverseDiscont:Face+3/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "ReverseDiscont:Face+3/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
    //std::cout << "ReverseDiscont:Face+3/2 (" << i << ") result:   " << result << std::endl;
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:ReverseDiscont:Face+3/2", result);

  // > Face+1/2 <
  // Discontinuity occurs at face i+1/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 2;
  pvalsExpResD[0] = 2;
  pvalsExpResD[1] = 5/3.0;
  pvalsExpResD[2] = 4/3.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[2-i], cw, 0, i, pval);

    //std::cout << "ReverseDiscont:Face+1/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "ReverseDiscont:Face+1/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
    //std::cout << "ReverseDiscont:Face+1/2 (" << i << ") result:   " << result << std::endl;
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:ReverseDiscont:Face+1/2", result);

  // > Face-1/2 <
  // Discontinuity occurs at face i-1/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 2;
  pvalsExpResD[0] = 2;
  pvalsExpResD[1] = 2;
  pvalsExpResD[2] = 13/6.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double pval = -1000; // we don't want test to pass because this is initialized to 0
    // note that using "i" to index valsD below takes advantage of knowledge about the stencils and their ordering
    WENO::ReconstPointValSten(&valsD[2-i], cw, 0, i, pval);

    //std::cout << "ReverseDiscont:Face-1/2 (" << i << ") Actual:   " << pval << std::endl;
    //std::cout << "ReverseDiscont:Face-1/2 (" << i << ") Expected: " << pvalsExpResD[i] << std::endl;
    result = result && (std::abs(pval - pvalsExpResD[i]) < eps);
    //std::cout << "ReverseDiscont:Face-1/2 (" << i << ") result:   " << result << std::endl;
  }
  serialUnitResults.UpdateResult("WENO:ReconstPointValSten:ReverseDiscont:Face-1/2", result);

}

void TestWENO_SmoothInd(ix::test::results &serialUnitResults) {
  // Tests for WENO smoothness indicator
  // Assumes formulation given in
  //  Jiang and Shu, J. Comput. Phys. 126 (1996) 202

  bool result;
  double eps = 1e-8;
  //std::cout << "WENO:SmoothInd" << std::endl;

  WENO::CoeffsWENO cw;

  // --- Equal ---
  // Test for correct result on each stencil if all values are equal
  // Since stencils are all of length 3 and all values are equal, we can pass
  //  the same array of values for each stencil
  int nGroups = 2;
  int nSten = 3;
  double vals[3] = {1, 1, 1};
  double siExpRes = 0;   // all results should be same

  result = true;
  for (int g=0; g<nGroups; g++) {
    for (int i=0; i<nSten; i++) {
      double si = -1000; // we don't want test to pass because this is initialized to 0
      WENO::SmoothInd(vals, cw, g, i, si);

      result = result && (std::abs(si - siExpRes) < eps);
    }
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:Equal", result);

  // --- Discontinuous ---
  // Test for correct result on each stencil if there is a discontinuity
  // Test for a discontinuity at each face
  // For this test vals will be long enough to cover all stencils
  // Using group 1 -- positive upwinding

  // > Face-3/2 <
  // Discontinuity occurs at face i-3/2, where central point is i
  double valsD[5] = {2, 1, 1, 1, 1};
  double siExpResD[3] = {4/3.0, 0, 0};

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[i], cw, 1, i, si);

    //std::cout << "Discontinuous:Face-3/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "Discontinuous:Face-3/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:Discontinuous:Face-3/2", result);

  // > Face-1/2 <
  // Discontinuity occurs at face i-1/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 1;
  valsD[3] = 1;
  valsD[4] = 1;
  siExpResD[0] = 10/3.0;
  siExpResD[1] = 4/3.0;
  siExpResD[2] = 0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[i], cw, 1, i, si);

    //std::cout << "Discontinuous:Face-1/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "Discontinuous:Face-1/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:Discontinuous:Face-1/2", result);

  // > Face+1/2 <
  // Discontinuity occurs at face i+1/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 1;
  valsD[4] = 1;
  siExpResD[0] = 0;
  siExpResD[1] = 4/3.0;
  siExpResD[2] = 10/3.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[i], cw, 1, i, si);

    //std::cout << "Discontinuous:Face+1/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "Discontinuous:Face+1/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:Discontinuous:Face+1/2", result);

  // > Face+3/2 <
  // Discontinuity occurs at face i+3/2, where central point is i
  valsD[0] = 2;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 1;
  siExpResD[0] = 0;
  siExpResD[1] = 0;
  siExpResD[2] = 4/3.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[i], cw, 1, i, si);

    //std::cout << "Discontinuous:Face+3/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "Discontinuous:Face+3/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:Discontinuous:Face+3/2", result);

  // --- ReverseDiscont ---
  // Test for correct result on each stencil if there is a discontinuity
  // Test for a discontinuity at each face
  // For this test vals will be long enough to cover all stencils
  // Using group 0 -- negative upwinding

  // > Face+5/2 <
  // Discontinuity occurs at face i+5/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 1;
  valsD[3] = 1;
  valsD[4] = 2;
  siExpResD[0] = 4/3.0;
  siExpResD[1] = 0;
  siExpResD[2] = 0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[2-i], cw, 0, i, si);

    //std::cout << "ReverseDiscont:Face+5/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "ReverseDiscont:Face+5/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:ReverseDiscont:Face+5/2", result);

  // > Face+3/2 <
  // Discontinuity occurs at face i+3/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 1;
  valsD[3] = 2;
  valsD[4] = 2;
  siExpResD[0] = 10/3.0;
  siExpResD[1] = 4/3.0;
  siExpResD[2] = 0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[2-i], cw, 0, i, si);

    //std::cout << "ReverseDiscont:Face+3/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "ReverseDiscont:Face+3/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:ReverseDiscont:Face+3/2", result);

  // > Face+1/2 <
  // Discontinuity occurs at face i+1/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 1;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 2;
  siExpResD[0] = 0;
  siExpResD[1] = 4/3.0;
  siExpResD[2] = 10/3.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[2-i], cw, 0, i, si);

    //std::cout << "ReverseDiscont:Face+1/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "ReverseDiscont:Face+1/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:ReverseDiscont:Face+1/2", result);

  // > Face-1/2 <
  // Discontinuity occurs at face i-1/2, where central point is i+1
  valsD[0] = 1;
  valsD[1] = 2;
  valsD[2] = 2;
  valsD[3] = 2;
  valsD[4] = 2;
  siExpResD[0] = 0;
  siExpResD[1] = 0;
  siExpResD[2] = 4/3.0;

  result = true;
  for (int i=0; i<nSten; i++) {
    double si = -1000; // we don't want test to pass because this is initialized to 0
    WENO::SmoothInd(&valsD[2-i], cw, 0, i, si);

    //std::cout << "Discontinuous:Face-1/2 (" << i << ") Actual:   " << si << std::endl;
    //std::cout << "Discontinuous:Face-1/2 (" << i << ") Expected: " << siExpResD[i] << std::endl;
    result = result && (std::abs(si - siExpResD[i]) < eps);
  }
  serialUnitResults.UpdateResult("WENO:SmoothInd:ReverseDiscont:Face-1/2", result);

}

void TestWENO_Project(ix::test::results &serialUnitResults) {
  // Tests characteristic projection functions

  bool result;
  double eps = 1e-8;

  // --- Mult ---
  // Test multiplication routines (Project)
  // First argument is 2D array stored in column-major order in a double[]
  // Second argument is 2D array stored in row-major order in a double[]
  // Result is 2D array stored in row-major order in a double[]

  // > MatMat <
  // Test matrix-matrix multiplication
  {
    double a[4] = {1, 3, 2, 4};
    double b[4] = {1, 2, 3, 4};
    double c[4] = {-1000, -1000, -1000, -1000};
    double cExp[4] = {7, 10, 15, 22};

    WENO::Project(2, 2, &a[0], &b[0], &c[0]);
    result = true;
    for (int i=0; i<4; i++) {
      result = result && (std::abs(c[i] - cExp[i]) < eps);
    }
    serialUnitResults.UpdateResult("WENO:Project:Mult:MatMat", result);
  }

  // > MatVec <
  // Test matrix-vector multiplication
  // Using same a matrix
  {
    double a[4] = {1, 3, 2, 4};
    double b[2] = {5, 6};
    double c[2] = {-1000, -1000};
    double cExp[2] = {17, 39};

    WENO::Project(2, 1, &a[0], &b[0], &c[0]);
    result = true;
    for (int i=0; i<2; i++) {
      result = result && (std::abs(c[i] - cExp[i]) < eps);
    }
    serialUnitResults.UpdateResult("WENO:Project:Mult:MatVec", result);
  }

  // --- Roe ---
  // Test Roe matrix routines
  // Assumes Euler equations with ideal gas EOS, no species transport

  double gamma = 1.4;

  // > 2D <
  // Test 2D matrix

  {
    int numDim = 2;
    double uL[4] = {1, 1, 1, 5};
    double uR[4] = {2, 4, 4, 16};
    double ujump[4];
    for (int i=0; i<numDim+2; i++) ujump[i] = uR[i] - uL[i];

    double tmat[16];
    double tinv[16];
    double lambda[16];

    // > Identity <
    // Test the identity Tmat*Tinv = I
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[2];
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
        } else {
          norm[i] = 0;
        }
      }

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[4];
      double tmp[4];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp1[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "ujump[" << i << "] = " << ujump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(ujump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:2D:Identity", result);
    }

    if (!result) return;

    // > RoeProp <
    // Test the identity A_hat [U] = [F], where A_hat = T*Lambda*Tinv
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[2];
      double fL[4];
      double fR[4];
      double fjump[4];
      fL[0] = 1;
      fL[numDim+1] = 6.6;
      fR[0] = 4;
      fR[numDim+1] = 38.4;
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
          fL[i+1] = 2.6;
          fR[i+1] = 11.2;
        } else {
          norm[i] = 0;
          fL[i+1] = 1;
          fR[i+1] = 8;
        }
      }
      for (int i=0; i<numDim+2; i++) fjump[i] = fR[i] - fL[i];

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[4];
      double tmp2[4];
      double tmp[4];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &lambda[0], &tmp1[0], &tmp2[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp2[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "fjump[" << i << "] = " << fjump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(fjump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:2D:RoeProp", result);
    }

    if (!result) return;
  }

  // > 3D <
  // Test 3D matrix

  {
    int numDim = 3;
    double uL[5] = {1, 1, 1, 1, 5.5};
    double uR[5] = {2, 4, 4, 4, 20};
    double ujump[5];
    for (int i=0; i<numDim+2; i++) ujump[i] = uR[i] - uL[i];

    double tmat[25];
    double tinv[25];
    double lambda[25];

    // > Identity <
    // Test the identity Tmat*Tinv = I
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[3];
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
        } else {
          norm[i] = 0;
        }
      }

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[5];
      double tmp[5];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp1[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "ujump[" << i << "] = " << ujump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(ujump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:3D:Identity", result);
    }

    if (!result) return;

    // > RoeProp <
    // Test the identity A_hat [U] = [F], where A_hat = T*Lambda*Tinv
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[3];
      double fL[5];
      double fR[5];
      double fjump[5];
      fL[0] = 1;
      fL[numDim+1] = 7.1;
      fR[0] = 4;
      fR[numDim+1] = 46.4;
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
          fL[i+1] = 2.6;
          fR[i+1] = 11.2;
        } else {
          norm[i] = 0;
          fL[i+1] = 1;
          fR[i+1] = 8;
        }
      }
      for (int i=0; i<numDim+2; i++) fjump[i] = fR[i] - fL[i];

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[5];
      double tmp2[5];
      double tmp[5];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &lambda[0], &tmp1[0], &tmp2[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp2[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "fjump[" << i << "] = " << fjump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(fjump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:3D:RoeProp", result);
    }
  }

  // --- RoeNeg ---
  // Test Roe matrix routines with negative background speed
  // Assumes Euler equations with ideal gas EOS, no species transport

  // > 2D <
  // Test 2D matrix

  {
    int numDim = 2;
    double uL[4] = {1, -1, -1, 5};
    double uR[4] = {2, -4, -4, 16};
    double ujump[4];
    for (int i=0; i<numDim+2; i++) ujump[i] = uR[i] - uL[i];

    double tmat[16];
    double tinv[16];
    double lambda[16];

    // > Identity <
    // Test the identity Tmat*Tinv = I
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[2];
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
        } else {
          norm[i] = 0;
        }
      }

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[4];
      double tmp[4];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp1[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "ujump[" << i << "] = " << ujump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(ujump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Neg:2D:Identity", result);
    }

    if (!result) return;

    // > RoeProp <
    // Test the identity A_hat [U] = [F], where A_hat = T*Lambda*Tinv
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[2];
      double fL[4];
      double fR[4];
      double fjump[4];
      fL[0] = -1;
      fL[numDim+1] = -6.6;
      fR[0] = -4;
      fR[numDim+1] = -38.4;
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
          fL[i+1] = 2.6;
          fR[i+1] = 11.2;
        } else {
          norm[i] = 0;
          fL[i+1] = 1;
          fR[i+1] = 8;
        }
      }
      for (int i=0; i<numDim+2; i++) fjump[i] = fR[i] - fL[i];

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[4];
      double tmp2[4];
      double tmp[4];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &lambda[0], &tmp1[0], &tmp2[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp2[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "fjump[" << i << "] = " << fjump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(fjump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Neg:2D:RoeProp", result);
    }

    if (!result) return;
  }

  // > 3D <
  // Test 3D matrix

  {
    int numDim = 3;
    double uL[5] = {1, -1, -1, -1, 5.5};
    double uR[5] = {2, -4, -4, -4, 20};
    double ujump[5];
    for (int i=0; i<numDim+2; i++) ujump[i] = uR[i] - uL[i];

    double tmat[25];
    double tinv[25];
    double lambda[25];

    // > Identity <
    // Test the identity Tmat*Tinv = I
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[3];
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
        } else {
          norm[i] = 0;
        }
      }

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[5];
      double tmp[5];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp1[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "ujump[" << i << "] = " << ujump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(ujump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Neg:3D:Identity", result);
    }

    if (!result) return;

    // > RoeProp <
    // Test the identity A_hat [U] = [F], where A_hat = T*Lambda*Tinv
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[3];
      double fL[5];
      double fR[5];
      double fjump[5];
      fL[0] = -1;
      fL[numDim+1] = -7.1;
      fR[0] = -4;
      fR[numDim+1] = -46.4;
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
          fL[i+1] = 2.6;
          fR[i+1] = 11.2;
        } else {
          norm[i] = 0;
          fL[i+1] = 1;
          fR[i+1] = 8;
        }
      }
      for (int i=0; i<numDim+2; i++) fjump[i] = fR[i] - fL[i];

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[5];
      double tmp2[5];
      double tmp[5];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &lambda[0], &tmp1[0], &tmp2[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp2[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "fjump[" << i << "] = " << fjump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(fjump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Neg:3D:RoeProp", result);
    }
  }

  // --- RoeAlt ---
  // Test Roe matrix routines with negative background speed
  // Assumes Euler equations with ideal gas EOS, no species transport

  // > 2D <
  // Test 2D matrix

  {
    int numDim = 2;
    double uL[4] = {1, 2, 2, 8};
    double uR[4] = {2, 6, 6, 26};
    double ujump[4];
    for (int i=0; i<numDim+2; i++) ujump[i] = uR[i] - uL[i];

    double tmat[16];
    double tinv[16];
    double lambda[16];

    // > Identity <
    // Test the identity Tmat*Tinv = I
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[2];
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
        } else {
          norm[i] = 0;
        }
      }

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[4];
      double tmp[4];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp1[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "ujump[" << i << "] = " << ujump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(ujump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Alt:2D:Identity", result);
    }

    if (!result) return;

    // > RoeProp <
    // Test the identity A_hat [U] = [F], where A_hat = T*Lambda*Tinv
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[2];
      double fL[4];
      double fR[4];
      double fjump[4];
      fL[0] = 2;
      fL[numDim+1] = 19.2;
      fR[0] = 6;
      fR[numDim+1] = 87.6;
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
          fL[i+1] = 5.6;
          fR[i+1] = 21.2;
        } else {
          norm[i] = 0;
          fL[i+1] = 4;
          fR[i+1] = 18;
        }
      }
      for (int i=0; i<numDim+2; i++) fjump[i] = fR[i] - fL[i];

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[4];
      double tmp2[4];
      double tmp[4];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &lambda[0], &tmp1[0], &tmp2[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp2[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "fjump[" << i << "] = " << fjump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(fjump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Alt:2D:RoeProp", result);
    }

    if (!result) return;
  }

  // > 3D <
  // Test 3D matrix

  {
    int numDim = 3;
    double uL[5] = {1, 2, 2, 2, 10};
    double uR[5] = {2, 6, 6, 6, 35};
    double ujump[5];
    for (int i=0; i<numDim+2; i++) ujump[i] = uR[i] - uL[i];

    double tmat[25];
    double tinv[25];
    double lambda[25];

    // > Identity <
    // Test the identity Tmat*Tinv = I
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[3];
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
        } else {
          norm[i] = 0;
        }
      }

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[5];
      double tmp[5];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp1[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "ujump[" << i << "] = " << ujump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(ujump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Alt:3D:Identity", result);
    }

    if (!result) return;

    // > RoeProp <
    // Test the identity A_hat [U] = [F], where A_hat = T*Lambda*Tinv
    for (int iDim=0; iDim<numDim && result; iDim++) {
      double norm[3];
      double fL[5];
      double fR[5];
      double fjump[5];
      fL[0] = 2;
      fL[numDim+1] = 23.2;
      fR[0] = 6;
      fR[numDim+1] = 114.6;
      for (int i=0; i<numDim; i++) {
        if (i == iDim) {
          norm[i] = 1;
          fL[i+1] = 5.6;
          fR[i+1] = 21.2;
        } else {
          norm[i] = 0;
          fL[i+1] = 4;
          fR[i+1] = 18;
        }
      }
      for (int i=0; i<numDim+2; i++) fjump[i] = fR[i] - fL[i];

      FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
        (&numDim, &uL[0], &uR[0], &gamma, &norm[0], &tmat[0], &tinv[0], &lambda[0]);

      double tmp1[5];
      double tmp2[5];
      double tmp[5];
      WENO::Project(numDim+2, 1, &tinv[0], &ujump[0], &tmp1[0]);
      WENO::Project(numDim+2, 1, &lambda[0], &tmp1[0], &tmp2[0]);
      WENO::Project(numDim+2, 1, &tmat[0], &tmp2[0], &tmp[0]);

      result = true;
      for (int i=0; i<numDim+2; i++) {
        //std::cout << "fjump[" << i << "] = " << fjump[i] << std::endl;
        //std::cout << "tmp[" << i << "] =  " << tmp[i] << std::endl;
        result = result && (std::abs(fjump[i] - tmp[i]) < eps);
      }
      serialUnitResults.UpdateResult("WENO:Project:Alt:3D:RoeProp", result);
    }
  }

}

void TestWENO_EntropyFix(ix::test::results &serialUnitResults) {
  // Tests for Roe entropy fix

  bool result = true;
  double eps = 1e-15;

  // --- Eta ---
  // Tests that entropy fix parameter eta is computed correctly

  int len = 4;
  int nCases = 3;
  double lambdaL[3][4] = {{1, 2, 2, 3}, {4, 2, 2, 0}, {1, -1, -1, -3}};
  double lambdaR[3][4] = {{1.5, 2, 2, 2.5}, {4, 1, 1, -2}, {0, -2, -2, -4}};
  double etaExp[3] = {0.25, 1, 0.5};
  double lambdaUseL[len*len];
  double lambdaUseR[len*len];

  for (int iCase=0; iCase<nCases && result; iCase++) {
    for (int i=0; i<len; i++) {
      for (int j=0; j<len; j++) {
        if (i == j) {
          lambdaUseL[i*len+j] = lambdaL[iCase][i];
          lambdaUseR[i*len+j] = lambdaR[iCase][i];
        } else {
          lambdaUseL[i*len+j] = 0;
          lambdaUseR[i*len+j] = 0;
        }
      }
    }

    double eta = WENO::EntropyFixEta(len, &lambdaUseL[0], &lambdaUseR[0]);
    if (std::abs(eta - etaExp[iCase]) > eps) {
      result = false;
      std::cout << "Case " << iCase << ", eta (actual/expected) = " << eta << " / " << etaExp[iCase] << std::endl;
    }
  }

  serialUnitResults.UpdateResult("WENO:EntropyFix:Eta", result);

}