/*@@
   @file      GenericRK.c
   @date      Sun May 26 03:47:15 2002
   @author    Ian Hawke
   @desc 
   This routine performs a generic Runge-Kutta type integration
   given the set of coefficients defined in the RKAlphaCoefficients
   and RKBetaCoefficients arrays. See the article by Shu referenced
   in the documentation for more details.
   @enddesc 
   @version   $Header$
 @@*/

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "ExternalVariables.h"

static const char *rcsid = "$Header$";

CCTK_FILEVERSION(CactusBase_MoL_GenericRK_c);

/********************************************************************
 *********************     Local Data Types   ***********************
 ********************************************************************/

/********************************************************************
 ********************* Local Routine Prototypes *********************
 ********************************************************************/

CCTK_INT AlphaIndex(CCTK_INT Step_Number, CCTK_INT Scratch_Level);

/********************************************************************
 ***************** Scheduled Routine Prototypes *********************
 ********************************************************************/

void MoL_GenericRKAdd(CCTK_ARGUMENTS);

/********************************************************************
 ********************* Other Routine Prototypes *********************
 ********************************************************************/

/********************************************************************
 *********************     Local Data   *****************************
 ********************************************************************/

/********************************************************************
 *********************     External Routines   **********************
 ********************************************************************/

 /*@@
   @routine    MoL_GenericRKAdd
   @date       Sun May 26 03:50:44 2002
   @author     Ian Hawke
   @desc 
   Performs a single step of a generic Runge-Kutta type time
   integration.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

void MoL_GenericRKAdd(CCTK_ARGUMENTS)
{

  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
    
  cGroupDynamicData arraydata;
  CCTK_INT groupindex, ierr;
  CCTK_INT arraytotalsize, arraydim;

  /* FIXME */

#ifdef MOLDOESCOMPLEX

  CCTK_COMPLEX Complex_alpha, Complex_beta, Complex_Delta_Time;
  CCTK_COMPLEX *UpdateComplexVar;
  CCTK_COMPLEX *RHSComplexVar;
  CCTK_COMPLEX *ScratchComplexVar;

#endif

  CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
  CCTK_INT totalsize;
  CCTK_REAL alpha, beta;
  CCTK_REAL *UpdateVar;
  CCTK_REAL *RHSVar;
  CCTK_REAL *ScratchVar;

  CCTK_INT arrayscratchlocation;

  totalsize = 1;
  for (arraydim = 0; arraydim < cctk_dim; arraydim++)
  {
    totalsize *= cctk_lsh[arraydim];
  }
  
  beta = RKBetaCoefficients[MoL_Intermediate_Steps - 
                           (*MoL_Intermediate_Step)];

  /* FIXME */

#ifdef MOLDOESCOMPLEX
  
  Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time), 0);
  Complex_beta = CCTK_Cmplx(beta, 0);

#endif

  /* Real GFs */

  for (var = 0; var < MoLNumEvolvedVariables; var++)
  {
    
    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
                                              EvolvedVariableIndex[var]);
    RHSVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
                                           RHSVariableIndex[var]);
/* #define MOLDEBUG */
#ifdef MOLDEBUG
    printf("In generic RK. Variable %d (%s). RHS %d (%s).\n",
           EvolvedVariableIndex[var],
           CCTK_VarName(EvolvedVariableIndex[var]),
           RHSVariableIndex[var],
           CCTK_VarName(RHSVariableIndex[var]));
#endif

    for (index = 0; index < totalsize; index++)
    {
      UpdateVar[index] = (*Original_Delta_Time) * beta * RHSVar[index];
#ifdef MOLDEBUG
      printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
             var, index, (*Original_Delta_Time), beta, RHSVar[index], 
             UpdateVar[index]);
