From 77424202a473eae863478cbd7c901a0c494d652d Mon Sep 17 00:00:00 2001 From: hawke Date: Thu, 10 Nov 2005 14:05:35 +0000 Subject: Zeroth attempt at doing excision with Cartoon. Now requires the Spacemask thorn from CactusEinstein. Initial tests suggest perfect second order convergence but not quite perfect fourth order convergence, at least at reasonable resolutions. Convergence is still better than 3rd order, though. Excision with the new spacemask has been added but not tested. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@97 eec4d7dc-71c2-46d6-addf-10296150bf52 --- src/Cartoon2DBC.c | 217 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 212 insertions(+), 5 deletions(-) (limited to 'src/Cartoon2DBC.c') 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; + } + } + +} -- cgit v1.2.3