MaxwellSolver.C

Go to the documentation of this file.

#include <iostream>
#include "OperatorKernels.H"
#include "MaxwellSolver.H"

#define DEBUG_MAXWELL_RHS_BFIELD        0
#define DEBUG_MAXWELL_RHS_DFIELD        0
#define DEBUG_MAXWELL_ComputeRecipEpsMu 0
#define DEBUG_MAXWELL_CONVERT_ETODFIELD 0
#define DEBUG_MAXWELL_CONVERT_DTOEFIELD 0
#define DEBUG_MAXWELL_CONVERT_BTOHFIELD 0
#define DEBUG_MAXWELL_CONVERT_HTOBFIELD 0



namespace Maxwell {


  int ComputeCurl(int numDim, size_t *numX, int numComponents, 
                  size_t numPoints, size_t *opInterval, int *stencilID,
                  plascom2::operators::sbp::base &inOperator,
                  double *recipDeltaX, double *inputField, double *outputField)
  {

      /* Check no. of dimensions - it must be 3. */
      if (numDim != 3) {
        return 1;
      }

      int comp1, comp2;
      int comp1plusOne, comp2plusOne;
      double *partDDx1 = new double [numPoints];
      double *partDDx2 = new double [numPoints];

      for (int compCurl=0; compCurl<numDim; compCurl++) {

        comp1 = (compCurl+1)%3;
        comp2 = (comp1+1)%3;

        comp1plusOne = comp1+1;
        comp2plusOne = comp2+1;

        size_t comp1xnumPoints = comp1*numPoints;
        size_t comp2xnumPoints = comp2*numPoints;

        FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
          (&numDim, numX, &numComponents, &numPoints,
           &comp1plusOne, opInterval, &(inOperator.numStencils), inOperator.stencilSizes,
           inOperator.stencilStarts, &(inOperator.numValues),
           inOperator.stencilWeights, inOperator.stencilOffsets,
           &stencilID[comp1xnumPoints],
           &inputField[comp2xnumPoints], partDDx1);
        
        FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
          (&numDim, numX, &numComponents, &numPoints,
           &comp2plusOne, opInterval, &(inOperator.numStencils), inOperator.stencilSizes,
           inOperator.stencilStarts, &(inOperator.numValues),
           inOperator.stencilWeights, inOperator.stencilOffsets,
           &stencilID[comp2xnumPoints],
           &inputField[comp1xnumPoints], partDDx2);

        size_t compCurlxnumPoints = compCurl*numPoints;

        for (size_t iPoint=0; iPoint<numPoints; iPoint++) {

          partDDx1[iPoint] *= recipDeltaX[comp1xnumPoints+iPoint];
          partDDx2[iPoint] *= recipDeltaX[comp2xnumPoints+iPoint];

          outputField[compCurlxnumPoints+iPoint] = partDDx1[iPoint] - partDDx2[iPoint];

        }

      }

      delete [] partDDx1;
      delete [] partDDx2;

      return 0;

  } // end ComputeCurl()


