TestMaxwellSolver.C

Go to the documentation of this file.

#include "Testing.H"
#include "Simulation.H"
#include "OperatorKernels.H"
#include "MaxwellSolver.H"
#include "RK4Advancer.H"
#include "EulerUtil.H"
#include "MaxwellUtil.H"
#include "MaxwellBC.H"
#include "PC2IO.H"

#define MACHINE_EPS     1e-12

#define EPS0            8.85418782e-12  // permittivity of free space [F/m]
#define MU0             1.25663706e-6   // permeability of free space [H/m]



void CreateLinearScalarField(int numDim, size_t *numX, double *deltaX, double *field, double *coeff)
{

  if (numDim == 3) {

    coeff[0] = 1.0;
    coeff[1] = 2.0;
    coeff[2] = 3.0;
    coeff[3] = 4.0;

    for (int k=0; k<numX[2]; k++) {
      for (int j=0; j<numX[1]; j++) {
        for (int i=0; i<numX[0]; i++) {
          field[i+j*numX[0]+k*numX[0]*numX[1]] = coeff[0]*(i*deltaX[0]) + coeff[1]*(j*deltaX[1]) + coeff[2]*(k*deltaX[2]) + coeff[3];
        }
      }
    }

  } else if (numDim == 2) {

    coeff[0] = 1.0;
    coeff[1] = 2.0;
    coeff[2] = 3.0;

    for (int j=0; j<numX[1]; j++) {
      for (int i=0; i<numX[0]; i++) {
        field[i+j*numX[0]] = coeff[0]*(i*deltaX[0]) + coeff[1]*(j*deltaX[1]) + coeff[2];
      }
    }

  } else if (numDim == 1) {

    coeff[0] = 1.0;
    coeff[1] = 2.0;

    for (int i=0; i<numX[0]; i++) {
      field[i] = coeff[0]*(i*deltaX[0]) + coeff[1];
    }

  } else {
    std::cout << "Invalid no. of dimensions. No. of dimensions = " << numDim << std::endl;
  }

}


void CreateConstantVectorField(int numDim, size_t *numX, double *value, double *inputField)
{

  if (numDim==3) {

    size_t numPoints = numX[0]*numX[1]*numX[2];

    for (int iDim=0; iDim<numDim; iDim++) {
      size_t iDimxnumPoints = iDim*numPoints;

      for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
        inputField[iDimxnumPoints+iPoint] = value[iDim];
      }
    }

  } else {
    std::cout << "Invalid no. of dimensions. Only 3-D is allowed." << std::endl;
  }

}


void CreateCurlFreeVectorFieldNonzeroDDx1(int numDim, size_t *numX, double *deltaX, double *inputField)
/******************************************************************************************************

                Create curl-free vector field with non-zero derivatives.
                Vector field is derived from the gradient of the potential:
                                phi(x,y,z) = coeff*x*y*z

 ******************************************************************************************************/
{

  if (numDim==3) {

    double coeff = 3.0;
    double X, Y, Z;

    size_t numPoints      = numX[0]*numX[1]*numX[2];
    size_t numX0xnumX1    = numX[0]*numX[1];
    size_t XYZIndex, XYIndex;

    // Calculate x-component
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex  = iZ*numX0xnumX1;
      Z         = iZ*deltaX[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex   = iY*numX[0];
        Y         = iY*deltaX[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          inputField[XYZIndex+XYIndex+iX] = coeff*Y*Z;
        }
      }
    }

    // Calculate y-component
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex  = iZ*numX0xnumX1;
      Z         = iZ*deltaX[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex    = iY*numX[0];

        for (size_t iX=0; iX<numX[0]; iX++) {
          X     = iX*deltaX[0];
          inputField[numPoints+XYZIndex+XYIndex+iX] = coeff*X*Z;
        }
      }
    }

    // Calculate z-component
    size_t twoxnumPoints = 2*numPoints;
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex   = iZ*numX0xnumX1;

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex = iY*numX[0];
        Y       = iY*deltaX[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          X     = iX*deltaX[0];
          inputField[twoxnumPoints+XYZIndex+XYIndex+iX] = coeff*X*Y;
        }
      }
    }

  } else {
    std::cout << "Invalid no. of dimensions. Only 3-D is allowed." << std::endl;
  }

}


void CreateCurlFreeVectorFieldNonzeroDDx2(int numDim, size_t *numX, double *deltaX, double *inputField)
/******************************************************************************************************

                Create curl-free vector field with non-zero derivatives.
                Vector field is derived from the gradient of the potential:
                                phi(x,y,z) = coeff*(x*y*z)^2

 ******************************************************************************************************/
{

  if (numDim==3) {

    double coeff = 3.0;
    double Y, Z;
    double Xsq, Ysq, Zsq;

    size_t numPoints      = numX[0]*numX[1]*numX[2];
    size_t numX0xnumX1    = numX[0]*numX[1];
    size_t XYZIndex, XYIndex;

    // Calculate squares of grid size in each dim.
    double deltaXsq[3];
    for (int iDim=0; iDim<numDim; iDim++) {
      deltaXsq[iDim] = deltaX[iDim]*deltaX[iDim];
    }

    // Calculate x-component
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex   = iZ*numX0xnumX1;
      Zsq        = iZ*iZ*deltaXsq[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex    = iY*numX[0];
        Ysq        = iY*iY*deltaXsq[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          inputField[XYZIndex+XYIndex+iX] = coeff*iX*deltaX[0]*Ysq*Zsq;
        }
      }
    }

    // Calculate y-component
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex   = iZ*numX0xnumX1;
      Zsq        = iZ*iZ*deltaXsq[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex    = iY*numX[0];
        Y          = iY*deltaX[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          Xsq = iX*iX*deltaXsq[0];
          inputField[numPoints+XYZIndex+XYIndex+iX] = coeff*Xsq*Y*Zsq;
        }
      }
    }

    // Calculate z-component
    size_t twoxnumPoints = 2*numPoints; 
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex   = iZ*numX0xnumX1;
      Z          = iZ*deltaX[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex    = iY*numX[0];
        Ysq        = iY*iY*deltaXsq[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          Xsq = iX*iX*deltaXsq[0];
          inputField[twoxnumPoints+XYZIndex+XYIndex+iX] = coeff*Xsq*Ysq*Z;
        }
      }
    }

  } else {
    std::cout << "Invalid no. of dimensions. Only 3-D is allowed." << std::endl; 
  }

}


void CreateConstCurlVectorField(int numDim, size_t *numX, double *deltaX, double *inputField, double *curlComp)
/**************************************************************************************************************

                Create a vector field with a constant curl everywhere.
                Vector field is given by:
                V(x,y,z) = (C00*x+C01*y+C02*z, C10*x+C11*y+C12*z, C20*x+C21*y+C22*z)

 **************************************************************************************************************/
{

  if (numDim==3) {

    double X, Y, Z;

    double C0[3] = { 3.0, -7.0, 1.0};
    double C1[3] = {-4.0,  9.0, 6.0};
    double C2[3] = { 5.0,  8.0, 2.0};

    size_t numPoints      = numX[0]*numX[1]*numX[2];
    size_t numX0xnumX1    = numX[0]*numX[1];
    size_t XYZIndex, XYIndex;

    // Calculate expected curl components
    curlComp[0] = C2[1]-C1[2];
    curlComp[1] = C0[2]-C2[0];
    curlComp[2] = C1[0]-C0[1];

    // Calculate x-component
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex  = iZ*numX0xnumX1;
      Z         = iZ*deltaX[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex   = iY*numX[0];
        Y         = iY*deltaX[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          X     = iX*deltaX[0];
          inputField[XYZIndex+XYIndex+iX] = C0[0]*X+C0[1]*Y+C0[2]*Z;
        }
      }
    }

    // Calculate y-component
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex  = iZ*numX0xnumX1;
      Z         = iZ*deltaX[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex   = iY*numX[0];
        Y         = iY*deltaX[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          X     = iX*deltaX[0];
          inputField[numPoints+XYZIndex+XYIndex+iX] = C1[0]*X+C1[1]*Y+C1[2]*Z;
        }
      }
    }

    // Calculate z-component
    size_t twoxnumPoints = 2*numPoints;
    for (size_t iZ=0; iZ<numX[2]; iZ++) {
      XYZIndex  = iZ*numX0xnumX1;
      Z         = iZ*deltaX[2];

      for (size_t iY=0; iY<numX[1]; iY++) {
        XYIndex = iY*numX[0];
        Y       = iY*deltaX[1];

        for (size_t iX=0; iX<numX[0]; iX++) {
          X     = iX*deltaX[0];
          inputField[twoxnumPoints+XYZIndex+XYIndex+iX] = C2[0]*X+C2[1]*Y+C2[2]*Z;
        }
      }
    }

  } else {
    std::cout << "Invalid no. of dimensions. Only 3-D is allowed." << std::endl;
  }

}


void CreateLinearIncreasingField(size_t numPoints, double coeff, double *inputField)
{

  for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
    inputField[iPoint] = coeff*(1.0+static_cast<double>(iPoint));
  }

}


