From 52a24adeead81157f5d205c6a94a3d4a8731c453 Mon Sep 17 00:00:00 2001 From: schnetter Date: Wed, 5 Jul 2006 19:05:46 +0000 Subject: Introduce a function SBP_determine_onesided_stencil which determines which faces should use one-sided stencils, depending on which boundaries are inter-processor boundaries, symmetry boundaries, and multi-patch boundaries. Use this function everywhere. Remove the previous mechinisms; some were not in all cases correct. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@75 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- src/GetScalProdDiag.F90 | 23 ++--------- src/Get_Coeff.F90 | 6 +++ src/call_derivs.c | 55 ++++++-------------------- src/call_derivs_name.c | 20 ++++++---- src/call_up_derivs.c | 61 ++++++++--------------------- src/dissipation.c | 28 ++++---------- src/get_coeffs.c | 101 ++++++++++++++++++++++++++---------------------- src/make.code.defn | 3 +- src/set_norm_mask.F90 | 18 ++++----- src/stencil.c | 66 +++++++++++++++++++++++++++++++ src/stencil.h | 8 ++++ 11 files changed, 196 insertions(+), 193 deletions(-) create mode 100644 src/stencil.c create mode 100644 src/stencil.h diff --git a/src/GetScalProdDiag.F90 b/src/GetScalProdDiag.F90 index e2a7e99..00e3dff 100644 --- a/src/GetScalProdDiag.F90 +++ b/src/GetScalProdDiag.F90 @@ -19,12 +19,7 @@ subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad ) CCTK_REAL, parameter :: zero = 0.0 integer, parameter :: wp = kind(zero) - CCTK_INT :: symtable, n_elements, nchar, pen_sym_handle, np - CCTK_INT, dimension(6) :: symbnd - CCTK_POINTER :: psym_name - character(len=256) :: symmetry_name - integer :: status - integer, dimension(6) :: bbox + integer, dimension(6) :: onesided CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 0.5_wp, 1.0_wp /) CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 17.0_wp/48.0_wp, & @@ -50,21 +45,11 @@ subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad ) call CCTK_WARN(0, 'Error: dir is outside the legal range') end if - call CCTK_GroupbboxGN ( status, cctkGH, 6, bbox, & - 'SummationByParts::normmask' ) - if ( status < 0 ) then - call CCTK_WARN(0,'Error: unable to get bounding box information') - end if - - symtable = SymmetryTableHandleForGrid ( cctkGH ) - call Util_TableGetIntArray ( n_elements, symtable, 6, & - symbnd, 'symmetry_handle' ) + call SBP_determine_onesided_stencil (cctkGH, onesided) - pen_sym_handle = SymmetryHandleOfName ( 'multipatch' ) - sigmad = 1.0_wp - if ( symbnd(dir+1) == pen_sym_handle .and. bbox(dir*2+1) == 1 ) then + if ( onesided(dir*2+1) == 1 ) then select case (order) case (2) sigmad(1:2) = bmask_2 @@ -76,7 +61,7 @@ subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad ) sigmad(1:8) = bmask_8 end select end if - if ( symbnd(dir+2) == pen_sym_handle .and. bbox(dir*2+2) == 1 ) then + if ( onesided(dir*2+2) == 1 ) then select case (order) case (2) sigmad(nsize:nsize-1:-1) = bmask_2 diff --git a/src/Get_Coeff.F90 b/src/Get_Coeff.F90 index 3c7a09f..60dab5b 100644 --- a/src/Get_Coeff.F90 +++ b/src/Get_Coeff.F90 @@ -22,6 +22,8 @@ CCTK_REAL function GetCoeff () GetCoeff = 13649.0_wp/43200.0_wp case (8) GetCoeff = 1498139.0_wp/5080320.0_wp + case default + call CCTK_WARN (0, "operator not implemented") end select else if ( CCTK_EQUALS(operator_type,'Minimal Bandwidth') ) then @@ -30,6 +32,8 @@ CCTK_REAL function GetCoeff () GetCoeff = 3.0_wp/11.0_wp case(6) GetCoeff = 30.0_wp/137.0_wp + case default + call CCTK_WARN (0, "operator not implemented") end select else select case (order) @@ -37,6 +41,8 @@ CCTK_REAL function GetCoeff () GetCoeff = 0.2388575707774486064323157210922003533466_wp case(6) GetCoeff = 0.2028105550720356346665604029847379994496_wp + case default + call CCTK_WARN (0, "operator not implemented") end select end if end if diff --git a/src/call_derivs.c b/src/call_derivs.c index 5dd4530..22d6527 100644 --- a/src/call_derivs.c +++ b/src/call_derivs.c @@ -5,28 +5,7 @@ #include "util_ErrorCodes.h" #include "util_Table.h" - - -/* Determine whether a boundary with the symmetry handle symbnd is a - "regular" symmetry boundary (where the SBP stencils should not be - modified), or an outer boundary (or a multi-block boundary), where - the SBP stencils need to be modified. */ -static int is_regular_symbnd (int const symbnd, int const pen_sym_handle) -{ - /* On an outer boundary (which is not a symmetry boundary), symbnd < - 0. */ - if (symbnd < 0) return 0; - - if (pen_sym_handle >= 0) { - /* If the symmetry boundary is a multi-block boundary, then symbnd - = pen_sym_handle. However, we can only check that if thorn - MultiPatch is active, i.e., if pen_sym_handle >= 0. */ - if (symbnd == pen_sym_handle) return 0; - } - - /* This is a "regular" symmetry boundary. */ - return 1; -} +#include "stencil.h" @@ -38,12 +17,10 @@ void DiffGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS - CCTK_INT ni, nj, nk, gsize, ic, loc_order,i; + CCTK_INT ni, nj, nk, gsize, loc_order; CCTK_REAL delta; - CCTK_INT ierr; - CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3]; - CCTK_INT symtable, pen_sym_handle; - CCTK_INT symbnd[6]; + CCTK_INT bb[2]; + int onesided[6]; int nelements; void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(const CCTK_REAL *var, @@ -151,37 +128,27 @@ void DiffGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, loc_order = order; } - symtable = SymmetryTableHandleForGrid (cctkGH); - if (symtable<0) { - CCTK_WARN(0,"symtable is out of bounds"); - } - - ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle"); - if (ierr!=6) { - CCTK_WARN(0,"Util_TableGetIntArray returned error"); - } - - pen_sym_handle = SymmetryHandleOfName ( "multipatch" ); + SBP_determine_onesided_stencil (cctkGH, onesided); switch(dir) { case 0: { delta = CCTK_DELTA_SPACE(0); - bb[0] = cctk_bbox[0] && ! is_regular_symbnd (symbnd[0], pen_sym_handle); - bb[1] = cctk_bbox[1] && ! is_regular_symbnd (symbnd[1], pen_sym_handle); + bb[0] = onesided[0]; + bb[1] = onesided[1]; gsize = cctk_nghostzones[0]; break; } case 1: { delta = CCTK_DELTA_SPACE(1); - bb[0] = cctk_bbox[2] && ! is_regular_symbnd (symbnd[2], pen_sym_handle); - bb[1] = cctk_bbox[3] && ! is_regular_symbnd (symbnd[3], pen_sym_handle); + bb[0] = onesided[2]; + bb[1] = onesided[3]; gsize = cctk_nghostzones[1]; break; } case 2: { delta = CCTK_DELTA_SPACE(2); - bb[0] = cctk_bbox[4] && ! is_regular_symbnd (symbnd[4], pen_sym_handle); - bb[1] = cctk_bbox[5] && ! is_regular_symbnd (symbnd[5], pen_sym_handle); + bb[0] = onesided[4]; + bb[1] = onesided[5]; gsize = cctk_nghostzones[2]; break; } diff --git a/src/call_derivs_name.c b/src/call_derivs_name.c index f7804b0..89e4b39 100644 --- a/src/call_derivs_name.c +++ b/src/call_derivs_name.c @@ -1,8 +1,12 @@ +#include + #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" -#include +#include "stencil.h" + + void DiffGf ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, const char *var_name, const char *dvar_name ) @@ -12,10 +16,10 @@ void DiffGf ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, DECLARE_CCTK_ARGUMENTS CCTK_REAL *var, *dvar; - CCTK_INT ni, nj, nk, gsize, ic; + CCTK_INT ni, nj, nk, gsize; CCTK_REAL delta; - CCTK_INT ierr; - CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3]; + CCTK_INT bb[2]; + int onesided[6]; void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(const CCTK_REAL *var, const CCTK_INT *ni, const CCTK_INT *nj, @@ -65,22 +69,24 @@ void DiffGf ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2]; + SBP_determine_onesided_stencil (cctkGH, onesided); + switch(dir) { case 0: { delta = CCTK_DELTA_SPACE(0); - bb[0] = cctk_bbox[0]; bb[1] = cctk_bbox[1]; + bb[0] = onesided[0]; bb[1] = onesided[1]; gsize = cctk_nghostzones[0]; break; } case 1: { delta = CCTK_DELTA_SPACE(1); - bb[0] = cctk_bbox[2]; bb[1] = cctk_bbox[3]; + bb[0] = onesided[2]; bb[1] = onesided[3]; gsize = cctk_nghostzones[1]; break; } case 2: { delta = CCTK_DELTA_SPACE(2); - bb[0] = cctk_bbox[4]; bb[1] = cctk_bbox[5]; + bb[0] = onesided[4]; bb[1] = onesided[5]; gsize = cctk_nghostzones[2]; break; } diff --git a/src/call_up_derivs.c b/src/call_up_derivs.c index e6f5afa..10607a1 100644 --- a/src/call_up_derivs.c +++ b/src/call_up_derivs.c @@ -1,3 +1,5 @@ +#include + #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" @@ -5,6 +7,10 @@ #include "util_ErrorCodes.h" #include "util_Table.h" +#include "stencil.h" + + + void DiffUpGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, const CCTK_REAL *var, CCTK_REAL *dvar, const CCTK_REAL *up, const CCTK_INT table_handle ) @@ -16,9 +22,8 @@ void DiffUpGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_INT ni, nj, nk, gsize, ic, loc_order; CCTK_REAL delta; CCTK_INT ierr; - CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3]; - CCTK_INT symtable, pen_sym_handle; - CCTK_INT symbnd[6]; + CCTK_INT lsh[3], bb[2], nghostzones[3]; + int onesided[6]; int nelements; void CCTK_FCALL CCTK_FNAME(up_deriv_gf_2_1)(const CCTK_REAL *var, @@ -76,61 +81,27 @@ void DiffUpGv ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, loc_order = order; } - symtable = SymmetryTableHandleForGrid (cctkGH); - if (symtable<0) { - CCTK_WARN(0,"symtable is out of bounds"); - } - - ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle"); - if (ierr!=6) { - CCTK_WARN(0,"Util_TableGetIntArray returned error"); - } - - pen_sym_handle = SymmetryHandleOfName ( "multipatch" ); - + SBP_determine_onesided_stencil (cctkGH, onesided); + switch(dir) { case 0: { delta = CCTK_DELTA_SPACE(0); - if ( symbnd[0] == pen_sym_handle && cctk_bbox[0] == 1 ) { - bb[0] = 1; - } else { - bb[0] = 0; - } - if ( symbnd[1] == pen_sym_handle && cctk_bbox[1] == 1 ) { - bb[1] = 1; - } else { - bb[1] = 0; - } + bb[0] = onesided[0]; + bb[1] = onesided[1]; gsize = cctk_nghostzones[0]; break; } case 1: { delta = CCTK_DELTA_SPACE(1); - if ( symbnd[2] == pen_sym_handle && cctk_bbox[2] == 1) { - bb[0] = 1; - } else { - bb[0] = 0; - } - if ( symbnd[3] == pen_sym_handle && cctk_bbox[3] == 1) { - bb[1] = 1; - } else { - bb[1] = 0; - } + bb[0] = onesided[2]; + bb[1] = onesided[3]; gsize = cctk_nghostzones[1]; break; } case 2: { delta = CCTK_DELTA_SPACE(2); - if ( symbnd[4] == pen_sym_handle && cctk_bbox[4] == 1) { - bb[0] = 1; - } else { - bb[0] = 0; - } - if ( symbnd[4] == pen_sym_handle && cctk_bbox[5] == 1) { - bb[1] = 1; - } else { - bb[1] = 0; - } + bb[0] = onesided[4]; + bb[1] = onesided[5]; gsize = cctk_nghostzones[2]; break; } diff --git a/src/dissipation.c b/src/dissipation.c index 8f85d32..bf6fc29 100644 --- a/src/dissipation.c +++ b/src/dissipation.c @@ -11,6 +11,8 @@ #include "util_ErrorCodes.h" #include "util_Table.h" +#include "stencil.h" + void CCTK_FCALL CCTK_FNAME(dissipation_4_3_opt) (const CCTK_REAL *var, const CCTK_INT *lsh, const CCTK_INT *gsh, @@ -155,13 +157,12 @@ apply (int const varindex, char const * const optstring, void * const arg) cGroup vardata, rhsdata; CCTK_REAL const * varptr; CCTK_REAL * rhsptr; - CCTK_INT ni, nj, nk, bbox[6], gsize[3]; + CCTK_INT gsize[3]; CCTK_REAL dx[3]; - int n; int d; int ierr; - CCTK_INT symtable, pen_sym_handle; - CCTK_INT symbnd[6]; + CCTK_INT bbox[6]; + int onesided[6]; CCTK_INT npatches, patch; assert (varindex >= 0); @@ -171,24 +172,9 @@ apply (int const varindex, char const * const optstring, void * const arg) gsize[d] = cctk_nghostzones[d]; } - symtable = SymmetryTableHandleForGrid (cctkGH); - if (symtable<0) { - CCTK_WARN(0,"symtable is out of bounds"); - } - - ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle"); - if (ierr!=6) { - CCTK_WARN(0,"Util_TableGetIntArray returned error"); - } - - pen_sym_handle = SymmetryHandleOfName ( "multipatch" ); - + SBP_determine_onesided_stencil (cctkGH, onesided); for (d=0; d<6; ++d) { - if ( symbnd[d] == pen_sym_handle && cctk_bbox[d] == 1 ) { - bbox[d] = 1; - } else { - bbox[d] = 0; - } + bbox[d] = onesided[d]; } rhsindex = MoLQueryEvolvedRHS (varindex); diff --git a/src/get_coeffs.c b/src/get_coeffs.c index 3d441a3..c2aa224 100644 --- a/src/get_coeffs.c +++ b/src/get_coeffs.c @@ -8,6 +8,55 @@ #include "util_ErrorCodes.h" #include "util_Table.h" +#include "stencil.h" + + + +void CCTK_FCALL CCTK_FNAME(set_coeff_2_1)(const CCTK_INT *nsize, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); +void CCTK_FCALL CCTK_FNAME(set_coeff_4_2)(const CCTK_INT *nsize, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); +void CCTK_FCALL CCTK_FNAME(set_coeff_6_3)(const CCTK_INT *nsize, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); +void CCTK_FCALL CCTK_FNAME(set_coeff_8_4)(const CCTK_INT *nsize, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); +void CCTK_FCALL CCTK_FNAME(set_coeff_4_3)(const CCTK_INT *nsize, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); +void CCTK_FCALL CCTK_FNAME(set_coeff_6_5)(const CCTK_INT *nsize, + const CCTK_INT *loc_order, + const CCTK_INT *bb, + const CCTK_INT *gsize, + CCTK_INT *imin, + CCTK_INT *imax, + CCTK_REAL *q); + + + void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, const CCTK_INT nsize, CCTK_INT *imin, CCTK_INT *imax, CCTK_REAL *q, const CCTK_INT table_handle ) @@ -19,49 +68,7 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, CCTK_INT ni, nj, nk, gsize, loc_order; int nelements; CCTK_INT bb[2]; - void CCTK_FCALL CCTK_FNAME(set_coeff_2_1)(const CCTK_INT *nsize, - const CCTK_INT *loc_order, - const CCTK_INT *bb, - const CCTK_INT *gsize, - CCTK_INT *imin, - CCTK_INT *imax, - CCTK_REAL *q); - void CCTK_FCALL CCTK_FNAME(set_coeff_4_2)(const CCTK_INT *nsize, - const CCTK_INT *loc_order, - const CCTK_INT *bb, - const CCTK_INT *gsize, - CCTK_INT *imin, - CCTK_INT *imax, - CCTK_REAL *q); - void CCTK_FCALL CCTK_FNAME(set_coeff_6_3)(const CCTK_INT *nsize, - const CCTK_INT *loc_order, - const CCTK_INT *bb, - const CCTK_INT *gsize, - CCTK_INT *imin, - CCTK_INT *imax, - CCTK_REAL *q); - void CCTK_FCALL CCTK_FNAME(set_coeff_8_4)(const CCTK_INT *nsize, - const CCTK_INT *loc_order, - const CCTK_INT *bb, - const CCTK_INT *gsize, - CCTK_INT *imin, - CCTK_INT *imax, - CCTK_REAL *q); - void CCTK_FCALL CCTK_FNAME(set_coeff_4_3)(const CCTK_INT *nsize, - const CCTK_INT *loc_order, - const CCTK_INT *bb, - const CCTK_INT *gsize, - CCTK_INT *imin, - CCTK_INT *imax, - CCTK_REAL *q); - void CCTK_FCALL CCTK_FNAME(set_coeff_6_5)(const CCTK_INT *nsize, - const CCTK_INT *loc_order, - const CCTK_INT *bb, - const CCTK_INT *gsize, - CCTK_INT *imin, - CCTK_INT *imax, - CCTK_REAL *q); - + int onesided[6]; ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2]; if ( table_handle >=0 ) { @@ -75,22 +82,24 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir, loc_order = order; } + SBP_determine_onesided_stencil (cctkGH, onesided); + switch(dir) { case 0: { assert(nsize==ni); - bb[0] = cctk_bbox[0]; bb[1] = cctk_bbox[1]; + bb[0] = onesided[0]; bb[1] = onesided[1]; gsize = cctk_nghostzones[0]; break; } case 1: { assert(nsize==nj); - bb[0] = cctk_bbox[2]; bb[1] = cctk_bbox[3]; + bb[0] = onesided[2]; bb[1] = onesided[3]; gsize = cctk_nghostzones[1]; break; } case 2: { assert(nsize==nk); - bb[0] = cctk_bbox[4]; bb[1] = cctk_bbox[5]; + bb[0] = onesided[4]; bb[1] = onesided[5]; gsize = cctk_nghostzones[2]; break; } diff --git a/src/make.code.defn b/src/make.code.defn index dd2255b..b9672ed 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -39,7 +39,8 @@ SRCS = call_derivs.c \ Coefficients_6_3.F90 \ Coefficients_8_4.F90 \ Coefficients_4_3.F90 \ - Coefficients_6_5.F90 + Coefficients_6_5.F90 \ + stencil.c # Subdirectories containing source files SUBDIRS = diff --git a/src/set_norm_mask.F90 b/src/set_norm_mask.F90 index c9b8ffb..40bae3e 100644 --- a/src/set_norm_mask.F90 +++ b/src/set_norm_mask.F90 @@ -17,8 +17,6 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS) integer, parameter :: wp = kind(zero) CCTK_INT :: symtable, n_elements, nchar, pen_sym_handle, np CCTK_INT, dimension(6) :: symbnd - CCTK_POINTER :: psym_name - character(len=256) :: symmetry_name CCTK_REAL, dimension(:), allocatable :: mask_x, mask_y, mask_z ! Note: The first number is twice the value from the paper, since Carpet @@ -55,7 +53,7 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS) pen_sym_handle = SymmetryHandleOfName ( 'multipatch' ) - if ( any ( symbnd == pen_sym_handle ) ) then + if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. any ( symbnd == pen_sym_handle )) ) then allocate ( mask_x(ni), mask_y(nj), mask_z(nk) ) mask_x = 1.0d0; mask_y = 1.0d0; mask_z = 1.0d0 @@ -79,22 +77,22 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS) end select end if - if ( symbnd(1) == pen_sym_handle ) then - mask_x(1:np) = bmask(1:np) + if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(1) == pen_sym_handle) ) then + mask_x(1:np) = bmask(1:np) end if - if ( symbnd(2) == pen_sym_handle ) then + if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(2) == pen_sym_handle) ) then mask_x(ni:ni-np+1:-1) = bmask(1:np) end if - if ( symbnd(3) == pen_sym_handle ) then + if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(3) == pen_sym_handle) ) then mask_y(1:np) = bmask(1:np) end if - if ( symbnd(4) == pen_sym_handle ) then + if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(4) == pen_sym_handle) ) then mask_y(nj:nj-np+1:-1) = bmask(1:np) end if - if ( symbnd(5) == pen_sym_handle ) then + if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(5) == pen_sym_handle) ) then mask_z(1:np) = bmask(1:np) end if - if ( symbnd(6) == pen_sym_handle ) then + if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(6) == pen_sym_handle) ) then mask_z(nk:nk-np+1:-1) = bmask(1:np) end if diff --git a/src/stencil.c b/src/stencil.c new file mode 100644 index 0000000..007a7c6 --- /dev/null +++ b/src/stencil.c @@ -0,0 +1,66 @@ +#include "cctk.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "stencil.h" + + + +/* Determine whether a boundary with the symmetry handle symbnd is a + "regular" symmetry boundary (where the SBP stencils should not be + modified), or an outer boundary (or a multi-block boundary), where + the SBP stencils need to be modified. */ + +void SBP_determine_onesided_stencil (const cGH * cctkGH, int * onesided) +{ + int symtable; + int pen_sym_handle; + CCTK_INT symbnd[6]; + int ierr; + int d; + + symtable = SymmetryTableHandleForGrid (cctkGH); + if (symtable<0) { + CCTK_WARN(0,"Cannot get symmetry table handle -- maybe thorn SymBase is not active?"); + } + + ierr = Util_TableGetIntArray (symtable, 6, symbnd, "symmetry_handle"); + if (ierr!=6) { + CCTK_WARN(0,"Cannot get symmetry handles"); + } + + pen_sym_handle = SymmetryHandleOfName ( "multipatch" ); + + for (d=0; d<6; ++d) { + if (! cctkGH->cctk_bbox[d]) { + /* This is an inter-processor boundary */ + onesided[d] = 0; + } else { + /* On an outer boundary (which is not a symmetry boundary), + symbnd < 0. */ + if (symbnd[d] < 0) { + onesided[d] = 1; + } else { + if (pen_sym_handle >= 0) { + /* If the symmetry boundary is a multi-block boundary, then + symbnd = pen_sym_handle. However, we can only check that + thorn MultiPatch is active, i.e., whether pen_sym_handle + >= 0. */ + if (symbnd[d] == pen_sym_handle) { + onesided[d] = 1; + } else { + onesided[d] = 0; + } + } else { + onesided[d] = 0; + } + } + } + } +} + +CCTK_FCALL void CCTK_FNAME (SBP_determine_onesided_stencil) (CCTK_POINTER_TO_CONST cctkGH, int * onesided) +{ + SBP_determine_onesided_stencil (cctkGH, onesided); +} diff --git a/src/stencil.h b/src/stencil.h new file mode 100644 index 0000000..7937e1e --- /dev/null +++ b/src/stencil.h @@ -0,0 +1,8 @@ +#ifndef STENCIL_H +#define STENCIL_H + +#include "cctk.h" + +void SBP_determine_onesided_stencil (const cGH * cctkGH, int * onesided); + +#endif /* #ifndef STENCIL_H */ -- cgit v1.2.3