diff options
Diffstat (limited to 'src/apply.c')
| -rw-r--r-- | src/apply.c | 155 |
1 files changed, 151 insertions, 4 deletions
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) { |