void TestApplyOperatorAppl(ix::test::results &serialUnitResults)
/****************************************************************

        Test if ApplyOperator is applied correctly.

 ****************************************************************/
{

  int interiorOrder = 2;
  plascom2::operators::sbp::base sbpXX;
  plascom2::operators::sbp::Initialize(sbpXX, interiorOrder);

  // Forticate the stencil starts
  for (int iStencil=0; iStencil<sbpXX.numStencils; iStencil++)
    sbpXX.stencilStarts[iStencil]++;

  // Initialize test boolean operator
  bool MaxwellApplyOperator = true;

  int numDim        = 3;                        // no. of dims.
  int numComponents = 1;                        // some variable - currently unused
  int boundaryDepth = (sbpXX.numStencils-1)/2;  // no. of pts. near boundaries that must use boundary stencils
  size_t numPoints  = 1;                        // total no. of pts.

  size_t *numX          = new size_t [numDim];   // no. of pts. in each dim.
  double *lengthX       = new double [numDim];   // length in each dim.
  double *deltaX        = new double [numDim];   // grid size in each dim.
  double *recipDeltaX   = new double [numDim];   // reciprocal of grid size in each dim.
  size_t *opInterval    = new size_t [2*numDim]; // start and end indices in each dim.

  double *inputData     = NULL; // input field
  double *opData        = NULL; // operator output field
  int *stencilID        = NULL; // stencil ID for each pt. in each dim.

  // Set and calculate domain parameters
  for (int iDim=0; iDim<numDim; iDim++) {
    numX[iDim]          = std::pow(2,static_cast<double>(iDim))*2*boundaryDepth+20;
    lengthX[iDim]       = 1.0;
    deltaX[iDim]        = lengthX[iDim]/(numX[iDim]-1);
    recipDeltaX[iDim]   = 1.0/deltaX[iDim];
    numPoints *= numX[iDim];

    opInterval[2*iDim]   = 1;
    opInterval[2*iDim+1] = numX[iDim];
  }

  inputData     = new double [numPoints];
  opData        = new double [numDim*numPoints];
  stencilID     = new int [numDim*numPoints];

  // Create connectivity for stencil - set stencil ID
  plascom2::operators::sbp::CreateStencilConnectivity(numDim, numX, opInterval, boundaryDepth, stencilID, true);

  // Create linear scalar field: f = coeff*[x 1]'
  double *coeff = new double [numDim+1];
  CreateLinearScalarField(numDim, numX, deltaX, inputData, coeff);

  // Apply ApplyOperator to obtain df
  for (int opDir=1; opDir<=numDim; opDir++) {
    FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)(&numDim, numX, &numComponents, &numPoints,
                                                               &opDir, opInterval, &(sbpXX.numStencils), sbpXX.stencilSizes,
                                                               sbpXX.stencilStarts, &(sbpXX.numValues),
                                                               sbpXX.stencilWeights, sbpXX.stencilOffsets,
                                                               &stencilID[(opDir-1)*numPoints], inputData, &opData[(opDir-1)*numPoints]);
  }

  // Run test
  for (int iDim=0; iDim<numDim && MaxwellApplyOperator; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (int iPoints=0; iPoints<numPoints && MaxwellApplyOperator; iPoints++) {
      if (std::fabs(recipDeltaX[iDim]*opData[iPoints+iDimxnumPoints]-coeff[iDim]) > MACHINE_EPS) {
        MaxwellApplyOperator = false;
//      break;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestApplyOperatorAppl", MaxwellApplyOperator);

  // Free memory
  delete [] numX;
  delete [] lengthX;
  delete [] deltaX;
  delete [] recipDeltaX;
  delete [] opInterval;

  delete [] inputData;
  delete [] opData;
  delete [] stencilID;

  delete [] coeff;

}


void TestCurlOperator(ix::test::results &serialUnitResults)
/***********************************************************

        Test curl operator function.

 ***********************************************************/
{

  int interiorOrder = 4;
  plascom2::operators::sbp::base sbpXX;
  plascom2::operators::sbp::Initialize(sbpXX, interiorOrder);

  // Forticate the stencil starts
  for (int iStencil=0; iStencil<sbpXX.numStencils; iStencil++)
    sbpXX.stencilStarts[iStencil]++;

  // Initialize Boolean operator
  int numTests = 3;
  std::vector<bool> MaxwellCurlOperator(numTests, true);

  int numDim        = 3;                        // no. of dims.
  int numComponents = 1;                        // some variable - currently unused
  int boundaryDepth = (sbpXX.numStencils-1)/2;  // no. of pts. near boundaries that must use boundary stencils
  size_t numPoints  = 1;                        // total no. of pts.

  size_t *numX          = new size_t [numDim];   // no. of pts. in each dim.
  double *lengthX       = new double [numDim];   // length in each dim.
  double *deltaX        = new double [numDim];   // grid size in each dim.
  double *recipDeltaX   = new double [numDim];   // reciprocal of grid size in each dim.
  size_t *opInterval    = new size_t [2*numDim]; // start and end indices in each dim.

  double *inputField    = NULL; // input field
  double *outputField   = NULL; // operator output field
  int *stencilID        = NULL; // stencil ID for each pt. in each dim.

  double *recipDeltaXField = NULL; // field of reciprocal of grid size in each dim.

  // Set and calculate domain parameters
  for (int iDim=0; iDim<numDim; iDim++) {
    numX[iDim]          = std::pow(2,static_cast<double>(iDim))*2*boundaryDepth+20;
    lengthX[iDim]       = 1.0;
    deltaX[iDim]        = lengthX[iDim]/(numX[iDim]-1);
    recipDeltaX[iDim]   = 1.0/deltaX[iDim];
    numPoints *= numX[iDim];

    opInterval[2*iDim]   = 1;
    opInterval[2*iDim+1] = numX[iDim];
  }

  recipDeltaXField = new double [numDim*numPoints];
  for (int iDim=0; iDim<numDim; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
      recipDeltaXField[iDimxnumPoints+iPoint] = recipDeltaX[iDim];
    }
  }

  // Create connectivity for stencil - set stencil ID
  stencilID     = new int [numDim*numPoints];
  plascom2::operators::sbp::CreateStencilConnectivity(numDim, numX, opInterval, boundaryDepth, stencilID, true);


  /* ======== CASE 1: Constant field ======== */
  inputField    = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate field
  double constComp[3] = {0.0, 1.0, 2.0};
  CreateConstantVectorField(numDim, numX, constComp, inputField);

  // Compute curl of field
  Maxwell::ComputeCurl(numDim, numX, numComponents,
                       numPoints, opInterval, stencilID,
                       sbpXX, recipDeltaXField, inputField, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellCurlOperator[0]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellCurlOperator[0]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellCurlOperator[0] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestCurlOperator:Case1", MaxwellCurlOperator[0]);

  delete [] inputField;
  delete [] outputField;


//  /* ======== CASE 1: Curl-free field (zero partial derivatives) ======== */
//  inputField = new double [numDim*numPoints];
//
//  // Populate field
//  CreateCurlFreeVectorFieldZeroDDx();
//
//  // Compute curl of field
//
//  // Run test
//
//  serialUnitResults.UpdateResult("Maxwell:TestCurlOperator:Case1", MaxwellCurlOperator[0]);
//
//  delete [] inputField;


  /* ======== CASE 2: Curl-free field (non-zero partial derivatives) ======== */
  inputField    = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate field
  CreateCurlFreeVectorFieldNonzeroDDx1(numDim, numX, deltaX, inputField);

  // Compute curl of field
  Maxwell::ComputeCurl(numDim, numX, numComponents,
                       numPoints, opInterval, stencilID,
                       sbpXX, recipDeltaXField, inputField, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellCurlOperator[1]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellCurlOperator[1]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellCurlOperator[1] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestCurlOperator:Case2", MaxwellCurlOperator[1]);

  delete [] inputField;
  delete [] outputField;


  /* ======== CASE 3: Constant-curl field ======== */
  double curlComp[3];
  inputField    = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate field
  CreateConstCurlVectorField(numDim, numX, deltaX, inputField, curlComp);

  // Compute curl of field
  Maxwell::ComputeCurl(numDim, numX, numComponents,
                       numPoints, opInterval, stencilID,
                       sbpXX, recipDeltaXField, inputField, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellCurlOperator[2]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellCurlOperator[2]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]-curlComp[iDim]) > MACHINE_EPS) {
        MaxwellCurlOperator[2] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestCurlOperator:Case3", MaxwellCurlOperator[2]);

  delete [] inputField;
  delete [] outputField;


  // Free memory
  delete [] numX;
  delete [] lengthX;
  delete [] deltaX;
  delete [] recipDeltaX;
  delete [] opInterval;

  delete [] stencilID;
  delete [] recipDeltaXField;

}


void TestMaxwellRHS_Bfield(ix::test::results &serialUnitResults)
/***************************************************************

        Test RHS function for B-field (Faraday's law):
                        dBdt = -curl(E)

 ***************************************************************/
{

  int interiorOrder = 4;
  plascom2::operators::sbp::base sbpXX;
  plascom2::operators::sbp::Initialize(sbpXX, interiorOrder);

  // Forticate the stencil starts
  for (int iStencil=0; iStencil<sbpXX.numStencils; iStencil++)
    sbpXX.stencilStarts[iStencil]++;

  // Initialize Boolean operator
  int numTests = 3;
  std::vector<bool> MaxwellRHS_Bfield(numTests, true);

  int numDim        = 3;                        // no. of dims.
  int numComponents = 1;                        // some variable - currently unused
  int boundaryDepth = (sbpXX.numStencils-1)/2;  // no. of pts. near boundaries that must use boundary stencils
  size_t numPoints  = 1;                        // total no. of pts.

  size_t *numX          = new size_t [numDim];   // no. of pts. in each dim.
  double *lengthX       = new double [numDim];   // length in each dim.
  double *deltaX        = new double [numDim];   // grid size in each dim.
  double *recipDeltaX   = new double [numDim];   // reciprocal of grid size in each dim.
  size_t *opInterval    = new size_t [2*numDim]; // start and end indices in each dim.

  double *inputField    = NULL; // input field
  double *outputField   = NULL; // operator output field
  int *stencilID        = NULL; // stencil ID for each pt. in each dim.

  double *recipDeltaXField = NULL; // field of reciprocal of grid size in each dim.

  // Set and calculate domain parameters
  for (int iDim=0; iDim<numDim; iDim++) {
    numX[iDim]          = std::pow(2,static_cast<double>(iDim))*2*boundaryDepth+20;
    lengthX[iDim]       = 1.0;
    deltaX[iDim]        = lengthX[iDim]/(numX[iDim]-1);
    recipDeltaX[iDim]   = 1.0/deltaX[iDim];
    numPoints *= numX[iDim];

    opInterval[2*iDim]   = 1;
    opInterval[2*iDim+1] = numX[iDim];
  }

  recipDeltaXField = new double [numDim*numPoints];
  for (int iDim=0; iDim<numDim; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
      recipDeltaXField[iDimxnumPoints+iPoint] = recipDeltaX[iDim];
    }
  }

  // Create connectivity for stencil - set stencil ID
  stencilID     = new int [numDim*numPoints];
  plascom2::operators::sbp::CreateStencilConnectivity(numDim, numX, opInterval, boundaryDepth, stencilID, true);


  /* ======== CASE 1: Constant field ======== */
  inputField    = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate field
  double constComp[3] = {0.0, 1.0, 2.0};
  CreateConstantVectorField(numDim, numX, constComp, inputField);

  // Compute RHS for B-field
  Maxwell::ComputeMaxwellRHS_Bfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Bfield[0]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Bfield[0]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellRHS_Bfield[0] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Bfield:Case1", MaxwellRHS_Bfield[0]);

  delete [] inputField;
  delete [] outputField;


//  /* ======== CASE 1: Curl-free field (zero partial derivatives) ======== */
//  inputField = new double [numDim*numPoints];
//
//  // Populate field
//  CreateCurlFreeVectorFieldZeroDDx();
//
//  // Compute curl of field
//
//  // Run test
//
//  serialUnitResults.UpdateResult("Maxwell:TestCurlOperator:Case1", MaxwellCurlOperator[0]);
//
//  delete [] inputField;


  /* ======== CASE 2: Curl-free field (non-zero partial derivatives) ======== */
  inputField    = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate field
  CreateCurlFreeVectorFieldNonzeroDDx1(numDim, numX, deltaX, inputField);

  // Compute RHS for B-field
  Maxwell::ComputeMaxwellRHS_Bfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Bfield[1]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Bfield[1]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellRHS_Bfield[1] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Bfield:Case2", MaxwellRHS_Bfield[1]);

  delete [] inputField;
  delete [] outputField;


  /* ======== CASE 3: Constant-curl field ======== */
  double curlComp[3];
  inputField    = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate field
  CreateConstCurlVectorField(numDim, numX, deltaX, inputField, curlComp);

  // Compute RHS for B-field
  Maxwell::ComputeMaxwellRHS_Bfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Bfield[2]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Bfield[2]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]+curlComp[iDim]) > MACHINE_EPS) {
        MaxwellRHS_Bfield[2] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Bfield:Case3", MaxwellRHS_Bfield[2]);

  delete [] inputField;
  delete [] outputField;


  // Free memory
  delete [] numX;
  delete [] lengthX;
  delete [] deltaX;
  delete [] recipDeltaX;
  delete [] opInterval;

  delete [] stencilID;
  delete [] recipDeltaXField;

}


