first set of the slicing

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinBase/CoordGauge/trunk@4 edbb7e70-9571-45d5-a481-0a560a9b4751

2 files changed, 408 insertions, 50 deletions

diff --git a/src/Slicing.c b/src/Slicing.c
index 4d41b55..a6443de 100644
--- a/src/Slicing.c
+++ b/src/Slicing.c
@@ -1,17 +1,21 @@
+static char *rcsid = "$Id$";
+
  /*@@
    @file      Slicing.c
    @date      Wed Jul 14 16:03:12 1999
    @author    Gerd Lanfermann
    @desc 
      A flexible calling structure to calling slicing functions:
-     * a slicing is registered with its name, with a flag to tell if its activated by
-       a parameter, with a function, which tells if the slicing should be used 
-       in the next interation. 
-     * The function is registered in a startup routine of the code providing the 
+     * a slicing is registered with its NAME, with a NUMBER derived from
+       the position in the parameter and indicating its priority, with a FUNCTION, 
+       which tells if the slicing should be used in the next iteration. 
+     * The slicing is REGISTERED in a startup routine of the thorn providing the 
        slicing.
-     * The param_active flag is set where  ? FIXME
-     * This structure is evaluted during a CCTK_PRESTEP, information can be queried 
-       by all slicings, which can take their turn or skip.
+     * This structure (HandledData) is created during STARTUP
+     * This strucure is initialized during CCTK_INITIAL
+     * This structure is evaluted during a CCTK_PRESTEP and an integer
+       active_slicing_handle defines which slicing is next. All slicings have to check
+       if they match this handle and only then take their turn.
    @enddesc 
  @@*/
 
@@ -19,18 +23,22 @@
 #include <stdlib.h>
 
 #include "cctk.h"
+#include "cctk_parameters.h"
+#include "FortranString.h"
 #include "StoreHandledData.h"
 #include "WarnLevel.h"
 #include "Slicing.h"
 #include "ErrorCodes.h"
 #include "Groups.h"
+#include "GroupsOnGH.h"
 #include "Misc.h"
 
-static char *rcsid = "$Id$";
+#define ESLIC_DEBUG
+
 
 /* Local data holding info on slicings..*/
 
-static int num_slics = 0;
+static int num_slicings = 0;
 static cHandledData *Eslicings = NULL;
 
  /*@@
@@ -38,9 +46,11 @@ static cHandledData *Eslicings = NULL;
    @date       Wed Jul 14 15:57:32 1999
    @author     Gerd Lanfermann
    @desc 
-   This routine registers a slicing in a structure whcih provides information
+   This C routine registers a slicing in a structure which provides information
    whether it is activated by parameter and which function is responsible 
    for evaluating whether the slicing should be used in the next iteartion.
+
+   Called at STARTUP by the thorn which provides the slicing.
    @enddesc 
    @calls     
    @calledby   
@@ -53,21 +63,26 @@ static cHandledData *Eslicings = NULL;
 int Einstein_RegisterSlicing(const char *slice_name)
 {
   int handle;
+  char *err;
 
   struct Einstein_slicing *new_slicing;
-  
+
   /*Check that this Slicing has not been registered, yet */
-  handle = Einstein_GetSlicing(slice_name);
+  handle = Util_GetHandle(Eslicings,slice_name,NULL);
   
   if (handle<0) {
-    
+  
     /* Allocate new slicing structure */
     new_slicing = (struct Einstein_slicing*) malloc(sizeof(struct Einstein_slicing));
 
     if (new_slicing) {
       
       /* store new_slicing in the StoreHandledData-array and returns the handle for it */
-      handle = CCTK_NewHandle(&Eslicings, slice_name, new_slicing);
+      handle = Util_NewHandle(&Eslicings, slice_name, new_slicing);
+    
+#ifdef ESLIC_DEBUG
+      printf("REGISTER_SLICING: name: >>%s<<  new handle %d \n",slice_name,handle);
+#endif
 
       /*Initialize the Einstein_slicing structure */
 
@@ -81,30 +96,40 @@ int Einstein_RegisterSlicing(const char *slice_name)
       
       /* Flag to indicate if this slicing activated by a parameter [0/1] */
       new_slicing->param_active = 0;
+
+      /*increase the counter for registered slicings */
+      num_slicings++;
       
     } else 
       {
 	/* Memory failure */
+	CCTK_WARN(0,"cannot allocate memory for new slicing");
 	handle = -2;  /* FIXME: some decent MACROS here */
       }
   }
   else