#endif
    }
    
    for (scratchstep = 0; 
         scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1;
         scratchstep++)
    {

      alphaindex = AlphaIndex(*MoL_Intermediate_Step, scratchstep);
      scratchindex = scratchstep - 1;

      alpha = RKAlphaCoefficients[alphaindex];

      if (scratchstep) 
      {
        /*
          The following would work if all drivers considered
          a vector group to have contiguous storage.
        */
        /*
          ScratchVar = &ScratchSpace[(var * MoL_Num_Scratch_Levels + 
                                   scratchindex) * totalsize];
        */
        ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
                                      CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
                                      + var 
                                      + MoL_Num_Evolved_Vars * scratchindex);
#ifdef MOLDEBUG
        printf("Reading from scratch space, initial address %ld index %d\n", 
               ScratchVar, (var * MoL_Num_Scratch_Levels + 
                            scratchindex) * totalsize);
#endif
      }
      else
      {
        ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 1, 
                                                  EvolvedVariableIndex[var]);
      }
      
      if ( (alpha > MoL_Tiny)||(alpha < -MoL_Tiny) )
      {
        for (index = 0; index < totalsize; index++)
        {
          UpdateVar[index] += alpha * ScratchVar[index];
#ifdef MOLDEBUG
          printf("Variable: %d. Index: %d. step: %d. "
                 "alpha: %f. Scratch: %f. q: %f.\n",
                 var, index, (*MoL_Intermediate_Step), alpha, 
                 ScratchVar[index], UpdateVar[index]);
#endif
        }
      }
      
    }
    
  }

  if (*MoL_Intermediate_Step > 1)
  {
    for (var = 0; var < MoLNumEvolvedVariables; var++)
    {
      UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
                                                EvolvedVariableIndex[var]);
      /*
        The following would work if all drivers considered
        a vector group to have contiguous storage.
      */
      /*
        ScratchVar = &ScratchSpace[(var * MoL_Num_Scratch_Levels + 
                                 MoL_Intermediate_Steps -
                                 (*MoL_Intermediate_Step)) * totalsize];
      */
      ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
                                    CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
                                    + var 
                                    + MoL_Num_Evolved_Vars * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)));
#ifdef MOLDEBUG
      printf("Writing to scratch space, initial address %ld, "
             "index %d, totalsize %d \n", 
             ScratchVar, (var * MoL_Num_Scratch_Levels + 
                          MoL_Intermediate_Steps -
                          (*MoL_Intermediate_Step)) * totalsize, totalsize);
#endif
      for (index = 0; index < totalsize; index++)
      {
        ScratchVar[index] = UpdateVar[index];
#ifdef MOLDEBUG
        printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n",
               var, index, (*MoL_Intermediate_Step), ScratchVar[index]);
        fflush(stdout);
#endif        
      }
    }
  }
  

  /* Complex GFs */

  /* FIXME */

#ifdef MOLDOESCOMPLEX

  for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
  {
    
    UpdateComplexVar = (CCTK_COMPLEX *)CCTK_VarDataPtrI(cctkGH, 0, 
                                                        EvolvedComplexVariableIndex[var]);
    RHSComplexVar = (CCTK_COMPLEX *)CCTK_VarDataPtrI(cctkGH, 0, 
                                                     RHSVariableIndex[var]);

    for (index = 0; index < totalsize; index++)
    {
      UpdateComplexVar[index] = CCTK_CmplxMul(Complex_Delta_Time,
                                              CCTK_CmplxMul(Complex_beta, 
                                                            RHSComplexVar[index])); 
    }
    
    for (scratchstep = 0; 
         scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1;
         scratchstep++)
    {

      alphaindex = AlphaIndex(*MoL_Intermediate_Step, scratchstep);
      scratchindex = scratchstep - 1;

      alpha = RKAlphaCoefficients[alphaindex];
      Complex_alpha = CCTK_Cmplx(alpha, 0);

      if (scratchstep) 
      {
        ScratchComplexVar = &ComplexScratchSpace[(var *
                                                  MoL_Num_Scratch_Levels +
                                                  scratchindex) * totalsize];
      }
      else
      {
        ScratchComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 1, 
                                                            EvolvedComplexVariableIndex[var]);
      }
      
      if ( (alpha > MoL_Tiny)||(alpha < -MoL_Tiny) )
      {
        for (index = 0; index < totalsize; index++)
        {
          UpdateComplexVar[index] =
            CCTK_CmplxAdd(UpdateComplexVar[index],
                          CCTK_CmplxMul(Complex_alpha,
                                        ScratchComplexVar[index]));
        }
      }
      
    }
    
  }

  if (*MoL_Intermediate_Step > 1)
  {
    for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
    {
      UpdateComplexVar = (CCTK_COMPLEX *)CCTK_VarDataPtrI(cctkGH, 0, 
                                                          EvolvedComplexVariableIndex[var]);
      ScratchComplexVar = &ComplexScratchSpace[(var * MoL_Num_Scratch_Levels + 
                                                MoL_Intermediate_Steps -
                                                (*MoL_Intermediate_Step)) * 
                                              totalsize];
      for (index = 0; index < totalsize; index++)
      {
        ScratchComplexVar[index] = UpdateComplexVar[index];
      }
    }
  }

  /* FIXME */