void TestMaxwellRHS_Dfield(ix::test::results &serialUnitResults)
/***************************************************************

        Test RHS function for D-field (Ampere's law):
                        dDdt = curl(H)-J

 ***************************************************************/
{

  int interiorOrder = 4;
  plascom2::operators::sbp::base sbpXX;
  plascom2::operators::sbp::Initialize(sbpXX, interiorOrder);

  // Forticate the stencil starts
  for (int iStencil=0; iStencil<sbpXX.numStencils; iStencil++)
    sbpXX.stencilStarts[iStencil]++;

  // Initialize Boolean operator
  int numTests = 5;
  std::vector<bool> MaxwellRHS_Dfield(numTests, true);

  int numDim        = 3;                        // no. of dims.
  int numComponents = 1;                        // some variable - currently unused
  int boundaryDepth = (sbpXX.numStencils-1)/2;  // no. of pts. near boundaries that must use boundary stencils
  size_t numPoints  = 1;                        // total no. of pts.

  size_t *numX          = new size_t [numDim];   // no. of pts. in each dim.
  double *lengthX       = new double [numDim];   // length in each dim.
  double *deltaX        = new double [numDim];   // grid size in each dim.
  double *recipDeltaX   = new double [numDim];   // reciprocal of grid size in each dim.
  size_t *opInterval    = new size_t [2*numDim]; // start and end indices in each dim.

  double *inputField    = NULL; // input field
  double *J             = NULL; // current density field
  double *outputField   = NULL; // operator output field
  int *stencilID        = NULL; // stencil ID for each pt. in each dim.

  double *recipDeltaXField = NULL; // field of reciprocal of grid size in each dim.

  // Set and calculate domain parameters
  for (int iDim=0; iDim<numDim; iDim++) {
    numX[iDim]          = std::pow(2,static_cast<double>(iDim))*2*boundaryDepth+20;
    lengthX[iDim]       = 1.0;
    deltaX[iDim]        = lengthX[iDim]/(numX[iDim]-1);
    recipDeltaX[iDim]   = 1.0/deltaX[iDim];
    numPoints *= numX[iDim];

    opInterval[2*iDim]   = 1;
    opInterval[2*iDim+1] = numX[iDim];
  }

  recipDeltaXField = new double [numDim*numPoints];
  for (int iDim=0; iDim<numDim; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
      recipDeltaXField[iDimxnumPoints+iPoint] = recipDeltaX[iDim];
    }
  }

  // Create connectivity for stencil - set stencil ID
  stencilID     = new int [numDim*numPoints];
  plascom2::operators::sbp::CreateStencilConnectivity(numDim, numX, opInterval, boundaryDepth, stencilID, true);


  /* ======== CASE 1: Constant field and J=0 ======== */
  inputField    = new double [numDim*numPoints];
  J             = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate fields
  double constComp1[3] = {0.0, 1.0, 2.0};
  CreateConstantVectorField(numDim, numX, constComp1, inputField);

  double constCompJ1[3] = {0.0, 0.0, 0.0};
  CreateConstantVectorField(numDim, numX, constCompJ1, J);

  // Compute RHS for D-field
  Maxwell::ComputeMaxwellRHS_Dfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, J, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Dfield[0]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Dfield[0]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellRHS_Dfield[0] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Dfield:Case1", MaxwellRHS_Dfield[0]);

  delete [] inputField;
  delete [] J;
  delete [] outputField;


