diff options
Diffstat (limited to 'src/setup_epsdis.c')
| -rw-r--r-- | src/setup_epsdis.c | 233 |
1 files changed, 233 insertions, 0 deletions
diff --git a/src/setup_epsdis.c b/src/setup_epsdis.c new file mode 100644 index 0000000..9a3b2cf --- /dev/null +++ b/src/setup_epsdis.c @@ -0,0 +1,233 @@ +/* $Header$ */ + +#include <assert.h> +#include <math.h> +#include <stdlib.h> +#include <stdio.h> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#define MAXDIM 3 + +void +setup_epsdis (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int ni,nj,nk; + int i,j,k,s; + int index,indexP; + int ierr; + int npts; + int *inds; + CCTK_REAL *xa,*ya,*za,*rads; + CCTK_REAL maxrad; + CCTK_REAL xmin,xmax, ymin,ymax, zmin,zmax; + CCTK_REAL odx,ody,odz; + CCTK_REAL radp; + const CCTK_INT MAXSURFNUM=100; /* XXX hard limit */ + CCTK_INT doBC[2*MAXDIM],symbnd[2*MAXDIM]; + CCTK_INT symtable; + + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING,"Setting up spatially varying dissipation at T=%g", + (double)cctk_time); + } + + ni=cctk_lsh[0]; + nj=cctk_lsh[1]; + nk=cctk_lsh[2]; + + for (i=0;i<ni*nj*nk;i++) { + epsdisA[i]=epsdis; + } + + + + if (extra_dissipation_at_outerbound) + { + symtable = SymmetryTableHandleForGrid (cctkGH); + if (symtable < 0) CCTK_WARN (1, "unable to get symmetry table"); + ierr=Util_TableGetIntArray(symtable, 6, symbnd, "symmetry_handle"); + if (ierr != 6) CCTK_WARN (1, "unable to get symmetry handle"); + for (i = 0; i < 6; i++) { + doBC[i] = cctk_bbox[i]!=0 && symbnd[i] < 0; + } + if(doBC[0]) { + for (k=0;k<nk;k++) { + for (j=0;j<nj;j++) { + for (i=0;i<outer_bound_npoints;i++) { + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + indexP= CCTK_GFINDEX3D(cctkGH,outer_bound_npoints,j,k); + epsdisA[index]=epsdis+epsdis*ob_slope*abs(x[index]-x[indexP]); + } + } + } + } + if(doBC[1]) { + for (k=0;k<nk;k++) { + for (j=0;j<nj;j++) { + for (i=ni;i>ni-outer_bound_npoints;i--) { + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + indexP= CCTK_GFINDEX3D(cctkGH,ni-outer_bound_npoints,j,k); + epsdisA[index]=epsdis+epsdis*ob_slope*abs(x[index]-x[indexP]); + } + } + } + } + + if(doBC[2]) { + for (k=0;k<nk;k++) { + for (j=0;j<outer_bound_npoints;j++) { + for (i=0;i<ni;i++) { + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + indexP= CCTK_GFINDEX3D(cctkGH,i,outer_bound_npoints,k); + epsdisA[index]=epsdis+epsdis*ob_slope*abs(y[index]-y[indexP]); + } + } + } + } + if(doBC[3]) { + for (k=0;k<nk;k++) { + for (j=nj;j>nj-outer_bound_npoints;j--) { + for (i=0;i<ni;i++) { + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + indexP= CCTK_GFINDEX3D(cctkGH,i,nj-outer_bound_npoints,k); + epsdisA[index]=epsdis+epsdis*ob_slope*abs(y[index]-y[indexP]); + } + } + } + } + + if(doBC[4]) { + for (k=0;k<outer_bound_npoints;k++) { + for (j=0;j<nj;j++) { + for (i=0;i<ni;i++) { + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + indexP= CCTK_GFINDEX3D(cctkGH,i,j,outer_bound_npoints); + epsdisA[index]=epsdis+epsdis*ob_slope*abs(z[index]-z[indexP]); + } + } + } + } + if(doBC[5]) { + for (k=nk;k>nk-outer_bound_npoints;k--) { + for (j=0;j<nj;j++) { + for (i=0;i<ni;i++) { + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + indexP= CCTK_GFINDEX3D(cctkGH,i,j,nk-outer_bound_npoints); + epsdisA[index]=epsdis+epsdis*ob_slope*abs(z[index]-z[indexP]); + } + } + } + } + } + + + + if (extra_dissipation_in_horizons && cctk_iteration%update_ah_every == 0) + { + if (verbose) { + CCTK_INFO("Linear Interpolation into AH surfaces"); + } + + for (s=0;s<MAXSURFNUM;s++) + { + if (surface_number[s]==-1 && horizon_number[s]==-1) { + continue; + } + + if (surface_number[s]<0 || horizon_number[s]<=0) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Invalid specification for horizon %d", s); + continue; + } + + if (! sf_valid[surface_number[s]]) { + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + "Invalid Surface: s=%d, sf_va=%d, surf_no=%d, ah_no=%d", + s,sf_valid[surface_number[s]],surface_number[s],horizon_number[s]); + } + continue; + } + + if (! CCTK_IsFunctionAliased ("HorizonRadiusInDirection")) + { + CCTK_WARN (1, "The aliased function \"HorizonRadiusInDirection\" must be defined when the parameter \"extra_dissipation_in_horizons\" is set and one of the sources is AHFinderDirect"); + continue; + } + + maxrad=sf_max_radius[surface_number[s]]; + odx=sf_origin_x[surface_number[s]]; + ody=sf_origin_y[surface_number[s]]; + odz=sf_origin_z[surface_number[s]]; + xmin=odx-maxrad; + xmax=odx+maxrad; + ymin=ody-maxrad; + ymax=ody+maxrad; + zmin=odz-maxrad; + zmax=odz+maxrad; + + npts=0; + for (i=0;i<ni*nj*nk;i++) { + if ( x[i]<=xmax&&x[i]>=xmin + &&y[i]<=ymax&&y[i]>=ymin + &&z[i]<=zmax&&z[i]>=zmin) + { + npts++; + } + } + + xa=(CCTK_REAL *) malloc(npts*sizeof(CCTK_REAL)); + ya=(CCTK_REAL *) malloc(npts*sizeof(CCTK_REAL)); + za=(CCTK_REAL *) malloc(npts*sizeof(CCTK_REAL)); + rads=(CCTK_REAL *) malloc(npts*sizeof(CCTK_REAL)); + inds=(CCTK_INT *) malloc(npts*sizeof(CCTK_INT)); + + for (i=0;i<npts;i++) { + rads[i]=0; + inds[i]=0; + } + + j=0; + for (i=0;i<ni*nj*nk;i++) { + if ( x[i]<=xmax&&x[i]>=xmin + &&y[i]<=ymax&&y[i]>=ymin + &&z[i]<=zmax&&z[i]>=zmin) + { + xa[j]=x[i]; + ya[j]=y[i]; + za[j]=z[i]; + inds[j]=i; + j++; + } + } + + ierr=HorizonRadiusInDirection(horizon_number[s], + npts, + xa, ya, za, rads); + assert(!ierr); + + for (i=0;i<npts;i++) { + radp=sqrt(xa[i]*xa[i]+ya[i]*ya[i]+za[i]*za[i]); + if (radp<=rads[i]) { + epsdisA[inds[i]]=epsdis+ ah_slope*(rads[i]-radp)*epsdis; + } + } + + free(xa); + free(ya); + free(za); + free(rads); + free(inds); + } + } +} |