  int ComputeMaxwellRHS_Bfield(int numDim, size_t *numX, int numComponents, size_t numPoints,
                               size_t *opInterval, int *stencilID,
                               plascom2::operators::sbp::base &inOperator,
                               double *recipDeltaX, double *Efield, double *dBdt)
  /******************************************************************************************

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

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

    if (ComputeCurl(numDim, numX, numComponents, numPoints, opInterval, stencilID,
                    inOperator, recipDeltaX, Efield, dBdt))
    {
      return 1;
    }

    double negOne = -1.0;

    for (int iDim=0; iDim<numDim; iDim++) {
      FC_MODULE(operators,xax,OPERATORS,XAX)(&numDim, &numPoints, numX, opInterval, &negOne, &dBdt[iDim*numPoints]);
    }

#if DEBUG_MAXWELL_RHS_BFIELD
    // FOR DEBUGGING PURPOSES
    for (int iDim=0; iDim<numDim; iDim++) {
      size_t iDimxnumPoints = iDim*numPoints;
      for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
        std::cout << "[Maxwell RHS B-field] dim=" << iDim << ", point=" << iPoint << ", Efield=" << Efield[iDimxnumPoints+iPoint] << ", dBdt=" << dBdt[iDimxnumPoints+iPoint] << std::endl;
      }
    }
#endif

    return 0;

  } // end ComputeMaxwellRHS_Bfield()


  int ComputeMaxwellRHS_Dfield(int numDim, size_t *numX, int numComponents, size_t numPoints,
                               size_t *opInterval, int *stencilID,
                               plascom2::operators::sbp::base &inOperator,
                               double *recipDeltaX, double *Hfield, double *J, double *dDdt)
  /******************************************************************************************

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

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

    if (ComputeCurl(numDim, numX, numComponents, numPoints, opInterval, stencilID,
                    inOperator, recipDeltaX, Hfield, dDdt))
    {
      return 1;
    }

    double negOne = -1.0;

    for (int iDim=0; iDim<numDim; iDim++) {
      FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)(&numDim, &numPoints, numX, opInterval, &negOne, &J[iDim*numPoints], &dDdt[iDim*numPoints]);
    }

#if DEBUG_MAXWELL_RHS_DFIELD
    // FOR DEBUGGING PURPOSES
    for (int iDim=0; iDim<numDim; iDim++) {
      size_t iDimxnumPoints = iDim*numPoints;
      for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
        std::cout << "[Maxwell RHS D-field] dim=" << iDim << ", point=" << iPoint << ", Hfield=" << Hfield[iDimxnumPoints+iPoint] << ", J=" << J[iDimxnumPoints+iPoint] << ", dDdt=" << dDdt[iDimxnumPoints+iPoint] << std::endl;
      }
    }
#endif

    return 0;

  } // end ComputeMaxwellRHS_Dfield()


  int ComputeRecipEpsMu(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                        double *epsMu, double *recipEpsMu)
  /************************************************************************************

        Compute reciprocals of epsilon and mu.

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

    if (numDim != 3) {
      return 1;
    }

    double one = 1.0;

    FC_MODULE(operators,yaxm1,OPERATORS,YAXM1)(&numDim, &numPoints, numX, opInterval, &one, &epsMu[0], &recipEpsMu[0]);
    FC_MODULE(operators,yaxm1,OPERATORS,YAXM1)(&numDim, &numPoints, numX, opInterval, &one, &epsMu[numPoints], &recipEpsMu[numPoints]);

#if DEBUG_MAXWELL_ComputeRecipEpsMu
    // FOR DEBUGGING PURPOSES
    for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
      std::cout << "[Maxwell Compute reciprocal eps, mu] point=" << iPoint << ", eps=" << epsMu[iPoint] << ", 1/eps=" << recipEpsMu[iPoint] << ", mu=" << epsMu[numPoints+iPoint] << ", 1/mu=" << recipEpsMu[numPoints+iPoint] << std::endl;
    }
#endif

    return 0;

  } // end ComputeRecipEpsMu()


  int ConvertEfieldtoDfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *epsMu, double *Efield, double *Dfield)
  /****************************************************************************************

        Calculate D = epsilon*E

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

    if (numDim != 3) {
      return 1;
    }

    for (int iDim=0; iDim<numDim; iDim++) {
      FC_MODULE(operators,zxy,OPERATORS,ZXY)(&numDim, &numPoints, numX, opInterval,
                                             &epsMu[0], &Efield[iDim*numPoints], &Dfield[iDim*numPoints]);
    }

#if DEBUG_MAXWELL_CONVERT_ETODFIELD
    // FOR DEBUGGING PURPOSES
    for (int iDim=0; iDim<numDim; iDim++) {
      size_t iDimxnumPoints = iDim*numPoints;
      for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
        std::cout << "[Maxwell Convert E->D] dim=" << iDim << ", point=" << iPoint << ", eps=" << epsMu[iPoint] << ", Efield=" << Efield[iDimxnumPoints+iPoint] << ", Dfield=" << Dfield[iDimxnumPoints+iPoint] << std::endl;
      }
    }
#endif

    return 0;

  } // end ConvertEfieldtoDfield()


  int ConvertDfieldtoEfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *recipEpsMu, double *Dfield, double *Efield)
  /******************************************************************************************

        Calculate E = D/epsilon

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

    if (numDim != 3) {
      return 1;
    }

    for (int iDim=0; iDim<numDim; iDim++) {
      FC_MODULE(operators,zxy,OPERATORS,ZXY)(&numDim, &numPoints, numX, opInterval,
                                             &recipEpsMu[0], &Dfield[iDim*numPoints], &Efield[iDim*numPoints]);
    }

#if DEBUG_MAXWELL_CONVERT_DTOEFIELD
    // FOR DEBUGGING PURPOSES
    for (int iDim=0; iDim<numDim; iDim++) {
      size_t iDimxnumPoints = iDim*numPoints;
      for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
        std::cout << "[Maxwell Convert D->E] dim=" << iDim << ", point=" << iPoint << ", 1/eps=" << recipEpsMu[iPoint] << ", Dfield=" << Dfield[iDimxnumPoints+iPoint] << ", Efield=" << Efield[iDimxnumPoints+iPoint] << std::endl;
      }
    }
#endif

    return 0;

  } // end ConvertDfieldtoEfield()


  int ConvertBfieldtoHfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *recipEpsMu, double *Bfield, double *Hfield)
  /******************************************************************************************

        Calculate H = B/mu

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

    if (numDim != 3) {
      return 1;
    }

    for (int iDim=0; iDim<numDim; iDim++) {
      FC_MODULE(operators,zxy,OPERATORS,ZXY)(&numDim, &numPoints, numX, opInterval,
                                             &recipEpsMu[numPoints], &Bfield[iDim*numPoints], &Hfield[iDim*numPoints]);
    }

#if DEBUG_MAXWELL_CONVERT_BTOHFIELD
    // FOR DEBUGGING PURPOSES
    for (int iDim=0; iDim<numDim; iDim++) {
      size_t iDimxnumPoints = iDim*numPoints;
      for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
        std::cout << "[Maxwell Convert B->H] dim=" << iDim << ", point=" << iPoint << ", 1/mu=" << recipEpsMu[numPoints+iPoint] << ", Bfield=" << Bfield[iDimxnumPoints+iPoint] << ", Hfield=" << Hfield[iDimxnumPoints+iPoint] << std::endl;
      }
    }
#endif

    return 0;

  } // end ConvertBfieldtoHfield()


  int ConvertHfieldtoBfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *epsMu, double *Hfield, double *Bfield)
  /******************************************************************************************

        Calculate B = mu*H

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

    if (numDim != 3) {
      return 1;
    }

    for (int iDim=0; iDim<numDim; iDim++) {
      FC_MODULE(operators,zxy,OPERATORS,ZXY)(&numDim, &numPoints, numX, opInterval,
                                             &epsMu[numPoints], &Hfield[iDim*numPoints], &Bfield[iDim*numPoints]);
    }

#if DEBUG_MAXWELL_CONVERT_HTOBFIELD
    // FOR DEBUGGING PURPOSES
    for (int iDim=0; iDim<numDim; iDim++) {
      size_t iDimxnumPoints = iDim*numPoints;
      for (size_t iPoint=0; iPoint<numPoints; iPoint++) {
        std::cout << "[Maxwell Convert H->B] dim=" << iDim << ", point=" << iPoint << ", mu=" << epsMu[numPoints+iPoint] << ", Hfield=" << Hfield[iDimxnumPoints+iPoint] << ", Bfield=" << Bfield[iDimxnumPoints+iPoint] << std::endl;
      }
    }
#endif

    return 0;

  } // end ConvertHfieldtoBfield()


} // end namespace Maxwell