MoL: add Multirate capabilities. This add three new multirate RK schemes to MoL.

Flags indicate whether it is time to execute slow RHS computation.
For instance, in the RK4-RK2 scheme, there are 4 substeps in total, but the RK2 RHS are only evaluated in the very first and in the very last step of the four substeps.

From: Christian Reisswig, minor changes by Roland Haas



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

1 files changed, 409 insertions, 0 deletions

diff --git a/src/ChangeType.c b/src/ChangeType.c
index e0bc333..9c1db4a 100644
--- a/src/ChangeType.c
+++ b/src/ChangeType.c
@@ -27,6 +27,7 @@ CCTK_FILEVERSION(CactusBase_MoL_ChangeType_c);
 #define MOL_EVOLVED_VARTYPE 1
 #define MOL_CONSTRAINED_VARTYPE 2
 #define MOL_SANDR_VARTYPE 3
+#define MOL_EVOLVEDSLOW_VARTYPE 4
 
 /********************************************************************
  ********************* Local Routine Prototypes *********************
@@ -42,6 +43,8 @@ CCTK_FILEVERSION(CactusBase_MoL_ChangeType_c);
 
 CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
 
+CCTK_INT MoL_ChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow);
+
 CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex);
 
 CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex);
@@ -95,6 +98,16 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -174,6 +187,46 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to evolved type from constrained.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariables < MoL_Num_Evolved_Vars))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; 
+             index++)
+        {
+          EvolvedVariableIndexSlow[index] = EvolvedVariableIndexSlow[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        EvolvedVariableIndex[MoLNumEvolvedVariables] = EvolvedIndex;
+        RHSVariableIndex[MoLNumEvolvedVariables] = RHSIndex;
+        MoLNumEvolvedVariables++;
+#ifdef MOLDEBUG
+        printf("Changing (constrained): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariables, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
     case MOL_CONSTRAINED_VARTYPE:
       
       {
@@ -268,6 +321,266 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
 }
 
  /*@@
+   @routine    MoL_ChangeToEvolvedSlow
+   @date       Thu May 30 16:45:30 2002
+   @author     Ian Hawke, Roland Haas
+   @desc 
+   Changes a variable to evolved slow type. Checks to see which type it was
+   before and does the bookkeeping on the index arrays.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+CCTK_INT MoL_ChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
+{
+
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_INT index, usedindex;
+  CCTK_INT vartype; /* See the defines at the top of file */
+  CCTK_INT timelevs;
+  
+  vartype = 0;
+  usedindex = -1;
+  
+  for (index = 0; (index < MoLNumEvolvedVariables)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVED_VARTYPE * 
+      (EvolvedVariableIndex[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
+  for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
+  {
+    vartype = MOL_CONSTRAINED_VARTYPE * 
+      (ConstrainedVariableIndex[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
+  for (index = 0; (index < MoLNumSandRVariables)&&(!vartype); index++)
+  {
+    vartype = MOL_SANDR_VARTYPE * 
+      (SandRVariableIndex[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+
+  switch (vartype)
+  {
+    
+    case MOL_UNKNOWN_VARTYPE:
+
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(EvolvedIndex);
+        if (timelevs < 1)
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved.", EvolvedIndex);
+          CCTK_WARN(0, "The index passed does not correspond to a GF.");
+        }
+        else if (timelevs == 1) 
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved.", EvolvedIndex);
+          CCTK_WARN(0, "The index passed only has a single timelevel.");
+        }
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved (RHS GF).", RHSIndex);
+          CCTK_WARN(0, "The index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (unknown): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariablesSlow, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    case MOL_EVOLVED_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved type from constrained.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumEvolvedVariables - 1; 
+             index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariables--;
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (constrained): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariablesSlow, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        RHSVariableIndex[usedindex] = RHSIndex;
+        break;
+      }
+
+    case MOL_CONSTRAINED_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved type from constrained.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumConstrainedVariables - 1; 
+             index++)
+        {
+          ConstrainedVariableIndex[index] = ConstrainedVariableIndex[index+1];
+        }
+        MoLNumConstrainedVariables--;
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (constrained): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariables, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    case MOL_SANDR_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved type from save and "
+                     "restore.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumSandRVariables - 1; index++)
+        {
+          SandRVariableIndex[index] = SandRVariableIndex[index+1];
+        }
+        MoLNumSandRVariables--;
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (SandR): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariablesSlow, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    default:
+      
+      {
+        CCTK_WARN(0, "Something is seriously wrong in ChangeType.c! "
+                  "Case out of range in switch statement.");
+      }
+      
+  }
+      
+  return 0;
+  
+}
+
+ /*@@
    @routine    MoL_ChangeToConstrained
    @date       Thu May 30 16:47:08 2002
    @author     Ian Hawke
@@ -321,6 +634,16 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == ConstrainedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -397,6 +720,34 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+          RHSVariableIndex[index] = RHSVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        if (timelevs > 1)
+        {
+          if ( !(MoLNumConstrainedVariables < MoL_Num_Constrained_Vars))
+          {
+            CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                       "Warning whilst trying to change variable "
+                       "index %i (%s) to constrained.", 
+                       ConstrainedIndex, CCTK_VarName(ConstrainedIndex));
+            CCTK_WARN(0, "When changing type there are more variables "
+                      "than the accumulator parameter "
+                      "MoL_Num_Constrained_Vars allows. Check that "
+                      "you are accumulating onto this parameter correctly");
+          }
+          ConstrainedVariableIndex[MoLNumConstrainedVariables] = ConstrainedIndex;
+          MoLNumConstrainedVariables++;
+        }
+        break;
+      }
+
     case MOL_CONSTRAINED_VARTYPE:
       
       {
@@ -489,6 +840,16 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == SandRIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -558,6 +919,32 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+
+        if (!(MoLNumSandRVariables < MoL_Num_SaveAndRestore_Vars))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to save and restore.", 
+                     SandRIndex, CCTK_VarName(SandRIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_SaveAndRestore_Vars allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+          RHSVariableIndex[index] = RHSVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        SandRVariableIndex[MoLNumSandRVariables] = SandRIndex;
+        MoLNumSandRVariables++;
+        break;
+      }
+
     case MOL_CONSTRAINED_VARTYPE:
       
       {
@@ -639,6 +1026,16 @@ CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == RemoveIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -680,6 +1077,18 @@ CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+          RHSVariableIndex[index] = RHSVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        break;
+      }
+
     case MOL_CONSTRAINED_VARTYPE:
       
       {