diff options
| -rw-r--r-- | src/SymInterp.c | 371 |
1 files changed, 302 insertions, 69 deletions
diff --git a/src/SymInterp.c b/src/SymInterp.c index 24668ab..44a3ddd 100644 --- a/src/SymInterp.c +++ b/src/SymInterp.c @@ -31,6 +31,22 @@ + +/* 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 { @@ -225,6 +241,35 @@ 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 @@ -252,6 +297,9 @@ CheckTensorType (struct tensor const * restrict const atensor) void Cartoon2D_CheckTensorTypes (CCTK_ARGUMENTS) { + int gi; + + /* Check all internal tensor type definitions */ CheckTensorType (&scalar); CheckTensorType (&vector); CheckTensorType (&tensor); @@ -259,6 +307,92 @@ Cartoon2D_CheckTensorTypes (CCTK_ARGUMENTS) CheckTensorType (&symmtensor3a); CheckTensorType (&symmtensor3b); CheckTensorType (&symmtensor4); + + /* Check tensor types of all groups */ + for (gi=0; gi<CCTK_NumGroups(); ++gi) { + + char tensortypealias[100]; + int numvars, firstvar; + int table; + int ierr; + + numvars = CCTK_NumVarsInGroupI(gi); + if (numvars == 0) continue; + assert (numvars>0); + firstvar = CCTK_FirstVarIndexI(gi); + assert (firstvar>=0); + table = CCTK_GroupTagsTableI(gi); + assert (table>=0); + + ierr = Util_TableGetString + (table, sizeof tensortypealias, tensortypealias, "tensortypealias"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "Tensor type alias not declared for group \"%s\" -- assuming a scalar", + groupname); + free (groupname); + strcpy (tensortypealias, ""); + } else if (ierr<0) { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in tensor type alias declaration for group \"%s\"", + groupname); + free (groupname); + } + + if (CCTK_EQUALS (tensortypealias, "")) { + /* do nothing */ + } else if (CCTK_EQUALS (tensortypealias, "scalar")) { + /* scalar */ + if (numvars != 1) { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING, + "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element", + groupname); + free (groupname); + } + } else if (CCTK_EQUALS (tensortypealias, "4scalar")) { + /* 4-scalar */ + if (numvars != 1) { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING, + "Group \"%s\" has the tensor type alias \"4scalar\", but contains more than 1 element", + groupname); + free (groupname); + } + } else if (CCTK_EQUALS (tensortypealias, "u") + || CCTK_EQUALS (tensortypealias, "d")) + { + /* vector */ + assert (numvars == DIM); + } else if (CCTK_EQUALS (tensortypealias, "4u") + || CCTK_EQUALS (tensortypealias, "4d")) + { + /* 4-vector */ + assert (numvars == DIM+1); + } else if (CCTK_EQUALS (tensortypealias, "uu_sym") + || CCTK_EQUALS (tensortypealias, "dd_sym")) { + /* symmetric tensor */ + assert (numvars == DIM*(DIM+1)/2); + } else if (CCTK_EQUALS (tensortypealias, "4uu_sym") + || CCTK_EQUALS (tensortypealias, "4dd_sym")) { + /* symmetric 4-tensor */ + assert (numvars == (DIM+1)*(DIM+2)/2); + } else { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Illegal tensor type alias for group \"%s\"", + groupname); + free (groupname); + } + + } } @@ -282,15 +416,19 @@ Cartoon2D_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]; - int * restrict * restrict thevars[3]; + 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; @@ -299,17 +437,19 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, 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 */ int iret; /* interpolator return value */ 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<DIM; ++d) { assert (N_interp_points == 0 || interp_coords[d]); } assert (N_output_arrays >= 0); @@ -334,7 +474,7 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Copy coordinates */ - for (d=0; d<3; ++d) { + for (d=0; d<DIM; ++d) { new_interp_coords[d] = malloc (N_interp_points * sizeof(CCTK_REAL)); assert (new_interp_coords[d]); } @@ -360,13 +500,6 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, for (m=0; m<N_output_arrays; ++m) { new_output_arrays[m] = malloc (N_interp_points * sizeof(CCTK_REAL)); assert (new_output_arrays[m]); - { - CCTK_INT mm; - for (mm = 0; mm < N_interp_points; ++mm) - { - ((CCTK_REAL*)new_output_arrays[m])[mm] = 0; - } - } } @@ -383,13 +516,27 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Free coordinates */ - for (d=0; d<3; ++d) { + for (d=0; d<DIM; ++d) { free (new_interp_coords[d]); } /* Find output variable indices */ + input_array_time_levels + = malloc (N_input_arrays * sizeof *input_array_time_levels); + assert (input_array_time_levels); + ierr = Util_TableGetIntArray + (param_table_handle, N_input_arrays, + input_array_time_levels, "input_array_time_levels"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + for (m=0; m<N_input_arrays; ++m) { + input_array_time_levels[m] = 0; /* time level is 0 */ + } + } else { + assert (ierr == N_input_arrays); + } + operand_indices = malloc (N_output_arrays * sizeof *operand_indices); assert (operand_indices); ierr = Util_TableGetIntArray @@ -408,14 +555,34 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, ierr = Util_TableGetIntArray (param_table_handle, N_output_arrays, operation_codes, "operation_codes"); if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - assert (N_output_arrays == N_input_arrays); + assert (N_output_arrays == N_input_arrays); /* why? */ for (m=0; m<N_output_arrays; ++m) { operation_codes[m] = 0; /* do not take derivatives */ } } else { assert (ierr == N_output_arrays); } + for (m=0; m<N_output_arrays; ++m) { + assert (operation_codes[m] >= 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<N_output_arrays; ++m) { + time_deriv_order[m] = 0; /* do not take time derivatives */ + } + } else { + assert (ierr == N_output_arrays); + } + for (m=0; m<N_output_arrays; ++m) { + assert (time_deriv_order[m] >= 0 + && time_deriv_order[m] <= MAX_TIME_DERIV_ORDER); + } + output_array_indices = malloc (N_output_arrays * sizeof *output_array_indices); assert (output_array_indices); @@ -429,14 +596,19 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Map Cactus variables to tensor objects */ - for (q=0; q<3; ++q) { - thetensor[q] = malloc (CCTK_NumVars() * sizeof *thetensor[q]); - assert (thetensor[q]); - thevars[q] = malloc (CCTK_NumVars() * sizeof *thevars[q]); - assert (thevars[q]); - for (n=0; n<CCTK_NumVars(); ++n) { - thetensor[q][n] = NULL; - thevars[q][n] = NULL; + 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) { + thetensor[p][q][tl] + = malloc (CCTK_NumVars() * sizeof *thetensor[p][q][tl]); + assert (thetensor[p][q][tl]); + thevars[p][q][tl] = malloc (CCTK_NumVars() * sizeof *thevars[p][q][tl]); + assert (thevars[p][q][tl]); + for (n=0; n<CCTK_NumVars(); ++n) { + thetensor[p][q][tl][n] = NULL; + thevars[p][q][tl][n] = NULL; + } + } } } @@ -458,9 +630,11 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, struct tensor const * tensortype; int basevar; int var; - int indices[10]; + int indices[MAX_RANK+1]; int oldrank; + int num_time_derivs; int num_derivs; + int time_level; /* Get some variable information */ vi = output_array_indices[m]; @@ -513,16 +687,57 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, var = 0; } else if ((CCTK_EQUALS (tensortypealias, "u")) ||(CCTK_EQUALS (tensortypealias, "d"))) { /* vector */ - assert (numvars == 3); + assert (numvars == DIM); tensortype = &vector; basevar = firstvar; var = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "4u") + || CCTK_EQUALS (tensortypealias, "4d")) + { + /* 4-vector */ + assert (numvars == DIM+1); + if (vi == firstvar) { + /* temporal component */ + int const off = 0; + tensortype = &scalar; + basevar = firstvar + off; + var = vi - basevar; + } else { + /* spatial components */ + int const off = 1; + tensortype = &vector; + basevar = firstvar + off; + var = vi - basevar; + } } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) { /* symmetric tensor */ - assert (numvars == 6); + assert (numvars == DIM*(DIM+1)/2); tensortype = &symmtensor; basevar = firstvar; var = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "4uu_sym") + || CCTK_EQUALS (tensortypealias, "4dd_sym")) { + /* symmetric 4-tensor */ + assert (numvars == (DIM+1)*(DIM+2)/2); + if (vi == firstvar) { + /* temporal-temporal component */ + int const off = 0; + tensortype = &scalar; + basevar = firstvar + off; + var = vi - basevar; + } else if (vi < firstvar+DIM+1) { + /* temporal-spatial components */ + int const off = 1; + tensortype = &vector; + basevar = firstvar + off; + var = vi - basevar; + } else { + /* spatial-spatial components */ + int const off = DIM+1; + tensortype = &symmtensor; + basevar = firstvar + off; + var = vi - basevar; + } } else { char * groupname = CCTK_GroupName(gi); assert (groupname); @@ -538,7 +753,7 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, int component; assert (var>=0 && var<tensortype->nvars); 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]<tensortype->dim); @@ -547,17 +762,18 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, assert (component == 0); oldrank = tensortype->rank; } + + /* 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. 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 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<=MAX_TIME_LEVELS); if (tensortype == &scalar) { switch (num_derivs) { @@ -587,7 +803,7 @@ Cartoon2D_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; r<tensortype->rank; ++r) { @@ -611,20 +827,28 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, } /* Create or cross-check the tensor object */ - if (! thetensor[num_derivs][basevar]) { - thetensor[num_derivs][basevar] = tensortype; - assert (! thevars[num_derivs][basevar]); - thevars[num_derivs][basevar] - = malloc (tensortype->nvars * sizeof *thevars[num_derivs][basevar]); - assert (thevars[num_derivs][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_time_derivs][num_derivs][time_level][basevar]); + assert (thevars[num_time_derivs][num_derivs][time_level][basevar]); for (i=0; i<tensortype->nvars; ++i) { - thevars[num_derivs][basevar][i] = -1; + thevars[num_time_derivs][num_derivs][time_level][basevar][i] = -1; } } - assert (thetensor[num_derivs][basevar] == tensortype); - assert (thevars[num_derivs][basevar]); - assert (thevars[num_derivs][basevar][var] == -1); - thevars[num_derivs][basevar][var] = m; + 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. */ + /* 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 */ @@ -633,30 +857,33 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, /* Loop over all output arrays */ for (m=0; m<N_output_arrays; ++m) { + int num_time_derivs; int num_derivs; + int time_level; int basevar; struct tensor const * restrict tensortype; int const * restrict vars; int var; - int indices[10]; + int indices[MAX_RANK+1]; - /* 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<=MAX_TIME_LEVELS); + /* Get the tensor type */ basevar = thebase[m]; assert (basevar>=0 && basevar<=CCTK_NumVars()); - tensortype = thetensor[num_derivs][basevar]; + tensortype = thetensor[num_time_derivs][num_derivs][time_level][basevar]; assert (tensortype); - vars = thevars[num_derivs][basevar]; + vars = thevars[num_time_derivs][num_derivs][time_level][basevar]; assert (vars); var = thevar[m]; assert (var>=0 && var<tensortype->nvars); @@ -664,7 +891,7 @@ Cartoon2D_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]<tensortype->dim); @@ -702,11 +929,11 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, int thecomponent; for (thecomponent=0; thecomponent<tensortype->ncomps; ++thecomponent) { - int theindices[10]; + int theindices[MAX_RANK+1]; int thevar; int theparity; int mm; - assert (tensortype->rank<10); + assert (tensortype->rank<MAX_RANK+1); { int r; int component = thecomponent; @@ -748,13 +975,19 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, } free (new_output_arrays); + free (input_array_time_levels); free (operand_indices); free (operation_codes); free (output_array_indices); - for (q=0; q<3; ++q) { - free ((void *) (thetensor[q])); - free ((void *) (thevars[q])); + + 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); free (thevar); |