/*@@
   @file      RK2.c
   @date      Sun May 26 04:13:45 2002
   @author    Ian Hawke
   @desc 
   A specialized second order Runge-Kutta time integrator. This is
   the integrator that Shu refers to as the optimal TVD second 
   order method (see reference in documentation). It is equivalent
   to Heun's predictor-corrector method, or the MacCormack method.
   @enddesc 
   @version   $Header$
 @@*/

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

#include "ExternalVariables.h"

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

CCTK_FILEVERSION(AlphaThorns_MoL_RK2_c);

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

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

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

void MoL_RK2Add(CCTK_ARGUMENTS);

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

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

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

 /*@@
   @routine    MoL_RK2Add
   @date       Sun May 26 04:17:23 2002
   @author     Ian Hawke
   @desc 
   Performs second order Runge-Kutta time integration.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

void MoL_RK2Add(CCTK_ARGUMENTS)
{

  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  CCTK_INT index, var;
  CCTK_INT totalsize;
  CCTK_REAL *OldVar;
  CCTK_REAL *UpdateVar;
  CCTK_REAL *RHSVar;

#ifdef MOLDEBUG
  printf("Inside RK2.\nStep %d.\nRefinement %d.\nTimestep %g.\nSpacestep %g.\nTime %g\n",
         MoL_Intermediate_Steps - *MoL_Intermediate_Step + 1,
         *cctk_levfac,
         cctk_delta_time,
         CCTK_DELTA_SPACE(0),
         cctk_time);
#endif  

  totalsize = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];

  switch (*MoL_Intermediate_Step)
  {
  
    case 2:
      {
        for (var = 0; var < MoLNumEvolvedVariables; var++)
        {
          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
                                                   EvolvedVariableIndex[var]);
          RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
                                                RHSVariableIndex[var]);
          
          for (index = 0; index < totalsize; index++)
          {
            UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index];
          }
        }
        break;
      }
    case 1:
      {
        for (var = 0; var < MoLNumEvolvedVariables; var++)
        {
          OldVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 1,
                                                EvolvedVariableIndex[var]);
          UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
                                                   EvolvedVariableIndex[var]);
          RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, 
                                                RHSVariableIndex[var]);
          
          for (index = 0; index < totalsize; index++)
          {
            UpdateVar[index] = 0.5 * (OldVar[index] + UpdateVar[index]) + 
                                      CCTK_DELTA_TIME * RHSVar[index];
          }
        }
        break;
      }
    default:
      {
        CCTK_WARN(0, "RK2 expects MoL_Intermediate_Step to be in [1,2]. This should be caught at ParamCheck - bug Ian!");
        break;
      }
      
  }

  return;

}

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