//  /* ======== CASE 1: Curl-free field (zero partial derivatives) ======== */
//  inputField = new double [numDim*numPoints];
//
//  // Populate field
//  CreateCurlFreeVectorFieldZeroDDx();
//
//  // Compute curl of field
//
//  // Run test
//
//  serialUnitResults.UpdateResult("Maxwell:TestCurlOperator:Case1", MaxwellCurlOperator[0]);
//
//  delete [] inputField;


  /* ======== CASE 2: Curl-free field (non-zero partial derivatives) and J=0 ======== */
  inputField    = new double [numDim*numPoints];
  J             = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate fields
  CreateCurlFreeVectorFieldNonzeroDDx1(numDim, numX, deltaX, inputField);

  double constCompJ2[3] = {0.0, 0.0, 0.0};
  CreateConstantVectorField(numDim, numX, constCompJ2, J);

  // Compute RHS for D-field
  Maxwell::ComputeMaxwellRHS_Dfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, J, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Dfield[1]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Dfield[1]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellRHS_Dfield[1] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Dfield:Case2", MaxwellRHS_Dfield[1]);

  delete [] inputField;
  delete [] J;
  delete [] outputField;


  /* ======== CASE 3: Constant-curl field and J=0 ======== */
  double curlComp3[3];
  inputField    = new double [numDim*numPoints];
  J             = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate fields
  CreateConstCurlVectorField(numDim, numX, deltaX, inputField, curlComp3);

  double constCompJ3[3] = {0.0, 0.0, 0.0};
  CreateConstantVectorField(numDim, numX, constCompJ3, J);

  // Compute RHS for D-field
  Maxwell::ComputeMaxwellRHS_Dfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, J, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Dfield[2]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Dfield[2]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]-curlComp3[iDim]) > MACHINE_EPS) {
        MaxwellRHS_Dfield[2] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Dfield:Case3", MaxwellRHS_Dfield[2]);

  delete [] inputField;
  delete [] J;
  delete [] outputField;


  /* ======== CASE 4: Curl-free field (non-zero partial derivatives) and J=J(x,y,z) ======== */
  inputField    = new double [numDim*numPoints];
  J             = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate fields
  CreateCurlFreeVectorFieldNonzeroDDx1(numDim, numX, deltaX, inputField);
  CreateCurlFreeVectorFieldNonzeroDDx1(numDim, numX, deltaX, J);

  // Compute RHS for D-field
  Maxwell::ComputeMaxwellRHS_Dfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, J, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Dfield[3]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Dfield[3]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]+J[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellRHS_Dfield[3] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Dfield:Case4", MaxwellRHS_Dfield[3]);

  delete [] inputField;
  delete [] J;
  delete [] outputField;


  /* ======== CASE 5: Constant-curl field and J=(c1,c2,c3) ======== */
  double curlComp5[3];
  inputField    = new double [numDim*numPoints];
  J             = new double [numDim*numPoints];
  outputField   = new double [numDim*numPoints];

  // Populate fields
  CreateConstCurlVectorField(numDim, numX, deltaX, inputField, curlComp5);

  double constCompJ5[3] = {-1.0, -2.0, 3.0};
  CreateConstantVectorField(numDim, numX, constCompJ5, J);

  // Compute RHS for D-field
  Maxwell::ComputeMaxwellRHS_Dfield(numDim, numX, numComponents,
                                    numPoints, opInterval, stencilID,
                                    sbpXX, recipDeltaXField, inputField, J, outputField);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellRHS_Dfield[4]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellRHS_Dfield[4]; iPoint++) {
      if (std::fabs(outputField[iDimxnumPoints+iPoint]-curlComp5[iDim]+J[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellRHS_Dfield[4] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS_Dfield:Case5", MaxwellRHS_Dfield[4]);

  delete [] inputField;
  delete [] J;
  delete [] outputField;


  // Free memory
  delete [] numX;
  delete [] lengthX;
  delete [] deltaX;
  delete [] recipDeltaX;
  delete [] opInterval;

  delete [] stencilID;
  delete [] recipDeltaXField;

}


void TestComputeRecipEpsMu(ix::test::results &serialUnitResults)
/****************************************************************

        Test function computing reciprocals
        of epsilon and mu.

 ****************************************************************/
{

  double eps0 = EPS0;
  double mu0  = MU0;

  int numDim        = 3;                        // no. of dims.
  size_t numPoints  = 1;                        // total no. of pts.
  size_t *numX          = new size_t [numDim];   // no. of pts. in each dim.
  size_t *opInterval    = new size_t [2*numDim]; // start and end indices in each dim.

  double *epsMu         = NULL;
  double *recipEpsMu    = NULL;

  // Initialize Boolean operator
  int numTests = 2;
  std::vector<bool> MaxwellComputeRecipEpsMu(numTests, true);

  // Set and calculate domain parameters
  for (int iDim=0; iDim<numDim; iDim++) {
    numX[iDim] = 3+iDim;
    numPoints *= numX[iDim];

    opInterval[2*iDim]   = 1;
    opInterval[2*iDim+1] = numX[iDim];
  }

  size_t twoxnumPoints =  2*numPoints;


  /* ======== CASE 1: Constant eps, mu ======== */
  epsMu         = new double [2*numPoints];
  recipEpsMu    = new double [2*numPoints];

  // Populate fields
  FC_MODULE(operators,assignmentxa,OPERATORS,ASSIGNMENTXA)(&numDim, &numPoints, numX, opInterval, &eps0, &epsMu[0]);
  FC_MODULE(operators,assignmentxa,OPERATORS,ASSIGNMENTXA)(&numDim, &numPoints, numX, opInterval, &mu0, &epsMu[numPoints]);

  // Compute reciprocals
  Maxwell::ComputeRecipEpsMu(numDim, numPoints, numX, opInterval, epsMu, recipEpsMu);

  // Execute test
  for (size_t iPoint=0; iPoint<twoxnumPoints && MaxwellComputeRecipEpsMu[0]; iPoint++) {
    if (std::fabs(epsMu[iPoint]*recipEpsMu[iPoint]-1.0) > MACHINE_EPS) {
      MaxwellComputeRecipEpsMu[0] = false;
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestComputeRecipEpsMu:Case1", MaxwellComputeRecipEpsMu[0]);

  delete [] epsMu;
  delete [] recipEpsMu;


  /* ======== CASE 2: Spatially-varying eps, mu ======== */
  epsMu         = new double [2*numPoints];
  recipEpsMu    = new double [2*numPoints];

  // Populate fields
  CreateLinearIncreasingField(numPoints, eps0, &epsMu[0]);
  CreateLinearIncreasingField(numPoints, mu0, &epsMu[numPoints]);

  // Compute reciprocals
  Maxwell::ComputeRecipEpsMu(numDim, numPoints, numX, opInterval, epsMu, recipEpsMu);

  // Execute test
  for (size_t iPoint=0; iPoint<twoxnumPoints && MaxwellComputeRecipEpsMu[1]; iPoint++) {
    if (std::fabs(epsMu[iPoint]*recipEpsMu[iPoint]-1.0) > MACHINE_EPS) {
      MaxwellComputeRecipEpsMu[1] = false;
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestComputeRecipEpsMu:Case2", MaxwellComputeRecipEpsMu[1]);

  delete [] epsMu;
  delete [] recipEpsMu;


  // Free memory
  delete [] numX;
  delete [] opInterval;

}


void TestConvertFields(ix::test::results &serialUnitResults)
/***********************************************************

        Test functions converting between fields:
        E->D, D->E, B->H and H->B

 ***********************************************************/
{

  double eps0 = EPS0;
  double mu0  = MU0;

  int numDim        = 3;                        // no. of dims.
  size_t numPoints  = 1;                        // total no. of pts.
  size_t *numX          = new size_t [numDim];   // no. of pts. in each dim.
  size_t *opInterval    = new size_t [2*numDim]; // start and end indices in each dim.

  double *epsMu         = NULL;
  double *recipEpsMu    = NULL;

  double *Efield_i      = NULL;
  double *Dfield        = NULL;
  double *Efield_f      = NULL;

  double *Bfield_i      = NULL;
  double *Hfield        = NULL;
  double *Bfield_f      = NULL;

  // Initialize Boolean operator
  int numTests = 12;
  std::vector<bool> MaxwellConvertFields(numTests, true);

  // Set and calculate domain parameters
  for (int iDim=0; iDim<numDim; iDim++) {
    numX[iDim] = 3+iDim;
    numPoints *= numX[iDim];

    opInterval[2*iDim]   = 1;
    opInterval[2*iDim+1] = numX[iDim];
  }


  /* ======== CASES 1-4: Constant eps, mu and fields ======== */
  epsMu         = new double [2*numPoints];
  recipEpsMu    = new double [2*numPoints];

  Efield_i      = new double [numDim*numPoints];
  Dfield        = new double [numDim*numPoints];
  Efield_f      = new double [numDim*numPoints];

  Bfield_i      = new double [numDim*numPoints];
  Hfield        = new double [numDim*numPoints];
  Bfield_f      = new double [numDim*numPoints];

  // Populate fields
  FC_MODULE(operators,assignmentxa,OPERATORS,ASSIGNMENTXA)(&numDim, &numPoints, numX, opInterval, &eps0, &epsMu[0]);
  FC_MODULE(operators,assignmentxa,OPERATORS,ASSIGNMENTXA)(&numDim, &numPoints, numX, opInterval, &mu0, &epsMu[numPoints]);

  Maxwell::ComputeRecipEpsMu(numDim, numPoints, numX, opInterval, epsMu, recipEpsMu);

  double constCompEi[3] = {1.2, 2.6, 3.8};
  CreateConstantVectorField(numDim, numX, constCompEi, Efield_i);

  double constCompBi[3] = {0.6, 4.4, 8.9};
  CreateConstantVectorField(numDim, numX, constCompBi, Bfield_i);

  /* == CASE 1: Test conversion E->D == */
  Maxwell::ConvertEfieldtoDfield(numDim, numPoints, numX, opInterval, epsMu, Efield_i, Dfield);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[0]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[0]; iPoint++) {
      if (std::fabs(Dfield[iDimxnumPoints+iPoint]-epsMu[iPoint]*Efield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[0] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case1", MaxwellConvertFields[0]);

  /* == CASE 2: Test conversion D->E == */
  Maxwell::ConvertDfieldtoEfield(numDim, numPoints, numX, opInterval, recipEpsMu, Dfield, Efield_f);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[1]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[1]; iPoint++) {
      if (std::fabs(Efield_f[iDimxnumPoints+iPoint]-recipEpsMu[iPoint]*Dfield[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[1] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case2", MaxwellConvertFields[1]);

  /* == CASE 1/2: Is Efield_f == Efield_i? == */
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[2]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[2]; iPoint++) {
      if (std::fabs(Efield_f[iDimxnumPoints+iPoint]-Efield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[2] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case1/2", MaxwellConvertFields[2]);

  /* == CASE 3: Test conversion B->H == */
  Maxwell::ConvertBfieldtoHfield(numDim, numPoints, numX, opInterval, recipEpsMu, Bfield_i, Hfield);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[3]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[3]; iPoint++) {
      if (std::fabs(Hfield[iDimxnumPoints+iPoint]-recipEpsMu[numPoints+iPoint]*Bfield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[3] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case3", MaxwellConvertFields[3]);

  /* == CASE 4: Test conversion H->B == */
  Maxwell::ConvertHfieldtoBfield(numDim, numPoints, numX, opInterval, epsMu, Hfield, Bfield_f);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[4]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[4]; iPoint++) {
      if (std::fabs(Bfield_f[iDimxnumPoints+iPoint]-epsMu[numPoints+iPoint]*Hfield[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[4] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case4", MaxwellConvertFields[4]);

  /* == CASE 3/4: Is Bfield_f == Bfield_i? == */
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[5]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[5]; iPoint++) {
      if (std::fabs(Bfield_f[iDimxnumPoints+iPoint]-Bfield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[5] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case3/4", MaxwellConvertFields[5]);

  delete [] epsMu;
  delete [] recipEpsMu;

  delete [] Efield_i;
  delete [] Dfield;
  delete [] Efield_f;

  delete [] Bfield_i;
  delete [] Hfield;
  delete [] Bfield_f;


  /* ======== CASES 5-8: Spatially-varying eps, mu and fields ======== */
  epsMu         = new double [2*numPoints];
  recipEpsMu    = new double [2*numPoints];

  Efield_i      = new double [numDim*numPoints];
  Dfield        = new double [numDim*numPoints];
  Efield_f      = new double [numDim*numPoints];

  Bfield_i      = new double [numDim*numPoints];
  Hfield        = new double [numDim*numPoints];
  Bfield_f      = new double [numDim*numPoints];

  // Populate fields
  CreateLinearIncreasingField(numPoints, eps0, &epsMu[0]);
  CreateLinearIncreasingField(numPoints, mu0, &epsMu[numPoints]);

  Maxwell::ComputeRecipEpsMu(numDim, numPoints, numX, opInterval, epsMu, recipEpsMu);

  double *deltaX = new double [numDim];
  for (int iDim=0; iDim<numDim; iDim++) {
    deltaX[iDim] = 0.1*(1.0+static_cast<double>(iDim));
  }
  CreateCurlFreeVectorFieldNonzeroDDx1(numDim, numX, deltaX, Efield_i);
  CreateCurlFreeVectorFieldNonzeroDDx1(numDim, numX, deltaX, Bfield_i);
  delete [] deltaX;

  /* == CASE 5: Test conversion E->D == */
  Maxwell::ConvertEfieldtoDfield(numDim, numPoints, numX, opInterval, epsMu, Efield_i, Dfield);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[6]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[6]; iPoint++) {
      if (std::fabs(Dfield[iDimxnumPoints+iPoint]-epsMu[iPoint]*Efield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[6] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case5", MaxwellConvertFields[6]);

  /* == CASE 6: Test conversion D->E == */
  Maxwell::ConvertDfieldtoEfield(numDim, numPoints, numX, opInterval, recipEpsMu, Dfield, Efield_f);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[7]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[7]; iPoint++) {
      if (std::fabs(Efield_f[iDimxnumPoints+iPoint]-recipEpsMu[iPoint]*Dfield[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[7] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case6", MaxwellConvertFields[7]);

  /* == CASE 5/6: Is Efield_f == Efield_i? == */
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[8]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[8]; iPoint++) {
      if (std::fabs(Efield_f[iDimxnumPoints+iPoint]-Efield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[8] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case5/6", MaxwellConvertFields[8]);

  /* == CASE 7: Test conversion B->H == */
  Maxwell::ConvertBfieldtoHfield(numDim, numPoints, numX, opInterval, recipEpsMu, Bfield_i, Hfield);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[9]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[9]; iPoint++) {
      if (std::fabs(Hfield[iDimxnumPoints+iPoint]-recipEpsMu[numPoints+iPoint]*Bfield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[9] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case7", MaxwellConvertFields[9]);

  /* == CASE 8: Test conversion H->B == */
  Maxwell::ConvertHfieldtoBfield(numDim, numPoints, numX, opInterval, epsMu, Hfield, Bfield_f);

  // Execute test
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[10]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[10]; iPoint++) {
      if (std::fabs(Bfield_f[iDimxnumPoints+iPoint]-epsMu[numPoints+iPoint]*Hfield[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[10] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case8", MaxwellConvertFields[10]);

  /* == CASE 7/8: Is Bfield_f == Bfield_i? == */
  for (int iDim=0; iDim<numDim && MaxwellConvertFields[11]; iDim++) {
    size_t iDimxnumPoints = iDim*numPoints;

    for (size_t iPoint=0; iPoint<numPoints && MaxwellConvertFields[11]; iPoint++) {
      if (std::fabs(Bfield_f[iDimxnumPoints+iPoint]-Bfield_i[iDimxnumPoints+iPoint]) > MACHINE_EPS) {
        MaxwellConvertFields[11] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestConvertFields:Case7/8", MaxwellConvertFields[11]);

  delete [] epsMu;
  delete [] recipEpsMu;

  delete [] Efield_i;
  delete [] Dfield;
  delete [] Efield_f;

  delete [] Bfield_i;
  delete [] Hfield;
  delete [] Bfield_f;

}


void TestDirichletBC(ix::test::results &serialUnitResults)
/**********************************************************

        Test Dirichlet boundary condition.

 **********************************************************/
{

  int numDim = 3;
  size_t numPointsBuffer, numPointsIndices;
  double *bufferVar = NULL;
  double *value     = NULL;

  // Initialize Boolean operator
  int numTests = 2;
  std::vector<bool> MaxwellDirichletBC(numTests, true);


  /* ==== CASE 1: Set 2nd half only ==== */
  numPointsBuffer  = 120; // = 4 x 5 x 6
  numPointsIndices = 60;
  std::vector<size_t> bufferIndices1(numPointsIndices, 0);
  bufferVar = new double [numDim*numPointsBuffer];
  value     = new double [numDim];

  // Initialize buffer variables to zero
  size_t numDimxnumPointsBuffer = numDim*numPointsBuffer;
  for(size_t iPoint=0; iPoint<numDimxnumPointsBuffer; iPoint++) {
    bufferVar[iPoint] = 0.0;
  }

  // Assign buffer indices
  for(size_t iPoint=0; iPoint<numPointsIndices; iPoint++) {
    bufferIndices1[iPoint] = 60+iPoint;
  }

  // Set values
  for(int iDim=0; iDim<numDim; iDim++) {
    value[iDim] = 1.0+static_cast<double>(iDim);
  }

  // Compute Dirichlet boundary condition
  for(int iDim=0; iDim<numDim; iDim++) {
    double *bufferPtr = bufferVar + iDim*numPointsBuffer;
    if(Maxwell::bc::ComputeDirichletBC(bufferIndices1, bufferPtr, value[iDim]))
      MaxwellDirichletBC[0] = false;
  }

  // Execute test
  for(int iDim=0; iDim<numDim; iDim++) {
    size_t iDimxnumPointsBuffer = iDim*numPointsBuffer;

    for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
      if(iPoint<60 && bufferVar[iPoint+iDimxnumPointsBuffer]!=0.0)
        MaxwellDirichletBC[0] = false;

      if(iPoint>=60 && bufferVar[iPoint+iDimxnumPointsBuffer]!=value[iDim])
        MaxwellDirichletBC[0] = false;
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestDirichletBC:Case1", MaxwellDirichletBC[0]);

  delete [] bufferVar;
  delete [] value;


  /* ==== CASE 2: Set every other element ==== */
  numPointsBuffer  = 120; // = 4 x 5 x 6
  numPointsIndices = 60;
  std::vector<size_t> bufferIndices2(numPointsIndices, 0);
  bufferVar = new double [numDim*numPointsBuffer];
  value     = new double [numDim];

  // Initialize buffer variables to zero
  numDimxnumPointsBuffer = numDim*numPointsBuffer;
  for(size_t iPoint=0; iPoint<numDimxnumPointsBuffer; iPoint++) {
    bufferVar[iPoint] = 0.0;
  }

  // Assign buffer indices
  for(size_t iPoint=0; iPoint<numPointsIndices; iPoint++) {
    bufferIndices2[iPoint] = 2*iPoint+1;
  }

  // Set values
  for(int iDim=0; iDim<numDim; iDim++) {
    value[iDim] = 1.0+static_cast<double>(iDim);
  }

  // Compute Dirichlet boundary condition
  for(int iDim=0; iDim<numDim; iDim++) {
    double *bufferPtr = bufferVar + iDim*numPointsBuffer;
    if(Maxwell::bc::ComputeDirichletBC(bufferIndices2, bufferPtr, value[iDim]))
      MaxwellDirichletBC[1] = false;
  }

  // Execute test
  for(int iDim=0; iDim<numDim; iDim++) {
    size_t iDimxnumPointsBuffer = iDim*numPointsBuffer;

    for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
      if(iPoint%2) {
        if(bufferVar[iPoint+iDimxnumPointsBuffer]!=value[iDim])
          MaxwellDirichletBC[1] = false;
      } else {
        if(bufferVar[iPoint+iDimxnumPointsBuffer]!=0.0)
          MaxwellDirichletBC[1] = false;
      }
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestDirichletBC:Case2", MaxwellDirichletBC[1]);

  delete [] bufferVar;
  delete [] value;

}


using pcpp::comm::CheckResult;

typedef simulation::grid::parallel_blockstructured                grid_t;
typedef simulation::state::base                                   state_t;
typedef simulation::grid::halo                                    halo_t;
typedef pcpp::IndexIntervalType                                   interval_t;
typedef plascom2::operators::sbp::base                            operator_t;
typedef pcpp::ParallelGlobalType                                  global_t;

void TestMaxwellRHS(ix::test::results &serialUnitResults)
/********************************************************

        Test individual components that form
        the RHS API for time integration.

 ********************************************************/
{

  int myThreadID = 0; // thread ID

  // basic types for the advection problem
  typedef plascom2::operators::sbp::base                     operator_t;
  typedef simulation::domain::base<grid_t,state_t>  domain_t;
  typedef Maxwell::rhs<grid_t,state_t,operator_t>            rhs_t;

  // -- Set up the RHS --
  rhs_t testRHS;

  // -- Set up the grid --
  grid_t testGrid;

  // -- Set up the halo --
  halo_t &testHalo(testGrid.Halo());

  int numDim = 3;
  int interiorOrder = 2;
  operator_t sbpXX;

  plascom2::operators::sbp::Initialize(sbpXX, interiorOrder);
  int boundaryDepth = (sbpXX.numStencils-1)/2;
  std::vector<int> haloDepths(2*numDim,boundaryDepth);
  std::vector<int> boundaryDepths(2*numDim,boundaryDepth);

  //   // Forticate the stencil starts
  //   for (int iStencil=0; iStencil<sbpXX.numStencils; iStencil++)
  //     sbpXX.stencilStarts[iStencil]++;

  std::vector<double> realGridExtent;
  realGridExtent.resize(6,0.0);
  realGridExtent[1] = 1.0;
  realGridExtent[3] = 1.0;
  realGridExtent[5] = 1.0;

  std::vector<size_t> numNodes;
  numNodes.resize(3);
  numNodes[0] = 5;
  numNodes[1] = 6;
  numNodes[2] = 9;
  testGrid.SetGridSizes(numNodes);

  std::vector<double> dX;
  dX.resize(3);
  for(int iDim = 0;iDim < 3;iDim++)
    dX[iDim] = (realGridExtent[2*iDim+1] - realGridExtent[2*iDim])/static_cast<double>(numNodes[iDim] - 1);

  testGrid.SetGridSpacings(dX);

  interval_t &partitionInterval(testGrid.PartitionInterval());
  partitionInterval.InitSimple(numNodes);
  
  // Set up buffers for periodic boundaries
  std::vector<size_t> extendGrid(2*numDim,0);
//  std::vector<bool> isPeriodic(numDim,true); // periodic true here
  std::vector<bool> isPeriodic(numDim,false); // periodic false here
  std::vector<int> periodicDirs(numDim,0);
  size_t numPointsPart = partitionInterval.NNodes();
  for(int iDim = 0;iDim < numDim; iDim++){
    if(isPeriodic[iDim]){
      extendGrid[2*iDim]   = haloDepths[2*iDim];
      extendGrid[2*iDim+1] = haloDepths[2*iDim+1];
    } else {
      if(partitionInterval[iDim].first != 0)
        extendGrid[2*iDim] = haloDepths[2*iDim];
      if(partitionInterval[iDim].second != (numNodes[iDim]-1))
        extendGrid[2*iDim + 1] = haloDepths[2*iDim+1];
    }
  }
  std::cout << "Dimension Extensions: (";
  pcpp::io::DumpContents(std::cout,extendGrid,",");
  std::cout << ")" << std::endl;
  testGrid.SetDimensionExtensions(extendGrid);
  testGrid.Finalize();
  testHalo.SetPeriodicDirs(periodicDirs);


  size_t numPointsBuffer = testGrid.BufferSize();
  const interval_t &partitionBufferInterval(testGrid.PartitionBufferInterval());

  // set up the state
  state_t paramState;
  state_t testState;

  Maxwell::util::SetupMaxwellStateAndParameters(testGrid, testState, paramState);
  testState.InitializeFieldHandles();

  // set up another state by copying
  state_t testState2(testState);

//  std::ostringstream OutStream;
//  testState.Report(OutStream);


  // Populate the state
  std::vector<double> constE(3, 1.0);
  std::vector<double> constB(3, 1.0);
  std::vector<double> constJ(3, 1.0);
  double constEps = EPS0*1e12;
  double constMu = MU0*1e6;

  Maxwell::util::InitializeMaxwellStateConstFields(testGrid, testState, &constE[0], &constB[0], 
                                                   &constJ[0], constEps, constMu, false);
  
  // Populate the parameters
  double CFLValue = 1.0;
  Maxwell::util::InitializeMaxwellParameters(testGrid, paramState, CFLValue);


  /* ======== Test RHS initialization routine ======== */
  bool MaxwellRHSInitialize = true;
  if(testRHS.Initialize(testGrid,testState,paramState,sbpXX))
  {
      MaxwellRHSInitialize = false;
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Initialize", MaxwellRHSInitialize);


  /* ======== Test SetParamState() function ======== */
  //   bool MaxwellRHSSetParamState = true;
  //   if(testRHS.SetParamState(paramState))
  //   {
  //     MaxwellRHSSetParamState = false;
  //   }
  //   serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:SetParamState", MaxwellRHSSetParamState);


  /* ======== Test RHS setup routine ======== */
  const std::vector<double *> &stateBuffers(testRHS.StateBuffers());
  const std::vector<double *> &rhsBuffers(testRHS.RHSBuffers());

  /**** CALL 1 ****/
  bool MaxwellRHSSetup = true;
  if(testRHS.SetupRHS(testGrid,testState,testState2,myThreadID))
  {
      MaxwellRHSSetup = false;
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Setup:Call1", MaxwellRHSSetup);

  // Check state buffers
  std::vector<std::string> stateFieldNames;
  stateFieldNames.push_back("Efield");
  stateFieldNames.push_back("Bfield");
  stateFieldNames.push_back("Dfield");
  stateFieldNames.push_back("Hfield");
  stateFieldNames.push_back("J");
  stateFieldNames.push_back("epsilonMu");
  
  bool MaxwellRHSSetupStateBuffers = true;
  if(stateBuffers.size() != stateFieldNames.size())
  {
      MaxwellRHSSetupStateBuffers = false;
  }

  std::vector<double *>::const_iterator stateBuffersIt = stateBuffers.begin();
  std::vector<std::string>::iterator stateFieldNamesIt = stateFieldNames.begin();
  while(stateBuffersIt != stateBuffers.end()) {
    double *stateBuffer = *stateBuffersIt++;
    std::string &stateFieldName(*stateFieldNamesIt++);

    if(stateFieldName.compare("Efield")==0 || stateFieldName.compare("Hfield")==0) {
      if(testState.GetFieldBuffer<double>(stateFieldName) != stateBuffer) {
        MaxwellRHSSetupStateBuffers = false;
      }
    }
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Setup:Call1:CheckStateBuffers", MaxwellRHSSetupStateBuffers);

  // Check RHS buffers
  std::vector<std::string> rhsFieldNames;
  rhsFieldNames.push_back("Efield");
  rhsFieldNames.push_back("Hfield");

  bool MaxwellRHSSetupRHSBuffers = true;
  if(rhsBuffers.size() != rhsFieldNames.size())
  {
      MaxwellRHSSetupRHSBuffers = false;
  }

  std::vector<double *>::const_iterator rhsBuffersIt = rhsBuffers.begin();
  std::vector<std::string>::iterator rhsFieldNamesIt = rhsFieldNames.begin();
  while(rhsBuffersIt != rhsBuffers.end()) {
    double *rhsBuffer = *rhsBuffersIt++;
    std::string &rhsFieldName(*rhsFieldNamesIt++);

    if(testState2.GetFieldBuffer<double>(rhsFieldName) != rhsBuffer)
      MaxwellRHSSetupRHSBuffers = false;
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Setup:Call1:CheckRHSBuffers", MaxwellRHSSetupRHSBuffers);

  /**** CALL 2 (switch input and output buffers) ****/
  bool MaxwellRHSSetup2 = true;
  if(testRHS.SetupRHS(testGrid,testState2,testState,myThreadID))
  {
      MaxwellRHSSetup2 = false;
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Setup:Call2", MaxwellRHSSetup2);
  
  // Check state buffers
  std::vector<std::string> stateFieldNames2;
  stateFieldNames2.push_back("Efield");
  stateFieldNames2.push_back("Bfield");
  stateFieldNames2.push_back("Dfield");
  stateFieldNames2.push_back("Hfield");
  stateFieldNames2.push_back("J");
  stateFieldNames2.push_back("epsilonMu");

  bool MaxwellRHSSetupStateBuffers2 = true;
  if(stateBuffers.size() != stateFieldNames2.size())
  {
    MaxwellRHSSetupStateBuffers2 = false;
  }

  std::vector<double *>::const_iterator stateBuffersIt2 = stateBuffers.begin();
  std::vector<std::string>::iterator stateFieldNamesIt2 = stateFieldNames2.begin();
  while(stateBuffersIt2 != stateBuffers.end()) {
    double *stateBuffer2 = *stateBuffersIt2++;
    std::string &stateFieldName2(*stateFieldNamesIt2++);

    if(stateFieldName2.compare("Efield")==0 || stateFieldName2.compare("Hfield")==0) {
      if(testState2.GetFieldBuffer<double>(stateFieldName2) != stateBuffer2) {
        MaxwellRHSSetupStateBuffers2 = false;
      }
    }
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Setup:Call2:CheckStateBuffers", MaxwellRHSSetupStateBuffers2);

  // Check RHS buffers
  std::vector<std::string> rhsFieldNames2;
  rhsFieldNames2.push_back("Efield");
  rhsFieldNames2.push_back("Hfield");

  bool MaxwellRHSSetupRHSBuffers2 = true;
  if(rhsBuffers.size() != rhsFieldNames2.size())
  {
    MaxwellRHSSetupRHSBuffers2 = false;
  }

  std::vector<double *>::const_iterator rhsBuffersIt2 = rhsBuffers.begin();
  std::vector<std::string>::iterator rhsFieldNamesIt2 = rhsFieldNames2.begin();
  while(rhsBuffersIt2 != rhsBuffers.end()) {
    double *rhsBuffer2 = *rhsBuffersIt2++;
    std::string &rhsFieldName2(*rhsFieldNamesIt2++);

    if(testState.GetFieldBuffer<double>(rhsFieldName2) != rhsBuffer2)
      MaxwellRHSSetupRHSBuffers2 = false;
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Setup:Call2:CheckRHSBuffers", MaxwellRHSSetupRHSBuffers2);


  /* ======== Test MaxwellRHS() function ======== */
  int numTests = 2;
  std::vector<bool> MaxwellRHSMaxwellRHS(numTests, true);

  /**** CASE 1: Constant epsilon, mu, E, B and J ****/

  // Set up RHS
  if(testRHS.SetupRHS(testGrid,testState,testState2,myThreadID))
  {
    std::cout << "SetupRHS failed." << std::endl;
    MaxwellRHSMaxwellRHS[0] = false;
  }

  std::vector<double *>::const_iterator rhsIt = rhsBuffers.begin();
  while(rhsIt != rhsBuffers.end()) {
    double *rhsBuffer = *rhsIt++;
    for(size_t iPoint = 0;iPoint < numDim*numPointsBuffer;iPoint++)
      rhsBuffer[iPoint] = 0.0;
  }

  // Compute Maxwell RHS
  if(testRHS.MaxwellRHS(myThreadID))
  {
    std::cout << "MaxwellRHS failed." << std::endl;
    MaxwellRHSMaxwellRHS[0] = false;
  }  

  // Execute test
  size_t nFailE = 0;
  size_t nFailH = 0;
  double *dE = rhsBuffers[0];
  double *dH = rhsBuffers[1];
  double *J     = testState.GetFieldBuffer<double>("J");
  double *epsMu = testState.GetFieldBuffer<double>("epsilonMu");
  for(int iDim=0; iDim<numDim; iDim++) {
    size_t iDimxnumPointsBuffer = iDim*numPointsBuffer;
    for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
      if(std::fabs(dH[iDimxnumPointsBuffer+iPoint]) > MACHINE_EPS) {
        nFailH++;
        MaxwellRHSMaxwellRHS[0] = false;
        std::cout << "H(" << iDim << "):" << dH[iDimxnumPointsBuffer+iPoint] << std::endl;
      }

      if(std::fabs(dE[iDimxnumPointsBuffer+iPoint]+J[iDimxnumPointsBuffer+iPoint]/epsMu[iPoint]) > MACHINE_EPS) {
        nFailE++;
        MaxwellRHSMaxwellRHS[0] = false;
        std::cout << "E(" << iDim << "):" << dE[iDimxnumPointsBuffer+iPoint] << std::endl;
      }
    }
    if(nFailE > 0) {
      std::cout << nFailE << "/" << numPointsBuffer << " E failures in dimension " << iDim << std::endl;
    }
    if(nFailH > 0) {
      std::cout << nFailH << "/" << numPointsBuffer << " H failures in dimension " << iDim << std::endl;
    }
    nFailE = 0;
    nFailH = 0;
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Execution:MaxwellRHS:Case1", MaxwellRHSMaxwellRHS[0]);

  /**** CASE 2: Constant epsilon, mu and spatially-varying E, B and J ****/

  // Re-populate fields
  double *E = testState.GetFieldBuffer<double>("Efield");
  double *D = testState.GetFieldBuffer<double>("Dfield");
  double curlCompE2[3];
  CreateConstCurlVectorField(numDim, &numNodes[0], &dX[0], &E[0], curlCompE2);

  double *B = testState.GetFieldBuffer<double>("Bfield");
  double *H = testState.GetFieldBuffer<double>("Hfield");
  double curlCompB2[3];
  CreateConstCurlVectorField(numDim, &numNodes[0], &dX[0], &B[0], curlCompB2);

  double curlCompJ2[3];
  CreateConstCurlVectorField(numDim, &numNodes[0], &dX[0], &J[0], curlCompJ2);

  size_t twoxnumPointsBuffer = 2*numPointsBuffer;
  for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
    D[iPoint]                     = epsMu[iPoint]*E[iPoint];
    D[iPoint+numPointsBuffer]     = epsMu[iPoint]*E[iPoint+numPointsBuffer];
    D[iPoint+twoxnumPointsBuffer] = epsMu[iPoint]*E[iPoint+twoxnumPointsBuffer];

    H[iPoint]                     = B[iPoint]/epsMu[iPoint+numPointsBuffer];
    H[iPoint+numPointsBuffer]     = B[iPoint+numPointsBuffer]/epsMu[iPoint+numPointsBuffer];
    H[iPoint+twoxnumPointsBuffer] = B[iPoint+twoxnumPointsBuffer]/epsMu[iPoint+numPointsBuffer];
  }

  // Set up RHS
  if(testRHS.SetupRHS(testGrid,testState,testState2,myThreadID))
  {
    MaxwellRHSMaxwellRHS[1] = false;
  }

  // Compute Maxwell RHS
  if(testRHS.MaxwellRHS(myThreadID))
  {
    MaxwellRHSMaxwellRHS[1] = false;
  }

  // Execute test
  for(int iDim=0; iDim<numDim; iDim++) {
    size_t iDimxnumPointsBuffer = iDim*numPointsBuffer;
    for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
      if(std::fabs(dH[iDimxnumPointsBuffer+iPoint]+curlCompE2[iDim]/epsMu[iPoint+numPointsBuffer]) > MACHINE_EPS)
        MaxwellRHSMaxwellRHS[1] = false;

      double epsxMu = epsMu[iPoint]*epsMu[numPointsBuffer+iPoint];
      if(std::fabs(dE[iDimxnumPointsBuffer+iPoint]+J[iDimxnumPointsBuffer+iPoint]/epsMu[iPoint]-curlCompB2[iDim]/epsxMu) > MACHINE_EPS)
        MaxwellRHSMaxwellRHS[1] = false;
    }
  }

  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Execution:MaxwellRHS:Case2", MaxwellRHSMaxwellRHS[1]);


  /* ======== Test RHS routine ======== */
  bool MaxwellRHSRHS = true;
  if(testRHS.RHS(myThreadID))
  {
    MaxwellRHSRHS = false;
  }
  serialUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:Execution:RHS", MaxwellRHSRHS);

}


void TestMaxwellRHSTimeIntegrate(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
/*********************************************************************************************************

        Test functions that set up the simulation for time integration,
        including initialization of the states.

 *********************************************************************************************************/
{

  /* ======== Set up global type ======== */
  global_t      testGlobal;

  const std::string testName("TestMaxwellRHSTimeIntegrate");
  testGlobal.Init(testName, testComm);
  testGlobal.StdOut("Entering "+testName+"\n");
 
  /* ======== Set up grid, halo and operator ======== */
  grid_t        testGrid;
  halo_t        &testHalo(testGrid.Halo());
  operator_t    testOperator;

  std::ostringstream messageStream;

  int opInteriorOrder = 4;
  std::vector<size_t> gridSizes(3, 21);
  testGrid.SetGridSizes(gridSizes);
  std::vector<double> physicalExtent(6,0.2);
  physicalExtent[0] = 0.0;
  physicalExtent[2] = 0.0;
  physicalExtent[4] = 0.0;
  testGrid.SetPhysicalExtent(physicalExtent);

  bool MaxwellRHSTimeIntegrateGridHaloOpSetup = true;
  if(euler::util::InitializeSimulationFixtures(testGrid, testOperator, opInteriorOrder, 
                                               testComm, &messageStream))
  {
    MaxwellRHSTimeIntegrateGridHaloOpSetup = false;
  }

  MaxwellRHSTimeIntegrateGridHaloOpSetup = CheckResult(MaxwellRHSTimeIntegrateGridHaloOpSetup, 
                                                       testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:GridHaloOpSetup", MaxwellRHSTimeIntegrateGridHaloOpSetup);
  
  testGlobal.StdOut(messageStream.str());
  pcpp::io::RenewStream(messageStream);

  /* ======== Set up states and parameters ======== */
  state_t paramState;
  state_t testState;

  bool MaxwellRHSTimeIntegrateStateParamsSetup = true;
  if(Maxwell::util::SetupMaxwellStateAndParameters(testGrid, testState, paramState))
  {
    MaxwellRHSTimeIntegrateStateParamsSetup = false;
  }
  testState.InitializeFieldHandles();
  MaxwellRHSTimeIntegrateStateParamsSetup = CheckResult(MaxwellRHSTimeIntegrateStateParamsSetup, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:StateParamsSetup", MaxwellRHSTimeIntegrateStateParamsSetup);

  testGlobal.StdOut("Creating states...\n");

  state_t testRHSState(testState);
  
  /* ======== Initialize states - constant fields ======== */
  size_t numPointsBuffer = testGrid.BufferSize();
//  const std::vector<size_t> &gridSizes(testGrid.GridSizes());
  const std::vector<size_t> &bufferSizes(testGrid.BufferSizes());
  const pcpp::IndexIntervalType &partitionInterval(testGrid.PartitionInterval());
  const pcpp::IndexIntervalType &partitionBufferInterval(testGrid.PartitionBufferInterval());
  int numDim = bufferSizes.size();

  // Get state buffers
  double *timePtr       = testState.GetFieldBuffer<double>("simTime");
  double *EfieldPtr     = testState.GetFieldBuffer<double>("Efield");
  double *BfieldPtr     = testState.GetFieldBuffer<double>("Bfield");
  double *DfieldPtr     = testState.GetFieldBuffer<double>("Dfield");
  double *HfieldPtr     = testState.GetFieldBuffer<double>("Hfield");
  double *JPtr          = testState.GetFieldBuffer<double>("J");
  double *epsilonMuPtr  = testState.GetFieldBuffer<double>("epsilonMu");

  /**** Initialize in entire buffer ****/
  std::vector<double> EfieldAll(3, 0.0);
  std::vector<double> BfieldAll(3, 0.0);
  std::vector<double> JAll(3, 0.0);
  double epsilonAll = -1.0;
  double muAll   = -1.0;

  testGlobal.StdOut("Initializing states...\n");
  bool MaxwellRHSTimeIntegrateStateInitEverywhere = true;
  if(Maxwell::util::InitializeMaxwellStateConstFields(testGrid, testState, &EfieldAll[0], &BfieldAll[0], &JAll[0], epsilonAll, muAll, true))
  {
    MaxwellRHSTimeIntegrateStateInitEverywhere = false;
  }
  MaxwellRHSTimeIntegrateStateInitEverywhere = CheckResult(MaxwellRHSTimeIntegrateStateInitEverywhere, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:StateInit:Everywhere", MaxwellRHSTimeIntegrateStateInitEverywhere);

  // Test initialization
  size_t twoxnumPointsBuffer = 2*numPointsBuffer;
  for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
    if(std::fabs(EfieldPtr[iPoint]                    -EfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(EfieldPtr[iPoint+numPointsBuffer]    -EfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(EfieldPtr[iPoint+twoxnumPointsBuffer]-EfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }

    if(std::fabs(BfieldPtr[iPoint]                    -BfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(BfieldPtr[iPoint+numPointsBuffer]    -BfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(BfieldPtr[iPoint+twoxnumPointsBuffer]-BfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }

    if(std::fabs(DfieldPtr[iPoint]                    -epsilonAll*EfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(DfieldPtr[iPoint+numPointsBuffer]    -epsilonAll*EfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(DfieldPtr[iPoint+twoxnumPointsBuffer]-epsilonAll*EfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }

    if(std::fabs(HfieldPtr[iPoint]                    -BfieldAll[0]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(HfieldPtr[iPoint+numPointsBuffer]    -BfieldAll[1]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(HfieldPtr[iPoint+twoxnumPointsBuffer]-BfieldAll[2]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }

    if(std::fabs(JPtr[iPoint]                    -JAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(JPtr[iPoint+numPointsBuffer]    -JAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(JPtr[iPoint+twoxnumPointsBuffer]-JAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }

    if(std::fabs(epsilonMuPtr[iPoint]-epsilonAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
    if(std::fabs(epsilonMuPtr[iPoint+numPointsBuffer]-muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitEverywhere = false; }
  }
  MaxwellRHSTimeIntegrateStateInitEverywhere = CheckResult(MaxwellRHSTimeIntegrateStateInitEverywhere, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:StateInit:Everywhere:CheckBuffers", MaxwellRHSTimeIntegrateStateInitEverywhere);

  /**** Initialize only in local partition ****/
  std::vector<double> Efield(3,  1.0); Efield[1] = -2.0; Efield[2] =  3.0;
  std::vector<double> Bfield(3, -1.0); Bfield[1] =  2.0; Bfield[2] = -3.0;
  std::vector<double> J(3, 4.0);       J[1]      = -5.0; J[2]      =  6.0;
  double epsilon = EPS0;
  double mu      = MU0;

  bool MaxwellRHSTimeIntegrateStateInitLocalPart = true;
  if(Maxwell::util::InitializeMaxwellStateConstFields(testGrid, testState, &Efield[0], &Bfield[0], &J[0], epsilon, mu, false))
  {
    MaxwellRHSTimeIntegrateStateInitLocalPart = false;
  }
  MaxwellRHSTimeIntegrateStateInitLocalPart = CheckResult(MaxwellRHSTimeIntegrateStateInitLocalPart, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:StateInit:LocalPartition", MaxwellRHSTimeIntegrateStateInitLocalPart);

  // Test initialization
  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 = partitionBufferInterval[2].first;
  size_t kEnd   = partitionBufferInterval[2].second;

  size_t nPlane = bufferSizes[0]*bufferSizes[1];

  for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
    size_t iIndex = static_cast<size_t>((iPoint%nPlane)%bufferSizes[0]);
    size_t jIndex = static_cast<size_t>(((iPoint%nPlane)-iIndex)/bufferSizes[0]);
    size_t kIndex = static_cast<size_t>((iPoint-jIndex*bufferSizes[0]-iIndex)/nPlane);

    if(iIndex<iStart || iIndex>iEnd || jIndex<jStart || jIndex>jEnd || kIndex<kStart || kIndex>kEnd) {
      // Halo buffers
      if(std::fabs(EfieldPtr[iPoint]                    -EfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(EfieldPtr[iPoint+numPointsBuffer]    -EfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(EfieldPtr[iPoint+twoxnumPointsBuffer]-EfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(BfieldPtr[iPoint]                    -BfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(BfieldPtr[iPoint+numPointsBuffer]    -BfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(BfieldPtr[iPoint+twoxnumPointsBuffer]-BfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(DfieldPtr[iPoint]                    -epsilonAll*EfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(DfieldPtr[iPoint+numPointsBuffer]    -epsilonAll*EfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(DfieldPtr[iPoint+twoxnumPointsBuffer]-epsilonAll*EfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(HfieldPtr[iPoint]                    -BfieldAll[0]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(HfieldPtr[iPoint+numPointsBuffer]    -BfieldAll[1]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(HfieldPtr[iPoint+twoxnumPointsBuffer]-BfieldAll[2]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(JPtr[iPoint]                    -JAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(JPtr[iPoint+numPointsBuffer]    -JAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(JPtr[iPoint+twoxnumPointsBuffer]-JAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(epsilonMuPtr[iPoint]-epsilonAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(epsilonMuPtr[iPoint+numPointsBuffer]-muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
    } else {
      // Partition buffers
      if(std::fabs(EfieldPtr[iPoint]                    -Efield[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(EfieldPtr[iPoint+numPointsBuffer]    -Efield[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(EfieldPtr[iPoint+twoxnumPointsBuffer]-Efield[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(BfieldPtr[iPoint]                    -Bfield[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(BfieldPtr[iPoint+numPointsBuffer]    -Bfield[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(BfieldPtr[iPoint+twoxnumPointsBuffer]-Bfield[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(DfieldPtr[iPoint]                    -epsilon*Efield[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(DfieldPtr[iPoint+numPointsBuffer]    -epsilon*Efield[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(DfieldPtr[iPoint+twoxnumPointsBuffer]-epsilon*Efield[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(HfieldPtr[iPoint]                    -Bfield[0]/mu) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(HfieldPtr[iPoint+numPointsBuffer]    -Bfield[1]/mu) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(HfieldPtr[iPoint+twoxnumPointsBuffer]-Bfield[2]/mu) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(JPtr[iPoint]                    -J[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(JPtr[iPoint+numPointsBuffer]    -J[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(JPtr[iPoint+twoxnumPointsBuffer]-J[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }

      if(std::fabs(epsilonMuPtr[iPoint]-epsilon) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
      if(std::fabs(epsilonMuPtr[iPoint+numPointsBuffer]-mu) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart = false; }
    }
  }
  MaxwellRHSTimeIntegrateStateInitLocalPart = CheckResult(MaxwellRHSTimeIntegrateStateInitLocalPart, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:StateInit:LocalPartition:CheckBuffers", MaxwellRHSTimeIntegrateStateInitLocalPart);

//  /* ======== Re-initialize states - using grid indices ======== */
//
//  /**** Initialize only in local partition ****/
//  bool MaxwellRHSTimeIntegrateStateInitLocalPart2 = true;
//  if(Maxwell::util::InitializeMaxwellStateGridIndices(testGrid, testState))
//  {
//    MaxwellRHSTimeIntegrateStateInitLocalPart2 = false;
//  }
//  MaxwellRHSTimeIntegrateStateInitLocalPart2 = CheckResult(MaxwellRHSTimeIntegrateStateInitLocalPart2, testComm);
//  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:StateInit:LocalPartition2", MaxwellRHSTimeIntegrateStateInitLocalPart2);
//
//  // Test initialization
//  size_t iStartG = partitionInterval[0].first;
//  size_t iEndG   = partitionInterval[0].second;
//  size_t jStartG = partitionInterval[1].first;
//  size_t jEndG   = partitionInterval[1].second;
//  size_t kStartG = partitionInterval[2].first;
//  size_t kEndG   = partitionInterval[2].second;
//
//  size_t nPlaneG = gridSizes[0]*gridSizes[1];
//
//  for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
//    size_t iIndex = static_cast<size_t>((iPoint%nPlane)%bufferSizes[0]);
//    size_t jIndex = static_cast<size_t>(((iPoint%nPlane)-iIndex)/bufferSizes[0]);
//    size_t kIndex = static_cast<size_t>((iPoint-jIndex*bufferSizes[0]-iIndex)/nPlane);
//
//    if(iIndex<iStart || iIndex>iEnd || jIndex<jStart || jIndex>jEnd || kIndex<kStart || kIndex>kEnd) {
//      // Halo buffers
//      if(std::fabs(EfieldPtr[iPoint]                  -EfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(EfieldPtr[iPoint+numPointsBuffer]    -EfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(EfieldPtr[iPoint+twoxnumPointsBuffer]-EfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(BfieldPtr[iPoint]                  -BfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(BfieldPtr[iPoint+numPointsBuffer]    -BfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(BfieldPtr[iPoint+twoxnumPointsBuffer]-BfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(DfieldPtr[iPoint]                  -epsilonAll*EfieldAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(DfieldPtr[iPoint+numPointsBuffer]    -epsilonAll*EfieldAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(DfieldPtr[iPoint+twoxnumPointsBuffer]-epsilonAll*EfieldAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(HfieldPtr[iPoint]                  -BfieldAll[0]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(HfieldPtr[iPoint+numPointsBuffer]    -BfieldAll[1]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(HfieldPtr[iPoint+twoxnumPointsBuffer]-BfieldAll[2]/muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(JPtr[iPoint]                          -JAll[0]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(JPtr[iPoint+numPointsBuffer]    -JAll[1]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(JPtr[iPoint+twoxnumPointsBuffer]-JAll[2]) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(epsilonMuPtr[iPoint]-epsilonAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(epsilonMuPtr[iPoint+numPointsBuffer]-muAll) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//    } else {
//      // Partition buffers
//      size_t gridIndex = (kStartG+kIndex-kStart)*nPlaneG + (jStartG+jIndex-jStart)*gridSizes[0] + (iStartG+iIndex-iStart+1);
//
//      if(std::fabs(EfieldPtr[iPoint]                  -gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(EfieldPtr[iPoint+numPointsBuffer]    -gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(EfieldPtr[iPoint+twoxnumPointsBuffer]-gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(BfieldPtr[iPoint]                  -gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(BfieldPtr[iPoint+numPointsBuffer]    -gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(BfieldPtr[iPoint+twoxnumPointsBuffer]-gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(DfieldPtr[iPoint]                  -gridIndex*gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(DfieldPtr[iPoint+numPointsBuffer]    -gridIndex*gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(DfieldPtr[iPoint+twoxnumPointsBuffer]-gridIndex*gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(HfieldPtr[iPoint]                  -1.0) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(HfieldPtr[iPoint+numPointsBuffer]    -1.0) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(HfieldPtr[iPoint+twoxnumPointsBuffer]-1.0) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(JPtr[iPoint]                          -gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(JPtr[iPoint+numPointsBuffer]    -gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(JPtr[iPoint+twoxnumPointsBuffer]-gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//
//      if(std::fabs(epsilonMuPtr[iPoint]-gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//      if(std::fabs(epsilonMuPtr[iPoint+numPointsBuffer]-gridIndex) > MACHINE_EPS) { MaxwellRHSTimeIntegrateStateInitLocalPart2 = false; }
//    }
//  }
//  MaxwellRHSTimeIntegrateStateInitLocalPart2 = CheckResult(MaxwellRHSTimeIntegrateStateInitLocalPart2, testComm);
//  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:StateInit:LocalPartition2:CheckBuffers", MaxwellRHSTimeIntegrateStateInitLocalPart2);

  /* ======== Initialize parameters ======== */
  double CFLValue = 1.0;

  bool MaxwellRHSTimeIntegrateParamsInit = true;
  if(Maxwell::util::InitializeMaxwellParameters(testGrid, paramState, CFLValue))
  {
    MaxwellRHSTimeIntegrateParamsInit = false;
  }
  MaxwellRHSTimeIntegrateParamsInit = CheckResult(MaxwellRHSTimeIntegrateParamsInit, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:ParamsInit", MaxwellRHSTimeIntegrateParamsInit);

  /* ======== Set up domain and RHS ======== */
  typedef plascom2::operators::sbp::base                     operator_t;
  typedef simulation::domain::base<grid_t,state_t>  domain_t;
  typedef Maxwell::rhs<grid_t,state_t,operator_t>            rhs_t;

  domain_t      testDomain;
  rhs_t         testRHS;

  testDomain.SetNumGrids(1);
  testDomain.SetGrid(0,testGrid);
  testDomain.SetGridState(0,testState);
  testDomain.SetGridParams(0,paramState);
  testDomain.PartitionDomain(0);
  testDomain.SetRHS(testRHS);

  bool MaxwellRHSTimeIntegrateInitialize = true;
  if(testRHS.Initialize(testGrid,testState,paramState,testOperator))
  {
    MaxwellRHSTimeIntegrateInitialize = false;
  }
  MaxwellRHSTimeIntegrateInitialize = CheckResult(MaxwellRHSTimeIntegrateInitialize, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:Initialize", MaxwellRHSTimeIntegrateInitialize);

  /* ======== Set up time advancer ======== */
  rk4advancer<domain_t> testAdvancer;
  testAdvancer.SetCommunication(true);
  testAdvancer.InitializeAdvancer(testDomain,testGlobal,messageStream);

  std::cout << "Initialization done." << std::endl;

  testGlobal.StdOut(messageStream.str());
  pcpp::io::RenewStream(messageStream);

  int myThreadID = 0;
  int numThreads = 1;
  int masterThread = 1;

#ifdef USE_OMP
  myThreadID = omp_get_thread_num();
  numThreads = omp_get_max_threads();
  masterThread = (myThreadID == 0);
#endif

  testAdvancer.SyncIntervals();
  //  testAdvancer.SyncStates();


  /* ======== Set up I/O ======== */
#ifdef ENABLE_HDF5
  fixtures::ConfigurationType   testConfig;
  std::ostringstream            outputFileNameOut;
  std::string                   outputFileNameRoot("TestMaxwellRHSTimeIntegrate");
  std::string                   domainName("TestDomain");
  std::string                   geometryName("TestGeom");
  std::string                   gridName("TestGrid");

  // Create directory to store test output files
  std::string   outputFileDirName("OutputFiles_TestMaxwellRHSTimeIntegrate");
  if(!ix::sys::FILEEXISTS(outputFileDirName)){
    if(ix::sys::CreateDirectory(outputFileDirName))
      testGlobal.ErrOut("ix::sys::CreateDirectory failed.\n");
  }
#endif

  /* ======== Advance in time ======== */
  size_t iStep;
  size_t numStepsMax    = 100;
  size_t numStepsStatus = 1;
  size_t numStepsIO     = 1;
  bool MaxwellRHSTimeIntegrateTimeAdvance = true;

  for(iStep=0; iStep<numStepsMax; iStep++) {

    if(!(iStep%numStepsStatus)) {
      if(masterThread) {
        messageStream << "MaxwellSolver:TestMaxwellRHS:TimeIntegration:TimeAdvance - Step = " 
                      << iStep << std::endl;
        testGlobal.StdOut(messageStream.str(),2);
        pcpp::io::RenewStream(messageStream);
      }
    }

    // Perform I/O operations at every numStepsIO
    if(!(iStep%numStepsIO)) {
#ifdef ENABLE_HDF5
      outputFileNameOut.str("");
      outputFileNameOut << outputFileDirName << "/" << outputFileNameRoot << "_tstep" << iStep << ".h5";

      std::string outputFileName(outputFileNameOut.str());
      if(plascom2::io::hdf5::OutputSingle(outputFileName,domainName,geometryName,gridName,
                                          testGrid,testState,paramState,testConfig,messageStream)) {
        testGlobal.ErrOut("Maxwell:TestMaxwellRHS:TimeIntegration:TimeAdvance - plascom2::io::hdf5::OutputSingle() failed.\n");
      }
      testGlobal.StdOut(messageStream.str(),3);
      pcpp::io::RenewStream(messageStream);
#else
      testGlobal.ErrOut("Maxwell:TestMaxwellRHS:TimeIntegration:TimeAdvance - WARNING! > HDF5 not enabled. Cannot perform I/O.\n");
#endif
    }
    
    testComm.Barrier();
    testGlobal.StdOut("Advancing domain.\n");
    
    if(testAdvancer.AdvanceDomain()) {
      MaxwellRHSTimeIntegrateTimeAdvance = false;
    }

    testComm.Barrier();
    testGlobal.StdOut("Advance done.\nComputing DV.");

    if(!MaxwellRHSTimeIntegrateTimeAdvance)
      break;

    if(testRHS.ComputeDV(myThreadID))
      testGlobal.ErrOut("Maxwell:TestMaxwellRHS:TimeIntegration:TimeAdvance - ComputeDV() failed.\n");

    testComm.Barrier();
    testGlobal.StdOut("Compute DV done.\n");

    if(iStep == 1)
      testAdvancer.SetTimers(true);

  }
  MaxwellRHSTimeIntegrateTimeAdvance = CheckResult(MaxwellRHSTimeIntegrateTimeAdvance, testComm);
  parallelUnitResults.UpdateResult("Maxwell:TestMaxwellRHS:TimeIntegration:TimeAdvance", MaxwellRHSTimeIntegrateTimeAdvance);

  testComm.Barrier();
  testGlobal.StdOut("Performing final I/O.\n");

  // Perform final I/O dump
  if(masterThread) {
#ifdef ENABLE_HDF5
      outputFileNameOut.str("");
      outputFileNameOut << outputFileDirName << "/" << outputFileNameRoot << "_tstep" << iStep << ".h5";

      std::string outputFileName(outputFileNameOut.str());
      if(plascom2::io::hdf5::OutputSingle(outputFileName,domainName,geometryName,gridName,
                                          testGrid,testState,paramState,testConfig,messageStream)) {
        testGlobal.ErrOut("Maxwell:TestMaxwellRHS:TimeIntegration:TimeAdvance - plascom2::io::hdf5::OutputSingle() failed.\n");
      }
      testGlobal.StdOut(messageStream.str(),3);
      pcpp::io::RenewStream(messageStream);
#else
      testGlobal.ErrOut("Maxwell:TestMaxwellRHS:TimeIntegration:TimeAdvance - WARNING! > HDF5 not enabled. Cannot perform I/O.\n");
#endif
  }

  if(masterThread) {
    messageStream << "Maxwell:TestMaxwellRHS:TimeIntegration:TimeAdvance - Ending at step " << iStep << std::endl;
    testGlobal.StdOut(messageStream.str(),1);
    pcpp::io::RenewStream(messageStream);
  }

}