initial commit.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Norms/trunk@2 f4913095-0e4f-0410-abea-a123d184f1f3

3 files changed, 355 insertions, 0 deletions

diff --git a/src/Compute_Norm.c b/src/Compute_Norm.c
new file mode 100644
index 0000000..083f498
--- /dev/null
+++ b/src/Compute_Norm.c
@@ -0,0 +1,184 @@
+/* XXX cleanup the include files ! */
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+
+static const char * rcsid = "$Header$";
+
+CCTK_FILEVERSION(Norms_Compute_Norm_c);
+
+
+void Norms_Compute_Norms (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  int vind; /* variable index */
+  CCTK_REAL *data;
+  int index,indexp, i,j,k;
+  int ierr;
+  int istart, iend, jstart, jend, kstart, kend;
+  int npoints;
+  CCTK_REAL hd, dx,dy,dz;
+  CCTK_REAL deriv;
+  int target_proc,reduction_handle;
+  CCTK_REAL result;
+
+  if (verbose>1)
+    CCTK_INFO("computing norms");
+
+  *norm=0.;
+  result=0.;
+
+  dx = CCTK_DELTA_SPACE(0);
+  dy = CCTK_DELTA_SPACE(1);
+  dz = CCTK_DELTA_SPACE(2);
+
+  /* XXX ghost size */
+  istart = 1;
+  jstart = 1;
+  kstart = 1;
+
+  iend = cctk_lsh[0]-1;
+  jend = cctk_lsh[1]-1;
+  kend = cctk_lsh[2]-1;
+
+  for (k=kstart;k<kend;k++) {
+    for (j=jstart;j<jend;j++) {
+      for (i=istart;i<iend;i++) {
+        index=CCTK_GFINDEX3D(cctkGH,i,j,k);
+        diff_term[index]=0.;
+      }
+    }
+  }
+
+  target_proc = 0;
+  reduction_handle = CCTK_ReductionArrayHandle("norm2");
+  if (reduction_handle < 0)
+       CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                        "can't get sum-reduction handle!");
+
+  if (verbose>2)
+    fprintf(stderr,"   nr1stvars %d nr2ndvars %d\n",
+                       *nr1stvars,*nr2ndvars);
+
+  /* L2-norm of 1st order variables */
+  for (vind=0; vind< *nr1stvars; vind++) {
+    ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1,
+                  CCTK_VARIABLE_REAL,&result,1,varindices_1st[vind]);
+    if (ierr<0)
+        CCTK_WARN(1,"reduction failed");
+    *norm=*norm+result*result;
+  }
+
+  if (verbose>2)
+    fprintf(stderr,"  after 1st order vars norm=%16.8g\n",*norm);
+
+  /* L2-norm of 2nd order variables */
+  for (vind=0; vind< *nr2ndvars; vind++) {
+    ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1,
+                  CCTK_VARIABLE_REAL,&result,1,varindices_2nd[vind]);
+    if (ierr<0)
+      CCTK_WARN(1,"reduction failed");
+    *norm=*norm+result*result;
+  }
+
+  if (verbose>2)
+    fprintf(stderr,"  after 2nd order vars norm=%16.8g\n",*norm);
+
+  /* L2-norm of D_{+i} terms */
+  /* d_x */
+  for (vind=0; vind< *nr2ndvars; vind++) {
+    data=CCTK_VarDataPtrI(cctkGH,0,varindices_2nd[vind]);
+    for (k=kstart;k<kend;k++) {
+      for (j=jstart;j<jend;j++) {
+        for (i=istart;i<iend;i++) {
+          index =CCTK_GFINDEX3D(cctkGH,i,j,k);
+          indexp=CCTK_GFINDEX3D(cctkGH,i+1,j,k);
+          diff_term[index]=(data[indexp]-data[index])/dx;
+        }
+      }
+    }
+    ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1,
+                  CCTK_VARIABLE_REAL,&result,1,
+                  CCTK_VarIndex("norms::diff_term"));
+    if (ierr<0)
+      CCTK_WARN(1,"reduction failed");
+    *norm=*norm+result*result;
+  }
+
+  if (verbose>2)
+    fprintf(stderr,"  after dx 2nd order terms norm=%16.8g\n",*norm);
+
+  /* d_y */
+  for (vind=0; vind< *nr2ndvars; vind++) {
+    data=CCTK_VarDataPtrI(cctkGH,0,varindices_2nd[vind]);
+    for (k=kstart;k<kend;k++) {
+      for (j=jstart;j<jend;j++) {
+        for (i=istart;i<iend;i++) {
+          index =CCTK_GFINDEX3D(cctkGH,i,j,k);
+          indexp=CCTK_GFINDEX3D(cctkGH,i,j+1,k);
+          diff_term[index]=(data[indexp]-data[index])/dy;
+        }
+      }
+    }
+    ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1,
+                  CCTK_VARIABLE_REAL,&result,1,
+                  CCTK_VarIndex("norms::diff_term"));
+    if (ierr<0)
+      CCTK_WARN(1,"reduction failed");
+    *norm=*norm+result*result;
+  }
+
+  if (verbose>2)
+    fprintf(stderr,"  after dy 2nd order terms norm=%16.8g\n",*norm);
+
+  /* d_z */
+  for (vind=0; vind< *nr2ndvars; vind++) {
+    data=CCTK_VarDataPtrI(cctkGH,0,varindices_2nd[vind]);
+    for (k=kstart;k<kend;k++) {
+      for (j=jstart;j<jend;j++) {
+        for (i=istart;i<iend;i++) {
+          index =CCTK_GFINDEX3D(cctkGH,i,j,k);
+          indexp=CCTK_GFINDEX3D(cctkGH,i,j,k+1);
+          diff_term[index]=(data[indexp]-data[index])/dz;
+        }
+      }
+    }
+    ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1,
+                  CCTK_VARIABLE_REAL,&result,1,
+                  CCTK_VarIndex("norms::diff_term"));
+    if (ierr<0)
+      CCTK_WARN(1,"reduction failed");
+    *norm=*norm+result*result;
+  }
+
+  if (verbose>2)
+    fprintf(stderr,"  after dz 2nd order terms norm=%16.8g\n",*norm);
+
+  if (verbose>2)
+    fprintf(stderr,"  after 2nd order vars (incl. D_{+i}): norm=%16.8g\n",
+                   *norm);
+  if (verbose>2)
+    fprintf(stderr,"final norm^2=%16.8g on this CPU\n",*norm);
+
+  *norm=sqrt(*norm);
+
+  /* multiply with h^d - see discussion above Eq. (12) in
+   * gr-qc/XXX */
+  hd=dx*dy*dz;
+  *norm=*norm*hd;
+
+  CCTK_VInfo(CCTK_THORNSTRING,"norm= %16.8g",*norm);
+
+}
diff --git a/src/Setup_Vars.c b/src/Setup_Vars.c
new file mode 100644
index 0000000..b09b594
--- /dev/null
+++ b/src/Setup_Vars.c
@@ -0,0 +1,163 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+
+static const char * rcsid = "$Header$";
+
+CCTK_FILEVERSION(Norms_Setup_Vars_c);
+
+
+struct norms_opts {
+  int active;
+  int vi; /* variable index */
+  int norm_type;
+};
+
+
+static void getopt (int const idx,
+                    const char * const optstring,
+                    void * const opts)
+{
+  struct norms_opts * norms_opts;
+  int table;
+  int cnt;
+  int ierr;
+  int norm_type;
+
+
+  assert (idx >= 0 && idx < CCTK_NumVars());
+  assert (opts);
+
+  norms_opts = &((struct norms_opts *)opts)[idx];
+
+  assert (! norms_opts->active);
+
+  if (optstring) {
+    assert (optstring);
+    table = Util_TableCreateFromString (optstring);
+    if (table < 0) {
+      char * fullname = CCTK_FullName (idx);
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The variable \"%s\" is ignored because it has an invalid option specification in the parameter \"Norms::gridfunctions...\"",
+                  fullname);
+      free (fullname);
+      return;
+    }
+    assert (table >= 0);
+  }
+
+  norms_opts->active = 1;
+  norms_opts->vi = idx;
+
+  if (optstring) {
+    /* XXX norm_type is not used currently !*/
+    cnt = Util_TableGetInt (table, &norm_type , "norm_type");
+
+    if (cnt < 0) {
+      norms_opts->norm_type = -100; /* XXX magic value */
+    } else {
+      norms_opts->norm_type = norm_type;
+    }
+
+    ierr = Util_TableDestroy (table);
+    assert (!ierr);
+  }
+
+  /* Don't compute norms if not output was requested */
+  if (norms_opts->active) {
+    if (norms_opts->norm_type==-100) {
+      norms_opts->active = 0;
+    }
+  }
+
+}
+
+
+
+
+
+void Norms_Setup_Vars (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  int nvars;
+  struct norms_opts * norms_opts_1st;
+  struct norms_opts * norms_opts_2nd;
+
+  int n,i;
+
+  int ierr;
+
+  /* init */
+  *nr1stvars=0;
+  *nr2ndvars=0;
+
+  if (verbose>6)
+    CCTK_INFO("here we are in Setup_vars");
+
+  if (verbose>0)
+    CCTK_VInfo(CCTK_THORNSTRING,"Starting Norms Computation at time %f",
+        cctkGH->cctk_time);
+  if (cctk_iteration % out_every != 0)
+  {
+    *do_nothing=1;
+    return;
+  }
+  else
+    *do_nothing=0;
+
+
+  nvars = CCTK_NumVars();
+  assert (nvars >= 0);
+  norms_opts_1st = malloc (nvars * sizeof *norms_opts_1st);
+  norms_opts_2nd = malloc (nvars * sizeof *norms_opts_2nd);
+  assert (norms_opts_1st);
+  assert (norms_opts_2nd);
+  for (n=0; n<nvars; ++n) {
+    norms_opts_1st[n].active = 0;
+    norms_opts_2nd[n].active = 0;
+  }
+  ierr = CCTK_TraverseString
+    (gridfunctions_1st, getopt, norms_opts_1st, CCTK_GROUP_OR_VAR);
+  assert (ierr >= 0);
+  ierr = CCTK_TraverseString
+    (gridfunctions_2nd, getopt, norms_opts_2nd, CCTK_GROUP_OR_VAR);
+  assert (ierr >= 0);
+
+  i=0;
+  if (verbose>0)
+    CCTK_INFO("We will compute norms for the following variables");
+  for (n=0; n<nvars; ++n) {
+    if (norms_opts_1st[n].active) {
+      varindices_1st[i]=norms_opts_1st[n].vi;
+      if (verbose>0) {
+          CCTK_VInfo(CCTK_THORNSTRING,"    %s (1st order var)",
+                        CCTK_FullName(varindices_1st[i]));
+      }
+      *nr1stvars=*nr1stvars+1;
+    }
+    if (norms_opts_2nd[n].active) {
+      varindices_2nd[i]=norms_opts_2nd[n].vi;
+      if (verbose>0) {
+          CCTK_VInfo(CCTK_THORNSTRING,"    %s (2nd order var)",
+                        CCTK_FullName(varindices_2nd[i]));
+      }
+      *nr2ndvars=*nr2ndvars+1;
+    }
+  }
+  if (verbose>2)
+    fprintf(stderr,"   nr1stvars %d nr2ndvars %d\n",*nr1stvars,*nr2ndvars);
+
+  free (norms_opts_1st);
+  free (norms_opts_2nd);
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..67bfb8e
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn Norms
+# $Header$
+
+# Source files in this directory
+SRCS = Setup_Vars.c Compute_Norm.c
+
+# Subdirectories containing source files
+SUBDIRS =