Courant condition implemented for timestep, as described in documentation.

This hasn't been very well tested yet.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Time/trunk@11 5633253d-7678-4964-a54d-f87795f8ee59

7 files changed, 98 insertions, 33 deletions

diff --git a/doc/documention.tex b/doc/documention.tex
index 4f32b81..316d514 100644
--- a/doc/documention.tex
+++ b/doc/documention.tex
@@ -10,19 +10,57 @@
 
 \section{Purpose}
 
-This, currently very simple, thorn provides a routine for calculating
-the timestep for an evolution based on the spatial Cartesian grid spacing.
-This thorn will be extended in the future to include features like
-a dynamic Courant condition. The thorn works for up to 3 spatial 
-dimensions.
-
+This thorn provides a routines for calculating
+the timestep for an evolution based on the spatial Cartesian grid spacing and
+a wave speed. 
 \section{Comments}
 
-The timestep, which is passed into the Cactus variable {\tt cctk\_delta\_time},
- is calculated by the algorithm
+There are currently two methods for calculating the timestep, either the
+simple courant method, or a dynamic courant method. The method is chosen
+using the keyword parameter {\tt time::method}
+The timestep,  is passed into the Cactus variable {\tt cctk\_delta\_time}.
+Both timesteps are based on the Courant condition, which states that for 
+numerical stability, the chosen timestep should satisfy
+$$
+\Delta t \le \mbox{min}(\Delta x^i)/\mbox{wave speed}/\sqrt(\mbox{dim})
+$$
+
+The two methods currently implemented are:
+\begin{itemize}
+
+\item {\tt time::method = "standard"} (this is the default). The timestep is calculated once
+   at the start of the run in the {\tt BASEGRID} timebin, and is then held static
+   throughout the run. The algorithm, which uses the parameter {\tt time::dtfac} is
 $$
 \Delta t = \mbox{\tt dtfac} * \mbox{min}(\Delta x^i)
 $$
+	Note that the parameter {\tt dtfac} should take into account the dimension
+        of the space being used, and the wave speed.
+
+\item {\tt time::method = "courant"}. The timestep is calculated dynamically at the
+	start of each iteration in the {\tt PRESTEP} timebin. The algorithm is
+$$
+\Delta t = \mbox{\tt courant\_fac} * \mbox{min}(\Delta x^i)/\mbox{maximum wavespeed}/\sqrt(\mbox{dim})
+$$
+
+
+\item {\tt time::method = "courant\_time"}. The timestep is calculated dynamically at the
+    start of each iteration in the {\tt PRESTEP} timebin. The algorithm is
+$$
+\Delta t = \mbox{\tt courant\_fac} * \mbox{minimum time to cross a zone}/\sqrt(\mbox{dim})
+$$
+
+\section{Technical Details}
+
+If a dynamic {\tt courant} condition is selected, a thorn must set the protected variable
+{\tt courant\_wave\_speed} for the maximum wave speed before Time sets the timestep. 
+
+If a dynamic {\tt courant_time} condition is selected, a thorn must set the protected variable
+{\tt courant\_time} for the minimum time for a wave to cross a zone
+ before Time sets the timestep. 
+
+\end{itemize}
+
 
 
 
diff --git a/interface.ccl b/interface.ccl
index 1009f1e..aea0edc 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -3,3 +3,16 @@
 
 implements: time
 
+protected:
+
+REAL speedvars type=SCALAR
+{
+  wave_speed
+} "Speed to use for Courant condition"
+
+private:
+
+REAL couranttemps type=SCALAR
+{
+  courant_dt
+} "Variable just for output"
\ No newline at end of file
diff --git a/param.ccl b/param.ccl
index e172a44..317be03 100644
--- a/param.ccl
+++ b/param.ccl
@@ -8,3 +8,23 @@ REAL dtfac "The standard timestep condition dt = dtfac*max(delta_space)"
  0:* :: "Probably only makes sense to be bigger than zero"
 } 0.5 
 
+REAL courant_fac "The courant timestep condition dt = courant_fac*max(delta_space)/speed/sqrt(dim)"
+{
+ 0:* :: "Probably only makes sense to be bigger than zero"
+} 0.9
+
+KEYWORD courant_method "Method for calculating timestep"
+{
+ "standard" :: "Courant condition at BASEGRID"
+ "courant"  :: "Courant condition at PRESTEP (using wavespeed)"
+ "courant_time" :: "Courant condition at PRESTEP (using min time)"
+} "standard"
+
+BOOLEAN courant_outonly "Only output courant timestep?"
+{
+} "no"
+
+INT outcourant_every "How often to output courant timestep"
+{
+ 0:* :: "Zero means no output"
+} 0
\ No newline at end of file
diff --git a/schedule.ccl b/schedule.ccl
index b8aca19..33a3133 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,10 +1,21 @@
 # Schedule definitions for thorn Time
 # $Header$
 
