This commit was generated by cvs2svn to compensate for changes in r2,

which included commits to RCS files with non-trunk default branches.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Noise/trunk@3 936e0995-0e4f-0410-aade-aa24bf7baf3d

9 files changed, 527 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..fa1c53a
--- /dev/null
+++ b/README
@@ -0,0 +1,9 @@
+Cactus Code Thorn IDRandom
+Authors: Denis Pollney <pollney@aei.mpg.de>
+CVS info: $Header$
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+This thorn puts down random initial data on a grid, for performing
+stability tests as proposed by Jeff Winicour.
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..37515db
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,5 @@
+# Interface definition for thorn IDBrBrPETSc
+# $Header$
+
+IMPLEMENTS: Noise
+
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..95006db
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,51 @@
+# Parameter definitions for thorn IDBrBrPETSc
+# $Header$
+
+#------------------------------------------------------------------------------
+# Private:
+#------------------------------------------------------------------------------
+private:
+
+BOOLEAN noisy_id "Add random noise to initial data"
+{
+} "no"
+
+BOOLEAN noisy_bc "Add random noise to initial data"
+{
+} "no"
+
+STRING noisy_id_vars "Initial data variables to modify with noise"
+{
+  .* :: "A regex which matches everything"
+} ""
+
+STRING noisy_bc_vars "Variables to modify with noise at boundary"
+{
+  .* :: "A regex which matches everything"
+} ""
+
+REAL amplitude "Maximum absolute value of random data"
+{
+  0: :: "Positive number"
+} 0.000001
+
+
+#------------------------------------------------------------------------------
+# parameters specific to bssn tests (potentially obsolete)
+#------------------------------------------------------------------------------
+BOOLEAN evolveK "Apply randomization to trK"
+{
+} "no"
+
+BOOLEAN evolveShift "Apply randomization to the shift components"
+{
+} "no"
+
+BOOLEAN evolveLapse "Apply randomization to the lapse"
+{
+} "no"
+
+BOOLEAN evolveGamma "Apply randomization to the Gamma components"
+{
+} "no"
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..eaa9123
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,16 @@
+# Schedule definitions for thorn IDBrBrPETSc
+# $Header$
+
+if (noisy_id) {
+  SCHEDULE id_noise AT CCTK_POSTINITIAL
+  {
+    LANG: C
+  } "Add noise to initial data"
+}
+
+if (noisy_bc) {
+  SCHEDULE bc_noise AT CCTK_POSTSTEP
+  {
+     LANG: C
+  } "Add noise to boundary condition"
+}
\ No newline at end of file
diff --git a/src/bc_noise.c b/src/bc_noise.c
new file mode 100644
index 0000000..0336cd8
--- /dev/null
+++ b/src/bc_noise.c
@@ -0,0 +1,270 @@
+
+#include "noise.h"
+
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_FortranString.h"
+
+#include "Symmetry.h"
+
+
+
+static int ApplyBndNoise (const cGH *GH,
+			  int stencil_dir,
+			  const int *stencil_alldirs,
+			  int dir,
+			  int first_var,
+			  int num_vars);
+
+int
+BndNoiseGI (const cGH *GH, const int *stencil, int gi)
+{
+  int first_vi, retval;
+
+  first_vi = CCTK_FirstVarIndexI (gi);
+  if (first_vi >= 0)
+    {
+      retval = ApplyBndNoise (GH, -1, stencil, 0, first_vi,
+			      CCTK_NumVarsInGroupI (gi));
+    }
+  else
+    {
+      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "Invalid group index %d in BndNoiseGI", gi);
+      retval = -1;
+    }
+
+  return (retval);
+}
+
+void
+CCTK_FCALL CCTK_FNAME (BndNoiseGI) (int *ierr, const cGH *GH,
+				    const int *stencil, const int *gi)
+{
+  *ierr = BndNoiseGI (GH, stencil, *gi);
+}
+
+
+int
+BndNoiseGN (const cGH *GH, const int *stencil, const char *gn)
+{
+  int gi, retval;
+
+  gi = CCTK_GroupIndex (gn);
+  if (gi >= 0)
+    {
+      retval = BndNoiseGI (GH, stencil, gi);
+    }
+  else
+    {
+      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "Invalid group name '%s' in BndNoiseGN", gn);
+      retval = -1;
+    }
+
+  return (retval);
+}
+
+void CCTK_FCALL CCTK_FNAME (BndNoiseGN)
+     (int *ierr,
+      const cGH *GH,
+      const int *stencil,
+      ONE_FORTSTRING_ARG)
+{
+  ONE_FORTSTRING_CREATE (gn)
+    *ierr = BndNoiseGN (GH, stencil, gn);
+  free (gn);
+}
+
+
+
+
+
+static int ApplyBndNoise (const cGH *GH,
+			  int stencil_dir,
+			  const int *stencil_alldirs,
+			  int dir,
+			  int first_var,
+			  int num_vars)
+{
+  int i, j, k;
+  int var, vtypesize, gindex, gdim, timelvl;
+  int doBC[2*MAXDIM], dstag[MAXDIM], lsh[MAXDIM], lssh[MAXDIM],
+    stencil[MAXDIM];
+  SymmetryGHex *sGHex;
+  CCTK_REAL amplitude;
+  CCTK_REAL offset;
+  int type;
+
+  stencil[0] = 1; stencil[1] = 1; stencil[2] = 1;
+
+  /* get the group index of the variables */
+  gindex = CCTK_GroupIndexFromVarI (first_var);
+
+  amplitude = *((CCTK_REAL*) CCTK_ParameterGet("amplitude", "Noise", &type))
+    / RAND_MAX;
+  offset = 0.5 * amplitude;
+
+  /* get the number of dimensions and the size of the variables' type */
+  gdim      = CCTK_GroupDimI (gindex);
+  vtypesize = CCTK_VarTypeSize (CCTK_VarTypeI (first_var));
+
+  /* make sure we can deal with this number of dimensions */
+  if (gdim > MAXDIM)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "ApplyBndNoise: Variable dimension of %d not supported", gdim);
+      return (-1);
+    }
+
+  /* check the direction parameter */
+  if (abs (dir) > gdim)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "ApplyBndNoise: direction %d greater than dimension %d",
+		  dir, gdim);
+      return (-2);
+    }
+
+  /* initialize arrays for variables with less dimensions than MAXDIM
+     so that we can use the INDEX_3D macro later on */
+  for (i = gdim; i < MAXDIM; i++)
+    {
+      lsh[i]  = 0;
+      lssh[i] = 1;
+    }
+
+  /* get the directional staggering of the group */
+  CCTK_GroupStaggerDirArrayGI (dstag, gdim, gindex);
+
+  /* get the current timelevel */
+  timelvl = 0;
+
+  /* see if we have a symmetry array */
+  sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry");
+
+
+  /* now loop over all variables */
+  for (var = first_var; var < first_var + num_vars; var++)
+    {
+      /* Apply condition if:
+	 + boundary is not a symmetry boundary (no symmetry or unset(=unsed))
+	 + boundary is a physical boundary
+	 + have enough grid points
+      */
+      memset (doBC, 1, sizeof (doBC));
+      if (sGHex)
+	{
+	  for (i = 0; i < 2 * gdim; i++)
+	    {
+	      doBC[i] = sGHex->GFSym[var][i] == GFSYM_NOSYM ||
+		sGHex->GFSym[var][i] == GFSYM_UNSET;
+	    }
+	}
+      for (i = 0; i < gdim; i++)
+	{
+	  lsh[i]       = GH->cctk_lsh[i];
+	  lssh[i]      = GH->cctk_lssh[CCTK_LSSH_IDX (dstag[i], i)];
+	  doBC[i*2]   &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2];
+	  doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1];
+	  if (dir != 0)
+	    {
+	      doBC[i*2]   &= (dir < 0 && (i + 1 == abs (dir)));
+	      doBC[i*2+1] &= (dir > 0 && (i + 1 == abs (dir)));
+	    }
+	}
+
+      /* now apply the boundaries face by face */
+      if (gdim > 0)
+	{
+#ifdef DEBUG_BOUNDARY
+	  if (doBC[0])
+	    {
+	      printf("Boundary: Applying lower x noise to boundary\n");
+	    }
+	  if (doBC[1])
+	    {
+	      printf("Boundary: Applying upper x noise to boundary\n");
+	    }
+#endif /* DEBUG_BOUNDARY */
+	  /* lower x */
+	  NOISY_BOUNDARY (doBC[0], stencil[0], lssh[1], lssh[2],
+			  i, j, k, rand);
+	  /* upper x */
+	  NOISY_BOUNDARY (doBC[1], stencil[0], lssh[1], lssh[2],
+			  lssh[0]-i-1, j, k, rand);
+
+	}
+      if (gdim > 1)
+
+	{
+#ifdef DEBUG_BOUNDARY
+	  if (doBC[2])
+	    {
+	      printf("Boundary: Applying lower y noise to boundary\n");
+	    }
+	  if (doBC[3])
+	    {
+	      printf("Boundary: Applying upper y noise to boundary\n");
+	    }
+#endif /* DEBUG_BOUNDARY */
+	  /* lower y */
+	  NOISY_BOUNDARY (doBC[2], lssh[0], stencil[1], lssh[2],
+			  i, j, k, rand);
+	  /* upper y */
+	  NOISY_BOUNDARY (doBC[3], lssh[0], stencil[1], lssh[2],
+			  i, lssh[1]-j-1, k, rand);
+	}
+      if (gdim > 2)
+	{
+#ifdef DEBUG_BOUNDARY
+	  if (doBC[4])
+	    {
+	      printf("Boundary: Applying lower z noise to boundary\n");
+	    }
+	  if (doBC[5])
+	    {
+	      printf("Boundary: Applying upper z noise to boundary\n");
+	    }
+#endif /* DEBUG_BOUNDARY */
+	  /* lower z */
+	  NOISY_BOUNDARY (doBC[4], lssh[0], lssh[1], stencil[2],
+			  i, j, k, rand);
+	  /* upper z */
+	  NOISY_BOUNDARY (doBC[5], lssh[0], lssh[1], stencil[2],
+			  i, j, lssh[2]-k-1, rand);
+	}
+    }
+
+  return(0);
+}
+
+
+void
+add_bc_noise_to_group (int idx, const char* optstring, void* cctkGH)
+{
+  int i, j, k, ijk;
+  CCTK_REAL* data;
+  cGH* GH = cctkGH;
+  int sw[3];
+
+  sw[0] = 1;
+  sw[1] = 1;
+  sw[2] = 1;
+
+  BndNoiseGI(GH, sw, idx);
+}
+
+
+void
+bc_noise(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+    DECLARE_CCTK_PARAMETERS
+
+    if (CCTK_TraverseString(noisy_bc_vars, add_bc_noise_to_group, cctkGH,
+			    CCTK_GROUP_OR_VAR) < 0)
+      {
+	CCTK_WARN (1, "Failed to parse 'Noise::noisy_bc_vars' parameter");
+      }
+}
diff --git a/src/id_noise.c b/src/id_noise.c
new file mode 100644
index 0000000..e94ec2b
--- /dev/null
+++ b/src/id_noise.c
@@ -0,0 +1,62 @@
+/*
+  IDRandom
+
+  Denis Pollney <pollney@aei.mpg.de>
+  14 February 2002
+
+  Put random data into the evolution variables for robust stability
+  testing a-la-winicour.
+
+  $Header$
+*/
+
+#include <stdlib.h>
+
+#include "noise.h"
+
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void
+add_noise_to_group (int idx, const char* optstring, void* cctkGH)
+{
+  int i, j, k, ijk;
+  CCTK_REAL* data;
+  CCTK_REAL amplitude;
+  CCTK_REAL offset;
+  cGH* GH = cctkGH;
+  int type;
+
+  data = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, idx);
+
+  amplitude = *((CCTK_REAL*) CCTK_ParameterGet("amplitude", "Noise", &type))
+    / RAND_MAX;
+  offset = 0.5 * amplitude;
+
+  for (k=1; k< GH->cctk_lsh[2]-1; ++k)
+    {
+      for (j=1; j< GH->cctk_lsh[1]-1; ++j)
+	{
+	  for (i=1; i< GH->cctk_lsh[0]-1; ++i)
+	    {
+	      ijk = CCTK_GFINDEX3D(GH, i, j, k);
+
+	      data[ijk] += RAND_VAL;
+      	    }
+	}
+    }
+}
+
+
+void
+id_noise(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  if (CCTK_TraverseString(noisy_id_vars, add_noise_to_group, cctkGH,
+			  CCTK_GROUP_OR_VAR) < 0)
+    {
+      CCTK_WARN (1, "Failed to parse 'IDRandom::noisy_id_vars' parameter");
+    }
+}
diff --git a/src/id_noise.h b/src/id_noise.h
new file mode 100644
index 0000000..6302306
--- /dev/null
+++ b/src/id_noise.h
@@ -0,0 +1,53 @@
+/* $Header$ */
+
+#ifndef NOISE_H
+#define NOISE_H
+
+#include <stdlib.h>
+#include "cctk.h"
+
+
+/* constants */
+
+#define TRUE 1
+#define FALSE 0
+
+#define MAXDIM  3
+#
+
+/* macros */
+
+#define RAND_VAL random()*amplitude - offset
+
+#define INDEX_3D(lsh, i, j, k)      ((i) + (lsh)[0]*((j) + (lsh)[1]*(k)))
+
+
+#define NOISY_BOUNDARY(doBC,                                                 \
+                        iend, jend, kend,                                     \
+                        ii, jj, kk,                                           \
+                        rand) \
+{                                                                             \
+if (doBC)                                        \
+  {                                              \
+    srandom(time(0)); \
+    for (k = 0; k < kend; k++)                   \
+    {                                                                         \
+      for (j = 0; j < jend; j++)                                              \
+      {                                                                       \
+        for (i = 0; i < iend; i++)                                            \
+        {                                                                     \
+          int _index = CCTK_GFINDEX3D(GH, ii, jj, kk);          \
+          CCTK_REAL* v= (CCTK_REAL *) GH->data[var][timelvl] + _index; \
+          v[0] += RAND_VAL;  \
+        }                                                  \
+      }                                                                 \
+    }                                                                   \
+  }                                                      \
+}
+
+
+
+#endif /* !define(NOISE_H) */
+
+
+
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..b4a36ec
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn IDRandom
+# $Header$
+
+# Source files in this directory
+SRCS = id_noise.c bc_noise.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/noise.h b/src/noise.h
new file mode 100644
index 0000000..6302306
--- /dev/null
+++ b/src/noise.h
@@ -0,0 +1,53 @@
+/* $Header$ */
+
+#ifndef NOISE_H
+#define NOISE_H
+
+#include <stdlib.h>
+#include "cctk.h"
+
+
+/* constants */
+
+#define TRUE 1
+#define FALSE 0
+
+#define MAXDIM  3
+#
+
+/* macros */
+
+#define RAND_VAL random()*amplitude - offset
+
+#define INDEX_3D(lsh, i, j, k)      ((i) + (lsh)[0]*((j) + (lsh)[1]*(k)))
+
+
+#define NOISY_BOUNDARY(doBC,                                                 \
+                        iend, jend, kend,                                     \
+                        ii, jj, kk,                                           \
+                        rand) \
+{                                                                             \
+if (doBC)                                        \
+  {                                              \
+    srandom(time(0)); \
+    for (k = 0; k < kend; k++)                   \
+    {                                                                         \
+      for (j = 0; j < jend; j++)                                              \
+      {                                                                       \
+        for (i = 0; i < iend; i++)                                            \
+        {                                                                     \
+          int _index = CCTK_GFINDEX3D(GH, ii, jj, kk);          \
+          CCTK_REAL* v= (CCTK_REAL *) GH->data[var][timelvl] + _index; \
+          v[0] += RAND_VAL;  \
+        }                                                  \
+      }                                                                 \
+    }                                                                   \
+  }                                                      \
+}
+
+
+
+#endif /* !define(NOISE_H) */
+
+
+