From 448b2720f5fe8dfb4bf4f3277c14d03e76adab8e Mon Sep 17 00:00:00 2001 From: schnetter Date: Fri, 28 May 2004 15:59:40 +0000 Subject: Provide symmetry interpolation. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@84 eec4d7dc-71c2-46d6-addf-10296150bf52 --- interface.ccl | 59 ++++- schedule.ccl | 10 +- src/Cartoon2D.h | 20 ++ src/RegisterSym.c | 10 + src/SymInterp.c | 758 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 2 +- 6 files changed, 848 insertions(+), 11 deletions(-) create mode 100644 src/SymInterp.c diff --git a/interface.ccl b/interface.ccl index dc60934..2047b2d 100644 --- a/interface.ccl +++ b/interface.ccl @@ -9,18 +9,59 @@ INCLUDES HEADER: Cartoon2D.h in Cartoon2D.h # Function aliases: CCTK_INT FUNCTION SymmetryRegister (CCTK_STRING IN sym_name) -USES FUNCTION SymmetryRegister +REQUIRES FUNCTION SymmetryRegister + +CCTK_INT FUNCTION \ + SymmetryRegisterGrid \ + (CCTK_POINTER IN cctkGH, \ + CCTK_INT IN sym_handle, \ + CCTK_INT IN ARRAY which_faces, \ # array [N_FACES] + CCTK_INT IN ARRAY symmetry_zone_width) # array [N_FACES] +REQUIRES FUNCTION SymmetryRegisterGrid + +CCTK_INT FUNCTION \ + SymmetryRegisterGridInterpolator \ + (CCTK_POINTER IN cctkGH, \ + CCTK_INT IN sym_handle, \ + CCTK_INT CCTK_FPOINTER IN symmetry_interpolate \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN N_dims, \ + CCTK_INT IN local_interp_handle, \ + CCTK_INT IN param_table_handle, \ + CCTK_INT IN coord_system_handle, \ + CCTK_INT IN N_interp_points, \ + CCTK_INT IN interp_coords_type, \ + CCTK_POINTER_TO_CONST ARRAY IN interp_coords, \ + CCTK_INT IN N_input_arrays, \ + CCTK_INT ARRAY IN input_array_indices, \ + CCTK_INT IN N_output_arrays, \ + CCTK_INT ARRAY IN output_array_types, \ + CCTK_POINTER ARRAY IN output_arrays, \ + CCTK_INT IN faces)) +REQUIRES FUNCTION SymmetryRegisterGridInterpolator + +CCTK_INT FUNCTION \ + SymmetryInterpolateFaces \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN N_dims, \ + CCTK_INT IN local_interp_handle, \ + CCTK_INT IN param_table_handle, \ + CCTK_INT IN coord_system_handle, \ + CCTK_INT IN N_interp_points, \ + CCTK_INT IN interp_coords_type, \ + CCTK_POINTER_TO_CONST ARRAY IN interp_coords, \ + CCTK_INT IN N_input_arrays, \ + CCTK_INT ARRAY IN input_array_indices, \ + CCTK_INT IN N_output_arrays, \ + CCTK_INT ARRAY IN output_array_types, \ + CCTK_POINTER ARRAY IN output_arrays, \ + CCTK_INT IN faces) +REQUIRES FUNCTION SymmetryInterpolateFaces + -CCTK_INT FUNCTION \ - SymmetryRegisterGrid \ - (CCTK_POINTER IN cctkGH, \ - CCTK_INT IN sym_handle, \ - CCTK_INT IN ARRAY which_faces, \ - CCTK_INT IN ARRAY symmetry_zone_width) -USES FUNCTION SymmetryRegisterGrid CCTK_INT FUNCTION Boundary_SelectedGVs(CCTK_POINTER_TO_CONST IN GH, \ CCTK_INT IN array_size, CCTK_INT ARRAY OUT var_indicies, \ CCTK_INT ARRAY OUT faces, CCTK_INT ARRAY OUT boundary_widths, \ CCTK_INT ARRAY OUT table_handles, CCTK_STRING IN bc_name) -USES FUNCTION Boundary_SelectedGVs +REQUIRES FUNCTION Boundary_SelectedGVs diff --git a/schedule.ccl b/schedule.ccl index 2902d30..f0f19a4 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -3,6 +3,12 @@ if (cartoon_active) { + schedule Cartoon2D_CheckTensorTypes at PARAMCHECK + { + LANG: C + OPTIONS: meta + } "Check tensor type definitions for consistency" + if (allow_grid_resize) { schedule Cartoon2D_SetGrid at CCTK_RECOVER_PARAMETERS @@ -14,16 +20,18 @@ if (cartoon_active) schedule Cartoon2D_RegisterSymmetries in SymmetryRegister { LANG: C + OPTIONS: global } "Register symmetry boundaries" schedule Cartoon2D_CheckParameters at CCTK_PARAMCHECK { LANG: C + OPTIONS: meta } "Check Cartoon2D parameters" # Apply cartoon boundary conditions after physical boundaries: schedule Cartoon_ApplyBoundaries in BoundaryConditions after Boundary_ApplyPhysicalBCs { LANG: C - } "Execute Cartoon boundary conditions" + } "Apply Cartoon boundary conditions" } diff --git a/src/Cartoon2D.h b/src/Cartoon2D.h index 476e985..b4fadc3 100644 --- a/src/Cartoon2D.h +++ b/src/Cartoon2D.h @@ -1,3 +1,7 @@ +/* $Header$ */ + +#include "cctk.h" + /* Function prototypes */ #ifdef __cplusplus extern "C" { @@ -5,6 +9,22 @@ extern "C" { int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi); +CCTK_INT +Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_INT const N_dims, + CCTK_INT const local_interp_handle, + CCTK_INT const param_table_handle, + CCTK_INT const coord_system_handle, + CCTK_INT const N_interp_points, + CCTK_INT const interp_coords_type, + CCTK_POINTER_TO_CONST const interp_coords[], + CCTK_INT const N_input_arrays, + CCTK_INT const input_array_indices[], + CCTK_INT const N_output_arrays, + CCTK_INT const output_array_types[], + CCTK_POINTER const output_arrays[], + CCTK_INT const faces); + #ifdef __cplusplus } #endif diff --git a/src/RegisterSym.c b/src/RegisterSym.c index 0574672..a5af8ab 100644 --- a/src/RegisterSym.c +++ b/src/RegisterSym.c @@ -7,10 +7,14 @@ @enddesc @@*/ +#include + #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" +#include "Cartoon2D.h" + static const char * const rcsid = "$Header$"; CCTK_FILEVERSION(BetaThorns_Cartoon2D_SetSym_c); @@ -48,4 +52,10 @@ void Cartoon2D_RegisterSymmetries (CCTK_ARGUMENTS) CCTK_WARN (0, "Could not register the Cartoon boundaries -- probably some other thorn has already registered the same boundary faces for a different symmetry"); } + ierr = SymmetryRegisterGridInterpolator + (cctkGH, handle, Cartoon2D_SymmetryInterpolate); + if (ierr < 0) { + CCTK_WARN (0, "Could not register the symmetry interpolator"); + } + } diff --git a/src/SymInterp.c b/src/SymInterp.c new file mode 100644 index 0000000..8904bcb --- /dev/null +++ b/src/SymInterp.c @@ -0,0 +1,758 @@ +/* $Header$ */ + +#include +#include +#include +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "Cartoon2D.h" + + + +/* + * NOTE: + * + * The contents of this file have been largely copied from the thorn + * AEIDevelopment/RotatingSymmetry90. Please do not invest much + * effort in changing this file here -- changes to its ancestor should + * be easy to port over. The only difference between this file and + * its ancestor are a few lines in the coordinate rotation routines. + * I hope that a generic tensor infrastructure will come up, and then + * this duplication will go away. + * + * Erik Schnetter, 2004-05-22 + */ + + + +/* A tensor description */ +struct tensor +{ + const char * name; /* description */ + int dim; + int rank; + int ncomps; /* dim^rank */ + int * restrict vars; /* map component to variable */ + int * restrict parity; /* parity for the above */ + int nvars; /* depends on symmetries */ + int * restrict comps; /* map variable to component */ +}; + + + +/* A scalar */ +static int scalar_vars[] = { + 0 +}; +static int scalar_parity[] = { + +1 +}; +static int scalar_comps[] = { + 0 +}; +static struct tensor scalar = { + "scalar", + 3, 0, 1, scalar_vars, scalar_parity, 1, scalar_comps +}; + + + +/* A vector */ +static int vector_vars[] = { + 0,1,2 +}; +static int vector_parity[] = { + +1,+1,+1 +}; +static int vector_comps[] = { + 0,1,2 +}; +static struct tensor vector = { + "vector", + 3, 1, 3, vector_vars, vector_parity, 3, vector_comps +}; + + + +/* A second rank tensor without symmetries */ +static int tensor_vars[] = { + 0,1,2, 3,4,5, 6,7,8 +}; +static int tensor_parity[] = { + +1,+1,+1, +1,+1,+1, +1,+1,+1 +}; +static int tensor_comps[] = { + 0,1,2, 3,4,5, 6,7,8 +}; +static struct tensor tensor = { + "tensor", + 3, 2, 9, tensor_vars, tensor_parity, 9, tensor_comps +}; + + + +/* A symmetric second rank tensor */ +static int symmtensor_vars[] = { + 0,1,2, 1,3,4, 2,4,5 +}; +static int symmtensor_parity[] = { + +1,+1,+1, +1,+1,+1, +1,+1,+1 +}; +static int symmtensor_comps[] = { + 0,1,2, 4,5, 8 +}; +static struct tensor symmtensor = { + "symmetric tensor T_(ij)", + 3, 2, 9, symmtensor_vars, symmtensor_parity, 6, symmtensor_comps +}; + + + +/* A third rank tensor with symmetry in the first two indices */ +static int symmtensor3a_vars[] = { + 0, 1, 2, + 3, 4, 5, + 6, 7, 8, + + 3, 4, 5, + 9,10,11, + 12,13,14, + + 6, 7, 8, + 12,13,14, + 15,16,17 +}; +static int symmtensor3a_parity[] = { + +1,+1,+1, + +1,+1,+1, + +1,+1,+1, + + +1,+1,+1, + +1,+1,+1, + +1,+1,+1, + + +1,+1,+1, + +1,+1,+1, + +1,+1,+1, +}; +static int symmtensor3a_comps[] = { + 0, 1, 2, + 3, 4, 5, + 6, 7, 8, + + 12,13,14, + 15,16,17, + + 24,25,26 +}; +static struct tensor symmtensor3a = { + "symmetric tensor T_(ij)k", + 3, 3, 27, symmtensor3a_vars, symmtensor3a_parity, 18, symmtensor3a_comps +}; + + + +/* A third rank tensor with symmetry in the last two indices */ +static int symmtensor3b_vars[] = { + 0, 1, 2, 1, 3, 4, 2, 4, 5, + 6, 7, 8, 7, 9,10, 8,10,11, + 12,13,14, 13,15,16, 14,16,17 +}; +static int symmtensor3b_parity[] = { + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1 +}; +static int symmtensor3b_comps[] = { + 0, 1, 2, 4, 5, 8, + 9,10,11, 13,14, 17, + 18,19,20, 22,23, 26 +}; +static struct tensor symmtensor3b = { + "symmetric tensor T_i(jk)", + 3, 3, 27, symmtensor3b_vars, symmtensor3b_parity, 18, symmtensor3b_comps +}; + + + +/* A fourth rank tensor with symmetries both in its first and last two + indices */ +static int symmtensor4_vars[] = { + 0, 1, 2, 1, 3, 4, 2, 4, 5, + 6, 7, 8, 7, 9,10, 8,10,11, + 12,13,14, 13,15,16, 14,16,17, + + 6, 7, 8, 7, 9,10, 8,10,11, + 18,19,20, 19,21,22, 20,22,23, + 24,25,26, 25,27,28, 26,28,29, + + 12,13,14, 13,15,16, 14,16,17, + 24,25,26, 25,27,28, 26,28,29, + 30,31,32, 31,33,34, 32,34,35 +}; +static int symmtensor4_parity[] = { + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1, + + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1, + + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1, + +1,+1,+1, +1,+1,+1, +1,+1,+1 +}; +static int symmtensor4_comps[] = { + 0, 1, 2, 4, 5, 8, + 9,10,11, 13,14, 17, + 18,19,20, 22,23, 26, + + 36,37,38, 40,41, 44, + 45,46,47, 49,50, 53, + + 72,73,74, 76,77, 80 +}; +static struct tensor symmtensor4 = { + "symmetric tensor T_(ij)(kl)", + 3, 4, 81, symmtensor4_vars, symmtensor4_parity, 36, symmtensor4_comps +}; + + + +/* Ensure that all tensor declarations are internally consistent */ +static void +CheckTensorType (struct tensor const * restrict const atensor) +{ + int i, n; + assert (atensor->name); + 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->nvars>=0 && atensor->nvars<=atensor->ncomps); + assert (atensor->vars); + for (i=0; incomps; ++i) { + assert (atensor->vars[i]>=0 && atensor->vars[i]nvars); + assert (abs(atensor->parity[i]) <= 1); + } + assert (atensor->comps); + for (n=0; nnvars; ++n) { + assert (atensor->comps[n]>=0 && atensor->comps[n]ncomps); + assert (atensor->vars[atensor->comps[n]] == n); + } +} + +void +Cartoon2D_CheckTensorTypes (CCTK_ARGUMENTS) +{ + CheckTensorType (&scalar); + CheckTensorType (&vector); + CheckTensorType (&tensor); + CheckTensorType (&symmtensor); + CheckTensorType (&symmtensor3a); + CheckTensorType (&symmtensor3b); + CheckTensorType (&symmtensor4); +} + + + +/* Symmetry interpolation */ +CCTK_INT +Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_INT const N_dims, + CCTK_INT const local_interp_handle, + CCTK_INT const param_table_handle, + CCTK_INT const coord_system_handle, + CCTK_INT const N_interp_points, + CCTK_INT const interp_coords_type, + CCTK_POINTER_TO_CONST const interp_coords[], + CCTK_INT const N_input_arrays, + CCTK_INT const input_array_indices[], + CCTK_INT const N_output_arrays, + CCTK_INT const output_array_types[], + CCTK_POINTER const output_arrays[], + CCTK_INT const faces) +{ + cGH const * restrict const cctkGH = cctkGH_; + + CCTK_POINTER new_interp_coords[3]; + CCTK_POINTER * restrict new_output_arrays; + CCTK_INT new_faces; + CCTK_INT * restrict operand_indices; + CCTK_INT * restrict operation_codes; + CCTK_INT * restrict output_array_indices; + + struct tensor const * restrict * restrict thetensor[3]; + int * restrict * restrict thevars[3]; + 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 r; /* rank */ + + int iret; /* interpolator return value */ + int ierr; + + /* Check arguments */ + assert (N_dims == 3); + assert (N_interp_points >= 0); + assert (interp_coords_type >= 0); + for (d=0; d<3; ++d) { + assert (N_interp_points == 0 || interp_coords[d]); + } + assert (N_output_arrays >= 0); + + /* Coordinates must be CCTK_REAL */ + assert (interp_coords_type == CCTK_VARIABLE_REAL); + for (m=0; m=0 && operand_indices[m]=0 + && output_array_indices[m]=0 && vi=0 && gi0); + firstvar = CCTK_FirstVarIndexI(gi); + assert (firstvar>=0); + table = CCTK_GroupTagsTableI(gi); + assert (table>=0); + + /* Get the tensor type alias */ + 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 (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Tensor type alias not declared for group \"%s\" -- assuming a scalar", + groupname); + free (groupname); + strcpy (tensortypealias, "scalar"); + } 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); + } + + /* Find the tensor type */ + tensortype = NULL; + basevar = -1; + var = -1; + if (CCTK_EQUALS (tensortypealias, "scalar")) { + /* scalar */ + if (numvars != 1) { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element", + groupname); + free (groupname); + } + tensortype = &scalar; + basevar = vi; + var = 0; + } else if (CCTK_EQUALS (tensortypealias, "u")) { + /* vector */ + assert (numvars == 3); + tensortype = &vector; + basevar = firstvar; + var = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) { + /* symmetric tensor */ + assert (numvars == 6); + tensortype = &symmtensor; + basevar = firstvar; + var = vi - firstvar; + } 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); + } + thebase[m] = basevar; + + /* Find my component */ + { + int component; + assert (var>=0 && varnvars); + component = tensortype->comps[var]; + assert (tensortype->rank<10); + for (r=tensortype->rank-1; r>=0; --r) { + indices[r] = component % tensortype->dim; + assert (indices[r]>=0 && indices[r]dim); + component /= tensortype->dim; + } + assert (component == 0); + 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); + } + + if (tensortype == &scalar) { + switch (num_derivs) { + case 0: break; + case 1: tensortype = &vector; break; + case 2: tensortype = &symmtensor; break; + default: assert (0); + } + } else if (tensortype == &vector) { + switch (num_derivs) { + case 0: break; + case 1: tensortype = &tensor; break; + case 2: tensortype = &symmtensor3b; break; + default: assert (0); + } + } else if (tensortype == &symmtensor) { + switch (num_derivs) { + case 0: break; + case 1: tensortype = &symmtensor3a; break; + case 2: tensortype = &symmtensor4; break; + default: assert (0); + } + } else { + assert (0); + } + + /* Find the additional indices */ + { + int code; + assert (tensortype->rank<10); + assert (tensortype->rank >= oldrank); + code = operation_codes[m]; + for (r=oldrank; rrank; ++r) { + const int thedir = code % 10 - 1; + code /= 10; + assert (thedir>=0 && thedirdim); + indices[r] = thedir; + } + } + + /* Re-calculate component */ + { + int component = 0; + for (r=0; rrank; ++r) { + component = component * tensortype->dim + indices[r]; + } + assert (component>=0 && componentncomps); + var = tensortype->vars[component]; + assert (var>=0 && varnvars); + thevar[m] = var; + } + + /* 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]); + for (i=0; invars; ++i) { + thevars[num_derivs][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; + + } /* for m */ + + + + /* Loop over all output arrays */ + for (m=0; m=0); + } + + /* Get the tensor type */ + basevar = thebase[m]; + assert (basevar>=0 && basevar<=CCTK_NumVars()); + tensortype = thetensor[num_derivs][basevar]; + assert (tensortype); + vars = thevars[num_derivs][basevar]; + assert (vars); + var = thevar[m]; + assert (var>=0 && varnvars); + + /* Transform into indices */ + { + int component = tensortype->comps[var]; + assert (tensortype->rank<10); + for (r=tensortype->rank-1; r>=0; --r) { + indices[r] = component % tensortype->dim; + assert (indices[r]>=0 && indices[r]dim); + component /= tensortype->dim; + } + } + + /* Loop over all grid points */ + for (n=0; nncomps; ++thecomponent) { + int theindices[10]; + int thevar; + int theparity; + int mm; + assert (tensortype->rank<10); + { + int r; + int component = thecomponent; + for (r=tensortype->rank-1; r>=0; --r) { + theindices[r] = component % tensortype->dim; + assert (theindices[r]>=0 && theindices[r]dim); + component /= tensortype->dim; + } + assert (component == 0); + } + theparity = tensortype->parity[thecomponent]; + assert (abs(theparity) <= 1); + thevar = tensortype->vars[thecomponent]; + assert (thevar>=0 && thevarnvars); + mm = vars[thevar]; + assert (mm>=0 && mmrank; ++r) { + val *= rot[indices[r]][theindices[r]]; + } + myval += val; + } + } + + ((CCTK_REAL *)output_arrays[m])[n] = myval; + + } /* for n */ + + } /* for m */ + + + + /* Free output variable descriptions */ + for (m=0; m