-    {
+    { 
       /* Extension handle already exists */
+      err =(char*)malloc(256*sizeof(char));
+      sprintf(err,"REGISTER_SLICING: new handle %s (# %d) allready in use \n",
+	      slice_name,handle);
+      CCTK_WARN(0,err);
+      if (err) free(err);
+
       handle = -1;
     }
-  
+
   return handle;
 }
 
- /*@@
-   @routine    Einstein_SetSliceActive
-   @date       Wed Jul 14 16:13:33 1999
+/*@@
+   @routine    Einstein_RegisterSlicing 
+   @date       Wed Jul 14 15:57:32 1999
    @author     Gerd Lanfermann
    @desc 
-     Sets the param_active flag to active (1) or inactive (0), reflecting that
-     it is specified as slicing in the parameter database, or that is not used
-     anymore.
+   This FORTRAN routine registers a slicing in a structure which provides information
+   whether it is activated by parameter and which function is responsible 
+   for evaluating whether the slicing should be used in the next iteartion.
    @enddesc 
    @calls     
    @calledby   
@@ -113,26 +138,16 @@ int Einstein_RegisterSlicing(const char *slice_name)
    @endhistory 
 
 @@*/
-
-int Einstein_SetSliceActive(int handle, int flag) 
-{
-  int return_code;
-  struct Einstein_slicing *slicing;
-
-  slicing = CCTK_GetHandledData(Eslicings, handle);
-  
-  if (slicing) {
-    slicing->param_active = flag;
-    return_code = 1;
-  }
-  else {
-    /* handle didn't return the data pointer */
-    return_code = 0;
-  }
   
-  return return_code;
+int FMODIFIER FORTRAN_NAME(Einstein_RegisterSlicing)(ONE_FORTSTRING_ARG)
+{
+  int handle;
+  ONE_FORTSTRING_CREATE(name)
+  handle = CCTK_VarIndex(name); 
+  free(name);
+  return(handle);
 }
-  
+ 
 
  /*@@
    @routine    Einstein_RegisterTimeToSlice
@@ -143,6 +158,8 @@ int Einstein_SetSliceActive(int handle, int flag)
       information whether the slicing is supposed to used during the next
       iteration. The user can register any function to check for iteration 
       conditions, grid function behavior, etc. 
+
+      This has to be in C currently. We cannot wrapp this easily (TOM).
    @enddesc 
    @calls     
    @calledby   The user
@@ -152,25 +169,344 @@ int Einstein_SetSliceActive(int handle, int flag)
 
 @@*/
 
-int Einstein_RegisterTimeToSlice(int handle, void (*func)(cGH *))
+int Einstein_RegisterTimeToSlice(int handle, int (*func)(cGH *))
 {
-  int return_code;
+  int return_code=1;
+  char *err;
   struct Einstein_slicing *slicing;
 
-  slicing = CCTK_GetHandledData(Eslicings, handle);
+  slicing = Util_GetHandledData(Eslicings, handle);
   
   if (slicing)  {
     slicing->timetoslice=func;
     return_code = 1;
   }
   else {
-    /* handle didn't return the data pointer */
+    /* handle didn't return the data pointer */ 
+    err=(char*)malloc(128*sizeof(char));
+    sprintf(err,"could not obtain slicing-structure for handle %d",handle);
+    CCTK_WARN(1,err);
+    if (err) free (err);
+
     return_code = 0;
-  }
+  }   
+
+#ifdef ESLIC_DEBUG
+    printf("REGISTER_SLICING: func OK for handle %d\n",handle);
+#endif
+
   return return_code;
 }
 
 
