Erik Schnetter's implementation of ICN with averaging, so the intermediate steps are always at t+dt.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@24 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b

3 files changed, 145 insertions, 0 deletions

diff --git a/src/ICN.c b/src/ICN.c
index 0a1cced..a5231a7 100644
--- a/src/ICN.c
+++ b/src/ICN.c
@@ -34,6 +34,7 @@ CCTK_FILEVERSION(AlphaThorns_MoL_ICN_c);
  ********************************************************************/
 
 void MoL_ICNAdd(CCTK_ARGUMENTS);
+void MoL_ICNAverage(CCTK_ARGUMENTS);
 
 /********************************************************************
  ********************* Other Routine Prototypes *********************
@@ -176,6 +177,132 @@ void MoL_ICNAdd(CCTK_ARGUMENTS)
 
 }
 
+ /*@@
+   @routine    MoL_ICNAverage
+   @date       Fri Jul 18 14:02:00 2003
+   @author     Ian Hawke, Erik Schnetter
+   @desc 
+   Averages between the current and the previous time level.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_ICNAverage(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  cGroupDynamicData arraydata;
+  CCTK_INT groupindex, ierr;
+  CCTK_INT arraytotalsize, arraydim;
+
+  CCTK_INT index, var;
+  CCTK_INT totalsize;
+  CCTK_REAL *OldVar;
+  CCTK_REAL *UpdateVar;
+  CCTK_REAL *RHSVar;
+  
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+  CCTK_COMPLEX *OldComplexVar;
+  CCTK_COMPLEX *UpdateComplexVar;
+  CCTK_COMPLEX *RHSComplexVar;
+  CCTK_COMPLEX Complex_Half = CCTK_Cmplx(0.5, 0);
+
+#endif
+
+#ifdef MOLDEBUG
+  printf("Inside ICN.\nProcessor %d.\nStep %d.\nRefinement %d.\nTimestep %g.\nSpacestep %g.\nTime %g\n",
+         CCTK_MyProc(cctkGH),
+         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];
+#ifdef MOLDEBUG
+  printf("MoL: the ICN routine says dt = %f.\n", CCTK_DELTA_TIME);
+#endif
+  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 * (UpdateVar[index] + OldVar[index]);
+    }
+  }
+  
+  for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+  {
+    OldVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 1,
+                                             EvolvedArrayVariableIndex[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] = 0.5 * (UpdateVar[index] + OldVar[index]);
+    }
+  }
+
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+  
+  for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
+  {
+    OldComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 1,
+                                                    EvolvedComplexVariableIndex[var]);
+    UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
+                                                       EvolvedComplexVariableIndex[var]);
+    RHSComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                    RHSComplexVariableIndex[var]);
+    
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateComplexVar[index] = CCTK_CmplxMul(Complex_Half, CCTK_CmplxAdd(UpdateComplexVar[index], OldComplexVar[index]));
+    }
+  }
+
+#endif
+  
+  return;
+
+}
+
 /********************************************************************
  *********************     Local Routines   *************************
  ********************************************************************/
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 2f843be..daf0d68 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -246,6 +246,12 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
     sprintf(infoline, "Iterative Crank Nicholson with %i iterations", 
             MoL_Intermediate_Steps);
   }
+  else if (CCTK_EQUALS(ODE_Method,"ICN-avg"))
+  {
+    sprintf(infoline,
+            "Averaging iterative Crank Nicholson with %i iterations", 
+            MoL_Intermediate_Steps);
+  }
   else 
   {
     CCTK_WARN(0, "ODE_Method not recognized!");
diff --git a/src/SetTime.c b/src/SetTime.c
index 927f5fa..b7634c9 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -80,6 +80,10 @@ int MoL_SetTime(CCTK_ARGUMENTS)
   {
     cctkGH->cctk_delta_time = 0.5*(*Original_Delta_Time);
   }
+  else if (CCTK_EQUALS(ODE_Method,"ICN-avg"))
+  {
+    cctkGH->cctk_delta_time = *Original_Delta_Time;
+  }
   else if (CCTK_EQUALS(ODE_Method,"Generic"))
   {
     beta = RKBetaCoefficients[0];
@@ -130,6 +134,10 @@ int MoL_ResetTime(CCTK_ARGUMENTS)
   {
     cctkGH->cctk_time = (*Original_Time)-0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac;
   }
+  else if (CCTK_EQUALS(ODE_Method,"ICN-avg"))
+  {
+    cctkGH->cctk_time = (*Original_Time);
+  }
   else if (CCTK_EQUALS(ODE_Method,"Generic"))
   {
     previous_times[0] = (*Original_Time) -
@@ -223,6 +231,10 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS)
       cctkGH->cctk_delta_time = 0.5*(*Original_Delta_Time);
     }
   }
+  else if (CCTK_EQUALS(ODE_Method,"ICN-avg"))
+  {
+    cctkGH->cctk_delta_time = (*Original_Delta_Time);
+  }
   else if (CCTK_EQUALS(ODE_Method,"Generic"))
   {
     cctkGH->cctk_delta_time = RKBetaCoefficients[MoL_Intermediate_Steps -