#endif

  /* Real arrays */

  arrayscratchlocation = 0;

  for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
  {
    
    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
                                              EvolvedArrayVariableIndex[var]);
    RHSVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
                                           RHSArrayVariableIndex[var]);
    
    groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
    ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
                                 &arraydata);
    if (ierr)
    {
      CCTK_VWarn(0, __LINE__, __FILE__, "MoL", 
                 "The driver does not return group information "
                 "for group '%s'.", 
                 CCTK_GroupName(groupindex));
    }
    arraytotalsize = 1;
    for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
    {
      arraytotalsize *= arraydata.lsh[arraydim];
    }

    for (index = 0; index < arraytotalsize; index++)
    {
      UpdateVar[index] = (*Original_Delta_Time) * beta * RHSVar[index];
    }
    
    for (scratchstep = 0; 
         scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1;
         scratchstep++)
    {

      alphaindex = AlphaIndex(*MoL_Intermediate_Step, scratchstep);
      scratchindex = scratchstep - 1;

      alpha = RKAlphaCoefficients[alphaindex];

      if (scratchstep) 
      {
        ScratchVar = &ArrayScratchSpace[scratchindex*
                                       (MoL_Max_Evolved_Array_Size+1) + 
                                        arrayscratchlocation];
#ifdef MOLDEBUG
        printf("Reading from scratch space, initial address %ld index %d\n", 
               ScratchVar, scratchindex*(MoL_Max_Evolved_Array_Size+1) + 
               arrayscratchlocation);
#endif
      }
      else
      {
        ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 1, 
                                                  EvolvedArrayVariableIndex[var]);
      }
      
      if ( (alpha > MoL_Tiny)||(alpha < -MoL_Tiny) )
      {
        for (index = 0; index < arraytotalsize; index++)
        {
          UpdateVar[index] += alpha * ScratchVar[index];
#ifdef MOLDEBUG
          printf("Variable: %d. Index: %d. step: %d. "
                 "alpha: %f. Scratch: %f. q: %f.\n",
                 var, index, (*MoL_Intermediate_Step), 
                 alpha, ScratchVar[index], UpdateVar[index]);
#endif
        }
      }
      
    }
   
    arrayscratchlocation += arraytotalsize;
 
  }

  arrayscratchlocation = 0;

  if (*MoL_Intermediate_Step > 1)
  {
    for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
    {
      UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
                                                EvolvedArrayVariableIndex[var]);
      ScratchVar = &ArrayScratchSpace[(MoL_Intermediate_Steps -
                                       (*MoL_Intermediate_Step)) * 
                                     (MoL_Max_Evolved_Array_Size+1) + 
                                     arrayscratchlocation];
      groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
      ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
                                   &arraydata);
      if (ierr)
      {
        CCTK_VWarn(0, __LINE__, __FILE__, "MoL", 
                   "The driver does not return group information "
                   "for group '%s'.", 
                   CCTK_GroupName(groupindex));
      }
      arraytotalsize = 1;
      for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
      {
        arraytotalsize *= arraydata.lsh[arraydim];
      }

#ifdef MOLDEBUG
      printf("Writing to scratch space, initial address %ld, index %d \n", 
             ScratchVar, (MoL_Intermediate_Steps -
                          (*MoL_Intermediate_Step)) * 
             (MoL_Max_Evolved_Array_Size+1) + 
             arrayscratchlocation);
#endif
      for (index = 0; index < arraytotalsize; index++)
      {
        ScratchVar[index] = UpdateVar[index];
#ifdef MOLDEBUG
        printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n",
               var, index, (*MoL_Intermediate_Step), ScratchVar[index]);
#endif        
      }
      arrayscratchlocation += arraytotalsize;
    }
  }
        
  return;
}

/********************************************************************
 *********************     Local Routines   *************************
 ********************************************************************/

CCTK_INT AlphaIndex(CCTK_INT Step_Number, CCTK_INT Scratch_Level)
{
  DECLARE_CCTK_PARAMETERS;
  
  return (MoL_Intermediate_Steps - Step_Number) * MoL_Intermediate_Steps + 
    Scratch_Level;
}