+/*@@
+   @routine    Einstein_ActivateSlicing
+   @date       Thu Jul 15 17:44:04 1999
+   @author     Gerd Lanfermann
+   @desc 
+     After the slicings have been registered by the thorns (in Startup .e.g.),
+     this routines checks which slicing is specified in the parameter database
+     and activates them. The activation flag takes a number, which reflect the order 
+     in which the keywords appear in the parameter: first = 1 = highest priority, 
+     sec. = 2 ,lower priority.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+  
+
+void Einstein_ActivateSlicing(cGH *GH) {
+  DECLARE_CCTK_PARAMETERS
+
+  struct Einstein_slicing *slic;
+     
+  int handle, priority, h;
+  char *err, *split_string;
+  char *item=NULL, *after=NULL;
+
+
+  /* we have a clearly defined slicing - no "mixed" */
+  if (!CCTK_Equals(slicing,"mixed")) 
+    { 
+      handle= Util_GetHandle(Eslicings,slicing,NULL);
+
+      if (handle<0) {	
+	err = (char*)malloc(256*sizeof(char));
+	sprintf(err,"%s%s%s%s",
+		"ERROR: cannot get handle for slicing: >>",
+		slicing,
+		"<< \nTypo? Is that slicing registered by its provider",
+		"with Einstein_RegisterSlicing?\n");
+	CCTK_WARN(1,err);
+	if (err) free(err);
+	return;
+      }
+
+      slic  = Util_GetHandledData(Eslicings,handle);
+      if (!slic) CCTK_WARN(0,"INTERNAL ERROR: cannot access handle structure!");
+      slic->param_active = 1;
+
+      /*Some error checking: */
+      if (slic->timetoslice!=NULL) {
+	err = (char*)malloc(256*sizeof(char));
+	sprintf(err,"%s%s%s%s",
+		"ERROR: You have specified a unique slicing (",slic->name,
+		") AND a condition-function is registered for that slicing.\n ",
+		"Impossible, what should I do inbetween ? Rather set slicing -> mixed");
+	CCTK_WARN(0,err);
+	if (err) free(err);
+      }
+    } 
+  /* We have mixed slicing: we check the mixed_slicing parameter */
+  else
+    {
+      /* for all slicings in the mix: e.g. mixed_slicing = "1+log maximal static" 
+	 1) get the handle for each slicing 
+	 2) set slicing->param_active, counting from 1 to n according to the
+	    position. The higher the number the lower the priority (0=inactive!!) 
+      */
+      priority    = 1;
+      split_string= mixed_slicing;
+
+      /* Here we take the string appart */
+      while (Util_SplitString(&item,&after,split_string," ")==0) {
+	slic = Util_GetHandledData(Eslicings,Util_GetHandle(Eslicings,item,NULL));
+
+	if (!slic) { 
+	  err=(char*)malloc(128*sizeof(char));
+	  sprintf(err,"ERROR: no registered slicing found for item >%s<\n",item);
+	  CCTK_WARN(0,err);
+	  if (err) free(err);
+	}
+
+	slic->param_active = priority;
+#ifdef  ESLIC_DEBUG
+	printf("ACTIVATE_SLICING: found %s priority %d \n",item, priority);
+#endif	
+	priority++;
+	split_string = after;
+      }
+      
+      if (item) free(item); 
+      if (after)free(after); 
+      if (split_string) free(split_string);
+    }
+
+  /* We initialize the active_slicing_handle at CACTUS_BASE 
+     for the first time. This is important as some thorns might depend on
+     slicing info quite early.
+     Problem: if s.th. uses this handle before we set it, its not initialized!
+     FIXEME: this routine should be registered at startup AFTER 
+     the Startup things. */
+    Einstein_SetNextSlicing(GH);
+
+
+#ifdef ESLIC_DEBUG   
+  if (CCTK_Equals(slicing_verbose,"yes")) 
+    { 
+      /* Some diagnostic output */
+      for (h=0;h<num_slicings;h++) {
+	slic = Util_GetHandledData(Eslicings,h);
+	if (!slic) CCTK_WARN(0,"INTERNAL ERROR: cannot get slicing for handle.");
+	
+	printf(" Einstein_ActivateSlicing:  >>%s<<   priority: %d   %s \n",
+	       slic->name,
+	       slic->param_active,
+	       (char*)((slic->timetoslice)?"FUNC":"nofunc "));
+      }
+    }
+#endif
+
+}
+ /*@@
+   @routine    Einstein_GetSlicingHandle
+   @date       Thu Jul 22 11:30:47 1999
+   @author     Gerd Lanfermann
+   @desc 
+     When called with the name of a (registered!!) slicing, 
+     it returns the associated handle (a number). This handle can be used
+     to compare to active_slicing_handle, which defines the slicing 
+     for the next iteration. 
+   @enddesc 
+   @calls   CCTK_GetHandle.c  
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int Einstein_GetSlicingHandle(const char *name) {
+  int handle;
+  char *err;
+  printf("NAME: %s\n",name);
+  
+  handle =  Util_GetHandle(Eslicings,name,NULL);
+  if (handle<0) {
+    err = (char*)malloc(256*sizeof(char));
+    sprintf(err,"ERROR: cannot get handle for slicing %s\n%s",name,
+	        "Typo ? Is this slicing registered by the thorn ? ");
+    CCTK_WARN(0,err);
+    if (err) free(err);
+  }
+  return(handle);
+}
+
+/*@@
+   @routine    Einstein_GetSlicingHandle-F
+   @date       Thu Jul 22 11:30:47 1999
+   @author     Gerd Lanfermann
+   @desc 
+     Fortran routine; error checking done in the C-routine;
+   @enddesc 
+   @calls   Einstein_GetSlicingHandle  
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int FMODIFIER FORTRAN_NAME(Einstein_GetSlicingHandle)(ONE_FORTSTRING_ARG) 
+{
+  int handle;
+  ONE_FORTSTRING_CREATE(name)
+  printf("NAME: %s\n",name);
+  handle = Util_GetHandle(Eslicings,name,NULL);
+  printf("HANDLE: %d\n",handle);
+  free(name);
+  return(handle);
+}
+
+
+
+/*@@
+   @routine    Einstein_SetNextSlicing
+   @date       Thu Jul 22 11:30:47 1999
+   @author     Gerd Lanfermann
+   @desc 
+   Sets the param_active flag to a number: 
+   0    indicates not active (e.g. not specified in the parameter file)
+   1..n priority of the slicing, in the order it is specified in the par file,
+        e.g. "1+log maximal static"
+
+	This is only called from CCTK_PRESTEP
+   @calls 
+   @calledby   
+   @history 
+   @endhistory 
+@@*/
+
+void Einstein_SetNextSlicing(cGH *GH) {
+  DECLARE_CCTK_PARAMETERS
+
+  struct    Einstein_slicing *slic;
+  CCTK_INT  *active_slicing_handle;
+  int i;
+  char *info;
+  
+  int h,doicare = 0;
+  int h_yes = num_slicings, 
+      h_no  = num_slicings,
+      h_egal= num_slicings;
+
+
+  /* Get our grid scalar pointer */
+  i = CCTK_VarIndex("Einstein::active_slicing_handle");
+  if (i<0) CCTK_WARN(0,"INTERNAL ERROR: cannot find grid scalar: active_slicing_handle");
+  active_slicing_handle = (int*)CCTK_VarDataPtrI(GH,0,i);
+
+
+  /* First check if we run mixed NOT: a single slicing specified,
+     that's easy: */
+  if (!CCTK_Equals(slicing,"mixed")) 
+    {
+
+      h = Util_GetHandle(Eslicings,slicing,NULL);
+      if (h<0) CCTK_WARN(0,"INTERNAL ERROR: only slicing not activ! \n");
+      
+      slic = (struct Einstein_slicing *)Util_GetHandledData(Eslicings,h);
+      if (!slic) 
+	CCTK_WARN(0,"INTERNAL ERROR: Slicing registry out of sync! No slicing found!");
+      
+      /* Set the active handle to the only slicing and return*/
+      /*$*active_slicing_handle = h;$*/
+      h_yes = h;
+
+    }
+  /* We run mixed: Loop over all registered handles, 
+     check if active, get the timetoslice function and decide which slicing
+     to use  */
+  else
+    {
+      for (h=0;h<num_slicings;h++) 
+	{
+#ifdef ESLIC_DEBUG	  
+	  printf("HANDLE LOOP: %d/%d\n",h,num_slicings);
+#endif	  
+	  slic = (struct Einstein_slicing *)Util_GetHandledData(Eslicings,h);
+	  if (!slic) 
+	    CCTK_WARN(0,"INTERNAL ERROR: Slicing registry out of sync! No slicing found!");
+	  
+	  /* we have two competing flags: the priority and the 
+	     timetoslice function. We pick two handles for
+	     1) highest priority (lowest param_active number) 
+	        for timetoslice = yes  (1)
+	     2) highest priority (lowest param_active number) 
+                for timetoslice = don't care (0)
+	     set the active_slicing_handle to one of the two handles in this order
+	  */
+#ifdef ESLIC_DEBUG
+	  printf(" h: %d pa: %d yes %d egal %d \n",h,slic->param_active,h_yes, h_egal);
+	  
+#endif
+	  if ((slic->param_active>0)) {
+	    if (slic->timetoslice!=NULL) {
+	      doicare = slic->timetoslice(GH);
+	    }
+	    if ((doicare== SLICING_YES)     &&(slic->param_active<h_yes))  h_yes = h; 
+	    if ((doicare== SLICING_DONTCARE)&&(slic->param_active<h_egal)) h_egal= h;
+	    if ((doicare== SLICING_NO)      &&(slic->param_active<h_no ))  h_no  = h;
+	    
+	  }
+	}
+    }
+  
+  /* First try to set the scalar to the highest priority that cares (h_yes), else
+     to the highest prio, that doesn;t care (h_egal). If that fails, we are in trouble. 
+  */
+
+  if (h_yes!=num_slicings) *active_slicing_handle = h_yes;
+  else 
+    if (h_egal!=num_slicings) *active_slicing_handle = h_egal;
+    else printf("ERROR: no active slicing found: hyes/hegal/hno: %d %d %d \n",
+		h_yes,h_no,h_egal);
+  
+  if (CCTK_Equals(slicing_verbose,"yes")) 
+    { 
+      info = (char*)malloc(256*sizeof(char));
+      slic = (struct  Einstein_slicing *)
+	Util_GetHandledData(Eslicings,*active_slicing_handle);
+      if (!slic) CCTK_WARN(0,"INTERNAL ERROR: Cannot find slicing handle");
+
+      sprintf(info,"next active_slicing_handle: %s  (handle# %d)",
+	      slic->name,*active_slicing_handle);
+      CCTK_INFO(info);
+      if (info) free(info);
+    }
+}
+
+
+
+
+
+   
+    
+    
 
 
 
