From a02eeb7ccbd75dc20ed00286e32e6d9fe4189a18 Mon Sep 17 00:00:00 2001 From: schnetter Date: Wed, 7 Dec 2005 00:44:45 +0000 Subject: Clean up the code: Replace several magic numbers by named constants DIM, MAX_RANK, MAX_TIME_LEVELS, and MAX_SPATIAL_DERIV_ORDER Handle the case where time derivatives are calculated while interpolating: Introduce a new constant MAX_TIME_DERIV_ORDER. Add a new rank to the variables thetensor and thevars. Take the option table key time_deriv_order into account. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/RotatingSymmetry90/trunk@34 c3c03602-0f4f-0410-b3fa-d2c81c8a7dc5 --- src/interpolate.c | 227 +++++++++++++++++++++++++++++++++++------------------- 1 file changed, 149 insertions(+), 78 deletions(-) diff --git a/src/interpolate.c b/src/interpolate.c index be75f3e..2bdb183 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -1,7 +1,6 @@ /* $Header$ */ #include -#include #include #include #include @@ -14,6 +13,21 @@ +/* This is pretty hard coded into all the tensor types and cannot be + changed easily. */ +#define DIM 3 /* spatial dimension */ + + + +/* These can be increased if necessary, but not all necessary tensor + types may be defined . */ +#define MAX_RANK 9 /* maximum tensor rank */ +#define MAX_TIME_LEVELS 3 /* maximum number of time levels */ +#define MAX_SPATIAL_DERIV_ORDER 2 /* maximum spatial derivative order */ +#define MAX_TIME_DERIV_ORDER 2 /* maximum time derivative order */ + + + /* A tensor description */ struct tensor { @@ -209,6 +223,36 @@ static struct tensor symmtensor4 = { +static int +ipow (int base, int expo) +{ + int res; + assert (expo >= 0); + res = expo & 1 ? base : 1; + while (expo >>= 1) + { + base *= base; + if (expo & 1) res *= base; + } + return res; +} + + + +static int +ilog (int res, int const base) +{ + int expo; + assert (base > 0); + assert (res >= 0); + for (expo = 0; res > 0; ++ expo) { + res /= base; + } + return expo; +} + + + /* Ensure that all tensor declarations are internally consistent */ static void CheckTensorType (struct tensor const * restrict const atensor) @@ -218,7 +262,7 @@ CheckTensorType (struct tensor const * restrict const atensor) assert (atensor->dim>=0); assert (atensor->rank>=0); assert (atensor->ncomps>=0); - assert (atensor->ncomps == floor(pow(atensor->dim, atensor->rank) + 0.5)); + assert (atensor->ncomps == ipow (atensor->dim, atensor->rank)); assert (atensor->nvars>=0 && atensor->nvars<=atensor->ncomps); assert (atensor->vars); for (i=0; incomps; ++i) { @@ -307,20 +351,20 @@ Rot90_CheckTensorTypes (CCTK_ARGUMENTS) || CCTK_EQUALS (tensortypealias, "d")) { /* vector */ - assert (numvars == 3); + assert (numvars == DIM); } else if (CCTK_EQUALS (tensortypealias, "4u") || CCTK_EQUALS (tensortypealias, "4d")) { /* 4-vector */ - assert (numvars == 4); + assert (numvars == DIM+1); } else if (CCTK_EQUALS (tensortypealias, "uu_sym") || CCTK_EQUALS (tensortypealias, "dd_sym")) { /* symmetric tensor */ - assert (numvars == 6); + assert (numvars == DIM*(DIM+1)/2); } else if (CCTK_EQUALS (tensortypealias, "4uu_sym") || CCTK_EQUALS (tensortypealias, "4dd_sym")) { /* symmetric 4-tensor */ - assert (numvars == 10); + assert (numvars == (DIM+1)*(DIM+2)/2); } else { char * groupname = CCTK_GroupName(gi); assert (groupname); @@ -354,25 +398,28 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, { cGH const * restrict const cctkGH = cctkGH_; - CCTK_POINTER new_interp_coords[3]; + CCTK_POINTER new_interp_coords[DIM]; CCTK_POINTER * restrict new_output_arrays; CCTK_INT new_faces; CCTK_INT * restrict input_array_time_levels; CCTK_INT * restrict operand_indices; CCTK_INT * restrict operation_codes; + CCTK_INT * restrict time_deriv_order; CCTK_INT * restrict output_array_indices; - struct tensor const * restrict * restrict thetensor[3][4]; - int * restrict * restrict thevars[3][4]; + struct tensor const * restrict * restrict thetensor + [MAX_TIME_DERIV_ORDER+1][MAX_SPATIAL_DERIV_ORDER+1][MAX_TIME_LEVELS+1]; + int * restrict * restrict thevars + [MAX_TIME_DERIV_ORDER+1][MAX_SPATIAL_DERIV_ORDER+1][MAX_TIME_LEVELS+1]; int * restrict thebase; int * restrict thevar; int m; /* output array */ int n; /* point */ int i; /* var */ - int f; /* face */ int d; /* dimension */ - int q; /* number of derivatives */ + int p; /* number of time derivatives */ + int q; /* number of spatial derivatives */ int tl; /* time level */ int r; /* rank */ @@ -380,10 +427,10 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, int ierr; /* Check arguments */ - assert (N_dims == 3); + assert (N_dims == DIM); assert (N_interp_points >= 0); assert (interp_coords_type >= 0); - for (d=0; d<3; ++d) { + for (d=0; d= 0); @@ -406,8 +453,8 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Copy coordinates */ - for (d=0; d<3; ++d) { - new_interp_coords[d] = (CCTK_REAL*)malloc (N_interp_points * sizeof(CCTK_REAL)); + for (d=0; d= 0 + && operation_codes[m] < ipow (10, MAX_SPATIAL_DERIV_ORDER)); + } + + time_deriv_order = malloc (N_output_arrays * sizeof *time_deriv_order); + assert (time_deriv_order); + ierr = Util_TableGetIntArray + (param_table_handle, N_output_arrays, time_deriv_order, "time_deriv_order"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + for (m=0; m= 0 + && time_deriv_order[m] <= MAX_TIME_DERIV_ORDER); + } output_array_indices = malloc (N_output_arrays * sizeof *output_array_indices); @@ -517,15 +584,18 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Map Cactus variables to tensor objects */ - for (q=0; q<3; ++q) { - for (tl=0; tl<4; ++tl) { - thetensor[q][tl] = malloc (CCTK_NumVars() * sizeof *thetensor[q][tl]); - assert (thetensor[q][tl]); - thevars[q][tl] = malloc (CCTK_NumVars() * sizeof *thevars[q][tl]); - assert (thevars[q][tl]); - for (n=0; n=0 && varnvars); component = tensortype->comps[var]; - assert (tensortype->rank<10); + assert (tensortype->rank>=0 && tensortype->rank<=MAX_RANK); for (r=tensortype->rank-1; r>=0; --r) { indices[r] = component % tensortype->dim; assert (indices[r]>=0 && indices[r]dim); @@ -696,20 +767,17 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, oldrank = tensortype->rank; } - /* Take operation code (i.e., derivatives) into account */ - { - int code = operation_codes[m]; - num_derivs - = code == 0 ? 0 - : code < 10 ? 1 - : code < 100 ? 2 - : -1; - assert (num_derivs>=0); - } + /* Take time derivative order (i.e., time derivatives) into account */ + num_time_derivs = time_deriv_order[m]; + assert (num_time_derivs>=0 && num_time_derivs<=MAX_TIME_DERIV_ORDER); + + /* Take operation code (i.e., spatial derivatives) into account */ + num_derivs = ilog (operation_codes[m], 10); + assert (num_derivs>=0 && num_derivs<=MAX_SPATIAL_DERIV_ORDER); /* Get the time level */ time_level = input_array_time_levels[operand_indices[m]]; - assert (time_level>=0 && time_level<4); + assert (time_level>=0 && time_level<=MAX_TIME_LEVELS); if (tensortype == &scalar) { switch (num_derivs) { @@ -739,7 +807,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Find the additional indices */ { int code; - assert (tensortype->rank<10); + assert (tensortype->rank>=0 && tensortype->rank<=MAX_RANK); assert (tensortype->rank >= oldrank); code = operation_codes[m]; for (r=oldrank; rrank; ++r) { @@ -763,25 +831,28 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, } /* Create or cross-check the tensor object */ - if (! thetensor[num_derivs][time_level][basevar]) { - thetensor[num_derivs][time_level][basevar] = tensortype; - assert (! thevars[num_derivs][time_level][basevar]); - thevars[num_derivs][time_level][basevar] + if (! thetensor[num_time_derivs][num_derivs][time_level][basevar]) { + thetensor[num_time_derivs][num_derivs][time_level][basevar] = tensortype; + assert (! thevars[num_time_derivs][num_derivs][time_level][basevar]); + thevars[num_time_derivs][num_derivs][time_level][basevar] = malloc (tensortype->nvars - * sizeof *thevars[num_derivs][time_level][basevar]); - assert (thevars[num_derivs][time_level][basevar]); + * sizeof *thevars[num_time_derivs][num_derivs][time_level][basevar]); + assert (thevars[num_time_derivs][num_derivs][time_level][basevar]); for (i=0; invars; ++i) { - thevars[num_derivs][time_level][basevar][i] = -1; + thevars[num_time_derivs][num_derivs][time_level][basevar][i] = -1; } } - assert (thetensor[num_derivs][time_level][basevar] == tensortype); - assert (thevars[num_derivs][time_level][basevar]); + assert (thetensor[num_time_derivs][num_derivs][time_level][basevar] + == tensortype); + assert (thevars[num_time_derivs][num_derivs][time_level][basevar]); /* This does not hold if the caller requests the same interpolation to be done into different output arrays. This may happen e.g. when CarpetInterp needs to differentiate in time. This is arguably a performance bug in CarpetInterp. */ - /* assert (thevars[num_derivs][time_level][basevar][var] == -1); */ - thevars[num_derivs][time_level][basevar][var] = m; + /* See whether this goes away now. */ + assert (thevars[num_time_derivs][num_derivs][time_level][basevar][var] + == -1); + thevars[num_time_derivs][num_derivs][time_level][basevar][var] = m; } /* for m */ @@ -790,35 +861,33 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Loop over all output arrays */ for (m=0; m=0); - } + /* Take time derivative order (i.e., time derivatives) into account */ + num_time_derivs = time_deriv_order[m]; + assert (num_time_derivs>=0 && num_time_derivs<=MAX_TIME_DERIV_ORDER); + + /* Take operation code (i.e., spatial derivatives) into account */ + num_derivs = ilog (operation_codes[m], 10); + assert (num_derivs>=0 && num_derivs<=MAX_SPATIAL_DERIV_ORDER); /* Get the time level */ time_level = input_array_time_levels[operand_indices[m]]; - assert (time_level>=0 && time_level<4); + assert (time_level>=0 && time_level<=MAX_TIME_LEVELS); /* Get the tensor type */ basevar = thebase[m]; assert (basevar>=0 && basevar<=CCTK_NumVars()); - tensortype = thetensor[num_derivs][time_level][basevar]; + tensortype = thetensor[num_time_derivs][num_derivs][time_level][basevar]; assert (tensortype); - vars = thevars[num_derivs][time_level][basevar]; + vars = thevars[num_time_derivs][num_derivs][time_level][basevar]; assert (vars); var = thevar[m]; assert (var>=0 && varnvars); @@ -826,7 +895,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Transform into indices */ { int component = tensortype->comps[var]; - assert (tensortype->rank<10); + assert (tensortype->rank>=0 && tensortype->rank<=MAX_RANK); for (r=tensortype->rank-1; r>=0; --r) { indices[r] = component % tensortype->dim; assert (indices[r]>=0 && indices[r]dim); @@ -840,7 +909,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, CCTK_REAL x = ((CCTK_REAL const *)interp_coords[0])[n]; CCTK_REAL y = ((CCTK_REAL const *)interp_coords[1])[n]; - int myindices[10]; + int myindices[MAX_RANK+1]; int parity; for (r=0; rrank; ++r) { @@ -916,10 +985,12 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, free (operation_codes); free (output_array_indices); - for (q=0; q<3; ++q) { - for (tl=0; tl<4; ++tl) { - free ((void *) (thetensor[q][tl])); - free ((void *) (thevars[q][tl])); + for (p=0; p<=MAX_TIME_DERIV_ORDER; ++p) { + for (q=0; q<=MAX_SPATIAL_DERIV_ORDER; ++q) { + for (tl=0; tl<=MAX_TIME_LEVELS; ++tl) { + free ((void *) (thetensor[p][q][tl])); + free ((void *) (thevars[p][q][tl])); + } } } free (thebase); -- cgit v1.2.3