-schedule Time_Simple at CCTK_BASEGRID after CartGrid3D
+if (CCTK_Equals(courant_method,"standard"))
 {
-  LANG: C
-} "Set timestep based on speed one Courant condition"
+  schedule Time_Simple at CCTK_BASEGRID after CartGrid3D
+  {
+     LANG: C
+  } "Set timestep based on Courant condition"
+}
 
+if (CCTK_Equals(courant_method,"courant") || CCTK_Equals(courant_method,"courant_time"))
+{
+  schedule Time_Courant at CCTK_PRESTEP 
+  {
+    LANG: C
+  } "Reset timestep each iteration"
+}
+  
 
 
diff --git a/src/Courant.c b/src/Courant.c
index 7bac44c..d0c36ef 100644
--- a/src/Courant.c
+++ b/src/Courant.c
@@ -26,9 +26,6 @@ void Time_Courant(CCTK_CARGUMENTS)
 
   int min_handle,ierr;
   
-  /* Get the handle for the MINIMUM operation */
-  /*$ierr = CCTK_ReductionArrayHandle("minimum");$*/
-
   /* Calculate the minimum grid spacing */
   if (cctk_dim>=1)
   {
@@ -52,15 +49,6 @@ void Time_Courant(CCTK_CARGUMENTS)
     CCTK_WARN(0,"Time Step not defined for greater than 4 dimensions");
   }
 
-  /* do the global minimum on local min_spacing to tmp and reassign that
-     to min_spacing */
-
-  /*$printf("Courant1: %d \n",min_spacing);
-  ierr = CCTK_ReduceLocalScalar(GH,-1,min_handle,min_spacing, tmp, CCTK_VAR_REAL);
-  printf("Courant2: %d \n",min_spacing,tmp);
-  min_spacing = tmp;
-  printf("Courant3: %d \n",min_spacing,tmp);$*/
-
   /* Calculate the courant timestep */
   courant_speed = *wave_speed;
   *courant_dt = courant_fac/courant_speed/sqrt((CCTK_REAL )cctk_dim);
@@ -85,7 +73,7 @@ void Time_Courant(CCTK_CARGUMENTS)
   }
   else
   {
-    cctkGH->cctk_delta_time = dtfac*min_spacing/courant_wave_speed;
+    cctkGH->cctk_delta_time = dtfac*min_spacing;
   }
 
 }
diff --git a/src/Simple.c b/src/Simple.c
index ac801e5..e09ad6a 100644
--- a/src/Simple.c
+++ b/src/Simple.c
@@ -13,8 +13,6 @@
 #include "cctk_arguments.h"
 #include "cctk_parameters.h"
 
-#include "cctk_WarnLevel.h"
-
 void Time_Simple(CCTK_CARGUMENTS)
 {
   DECLARE_CCTK_PARAMETERS
@@ -23,6 +21,7 @@ void Time_Simple(CCTK_CARGUMENTS)
   CCTK_REAL min_spacing;
   char *message;
 
+  /* Calculate the minimum grid spacing */
   if (cctk_dim>=1)
   {
     min_spacing = cctk_delta_space[0];
@@ -42,15 +41,11 @@ void Time_Simple(CCTK_CARGUMENTS)
 
   if (cctk_dim>=4)
   {
-    CCTK_WARN(0,"Time Step now defined for greater than 4 dimensions");
+    CCTK_WARN(0,"Time Step not defined for greater than 4 dimensions");
   }
 
+  /* Calculate the timestep */
   cctkGH->cctk_delta_time = dtfac*min_spacing;
-  
-  message = (char *)malloc(1024*sizeof(char));
-  sprintf(message,"Time step set to %f",cctkGH->cctk_delta_time);
-  CCTK_INFO(message);
-  free(message);
 
 }
 
diff --git a/src/make.code.defn b/src/make.code.defn
index 8ebd1fb..a37cea0 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Simple.c 
+SRCS = Simple.c Courant.c
 
 # Subdirectories containing source files
 SUBDIRS =