diff options
| -rw-r--r-- | configuration.ccl | 4 | ||||
| -rw-r--r-- | interface.ccl | 13 | ||||
| -rw-r--r-- | param.ccl | 28 | ||||
| -rw-r--r-- | schedule.ccl | 8 | ||||
| -rw-r--r-- | src/ApplyCartoon.c | 2 | ||||
| -rw-r--r-- | src/Cartoon2DBC.c | 217 | ||||
| -rw-r--r-- | src/SymInterp.c | 7 | ||||
| -rw-r--r-- | src/interpolate.c | 79 |
8 files changed, 342 insertions, 16 deletions
diff --git a/configuration.ccl b/configuration.ccl new file mode 100644 index 0000000..6582098 --- /dev/null +++ b/configuration.ccl @@ -0,0 +1,4 @@ +# Configuration definition for thorn Cartoon2D +# $Header$ + +REQUIRES SpaceMask diff --git a/interface.ccl b/interface.ccl index 6d97591..b432be3 100644 --- a/interface.ccl +++ b/interface.ccl @@ -6,6 +6,8 @@ implements: cartoon2d INCLUDES HEADER: Cartoon2D_tensors.h in Cartoon2D_tensors.h INCLUDES HEADER: Cartoon2D.h in Cartoon2D.h +USES INCLUDE: SpaceMask.h + # Function aliases: CCTK_INT FUNCTION SymmetryRegister (CCTK_STRING IN sym_name) @@ -67,3 +69,14 @@ CCTK_INT FUNCTION Boundary_SelectedGVs(CCTK_POINTER_TO_CONST IN GH, \ CCTK_INT ARRAY OUT faces, CCTK_INT ARRAY OUT boundary_widths, \ CCTK_INT ARRAY OUT table_handles, CCTK_STRING IN bc_name) REQUIRES FUNCTION Boundary_SelectedGVs + +private: + +int excision_variables type=SCALAR +{ + excision_active + old_mask_vi + new_mask_vi + new_excision_field + new_excision_descriptor +} "Internal variables to store information about excision" @@ -28,3 +28,31 @@ INT eno_order "The interpolation order applied to the ENO interpolator" BOOLEAN allow_grid_resize "Allow grid to be resized in a cartoon-compatible way" { } "no" + +BOOLEAN old_excision "Are we doing excision based on the old style mask?" +{ +} "no" + +BOOLEAN new_excision "Are we doing excision based on the new style mask?" +{ +} "no" + +STRING old_style_excision_var "The variable to be checked for old style excision" +{ + ".*" :: "Expected to be \'Spacemask::emask\'" +} "" + +STRING new_style_excision_var "The variable to be checked for new style excision" +{ + ".*" :: "Expected to be \'Spacemask::space_mask\'" +} "" + +STRING new_mask_field_name "The name of the field that describes excision for the new mask" +{ + ".*" :: "Could be anything" +} "" + +STRING new_mask_excised_name "The name of the descriptor that says the point is excised for the new mask" +{ + ".*" :: "Could be anything" +} "" diff --git a/schedule.ccl b/schedule.ccl index f0f19a4..55719b8 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -3,12 +3,20 @@ if (cartoon_active) { + STORAGE: excision_variables + schedule Cartoon2D_CheckTensorTypes at PARAMCHECK { LANG: C OPTIONS: meta } "Check tensor type definitions for consistency" + schedule Cartoon2D_InitExcisionVars at BASEGRID + { + LANG: C + OPTIONS: global + } "Initialize the excision variables" + if (allow_grid_resize) { schedule Cartoon2D_SetGrid at CCTK_RECOVER_PARAMETERS diff --git a/src/ApplyCartoon.c b/src/ApplyCartoon.c index 87680cd..b1417d6 100644 --- a/src/ApplyCartoon.c +++ b/src/ApplyCartoon.c @@ -178,7 +178,7 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) { if (CCTK_Equals(prolongtype, "None")) { - prolongmethod = PROLONG_LAGRANGE; /* But why? */ + prolongmethod = PROLONG_NONE; /* But why? */ } else if (CCTK_Equals(prolongtype, "Lagrange")) { diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c index f74c30f..7653424 100644 --- a/src/Cartoon2DBC.c +++ b/src/Cartoon2DBC.c @@ -24,6 +24,7 @@ #include "cctk_FortranString.h" #include "Cartoon2D_tensors.h" +#include "SpaceMask.h" static const char *rcsid = "$Id$"; @@ -38,9 +39,11 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, CCTK_REAL *f, int i, int di, int ijk, CCTK_REAL x); CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], CCTK_REAL x); + CCTK_REAL y[], CCTK_REAL x, + CCTK_INT excised_points[]); CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], CCTK_REAL x); + CCTK_REAL y[], CCTK_REAL x, + CCTK_INT excised_points[]); void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)(int *retval, const cGH **GH, int *tensortype, int *vi); void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)(int *retval, const cGH **GH, @@ -49,6 +52,9 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN) (int *retval, const cGH **GH, int *tensortype, ONE_FORTSTRING_ARG); +void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS); +void SetupExcisionVars(CCTK_ARGUMENTS); + /* set boundaries of a grid tensor assuming axisymmetry - handles lower boundary in x - does not touch other boundaries @@ -90,7 +96,6 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi) CCTK_REAL f[100], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz; CCTK_REAL *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz; - gi = CCTK_GroupIndexFromVarI (vi); assert (gi >= 0); n_vars = CCTK_NumVarsInGroupI (gi); @@ -342,6 +347,59 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f, int lnx; int n, offset; + CCTK_REAL *old_mask; + CCTK_INT *new_mask; + + CCTK_INT excised_points[6]; + + if (*excision_active < 0) + { + SetupExcisionVars((cGH*)cctkGH); + } + + if (old_excision) + { + old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi); + } + else + { + old_mask = NULL; + } + + if (new_excision) + { + new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi); + } + else + { + new_mask = NULL; + } + + /* + If this boundary point is excised then we return + without doing anything + */ + + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + if (old_mask[ijk] < 0.75) + { + return f[ijk]; + } + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + if (SpaceMask_CheckStateBits(new_mask, (ijk), + *new_excision_field, + *new_excision_descriptor)) + { + return f[ijk]; + } + } + } + lnx = cctkGH->cctk_lsh[0]; dx0 = CCTK_DELTA_SPACE(0); @@ -369,10 +427,25 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f, x0 = lx0 + dx0 * (i-offset); for (n = 0; n <= order; n++) { y[n] = f[ijk-offset+n]; + excised_points[n] = 0; + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 0 : 1; + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + excised_points[n] |= + (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n), + *new_excision_field, + *new_excision_descriptor)); + } + } } /* call interpolator */ - r = interpolate_local(order, x0, dx0, y, x); + r = interpolate_local(order, x0, dx0, y, x, excised_points); return r; } @@ -391,6 +464,59 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, int lnx; int n, offset; int tmp_order = 2; + + CCTK_REAL *old_mask; + CCTK_INT *new_mask; + + CCTK_INT excised_points[11]; + + if (*excision_active < 0) + { + SetupExcisionVars((cGH*)cctkGH); + } + + if (old_excision) + { + old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi); + } + else + { + old_mask = NULL; + } + + if (new_excision) + { + new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi); + } + else + { + new_mask = NULL; + } + + /* + If this boundary point is excised then we return + without doing anything + */ + + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + if (old_mask[ijk] < 0.75) + { + return f[ijk]; + } + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + if (SpaceMask_CheckStateBits(new_mask, (ijk), + *new_excision_field, + *new_excision_descriptor)) + { + return f[ijk]; + } + } + } lnx = cctkGH->cctk_lsh[0]; @@ -424,11 +550,92 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, for (n = 0; n < 2 * eno_order + 1; ++n) { y[n] = f[ijk-offset+n]; + excised_points[n] = 0; + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 1 : 0; + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + excised_points[n] |= + (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n), + *new_excision_field, + *new_excision_descriptor)); + } + } } /* call interpolator */ /* r = interpolate_eno(tmp_order, x0, dx0, y, x); */ - r = interpolate_eno(eno_order, x0, dx0, y, x); + r = interpolate_eno(eno_order, x0, dx0, y, x, excised_points); return r; } + +void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + *excision_active = -1; + *old_mask_vi = -1; + *new_mask_vi = -1; + *new_excision_field = -1; + *new_excision_descriptor = -1; + +} + +void SetupExcisionVars(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (old_excision) + { + *old_mask_vi = CCTK_VarIndex(old_style_excision_var); + } + else + { + *old_mask_vi = -1; + } + if (new_excision) + { + *new_mask_vi = CCTK_VarIndex(new_style_excision_var); + } + else + { + *new_mask_vi = -1; + } + + if (*old_mask_vi > -1) + { + if (*new_mask_vi > -1) + { + *excision_active = 3; + *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name); + *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name, + new_mask_excised_name); + } + else + { + *excision_active = 1; + } + } + else + { + if (*new_mask_vi > -1) + { + *excision_active = 2; + *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name); + *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name, + new_mask_excised_name); + } + else + { + *excision_active = 0; + } + } + +} diff --git a/src/SymInterp.c b/src/SymInterp.c index 0efc291..24668ab 100644 --- a/src/SymInterp.c +++ b/src/SymInterp.c @@ -360,6 +360,13 @@ 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; + } + } } diff --git a/src/interpolate.c b/src/interpolate.c index 1d774a2..e207cb2 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -23,11 +23,13 @@ CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x); + CCTK_REAL x, + CCTK_INT excised_points[]); CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x); + CCTK_REAL x, + CCTK_INT excised_points[]); /* prototypes for private functions defined in this file */ static CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, @@ -75,14 +77,70 @@ static CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x) + CCTK_REAL x, + CCTK_INT excised_points[]) { - switch (order) { - case 1: return interpolate_local_order1(x0, dx, y, x); - case 2: return interpolate_local_order2(x0, dx, y, x); - case 3: return interpolate_local_order3(x0, dx, y, x); - case 4: return interpolate_local_order4(x0, dx, y, x); - case 5: return interpolate_local_order5(x0, dx, y, x); + + /* + We assume that the excised region has one of the following 4 forms: + + o o o o o (No excision) + x o o o o (left excision) + o o o o x (right excision) + x o o o x (left and right excision) + + Anything else (e.g., x o x o x) will go horribly wrong. + */ + + int start, end, available_order; + + start = 0; + end = order; + + while ((start < order + 1)&&(excised_points[start])) + { + ++start; + } + while ((end > -1)&&(excised_points[end])) + { + --end; + } + + if (end < start) + { + /* The whole block is excised */ + return y[0]; + } + + /* The possible interpolation order that can be used */ + + available_order = end - start; + + if (available_order == 0) + { + /* There is only one non-excised point */ + return y[start]; + } + + if (dx * start + x0 > x) + { + /* We would have to extrapolate to the left */ + return y[start]; + } + + if (dx * end + x0 < x) + { + /* We would have to extrapolate to the right */ + return y[end]; + } + + switch (available_order) + { + case 1: return interpolate_local_order1(x0 + start * dx, dx, &y[start], x); + case 2: return interpolate_local_order2(x0 + start * dx, dx, &y[start], x); + case 3: return interpolate_local_order3(x0 + start * dx, dx, &y[start], x); + case 4: return interpolate_local_order4(x0 + start * dx, dx, &y[start], x); + case 5: return interpolate_local_order5(x0 + start * dx, dx, &y[start], x); default: CCTK_WARN(0, "Interpolation order not supported"); } } @@ -218,7 +276,8 @@ return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))) ); CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x) + CCTK_REAL x, + CCTK_INT excised_points[]) { /* CCTK_REAL diff[4]; */ |