Add support for evolving complex GFs and (real and complex) GAs.

Only works with ICN or RK2 for now - in fact this commit may break the generic RK methods temporarily.

Note the documentation isn't quite right - there's no longer a seperate function for each different type...


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

1 files changed, 184 insertions, 1 deletions

diff --git a/src/GenericRK.c b/src/GenericRK.c
index 939faa3..d1e3bc8 100644
--- a/src/GenericRK.c
+++ b/src/GenericRK.c
@@ -70,7 +70,16 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
 
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
-  
+    
+  cGroupDynamicData arraydata;
+  CCTK_INT groupindex, ierr;
+  CCTK_INT arraytotalsize, arraydim;
+
+  CCTK_COMPLEX Complex_alpha, Complex_beta, Complex_Delta_Time;
+  CCTK_COMPLEX *UpdateComplexVar;
+  CCTK_COMPLEX *RHSComplexVar;
+  CCTK_COMPLEX *ScratchComplexVar;
+
   CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
   CCTK_INT totalsize;
   CCTK_REAL alpha, beta;
@@ -80,8 +89,13 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
 
   totalsize = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
   
+  Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time), 0);
+
   beta = RKBetaCoefficients[MoL_Intermediate_Steps - 
                            (*MoL_Intermediate_Step)];
+  Complex_beta = CCTK_Cmplx(beta, 0);
+
+  /* Real GFs */
 
   for (var = 0; var < MoLNumEvolvedVariables; var++)
   {
@@ -176,6 +190,175 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
       }
     }
   }
+
+  /* Complex GFs */
+
+  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];
+      }
+    }
+  }
+
+  /* Real arrays */
+
+  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 = &ScratchSpace[(var * MoL_Num_Scratch_Levels + 
+                                   scratchindex) * totalsize];
+#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]);
+      ScratchVar = &ScratchSpace[(var * MoL_Num_Scratch_Levels + 
+                                 MoL_Intermediate_Steps -
+                                 (*MoL_Intermediate_Step)) * totalsize];
+#ifdef MOLDEBUG
+      printf("Writing to scratch space, initial address %ld, index %d \n", 
+             ScratchVar, (var * MoL_Num_Scratch_Levels + 
+                          MoL_Intermediate_Steps -
+                          (*MoL_Intermediate_Step)) * 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]);
+#endif        
+      }
+    }
+  }
         
   return;
 }