The Method of Lines thorn (version 2 - see below).

MoL provides generic integration methods for multiple thorns simultaneously. By providing a layer between the driver and evolution thorns, this should mean that some technical issues to do with mesh refinement can be ignored. It also allows you to choose different evolution methods (in time, at least). But the primary purpose is to unambiguously evolve models in different thorns at the same time.

This is version 2 - the one that will work with mesh refinement. It's a straight copy of HawkeCVS/Public/CactusMoL2/MoL2 and I haven't checked that it will work "as is" with the new name (one of the reasons it goes into Alpha).

At the moment the only evolution method guaranteed to work with mesh refinement is ICN, although RK2 should. The generic RK methods will do something that will be subtly wrong...

Note that the "old" way of registering variables (through the functions declared in header files) is still there. As soon as function aliasing settles down, this will be removed.

Also to be done:

Better documentation
Tidy up the code (especially the debugging statements)
Optimize
Add various useful time evolution methods.


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

1 files changed, 599 insertions, 0 deletions

diff --git a/src/Registration.c b/src/Registration.c
new file mode 100644
index 0000000..7152304
--- /dev/null
+++ b/src/Registration.c
@@ -0,0 +1,599 @@
+ /*@@
+   @file      Registration.c
+   @date      Thu May 30 11:21:44 2002
+   @author    Ian Hawke
+   @desc 
+   The external functions called (via function aliasing) by physics
+   thorns to tell MoL that they want these GFs to be treated as a
+   given type.
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "ExternalVariables.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusMoL2_MoL2_Registration_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
+
+CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex);
+
+CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex);
+
+CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex, 
+                                 CCTK_INT RHSGroupIndex);
+
+CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex);
+
+CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex);
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RegisterEvolved
+   @date       Thu May 30 11:36:59 2002
+   @author     Ian Hawke
+   @desc 
+   Given the index of the GF to be evolved and the RHS GF, it stores 
+   the indexes for later use together with various error checking.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
+{
+
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_INT ierr, index, varused, numtimelevs1, numtimelevs2;
+  
+  /*
+  printf("Arrived in MoLRegisterEvolved \n");
+  printf("The indexes are %d and %d.\n",EvolvedIndex, RHSIndex);
+  printf("These correspond to variables %s and %s.\n",
+         CCTK_VarName(EvolvedIndex),CCTK_VarName(RHSIndex));
+  printf("The pointer to EvolvedVariableIndex: %p\n",
+         EvolvedVariableIndex);
+  */
+
+  numtimelevs1 = CCTK_NumTimeLevelsFromVarI(EvolvedIndex);
+  numtimelevs2 = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+  
+  if ( (numtimelevs1 < 0) || (numtimelevs2 < 0) ) 
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning for variable index %i", EvolvedIndex); 
+    CCTK_WARN(0, "The index passed does not correspond to a GF.");
+  }
+
+  if (numtimelevs1 < 2)
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning for variable index %i name %s", EvolvedIndex, CCTK_VarName(EvolvedIndex));
+    CCTK_WARN(0, "The GF passed only has one timelevel. It must have at least two.");
+  }
+
+  varused = 0;
+
+  for (index = 0; (index < MoLNumEvolvedVariables)&&(!varused); index++)
+  {
+    varused = (EvolvedIndex == EvolvedVariableIndex[index]);
+    /*
+    printf("Registering %d. Checking index %d which is %d\n",EvolvedIndex,
+           index,EvolvedVariableIndex[index]);
+    */
+  }
+
+  if (varused)
+  {
+
+    CCTK_VWarn(2,__LINE__,__FILE__,"MoL",
+               "The GF %s has already been registered as an evolved variable with RHS GF %s. The attempt to register with RHS GF %s will be ignored",
+               CCTK_VarName(EvolvedIndex),
+               CCTK_VarName(RHSVariableIndex[index-1]),
+               CCTK_VarName(RHSIndex));
+
+  }
+  else
+  {
+  
+    EvolvedVariableIndex[MoLNumEvolvedVariables] = EvolvedIndex;
+    RHSVariableIndex[MoLNumEvolvedVariables] = RHSIndex;
+  
+    MoLNumEvolvedVariables++;
+    
+    if (MoLNumEvolvedVariables > MoL_Num_Evolved_Vars)
+    {
+      CCTK_WARN(0,"You have tried to register more evolved variables than the accumulator parameter MoL_Num_Evolved_Variables allows. Check that you are accumulating onto this parameter correctly");
+    }
+
+  }
+
+  varused = 0;
+
+  for (index = 0; (index < MoLNumConstrainedVariables)&&(!varused); index++)
+  {
+    varused = (EvolvedIndex == ConstrainedVariableIndex[index]);
+  }
+  
+  if (varused)
+  {
+    for (index = varused; index < MoLNumConstrainedVariables-1; index++)
+    {
+      ConstrainedVariableIndex[index] = ConstrainedVariableIndex[index+1];
+    }
+    MoLNumConstrainedVariables--;
+  }
+    
+  varused = 0;
+
+  for (index = 0; (index < MoLNumSandRVariables)&&(!varused); index++)
+  {
+    varused = (EvolvedIndex == SandRVariableIndex[index]);
+  }
+  
+  if (varused)
+  {
+    for (index = varused; index < MoLNumSandRVariables-1; index++)
+    {
+      SandRVariableIndex[index] = SandRVariableIndex[index+1];
+    }
+    MoLNumSandRVariables--;
+  }
+  
+  return 0;
+
+}
+
+ /*@@
+   @routine    MoL_RegisterConstrained
+   @date       Thu May 30 12:35:58 2002
+   @author     Ian Hawke
+   @desc 
+   Given the index of the GF, register it as a constrained variable.
+   If there's only one timelevel then ignore it as there will be no
+   rotation and so MoL doesn't have to do anything.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex)
+{
+
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT numtimelevs, varused, evolved, index;
+  
+  numtimelevs = CCTK_NumTimeLevelsFromVarI(ConstrainedIndex);
+
+  if (numtimelevs < 1) {
+
+    CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning for constrained variable index %i", ConstrainedIndex);
+    CCTK_WARN(0, "The index passed does not correspond to a GF.");
+
+  }
+  else if (numtimelevs > 1) {
+    
+    varused = 0;
+    
+    for (evolved = 0; (evolved < MoLNumEvolvedVariables)&&(!varused); evolved++)
+    {
+      varused = (EvolvedVariableIndex[evolved] == ConstrainedIndex);
+    }
+    
+    for (evolved = 0; (evolved < MoLNumConstrainedVariables)&&(!varused); evolved++)
+    {
+      varused = (ConstrainedVariableIndex[evolved] == ConstrainedIndex);
+    }
+    
+    if (!varused)
+    {
+      ConstrainedVariableIndex[MoLNumConstrainedVariables] = ConstrainedIndex;
+      MoLNumConstrainedVariables++;
+
+      if (MoLNumConstrainedVariables > MoL_Num_Constrained_Vars)
+      {
+        CCTK_WARN(0,"You have tried to register more evolved variables than the accumulator parameter MoL_Num_Evolved_Variables allows. Check that you are accumulating onto this parameter correctly");
+      }
+
+    }
+    
+    varused = 0;
+    
+    for (evolved = 0; (evolved < MoLNumSandRVariables)&&(!varused); evolved++)
+    {
+      varused = (SandRVariableIndex[evolved] == ConstrainedIndex);
+    }
+    
+    if (varused)
+    {
+      for (index = evolved; index < MoLNumSandRVariables-1; index++)
+      {
+        SandRVariableIndex[index] = SandRVariableIndex[index+1];
+      }
+      MoLNumSandRVariables--;
+    }
+    
+  }
+  else
+  {
+    
+    CCTK_VInfo(CCTK_THORNSTRING,
+               "MoL will not treat variable %s as a constrained variable at it has only one timelevel. This should not cause problems with the evolution.\n", 
+               CCTK_VarName(ConstrainedIndex));
+
+  }
+  
+  return 0;
+  
+}
+
+ /*@@
+   @routine    MoL_RegisterSaveAndRestore
+   @date       Thu May 30 12:37:40 2002
+   @author     Ian Hawke
+   @desc 
+   Given a GF index store it for later use as a save and restore type.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex)
+{
+
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_INT numtimelevs, varused, evolved;
+  
+  numtimelevs = CCTK_NumTimeLevelsFromVarI(SandRIndex);
+
+  if (numtimelevs < 1) {
+
+    CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning for save and restore variable index %i", SandRIndex);
+    CCTK_WARN(0, "The index passed does not correspond to a GF.");
+
+  }
+  else if (numtimelevs > 1) {
+    
+    varused = 0;
+    
+    for (evolved = 0; (evolved < MoLNumEvolvedVariables)&&(!varused); evolved++)
+    {
+      varused = (EvolvedVariableIndex[evolved] == SandRIndex);
+    }
+    
+    for (evolved = 0; (evolved < MoLNumConstrainedVariables)&&(!varused); evolved++)
+    {
+      varused = (ConstrainedVariableIndex[evolved] == SandRIndex);
+    }
+    
+    for (evolved = 0; (evolved < MoLNumSandRVariables)&&(!varused); evolved++)
+    {
+      varused = (SandRVariableIndex[evolved] == SandRIndex);
+    }
+
+    if (!varused)
+    {
+      SandRVariableIndex[MoLNumSandRVariables] = SandRIndex;
+      MoLNumSandRVariables++;
+
+      if (MoLNumSandRVariables > MoL_Num_SaveAndRestore_Vars)
+      {
+        CCTK_WARN(0,"You have tried to register more evolved variables than the accumulator parameter MoL_Num_Evolved_Variables allows. Check that you are accumulating onto this parameter correctly");
+      }
+
+    }
+        
+  }
+  else
+  {
+    
+    CCTK_VInfo(CCTK_THORNSTRING,
+               "MoL will not treat variable %s as a save and restore variable at it has only one timelevel. This should not cause problems with the evolution.\n", 
+               CCTK_VarName(SandRIndex));
+
+  }
+  
+  return 0;
+  
+}
+
+CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex, 
+                                 CCTK_INT RHSGroupIndex)
+{
+
+  CCTK_INT EvolvedGroupFirstVar, RHSGroupFirstVar, GroupNumVars, retval;
+  CCTK_INT EvolvedVar, RHSVar;
+  
+  EvolvedGroupFirstVar = CCTK_FirstVarIndexI(EvolvedGroupIndex);
+  if (EvolvedGroupFirstVar < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Evolved group index %i is not a real group index.",
+              EvolvedGroupIndex);
+  }
+  
+  RHSGroupFirstVar = CCTK_FirstVarIndexI(RHSGroupIndex);
+  if (RHSGroupFirstVar < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "RHS group index %d is not a real group index.",
+              RHSGroupIndex);
+  }
+
+  GroupNumVars = CCTK_NumVarsInGroupI(EvolvedGroupIndex);
+  if (CCTK_NumVarsInGroupI(RHSGroupIndex) != GroupNumVars)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "There are a different number of variables in evolved group %d and RHS group %d.", EvolvedGroupIndex, RHSGroupIndex);
+  }
+  
+  retval = 0;
+  
+  for (EvolvedVar = EvolvedGroupFirstVar, RHSVar = RHSGroupFirstVar;
+       EvolvedVar < EvolvedGroupFirstVar + GroupNumVars;
+       EvolvedVar++, RHSVar++)
+  {
+    retval += MoLRegisterEvolved(EvolvedVar, RHSVar);
+  }
+  
+  return retval;
+}
+
+
+CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex)
+{
+  
+  CCTK_INT ConstrainedGroupFirstVar, GroupNumVars, retval;
+  CCTK_INT ConstrainedVar;
+  
+  ConstrainedGroupFirstVar = CCTK_FirstVarIndexI(ConstrainedGroupIndex);
+  if (ConstrainedGroupFirstVar < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Constrained group index %i is not a real group index.",
+              ConstrainedGroupIndex);
+  }
+  
+  GroupNumVars = CCTK_NumVarsInGroupI(ConstrainedGroupIndex);
+  
+  retval = 0;
+  
+  for (ConstrainedVar = ConstrainedGroupFirstVar;
+       ConstrainedVar < ConstrainedGroupFirstVar + GroupNumVars;
+       ConstrainedVar++)
+  {
+    retval += MoLRegisterConstrained(ConstrainedVar);
+  }
+  
+  return retval;
+}
+
+CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex)
+{
+  
+  CCTK_INT SandRGroupFirstVar, GroupNumVars, retval;
+  CCTK_INT SandRVar;
+  
+  SandRGroupFirstVar = CCTK_FirstVarIndexI(SandRGroupIndex);
+  if (SandRGroupFirstVar < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Save and Restore group index %i is not a real group index.",
+              SandRGroupIndex);
+  }
+  
+  GroupNumVars = CCTK_NumVarsInGroupI(SandRGroupIndex);
+  
+  retval = 0;
+  
+  for (SandRVar = SandRGroupFirstVar;
+       SandRVar < SandRGroupFirstVar + GroupNumVars;
+       SandRVar++)
+  {
+    retval += MoLRegisterSaveAndRestore(SandRVar);
+  }
+  
+  return retval;
+}
+
+/*
+  Old function names. Just calls the new version. 
+  Included for compatibility
+*/
+
+CCTK_INT MoL_RegisterVar(CCTK_INT molvarindex,CCTK_INT molrhsvarindex)
+{
+  return MoL_RegisterEvolved(molvarindex, molrhsvarindex);
+}
+
+CCTK_INT MoL_RegisterPrimitive(CCTK_INT primitiveindex)
+{
+  return MoL_RegisterConstrained(primitiveindex);
+}
+
+CCTK_INT MoL_RegisterDepends(CCTK_INT dependsindex)
+{
+  return MoL_RegisterSaveAndRestore(dependsindex);
+}
+
+CCTK_INT MoL_RegisterVarGroup(CCTK_INT groupindex,CCTK_INT rhsgroupindex)
+{
+  return MoL_RegisterEvolvedGroup(groupindex, rhsgroupindex);
+}
+
+CCTK_INT MoL_RegisterPrimitiveGroup(CCTK_INT groupindex)
+{
+  return MoL_RegisterConstrainedGroup(groupindex);
+}
+
+CCTK_INT MoL_RegisterDependsGroup(CCTK_INT groupindex)
+{
+  return MoL_RegisterConstrainedGroup(groupindex);
+}
+
+CCTK_INT MoL_ChangeVarToEvolved(CCTK_INT varindex, CCTK_INT rhsindex)
+{
+  return MoL_ChangeToEvolved(varindex, rhsindex);
+}
+
+CCTK_INT MoL_ChangeVarToDependent(CCTK_INT dependsindex)
+{
+  return MoL_ChangeToSaveAndRestore(dependsindex);
+}
+
+CCTK_INT MoL_ChangeVarToPrimitive(CCTK_INT primitiveindex)
+{
+  return MoL_ChangeToConstrained(primitiveindex);
+}
+
+CCTK_INT MoL_ChangeVarToNone(CCTK_INT removeindex)
+{
+  return MoL_ChangeToNone(removeindex);
+}
+
+/* 
+  Fortran wrappers for the above functions. 
+  Should be replaced by using function aliasing eventually.
+*/
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterEvolved)(int *ierr, 
+                                            CCTK_INT *EvolvedIndex,
+                                            CCTK_INT *RHSIndex)
+{
+  *ierr = MoL_RegisterEvolved(*EvolvedIndex, *RHSIndex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterConstrained)(int *ierr, 
+                                            CCTK_INT *EvolvedIndex)
+{
+  *ierr = MoL_RegisterConstrained(*EvolvedIndex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterSaveAndRestore)(int *ierr, 
+                                            CCTK_INT *EvolvedIndex)
+{
+  *ierr = MoL_RegisterSaveAndRestore(*EvolvedIndex);
+  return;
+}
+
+/* And for MoL compatibility... */
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterVar)(int *ierr, 
+                                            CCTK_INT *molvarindex,
+                                            CCTK_INT *molrhsvarindex)
+{
+  *ierr = MoL_RegisterVar(*molvarindex, *molrhsvarindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterDepends)(int *ierr, 
+                                                CCTK_INT *moldependsindex)
+{
+  *ierr = MoL_RegisterDepends(*moldependsindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterPrimitive)(int *ierr, 
+                                                CCTK_INT *molprimitiveindex)
+{
+  *ierr = MoL_RegisterPrimitive(*molprimitiveindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterVarGroup)(int *ierr, 
+                                                 CCTK_INT *groupindex,
+                                                 CCTK_INT *rhsgroupindex)
+{
+  *ierr = MoL_RegisterVarGroup(*groupindex, *rhsgroupindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterPrimitiveGroup)(int *ierr, 
+                                                       CCTK_INT *groupindex)
+{
+  *ierr = MoL_RegisterPrimitiveGroup(*groupindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_RegisterDependsGroup)(int *ierr, 
+                                                       CCTK_INT *groupindex)
+{
+  *ierr = MoL_RegisterDependsGroup(*groupindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_ChangeVarToEvolved)(int *ierr, 
+                                                   CCTK_INT *varindex,
+                                                   CCTK_INT *rhsindex)
+{
+  *ierr = MoL_ChangeVarToEvolved(*varindex, *rhsindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_ChangeVarToDependent)(int *ierr, 
+                                                     CCTK_INT *dependsindex)
+{
+  *ierr = MoL_ChangeVarToDependent(*dependsindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_ChangeVarToPrimitive)(int *ierr, 
+                                                     CCTK_INT *primitiveindex)
+{
+  *ierr = MoL_ChangeVarToPrimitive(*primitiveindex);
+  return;
+}
+
+void CCTK_FCALL CCTK_FNAME(MoL_ChangeVarToNone)(int *ierr, 
+                                                CCTK_INT *removeindex)
+{
+  *ierr = MoL_ChangeVarToNone(*removeindex);
+  return;
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/