From 54620b9361553333aeec9f43a20ab87d5ac1dcbf Mon Sep 17 00:00:00 2001 From: schnetter Date: Thu, 19 Oct 2006 00:27:03 +0000 Subject: Add tensor types to handle the Weyl scalars. Check that the boundary specification is consistent between this thorn and CoordBase, when CoordBase is used. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/ReflectionSymmetry/trunk@17 082bdb00-0f4f-0410-b49e-b1835e5f2039 --- interface.ccl | 10 ++++ src/apply.c | 155 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- src/interpolate.c | 34 ++++++++++-- 3 files changed, 190 insertions(+), 9 deletions(-) diff --git a/interface.ccl b/interface.ccl index ae0af7d..d1c90c4 100644 --- a/interface.ccl +++ b/interface.ccl @@ -68,3 +68,13 @@ CCTK_INT FUNCTION Boundary_SelectedGVs \ CCTK_INT ARRAY OUT table_handles, \ CCTK_STRING IN bc_name) REQUIRES FUNCTION Boundary_SelectedGVs + + + +CCTK_INT FUNCTION GetBoundarySpecification \ + (CCTK_INT IN size, \ + CCTK_INT OUT ARRAY nboundaryzones, \ + CCTK_INT OUT ARRAY is_internal, \ + CCTK_INT OUT ARRAY is_staggered, \ + CCTK_INT OUT ARRAY shiftout) +USES FUNCTION GetBoundarySpecification diff --git a/src/apply.c b/src/apply.c index 6170234..50854d0 100644 --- a/src/apply.c +++ b/src/apply.c @@ -146,7 +146,7 @@ BndReflectVI (cGH const * restrict const cctkGH, int table; char tensortypealias[1000]; - enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR }; + enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL }; enum tensortype ttype; CCTK_INT tensorparity; int tcomponent; @@ -263,9 +263,8 @@ BndReflectVI (cGH const * restrict const cctkGH, /* scalar */ ttype = SCALAR; tcomponent = 0; - } - else if (CCTK_EQUALS (tensortypealias, "u") - || CCTK_EQUALS (tensortypealias, "d")) + } else if (CCTK_EQUALS (tensortypealias, "u") + || CCTK_EQUALS (tensortypealias, "d")) { /* vector */ assert (numvars == 3); @@ -277,6 +276,13 @@ BndReflectVI (cGH const * restrict const cctkGH, assert (numvars == 6); ttype = TENSOR; tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) { + /* Weyl scalars, stored as 10 real values. NOTE: This assumes + that Psi_0 comes first, which is NOT the default with + PsiKadelia. */ + assert (numvars == 10); + ttype = WEYLSCALARS_REAL; + tcomponent = vi - firstvar; } else { char * groupname = CCTK_GroupName(gi); assert (groupname); @@ -297,6 +303,9 @@ BndReflectVI (cGH const * restrict const cctkGH, case TENSOR: assert (tcomponent>=0 && tcomponent<6); break; + case WEYLSCALARS_REAL: + assert (tcomponent>=0 && tcomponent<10); + break; default: assert (0); } @@ -366,6 +375,22 @@ BndReflectVI (cGH const * restrict const cctkGH, default: assert (0); } break; + case WEYLSCALARS_REAL: + { + static int const weylparities[10][3] = + {{+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}}; + parity *= weylparities[tcomponent][dir]; + } + break; default: assert (0); } @@ -563,6 +588,126 @@ BndReflectVI (cGH const * restrict const cctkGH, +/* When CoordBase is used to specify the location of the boundary + points, then ensure that the CoordBase parameters and this thorn's + parameters are consistent. */ +static +void +CheckBoundaryParameters (cGH const * restrict const cctkGH, + int const vi, + int const * restrict const stencil) +{ + DECLARE_CCTK_PARAMETERS; + + static int did_check = 0; + + int type; /* Parameter type */ + void const * ptr; /* Pointer to parameter value */ + char const * coordtype; /* CartGrid3D::type */ + + int dim; /* Number of dimensions of vi */ + + CCTK_INT * restrict nboundaryzones; /* CoordBase boundary location */ + CCTK_INT * restrict is_internal; + CCTK_INT * restrict is_staggered; + CCTK_INT * restrict shiftout; + + int do_reflection[6]; /* This thorn's parameters */ + int do_stagger[6]; + + int d; + int ierr; + + + + /* Check only once to save time */ + if (did_check) return; + + /* Check only for grid functions */ + if (CCTK_GroupTypeFromVarI (vi) != CCTK_GF) return; + + did_check = 1; + + /* Check whether CartGrid3D is active */ + if (! CCTK_IsThornActive ("CartGrid3D")) return; + + /* Check whether CoordBase is used */ + ptr = CCTK_ParameterGet ("type", "CartGrid3D", & type); + assert (ptr != 0); + assert (type == PARAMETER_KEYWORD); + coordtype = * (char const * *) ptr; + if (! CCTK_EQUALS (coordtype, "coordbase")) return; + + /* Get the boundary specification */ + dim = CCTK_GroupDimFromVarI (vi); + assert (dim >= 0); + nboundaryzones = malloc (2*dim * sizeof *nboundaryzones); + is_internal = malloc (2*dim * sizeof *is_internal); + is_staggered = malloc (2*dim * sizeof *is_staggered); + shiftout = malloc (2*dim * sizeof *shiftout); + ierr = GetBoundarySpecification + (2*dim, nboundaryzones, is_internal, is_staggered, shiftout); + assert (! ierr); + + /* Reflection symmetry information */ + assert (dim == 3); + do_reflection[0] = reflection_x; + do_reflection[1] = reflection_upper_x; + do_reflection[2] = reflection_y; + do_reflection[3] = reflection_upper_y; + do_reflection[4] = reflection_z; + do_reflection[5] = reflection_upper_z; + + do_stagger[0] = avoid_origin_x; + do_stagger[1] = avoid_origin_upper_x; + do_stagger[2] = avoid_origin_y; + do_stagger[3] = avoid_origin_upper_y; + do_stagger[4] = avoid_origin_z; + do_stagger[5] = avoid_origin_upper_z; + + /* Check the boundary sizes */ + for (d=0; d<6; ++d) { + if (do_reflection[d]) { + char const dir = "xyz"[d/2]; + char const * const face = (d%2==0) ? "lower" : "upper"; + if (stencil[d/2] != nboundaryzones[d]) { + CCTK_VWarn (CCTK_WARN_ABORT, + __LINE__, __FILE__, CCTK_THORNSTRING, + "The %s %c face is a symmetry boundary. Since there are %d ghost zones in the %c direction, the corresponding CoordBase boundary width must also be %d. The boundary width is currently %d.", + face, dir, + stencil[d/2], dir, stencil[d/2], + (int) nboundaryzones[d]); + } + if (is_internal[d]) { + CCTK_VWarn (CCTK_WARN_ABORT, + __LINE__, __FILE__, CCTK_THORNSTRING, + "The %s %c face is a symmetry boundary. The corresponding CoordBase boundary must not be internal.", + face, dir); + } + if (do_stagger[d] != is_staggered[d]) { + CCTK_VWarn (CCTK_WARN_ABORT, + __LINE__, __FILE__, CCTK_THORNSTRING, + "The %s %c face is a symmetry boundary. The symmetry condition and the corresponding CoordBase boundary must either be both staggered or both not staggered.", + face, dir); + } + if ((do_stagger[d] ? 0 : 1) != shiftout[d]) { + CCTK_VWarn (CCTK_WARN_ABORT, + __LINE__, __FILE__, CCTK_THORNSTRING, + "The %s %c face is a symmetry boundary. If the symmetry condition is staggered, then the corresponding CoordBase shiftout must be 0; otherwise it must be 1.", + face, dir); + } + } + } + + /* Free memory */ + free (nboundaryzones); + free (is_internal); + free (is_staggered); + free (shiftout); +} + + + void ReflectionSymmetry_Apply (CCTK_ARGUMENTS) { @@ -655,6 +800,8 @@ ReflectionSymmetry_Apply (CCTK_ARGUMENTS) CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI"); } + CheckBoundaryParameters (cctkGH, vi, stencil); + ierr = BndReflectVI (cctkGH, stencil, vi); if (ierr) { diff --git a/src/interpolate.c b/src/interpolate.c index 884abe6..7c8831e 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -178,7 +178,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, int table; char tensortypealias[1000]; - enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR }; + enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL }; enum tensortype ttype; CCTK_INT tensorparity; int tcomponent; @@ -224,8 +224,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, ttype = UNKNOWN; tcomponent = 0; - if (CCTK_EQUALS (tensortypealias, "scalar")) - { + if (CCTK_EQUALS (tensortypealias, "scalar")) { /* scalar */ if (numvars != 1) { groupname = CCTK_GroupName(gi); @@ -239,8 +238,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, tcomponent = 0; } else if (CCTK_EQUALS (tensortypealias, "u") - || CCTK_EQUALS (tensortypealias, "d")) - { + || CCTK_EQUALS (tensortypealias, "d")) { /* vector */ assert (numvars == 3); ttype = VECTOR; @@ -251,6 +249,11 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, assert (numvars == 6); ttype = TENSOR; tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) { + /* Weyl scalars, stored as 10 real numbers */ + assert (numvars == 10); + ttype = WEYLSCALARS_REAL; + tcomponent = vi - firstvar; } else { char * groupname = CCTK_GroupName(gi); assert (groupname); @@ -271,6 +274,9 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, case TENSOR: assert (tcomponent>=0 && tcomponent<6); break; + case WEYLSCALARS_REAL: + assert (tcomponent>=0 && tcomponent<10); + break; default: assert (0); } @@ -311,6 +317,24 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, default: assert (0); } break; + case WEYLSCALARS_REAL: + { + static int const weylparities[10][3] = + {{+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}}; + for (dir=0; dir<3; ++dir) { + parities[dir] = weylparities[tcomponent][dir]; + } + } + break; default: assert (0); } -- cgit v1.2.3