diff options
Diffstat (limited to 'src/bc_noise.c')
| -rw-r--r-- | src/bc_noise.c | 270 |
1 files changed, 270 insertions, 0 deletions
diff --git a/src/bc_noise.c b/src/bc_noise.c new file mode 100644 index 0000000..0336cd8 --- /dev/null +++ b/src/bc_noise.c @@ -0,0 +1,270 @@ + +#include "noise.h" + +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_FortranString.h" + +#include "Symmetry.h" + + + +static int ApplyBndNoise (const cGH *GH, + int stencil_dir, + const int *stencil_alldirs, + int dir, + int first_var, + int num_vars); + +int +BndNoiseGI (const cGH *GH, const int *stencil, int gi) +{ + int first_vi, retval; + + first_vi = CCTK_FirstVarIndexI (gi); + if (first_vi >= 0) + { + retval = ApplyBndNoise (GH, -1, stencil, 0, first_vi, + CCTK_NumVarsInGroupI (gi)); + } + else + { + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "Invalid group index %d in BndNoiseGI", gi); + retval = -1; + } + + return (retval); +} + +void +CCTK_FCALL CCTK_FNAME (BndNoiseGI) (int *ierr, const cGH *GH, + const int *stencil, const int *gi) +{ + *ierr = BndNoiseGI (GH, stencil, *gi); +} + + +int +BndNoiseGN (const cGH *GH, const int *stencil, const char *gn) +{ + int gi, retval; + + gi = CCTK_GroupIndex (gn); + if (gi >= 0) + { + retval = BndNoiseGI (GH, stencil, gi); + } + else + { + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "Invalid group name '%s' in BndNoiseGN", gn); + retval = -1; + } + + return (retval); +} + +void CCTK_FCALL CCTK_FNAME (BndNoiseGN) + (int *ierr, + const cGH *GH, + const int *stencil, + ONE_FORTSTRING_ARG) +{ + ONE_FORTSTRING_CREATE (gn) + *ierr = BndNoiseGN (GH, stencil, gn); + free (gn); +} + + + + + +static int ApplyBndNoise (const cGH *GH, + int stencil_dir, + const int *stencil_alldirs, + int dir, + int first_var, + int num_vars) +{ + int i, j, k; + int var, vtypesize, gindex, gdim, timelvl; + int doBC[2*MAXDIM], dstag[MAXDIM], lsh[MAXDIM], lssh[MAXDIM], + stencil[MAXDIM]; + SymmetryGHex *sGHex; + CCTK_REAL amplitude; + CCTK_REAL offset; + int type; + + stencil[0] = 1; stencil[1] = 1; stencil[2] = 1; + + /* get the group index of the variables */ + gindex = CCTK_GroupIndexFromVarI (first_var); + + amplitude = *((CCTK_REAL*) CCTK_ParameterGet("amplitude", "Noise", &type)) + / RAND_MAX; + offset = 0.5 * amplitude; + + /* get the number of dimensions and the size of the variables' type */ + gdim = CCTK_GroupDimI (gindex); + vtypesize = CCTK_VarTypeSize (CCTK_VarTypeI (first_var)); + + /* make sure we can deal with this number of dimensions */ + if (gdim > MAXDIM) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "ApplyBndNoise: Variable dimension of %d not supported", gdim); + return (-1); + } + + /* check the direction parameter */ + if (abs (dir) > gdim) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "ApplyBndNoise: direction %d greater than dimension %d", + dir, gdim); + return (-2); + } + + /* initialize arrays for variables with less dimensions than MAXDIM + so that we can use the INDEX_3D macro later on */ + for (i = gdim; i < MAXDIM; i++) + { + lsh[i] = 0; + lssh[i] = 1; + } + + /* get the directional staggering of the group */ + CCTK_GroupStaggerDirArrayGI (dstag, gdim, gindex); + + /* get the current timelevel */ + timelvl = 0; + + /* see if we have a symmetry array */ + sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry"); + + + /* now loop over all variables */ + for (var = first_var; var < first_var + num_vars; var++) + { + /* Apply condition if: + + boundary is not a symmetry boundary (no symmetry or unset(=unsed)) + + boundary is a physical boundary + + have enough grid points + */ + memset (doBC, 1, sizeof (doBC)); + if (sGHex) + { + for (i = 0; i < 2 * gdim; i++) + { + doBC[i] = sGHex->GFSym[var][i] == GFSYM_NOSYM || + sGHex->GFSym[var][i] == GFSYM_UNSET; + } + } + for (i = 0; i < gdim; i++) + { + lsh[i] = GH->cctk_lsh[i]; + lssh[i] = GH->cctk_lssh[CCTK_LSSH_IDX (dstag[i], i)]; + doBC[i*2] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2]; + doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1]; + if (dir != 0) + { + doBC[i*2] &= (dir < 0 && (i + 1 == abs (dir))); + doBC[i*2+1] &= (dir > 0 && (i + 1 == abs (dir))); + } + } + + /* now apply the boundaries face by face */ + if (gdim > 0) + { +#ifdef DEBUG_BOUNDARY + if (doBC[0]) + { + printf("Boundary: Applying lower x noise to boundary\n"); + } + if (doBC[1]) + { + printf("Boundary: Applying upper x noise to boundary\n"); + } +#endif /* DEBUG_BOUNDARY */ + /* lower x */ + NOISY_BOUNDARY (doBC[0], stencil[0], lssh[1], lssh[2], + i, j, k, rand); + /* upper x */ + NOISY_BOUNDARY (doBC[1], stencil[0], lssh[1], lssh[2], + lssh[0]-i-1, j, k, rand); + + } + if (gdim > 1) + + { +#ifdef DEBUG_BOUNDARY + if (doBC[2]) + { + printf("Boundary: Applying lower y noise to boundary\n"); + } + if (doBC[3]) + { + printf("Boundary: Applying upper y noise to boundary\n"); + } +#endif /* DEBUG_BOUNDARY */ + /* lower y */ + NOISY_BOUNDARY (doBC[2], lssh[0], stencil[1], lssh[2], + i, j, k, rand); + /* upper y */ + NOISY_BOUNDARY (doBC[3], lssh[0], stencil[1], lssh[2], + i, lssh[1]-j-1, k, rand); + } + if (gdim > 2) + { +#ifdef DEBUG_BOUNDARY + if (doBC[4]) + { + printf("Boundary: Applying lower z noise to boundary\n"); + } + if (doBC[5]) + { + printf("Boundary: Applying upper z noise to boundary\n"); + } +#endif /* DEBUG_BOUNDARY */ + /* lower z */ + NOISY_BOUNDARY (doBC[4], lssh[0], lssh[1], stencil[2], + i, j, k, rand); + /* upper z */ + NOISY_BOUNDARY (doBC[5], lssh[0], lssh[1], stencil[2], + i, j, lssh[2]-k-1, rand); + } + } + + return(0); +} + + +void +add_bc_noise_to_group (int idx, const char* optstring, void* cctkGH) +{ + int i, j, k, ijk; + CCTK_REAL* data; + cGH* GH = cctkGH; + int sw[3]; + + sw[0] = 1; + sw[1] = 1; + sw[2] = 1; + + BndNoiseGI(GH, sw, idx); +} + + +void +bc_noise(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + if (CCTK_TraverseString(noisy_bc_vars, add_bc_noise_to_group, cctkGH, + CCTK_GROUP_OR_VAR) < 0) + { + CCTK_WARN (1, "Failed to parse 'Noise::noisy_bc_vars' parameter"); + } +} |