diff --git a/src/Slicing.h b/src/Slicing.h
index d7031db..2769422 100644
--- a/src/Slicing.h
+++ b/src/Slicing.h
@@ -4,6 +4,16 @@
    @author    Gerd Lanfermann
    @desc 
    Prototypes and structure for Einstein_slicing.
+
+   Einstein_slicing holds:
+   1) the name of the slicing as specified in the parameter 
+      slicing / mixed_slicing
+   2) the param_active number, which also indicates priority
+      when running mixed
+   3) a function pointer: a "time-to-slice" function can
+      be registered which is  called and which can be used
+      to do fancy conditionals, like use slicing XYZ if a 
+      grid functions has collapsed in some region.
    @enddesc 
  @@*/
 
@@ -12,19 +22,30 @@
 extern "C" {
 #endif
 
+  /* Return values for the time-to-slice fucntion */
+#define SLICING_YES       1
+#define SLICING_DONTCARE  0
+#define SLICING_NO       -1
+
 
-struct Einstein_slicing
-{ 
+
+struct Einstein_slicing { 
   char *name;
   int   param_active;
-  void  (*timetoslice)(cGH*);
+  int  (*timetoslice)(cGH*);
 };
 
 int Einstein_RegisterSlicing(const char *name);
 
-int Einstein_SetSliceActive(int handle,int flag);
+int Einstein_SetSlicingPriority(int handle,int flag);
+
+int Einstein_RegisterTimeToSlice(int handle,int (*func)(cGH *));
+
+int Einstein_GetSlicingHandle(const char *name);
+
+void Einstein_SetNextSlicing(cGH *);
+
 
-int Einstein_RegisterTimeToSlice(int handle,void (*func)(cGH *));
 
 #ifdef __cplusplus 
 }
@@ -33,3 +54,4 @@ int Einstein_RegisterTimeToSlice(int handle,void (*func)(cGH *));
 
 
 
+