From 5766861f7d6bbe2e106b20f479c447855f5a8448 Mon Sep 17 00:00:00 2001 From: allen Date: Wed, 24 Jan 2001 09:56:16 +0000 Subject: Adding directional boundary conditions for everything apart from Robin. See thorn documentation for details Next to come, boundary conditions for grid functions which aren't CCTK_REALs. git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Boundary/trunk@131 6a38eb6e-646e-4a02-a296-d141613ad6c4 --- doc/documentation.tex | 129 +++++++-- src/Boundary.h | 71 +++++ src/CopyBoundary.c | 529 ++++++++++++++++++++++++++++-------- src/FlatBoundary.c | 505 ++++++++++++++++++++++++++-------- src/RadiationBoundary.c | 702 +++++++++++++++++++++++++++++++++++------------- src/ScalarBoundary.c | 22 +- 6 files changed, 1517 insertions(+), 441 deletions(-) diff --git a/doc/documentation.tex b/doc/documentation.tex index d501260..c1c2176 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -177,7 +177,7 @@ integer \> group\_index\\ \end{tabbing} } -\subsection*{Arguments} +\subsubsection*{Arguments} \begin{Lentry} \item[{\tt ierr}] Return value, negative value indicates the boundary condition was not successfully applied @@ -201,6 +201,8 @@ boundary value of phi {\tt phi(nx,j,k)}, on the positive x-boundary will be copied from {\tt phi(nx-1,j,k)}. \subsubsection*{Calling from C:} + +{\bf All Coordinate Directions:} \begin{verbatim} int ierr = BndFlatVN(cGH *cctkGH, int *stencil_size, char *variable_name) int ierr = BndFlatGN(cGH *cctkGH, int *stencil_size, char *group_name) @@ -208,30 +210,62 @@ int ierr = BndFlatVI(cGH *cctkGH, int *stencil_size, int variable_index) int ierr = BndFlatGI(cGH *cctkGH, int *stencil_size, int group_index) \end{verbatim} +{\bf Individual Coordinate Directions:} +\begin{verbatim} +int ierr = BndFlatVN(cGH *cctkGH, int stencil, int dir, char *variable_name) +int ierr = BndFlatGN(cGH *cctkGH, int stencil, int dir, char *group_name) +int ierr = BndFlatVI(cGH *cctkGH, int stencil, int dir, int variable_index) +int ierr = BndFlatGI(cGH *cctkGH, int stencil, int dir, int group_index) +\end{verbatim} + \subsubsection*{Calling from Fortran:} +{\bf All Coordinate Directions:} +\begin{verbatim} +call BndFlatVN(ierr, cctkGH, stencil_array, variable_name) +call BndFlatGN(ierr, cctkGH, stencil_array, group_name) +call BndFlatVI(ierr, cctkGH, stencil_array, variable_index) +call BndFlatGI(ierr, cctkGH, stencil_array, group_index) +\end{verbatim} + +{\bf Individual Coordinate Directions:} \begin{verbatim} -call BndFlatVN(ierr, cctkGH, stencil_size, variable_name) -call BndFlatGN(ierr, cctkGH, stencil_size, group_name) -call BndFlatVI(ierr, cctkGH, stencil_size, variable_index) -call BndFlatGI(ierr, cctkGH, stencil_size, group_index) +call BndFlatVN(ierr, cctkGH, stencil, dir, variable_name) +call BndFlatGN(ierr, cctkGH, stencil, dir, group_name) +call BndFlatVI(ierr, cctkGH, stencil, dir, variable_index) +call BndFlatGI(ierr, cctkGH, stencil, dir, group_index) \end{verbatim} where +{\tt +\begin{tabbing} +character*(*) \= variable\_name\=\kill +integer \> ierr \\ +CCTK\_POINTER \> cctkGH\\ +integer \> dir\\ +integer \> stencil\\ +integer \> stencil\_array(dim)\\ +character*(*) \> variable\_name\\ +character*(*) \> group\_name\\ +integer \> variable\_index\\ +integer \> group\_index\\ +\end{tabbing} +} + +\subsubsection*{Arguments} \begin{Lentry} -\item[{\tt integer ierr}] return value, operation failed when return -value {\em negative} -\item[{\tt CCTK\_POINTER cctkGH}] grid hierarchy pointer -\item[{\tt CCTK\_REAL var0}] scalar value to apply -\item[{\tt integer stencil\_size(dim)}] array of size {\tt dim} -(dimension of the gridfunction). To how many points from the outer -boundary to apply the boundary condition. -\item[{\tt character*(*) variable\_name}] Name of the variable. -\item[{\tt character*(*) group\_name}] Name of the group. -\item[{\tt integer variable\_index}] Variable index. -\item[{\tt integer group\_index}] Group index. +\item[{\tt ierr}] Return value, negative value indicates the +boundary condition was not successfully applied +\item[{\tt cctkGH}] Grid hierarchy pointer +\item[{\tt dir}] Coordinate direction in which to apply boundary condition +\item[{\tt stencil\_size}] Array with dimension of the grid function, containing the stencil width to apply the boundary at +\item[{\tt variable\_name}] Name of the variable +\item[{\tt group\_name}] Name of the group +\item[{\tt variable\_index}] Variable index +\item[{\tt group\_index}] Group index \end{Lentry} + \subsection{Radiation Boundary Condition} This is a two level scheme. Grid functions are given for the current time @@ -293,22 +327,42 @@ the region where the characteristic speed is constant. Notice that this speed does not have to be 1. \subsubsection*{Calling from C:} + +{\bf All Coordinate Directions:} \begin{verbatim} int ierr = BndRadiativeVN(cGH *cctkGH, int *stencil_size, - CCTK_REAL speed, CCTK_REAL limit, + CCTK_REAL limit, CCTK_REAL speed, char *variable_name, char *variable_name_past) int ierr = BndRadiativeGN(cGH *cctkGH, int *stencil_size, - CCTK_REAL speed, CCTK_REAL limit, + CCTK_REAL limit, CCTK_REAL speed, char *group_name, char *group_name_past) int ierr = BndRadiativeVI(cGH *cctkGH, int *stencil_size, - CCTK_REAL speed, CCTK_REAL limit, + CCTK_REAL limit, CCTK_REAL speed, int variable_index, int variable_index_past) int ierr = BndRadiativeGI(cGH *cctkGH, int *stencil_size, - CCTK_REAL speed, CCTK_REAL limit, + CCTK_REAL limit, CCTK_REAL speed, + int group_index, int group_index_past) +\end{verbatim} + +{\bf Individual Coordinate Directions:} +\begin{verbatim} +int ierr = BndRadiativeVN(cGH *cctkGH, int stencil, int dir, + CCTK_REAL limit, CCTK_REAL speed, + char *variable_name, char *variable_name_past) +int ierr = BndRadiativeGN(cGH *cctkGH, int *stencil, int dir, + CCTK_REAL limit, CCTK_REAL speed, + char *group_name, char *group_name_past) +int ierr = BndRadiativeVI(cGH *cctkGH, int *stencil, int dir, + CCTK_REAL limit, CCTK_REAL speed, + int variable_index, int variable_index_past) +int ierr = BndRadiativeGI(cGH *cctkGH, int *stencil, int dir, + CCTK_REAL limit, CCTK_REAL speed, int group_index, int group_index_past) \end{verbatim} \subsubsection*{Calling from Fortran:} + +{\bf All Coordinate Directions:} \begin{verbatim} call BndRadiativeVN(ierr, cctkGH, stencil_size, speed, limit, variable_name, variable_name_past) @@ -319,6 +373,39 @@ call BndRadiativeVN(ierr, cctkGH, stencil_size, speed, limit, call BndRadiativeVN(ierr, cctkGH, stencil_size, speed, limit, group_index, group_index_past) \end{verbatim} + +{\bf Individual Coordinate Directions:} +\begin{verbatim} +call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, + variable_name, variable_name_past) +call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, + group_name, group_name_past) +call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, + variable_index, variable_index_past) +call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, + group_index, group_index_past) +\end{verbatim} +where +{\tt +\begin{tabbing} +character*(*) \= variable\_name\=\kill +integer \> ierr \\ +CCTK\_POINTER \> cctkGH\\ +integer \> dir\\ +integer \> stencil\\ +integer \> stencil\_array(dim)\\ +character*(*) \> variable\_name\\ +character*(*) \> group\_name\\ +integer \> variable\_index\\ +integer \> group\_index\\ +CCTK_REAL\>speed\\ +CCTK_REAL\>limit\\ +\end{tabbing} +} + + + + where \begin{Lentry} \item[{\tt integer ierr}] return value, operation failed when return @@ -329,7 +416,7 @@ value {\em negative} boundary to apply the boundary condition. \item[{\tt CCTK\_REAL speed}] wave speed used for boundary condition ($v$). -\item[{\tt CCTK\_REAL limit}] radation constant +\item[{\tt CCTK\_REAL limit}] asymptotic value of function at infinity \item[{\tt character*(*) variable\_name}] the name of the grid function to which the boundary condition will be applied diff --git a/src/Boundary.h b/src/Boundary.h index a07ebe9..acddeb5 100644 --- a/src/Boundary.h +++ b/src/Boundary.h @@ -57,6 +57,27 @@ int BndScalarVN(cGH *GH, /* Copying boundaries */ +int BndCopyDirGI(cGH *GH, + int stencil_size, + int dir, + int gi_to, + int gi_from); +int BndCopyDirGN(cGH *GH, + int stencil_size, + int dir, + const char *gn_to, + const char *gn_from); +int BndCopyDirVI(cGH *GH, + int stencil_size, + int dir, + int vi_to, + int vi_from); +int BndCopyDirVN(cGH *GH, + int stencil_size, + int dir, + const char *vn_to, + const char *vn_from); + int BndCopyGI(cGH *GH, int *stencil, int gi_to, @@ -76,6 +97,39 @@ int BndCopyVN(cGH *GH, /* Radiative boundaries */ +int BndRadiativeDirGI(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL v0, + int gi, + int gi_p); + +int BndRadiativeDirGN(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL v0, + const char *gn, + const char *gn_p); + +int BndRadiativeDirVI(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL v0, + int vi, + int vi_p); + +int BndRadiativeDirVN(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL v0, + const char *vn, + const char *vn_p); + + int BndRadiativeGI(cGH *GH, int *sw, CCTK_REAL var0, @@ -132,6 +186,23 @@ int BndRobinVN(cGH *GH, /* Flat boundaries */ +int BndFlatDirGI(cGH *GH, + int stencil, + int dir, + int gi); +int BndFlatDirGN(cGH *GH, + int stencil, + int dir, + const char *gn); +int BndFlatDirVI(cGH *GH, + int stencil, + int dir, + int vi); +int BndFlatDirVN(cGH *GH, + int stencil, + int dir, + const char *vn); + int BndFlatGI(cGH *GH, int *stencil, int gi); diff --git a/src/CopyBoundary.c b/src/CopyBoundary.c index 93ebcb3..117b742 100644 --- a/src/CopyBoundary.c +++ b/src/CopyBoundary.c @@ -25,34 +25,33 @@ /* Internal routine prototypes */ static int BndApplyCopy3Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, int *stencil_size, CCTK_REAL *var, CCTK_REAL *varp); static int BndApplyCopy2Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, int *stencil_size, CCTK_REAL *var, CCTK_REAL *varp); static int BndApplyCopy1Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, int *stencil_size, CCTK_REAL *var, CCTK_REAL *varp); static int ApplyBndCopy(cGH *GH, - int *stencil, + int stencil, + int *stencil_array, int first_var, int second_var, - int num_vars); + int num_vars, + int dir); /* Local variables */ @@ -62,6 +61,43 @@ static int ApplyBndCopy(cGH *GH, ******************** External Routines ************************ ********************************************************************/ +/*@@ + @routine BndCopyDirVI + @date Sat Mar 20 + @author Gabrielle Allen + @desc + Apply copy boundary routines by var index in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int BndCopyDirVI(cGH *GH, + int stencil_size, + int dir, + int vi_to, + int vi_from) +{ + return ApplyBndCopy(GH,stencil_size,NULL,vi_to,vi_from,1,dir); +} + +void CCTK_FCALL CCTK_FNAME(BndCopyDirVI) + (int *ierr, + cGH *GH, + int *stencil_size, + int *dir, + int *vi_to, + int *vi_from) +{ + *ierr = BndCopyDirVI(GH, *stencil_size, *dir, *vi_to, *vi_from); +} + + + /*@@ @routine BndCopyVI @date Thu Mar 2 11:02:10 2000 @@ -82,21 +118,54 @@ int BndCopyVI(cGH *GH, int vi_to, int vi_from) { - return ApplyBndCopy(GH,stencil,vi_to,vi_from,1); + return ApplyBndCopy(GH,-1,stencil,vi_to,vi_from,1,0); } -void CCTK_FCALL CCTK_FNAME(BndCopyVI) +/* ====================================================== */ + + /*@@ + @routine BndCopyDirGI + @date Sat Jan 20 2001 + @author Gabrielle Allen + @desc + Apply copy boundaries by group index in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int BndCopyDirGI(cGH *GH, + int stencil_size, + int dir, + int gi_to, + int gi_from) +{ + int first_vi_to,first_vi_from; + int numvars; + + first_vi_to = CCTK_FirstVarIndexI(gi_to); + first_vi_from= CCTK_FirstVarIndexI(gi_from); + numvars = CCTK_NumVarsInGroupI(gi_to); + + return ApplyBndCopy(GH,stencil_size,NULL,first_vi_to,first_vi_from,numvars,dir); +} + +void CCTK_FCALL CCTK_FNAME(BndCopyDirGI) (int *ierr, cGH *GH, - int *stencil, - int *vi_to, - int *vi_from) + int *stencil_size, + int *dir, + int *gi_to, + int *gi_from) { - *ierr = BndCopyVI(GH, stencil, *vi_to, *vi_from); + *ierr = BndCopyDirGI(GH, *stencil_size, *dir, *gi_to, *gi_from); + return; } - - /*@@ @routine BndCopyGI @date Thu Mar 2 11:07:11 2000 @@ -124,7 +193,7 @@ int BndCopyGI(cGH *GH, first_vi_from= CCTK_FirstVarIndexI(gi_from); numvars = CCTK_NumVarsInGroupI(gi_to); - return ApplyBndCopy(GH,stencil,first_vi_to,first_vi_from,numvars); + return ApplyBndCopy(GH,-1,stencil,first_vi_to,first_vi_from,numvars,0); } void CCTK_FCALL CCTK_FNAME(BndCopyGI) @@ -137,7 +206,62 @@ void CCTK_FCALL CCTK_FNAME(BndCopyGI) *ierr = BndCopyGI(GH, stencil, *gi_to, *gi_from); } +/* ======================================================= */ +/*@@ + @routine BndCopyDirGN + @date Sat Jan 20 2001 + @author Gabrielle Allen + @desc + Apply copy boundary routines by group name in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int BndCopyDirGN(cGH *GH, + int stencil_size, + int dir, + const char *gn_to, + const char *gn_from) +{ + int retval; + int gi_to, gi_from; + + gi_to = CCTK_GroupIndex(gn_to); + gi_from = CCTK_GroupIndex(gn_from); + + if ((gi_to>-1)&&(gi_from>-1)) + retval = BndCopyDirGI(GH, stencil_size, dir, gi_to, gi_from); + else + { + CCTK_VWarn(2,__LINE__,__FILE__,"Boundary", + "BndCopyDirGN: Cannot find group indices for %s, %s ", + gn_to,gn_from); + retval = -1; + } + + return retval; + +} + +void CCTK_FCALL CCTK_FNAME(BndCopyDirGN) + (int *ierr, + cGH *GH, + int *stencil_size, + int *dir, + TWO_FORTSTRINGS_ARGS) +{ + TWO_FORTSTRINGS_CREATE(gn_to, gn_from) + *ierr = BndCopyDirGN(GH, *stencil_size, *dir, gn_to, gn_from); + free(gn_to); + free(gn_from); + return; +} /*@@ @routine BndCopyGN @@ -189,9 +313,64 @@ void CCTK_FCALL CCTK_FNAME(BndCopyGN) *ierr = BndCopyGN(GH, stencil, gn_to, gn_from); free(gn_to); free(gn_from); + return; } +/* ================================================================ */ + +/*@@ + @routine BndCopyDirVN + @date Sat Jan 20 2001 + @author Gabrielle Allen + @desc + Apply copy boundary routines by variable name in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ +int BndCopyDirVN(cGH *GH, + int stencil_size, + int dir, + const char *vn_to, + const char *vn_from) +{ + int vi_to, vi_from; + int retval; + + vi_to = CCTK_VarIndex(vn_to); + vi_from = CCTK_VarIndex(vn_from); + + if ((vi_to>-1)&&(vi_from>-1)) + retval = BndCopyDirVI(GH, stencil_size, dir, vi_to, vi_from); + else + { + CCTK_VWarn(2,__LINE__,__FILE__,"Boundary", + "BndCopyDirVN: Cannot find variable indices for %s, %s", + vn_to,vn_from); + retval = -1; + } + + return retval; +} + +void CCTK_FCALL CCTK_FNAME(BndCopyDirVN) + (int *ierr, + cGH *GH, + int *stencil_size, + int *dir, + TWO_FORTSTRINGS_ARGS) +{ + TWO_FORTSTRINGS_CREATE(vn_to, vn_from) + *ierr = BndCopyDirVN(GH, *stencil_size, *dir, vn_to, vn_from); + free(vn_to); + free(vn_from); + return; +} /*@@ @routine BndCopyVN @@ -267,117 +446,214 @@ void CCTK_FCALL CCTK_FNAME(BndCopyVN) @@*/ static int ApplyBndCopy(cGH *GH, - int *stencil, + int stencil, + int *stencil_array, int first_var, int second_var, - int num_vars) + int num_vars, + int dir) { int symmetry_handle; /* handle for the optional symmetry structure */ int vi, vi2, gi, dim; int idim; + int retval=0; + int dirloop1,dirloop2; + int doBC,count; int berr,ierr; - int *doBC, *dstag, *lssh; - SymmetryGHex *sGHex; + int *dstag; + int *lssh; + int *sw; int timelevel; /* timelevel that condition applied on */ + SymmetryGHex *sGHex; + +#ifdef DEBUG_BOUNDARY + printf("Input arguments for ApplyBndCopy:\n"); + printf("GH = %x\n",GH); + printf("stencil = %d\n",stencil); + printf("stencil_array = %x\n",stencil_array); + printf("first_var = %d\n",first_var); + printf("second_var = %d\n",second_var); + printf("num_vars = %d\n",num_vars); + printf("dir = %d\n",dir); + printf("-----------------------------------\n"); +#endif + + /* Get group dimension */ + dim = CCTK_GroupDimFromVarI(first_var); - /* See if we have a symmetry array */ - symmetry_handle = CCTK_GHExtensionHandle("Symmetry"); - if (symmetry_handle < 0) + if (dim<1 || dim>3) { - sGHex = NULL; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndCopy: Invalid grid function dimension %d", + dim); + retval=-1; } else { - sGHex = (SymmetryGHex*)GH->extensions[symmetry_handle]; - } - /* Get group dimension */ - dim = CCTK_GroupDimFromVarI(first_var); + /* See if we have a symmetry array */ + symmetry_handle = CCTK_GHExtensionHandle("Symmetry"); + if (symmetry_handle < 0) + { + sGHex = NULL; + } + else + { + sGHex = (SymmetryGHex*)GH->extensions[symmetry_handle]; + } + + /* Set up stencil width array if needed */ + if (stencil_array && dir == 0) + { + sw = stencil_array; + } + else if (dir==0) + { + CCTK_WARN(0,"ApplyBndCopy: Direction 0 invalid for directional BC"); + } + else if (dir>dim||-dir>dim) + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndCopy: Direction %d greater than dimension %d", + dir,dim); + } + else + { + int i; + /* Should be single direction */ + sw = (int *)malloc(dim*sizeof(int)); + if (sw) + { + for (i=0;i 0) + { + sw[dir-1] = stencil; + } + else + { + sw[-dir-1] = stencil; + } + } + else + { + CCTK_WARN(0,"ApplyBndCopy: Malloc sw failed"); + } + } - /* allocate arrays */ - doBC = (int *)malloc((2*dim)*sizeof(int)); - dstag = (int *)malloc(dim*sizeof(int)); - lssh = (int *)malloc(dim*sizeof(int)); + /* allocate arrays */ + dstag = (int *)malloc(dim*sizeof(int)); + lssh = (int *)malloc(dim*sizeof(int)); - /* get the directional staggering of the group */ - gi = CCTK_GroupIndexFromVarI(first_var); - ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi); + /* get the directional staggering of the group */ + gi = CCTK_GroupIndexFromVarI(first_var); + ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi); - /* get the current timelevel */ - timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1; - /* if (timelevel > 0) - timelevel--;*/ + /* get the current timelevel */ + timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1; - for (vi=first_var; vi0) { - for (idim=0;idimGFSym[vi][idim*2]==GFSYM_NOSYM)|| - (sGHex->GFSym[vi][idim*2]==GFSYM_UNSET)) && - GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]); - doBC[idim*2+1] = (((sGHex->GFSym[vi][idim*2+1]==GFSYM_NOSYM)|| - (sGHex->GFSym[vi][idim*2+1]==GFSYM_UNSET)) && - GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]); - lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)]; - } + dirloop1 = 2*(dir-1)+1; + dirloop2 = dirloop1+1; + } + else if (dir<0) + { + dirloop1 = 2*(-dir-1); + dirloop2 = dirloop1+1; } else { - for (idim=0;idimcctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]); - doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]); - lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)]; + lssh[idim] = (int)GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)]; } - } + + for (count=dirloop1;countGFSym[vi][count]==GFSYM_NOSYM) || + (sGHex->GFSym[vi][count]==GFSYM_UNSET) ) && + GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? + count : -1; + } + else + { + doBC = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? + count : -1 ; + } - switch (dim) - { - case 1: berr = BndApplyCopy1Di(GH, - dim, + if (doBC >=0) + { + switch (dim) + { + case 1: + berr = BndApplyCopy1Di(GH, doBC, lssh, - stencil, + sw, GH->data[vi][timelevel], - GH->data[vi2][timelevel]); break; - case 2: berr = BndApplyCopy2Di(GH, - dim, + GH->data[vi2][timelevel]); + break; + case 2: + berr = BndApplyCopy2Di(GH, doBC, lssh, - stencil, + sw, GH->data[vi][timelevel], - GH->data[vi2][timelevel]); break; - case 3: berr = BndApplyCopy3Di(GH, - dim, + GH->data[vi2][timelevel]); + break; + case 3: + berr = BndApplyCopy3Di(GH, doBC, lssh, - stencil, + sw, GH->data[vi][timelevel], - GH->data[vi2][timelevel]); break; - default : berr = -1; - CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "ApplyBndCopy: No BC for dim>3: grid variables '%s'/'%s'", - CCTK_VarName(vi),CCTK_VarName(vi2)); + GH->data[vi2][timelevel]); + break; + default : + berr = -1; + CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", + "No copy boundary for dim>3: grid variable '%s'", + CCTK_VarName(vi)); + } + } + } + + /* retval = (berr>-1) ? 0 : -1;*/ + } + + free(dstag); + free(lssh); + + if (!stencil_array) + { + free(sw); } - berr = (berr>-1) ? 0 : -1; } - - free(dstag); - free(doBC); - free(lssh); - - return(ierr); + + return retval; } @@ -391,7 +667,7 @@ static int ApplyBndCopy(cGH *GH, for the 3d case @enddesc @calls - @calledby BndCopyGI, BndCopyVN, BndCopyGN, BndCopyVI, + @calledby @history @endhistory @@ -399,8 +675,7 @@ static int ApplyBndCopy(cGH *GH, @@*/ static int BndApplyCopy3Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, int *stencil_size, CCTK_REAL *var, @@ -409,8 +684,11 @@ static int BndApplyCopy3Di(cGH *GH, int i,j,k,sw,lin; /* lower x */ - if (doBC[0] == 1) + if (doBC == 0) { +#ifdef DEBUG_BOUNDARY + printf("Applying copy boundary (lower x)\n"); +#endif for (k=0;k #include @@ -28,30 +28,29 @@ /* Internal routine prototypes */ static int BndApplyFlat3Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, - int *stencil, + int *stencilarray, CCTK_REAL *var); static int BndApplyFlat2Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, - int *stencil, + int *stencil_array, CCTK_REAL *var); static int BndApplyFlat1Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, - int *stencil, + int *stencil_array, CCTK_REAL *var); static int ApplyBndFlat(cGH *GH, - int *stencil, + int stencil, + int *stencil_array, int first_var, - int num_vars); + int num_vars, + int dir); /* Local variables */ @@ -61,6 +60,39 @@ static int ApplyBndFlat(cGH *GH, ******************** External Routines ************************ ********************************************************************/ +/*@@ + @routine BndFlatDirGI + @date Sun Jan 21 2001 + @author Gabrielle Allen + @desc + Apply flat boundary conditions by group index in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int BndFlatDirGI(cGH *GH, int stencil, int dir, int gi) +{ + int numvars, first_vi; + + first_vi = CCTK_FirstVarIndexI(gi); + numvars = CCTK_NumVarsInGroupI(gi); + + return ApplyBndFlat(GH,stencil,NULL,first_vi,numvars,dir); +} + +void CCTK_FCALL CCTK_FNAME(BndFlatDirGI) + (int *ierr, cGH *GH, int *stencil, int *dir, int *gi) +{ + *ierr = BndFlatDirGI(GH, *stencil, *dir, *gi); + return; +} + + /*@@ @routine BndFlatGI @date Thu Mar 2 11:11:40 2000 @@ -83,7 +115,7 @@ int BndFlatGI(cGH *GH, int *stencil, int gi) first_vi = CCTK_FirstVarIndexI(gi); numvars = CCTK_NumVarsInGroupI(gi); - return ApplyBndFlat(GH,stencil,first_vi,numvars); + return ApplyBndFlat(GH,-1,stencil,first_vi,numvars,0); } void CCTK_FCALL CCTK_FNAME(BndFlatGI) @@ -94,6 +126,51 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGI) } + /* ********************************************************************/ + + +/*@@ + @routine BndFlatDirGN + @date Sun Jan 21 2001 + @author Gabrielle Allen + @desc + Apply flat boundary conditions by group name in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int BndFlatDirGN(cGH *GH, int stencil, int dir, const char *gn) +{ + int retval; + int gi; + gi = CCTK_GroupIndex(gn); + + if (gi>-1) + { + retval=BndFlatDirGI(GH, stencil, dir, gi); + } + else + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "BndFlatDirGN: Grid variable %s not found",gn); + retval=-1; + } + return retval; +} + +void CCTK_FCALL CCTK_FNAME(BndFlatDirGN) + (int *ierr, cGH *GH, int *stencil, int *dir, ONE_FORTSTRING_ARG) +{ + ONE_FORTSTRING_CREATE(gn) + *ierr = BndFlatDirGN(GH, *stencil, *dir, gn); + free(gn); + return; +} /*@@ @@ -140,6 +217,37 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGN) } + /* ********************************************************************/ + + +/*@@ + @routine BndFlatDirVI + @date Sun Jan 21 2001 + @author Gabrielle Allen + @desc + Apply flat boundary conditions by variable index in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int BndFlatDirVI(cGH *GH, int stencil, int dir, int vi) +{ + return ApplyBndFlat(GH,stencil,NULL,vi,1,dir); +} + +void CCTK_FCALL CCTK_FNAME(BndFlatDirVI) + (int *ierr, cGH *GH, int *stencil, int *dir, int *vi) +{ + *ierr = BndFlatDirVI(GH, *stencil, *dir, *vi); + return; +} + + /*@@ @routine BndFlatVI @@ -158,7 +266,7 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGN) int BndFlatVI(cGH *GH, int *stencil, int vi) { - return ApplyBndFlat(GH,stencil,vi,1); + return ApplyBndFlat(GH,-1,stencil,vi,1,0); } void CCTK_FCALL CCTK_FNAME(BndFlatVI) @@ -169,6 +277,54 @@ void CCTK_FCALL CCTK_FNAME(BndFlatVI) } + /* ********************************************************************/ + + +/*@@ + @routine BndFlatDirVN + @date Sun Jan 21 2001 + @author Gabrielle Allen + @desc + Apply flat boundary conditions by variable name in given direction + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +int BndFlatDirVN(cGH *GH, int stencil, int dir, const char *vn) +{ + int vi; + int retval; + + vi = CCTK_VarIndex(vn); + + if (vi>-1) + { + retval = BndFlatDirVI(GH, stencil, dir, vi); + } + else + { + retval = -1; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "BndFlatDirVN: Grid variable %s not found",vn); + } + return retval; +} + +void CCTK_FCALL CCTK_FNAME(BndFlatDirVN) + (int *ierr, cGH *GH, int *stencil, int *dir, ONE_FORTSTRING_ARG) +{ + ONE_FORTSTRING_CREATE(vn) + *ierr = BndFlatDirVN(GH, *stencil, *dir, vn); + free(vn); + return; +} + + /*@@ @routine BndFlatVN @@ -236,118 +392,208 @@ void CCTK_FCALL CCTK_FNAME(BndFlatVN) @@*/ int ApplyBndFlat(cGH *GH, - int *stencil, + int stencil, + int *stencil_array, int first_var, - int num_vars) + int num_vars, + int dir) { int symmetry_handle; /* handle for the optional symmetry structure */ int vi, gi, dim; int idim; - int berr=0,ierr=0; - int *doBC, *dstag, *lssh; + int berr=0,ierr=0; + int retval=0; + int dirloop1,dirloop2; + int doBC; + int count; + int *dstag; + int *sw; + int *lssh; int timelevel; /* timelevel that condition applied on */ SymmetryGHex *sGHex; - /* See if we have a symmetry array */ - symmetry_handle = CCTK_GHExtensionHandle("Symmetry"); - if (symmetry_handle < 0) +#ifdef DEBUG_BOUNDARY + printf("Input arguments for ApplyBndFlat:\n"); + printf("GH = %x\n",GH); + printf("stencil = %d\n",stencil); + printf("stencil_array = %x\n",stencil_array); + printf("first_var = %d\n",first_var); + printf("num_vars = %d\n",num_vars); + printf("dir = %d\n",dir); + printf("-----------------------------------\n"); +#endif + + /* Get group dimension */ + dim = CCTK_GroupDimFromVarI(first_var); + + if (dim<1 || dim>3) { - sGHex = NULL; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndScalar: Invalid grid function dimension %d", + dim); + retval=-1; } else { - sGHex = (SymmetryGHex*)GH->extensions[symmetry_handle]; - } + /* See if we have a symmetry array */ + symmetry_handle = CCTK_GHExtensionHandle("Symmetry"); + if (symmetry_handle < 0) + { + sGHex = NULL; + } + else + { + sGHex = (SymmetryGHex*)GH->extensions[symmetry_handle]; + } - /* Get group dimension */ - dim = CCTK_GroupDimFromVarI(first_var); + /* Set up stencil width array if needed */ + if (stencil_array && dir == 0) + { + sw = stencil_array; + } + else if (dir==0) + { + CCTK_WARN(0,"ApplyBndScalar: Direction 0 invalid for directional BC"); + } + else if (dir>dim||-dir>dim) + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndScalar: Direction %d greater than dimension %d", + dir,dim); + } + else + { + int i; + /* Should be single direction */ + sw = (int *)malloc(dim*sizeof(int)); + if (sw) + { + for (i=0;i 0) + { + sw[dir-1] = stencil; + } + else + { + sw[-dir-1] = stencil; + } + } + else + { + CCTK_WARN(0,"ApplyBndScalar: Malloc sw failed"); + } + } - /* allocate arrays */ - doBC = (int *)malloc((2*dim)*sizeof(int)); - dstag = (int *)malloc(dim*sizeof(int)); - lssh = (int *)malloc(dim*sizeof(int)); - - /* get the directional staggering of the group */ - gi = CCTK_GroupIndexFromVarI(first_var); - ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi); - if (ierr<0) CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "Cannot get stagger information for group '%s' (GF '$s')\n", - CCTK_GroupName(gi), CCTK_VarName(first_var)); - + /* allocate arrays */ + dstag = (int *)malloc(dim*sizeof(int)); + lssh = (int *)malloc(dim*sizeof(int)); - /* get the current timelevel */ - timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1; - /* if (timelevel > 0) - timelevel--;*/ + /* get the directional staggering of the group */ + gi = CCTK_GroupIndexFromVarI(first_var); + ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi); - for (vi=first_var; vi0) { - for (idim=0;idimGFSym[vi][idim*2]==GFSYM_NOSYM)|| - (sGHex->GFSym[vi][idim*2]==GFSYM_UNSET)) && - GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]); - doBC[idim*2+1] = (((sGHex->GFSym[vi][idim*2+1]==GFSYM_NOSYM)|| - (sGHex->GFSym[vi][idim*2+1]==GFSYM_UNSET)) && - GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]); - lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)]; - } + dirloop1 = 2*(dir-1)+1; + dirloop2 = dirloop1+1; + } + else if (dir<0) + { + dirloop1 = 2*(-dir-1); + dirloop2 = dirloop1+1; } else { - for (idim=0;idimcctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]); - doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]); - lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)]; + lssh[idim] = (int)GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)]; } + + for (count=dirloop1;countGFSym[vi][count]==GFSYM_NOSYM) || + (sGHex->GFSym[vi][count]==GFSYM_UNSET) ) && + GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? + count : -1; + } + else + { + doBC = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? + count : -1 ; + } + + if (doBC >=0) + { + switch (dim) + { + case 1: + berr = BndApplyFlat1Di(GH, + doBC, + lssh, + sw, + GH->data[vi][timelevel]); + break; + case 2: + berr = BndApplyFlat2Di(GH, + doBC, + lssh, + sw, + GH->data[vi][timelevel]); + break; + case 3: + berr = BndApplyFlat3Di(GH, + doBC, + lssh, + sw, + GH->data[vi][timelevel]); + break; + default : + berr = -1; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndFlat: No BC for dim>3: grid variable %s", + CCTK_VarName(vi)); + } + } + } + /* retval = (berr>-1) ? 0 : -1;*/ } + + free(dstag); + free(lssh); - switch (dim) + if (!stencil_array) { - case 1: berr = BndApplyFlat1Di(GH, - dim, - doBC, - lssh, - stencil, - GH->data[vi][timelevel]); break; - case 2: berr = BndApplyFlat2Di(GH, - dim, - doBC, - lssh, - stencil, - GH->data[vi][timelevel]); break; - case 3: berr = BndApplyFlat3Di(GH, - dim, - doBC, - lssh, - stencil, - GH->data[vi][timelevel]); break; - default : berr = -1; - CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "ApplyBndFlat: No BC for dim>3: grid variable '%s'", - CCTK_VarName(vi)); + free(sw); } - berr = (berr>-1) ? 0 : -1; } - - free(dstag); - free(doBC); - free(lssh); - - return(berr); + + return retval; } + /*@@ @routine BndApplyFlat3Di @date Thu Mar 2 11:14:08 2000 @@ -364,8 +610,7 @@ int ApplyBndFlat(cGH *GH, @@*/ int BndApplyFlat3Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, int *stencil, CCTK_REAL *var) @@ -373,8 +618,11 @@ int BndApplyFlat3Di(cGH *GH, int i,j,k,sw; /* lower x */ - if (doBC[0] == 1) + if (doBC == 0) { +#ifdef DEBUG_BOUNDARY + printf("Applying flat boundary (lower x)\n"); +#endif for (k=0;kcctk_lsh[2];k++) { for (j=0;jcctk_lsh[1];j++) @@ -389,8 +637,11 @@ int BndApplyFlat3Di(cGH *GH, } /* upper x */ - if (doBC[1] == 1) + if (doBC == 1) { +#ifdef DEBUG_BOUNDARY + printf("Applying flat boundary (upper x)\n"); +#endif for (k=0;kcctk_lsh[2];k++) { for (j=0;jcctk_lsh[1];j++) @@ -405,8 +656,11 @@ int BndApplyFlat3Di(cGH *GH, } /*lower y */ - if (doBC[2] == 1) + if (doBC == 2) { +#ifdef DEBUG_BOUNDARY + printf("Applying flat boundary (lower y)\n"); +#endif for (k=0;kcctk_lsh[2];k++) { for (i=0;icctk_lsh[0];i++) @@ -421,8 +675,11 @@ int BndApplyFlat3Di(cGH *GH, } /* upper y */ - if (doBC[3] == 1) + if (doBC == 3) { +#ifdef DEBUG_BOUNDARY + printf("Applying flat boundary (upper y)\n"); +#endif for (k=0;kcctk_lsh[2];k++) { for (i=0;icctk_lsh[0];i++) @@ -437,8 +694,11 @@ int BndApplyFlat3Di(cGH *GH, } /* lower z */ - if (doBC[4] == 1) + if (doBC == 4) { +#ifdef DEBUG_BOUNDARY + printf("Applying flat boundary (lower z)\n"); +#endif for (j=0;jcctk_lsh[1];j++) { for (i=0;icctk_lsh[0];i++) @@ -453,8 +713,11 @@ int BndApplyFlat3Di(cGH *GH, } /* upper z */ - if (doBC[5] == 1) + if (doBC == 5) { +#ifdef DEBUG_BOUNDARY + printf("Applying flat boundary (upper z)\n"); +#endif for (j=0;jcctk_lsh[1];j++) { for (i=0;icctk_lsh[0];i++) @@ -487,8 +750,7 @@ int BndApplyFlat3Di(cGH *GH, @@*/ int BndApplyFlat2Di(cGH *GH, - int gdim, - int *doBC, + int doBC, int *lssh, int *stencil, CCTK_REAL *var) @@ -496,8 +758,11 @@ int BndApplyFlat2Di(cGH *GH, int i,j,sw; /* lower x */ - if (doBC[0] == 1) + if (doBC == 0) { +#ifdef DEBUG_BOUNDARY + printf("Applying flat boundary (lower x)\n"); +#endif for (j=0;jcctk_lsh[1];j++) { for (sw=0;swcctk_lsh[1];j++) { for (sw=0;swcctk_lsh[0];i++) { for (sw=0;swcctk_lsh[0];i++) { for (sw=0;swcctk_lsh[0]-sw-1)] diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c index dcc853f..68ab210 100644 --- a/src/RadiationBoundary.c +++ b/src/RadiationBoundary.c @@ -8,7 +8,7 @@ @version $Header$ @@*/ -/*#define DEBUG_BOUND*/ +/*#define DEBUG_BOUNDARY*/ #include #include @@ -34,10 +34,9 @@ CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundaryWrappers_c) /* Internal routine prototypes */ static int BndApplyRadiative3Di(cGH *GH, - int gdim, - int *sw, - int *doBC, + int doBC, int *lssh, + int *stencil_size, CCTK_REAL *dxyz, CCTK_REAL dt, CCTK_REAL *var_n, @@ -50,12 +49,14 @@ static int BndApplyRadiative3Di(cGH *GH, CCTK_REAL v0) ; static int ApplyBndRadiative(cGH *GH, - int *sw, + int stencil, + int *stencil_array, CCTK_REAL var0, CCTK_REAL v0, int var_current, int var_previous, - int num_vars); + int num_vars, + int dir); /* Local variables */ @@ -64,6 +65,46 @@ static int ApplyBndRadiative(cGH *GH, ******************** External Routines ************************ ********************************************************************/ +/*@@ + @routine BndRadiativeDirGI + @date Sat Jan 20 2001 + @author Gabrielle Allen + @desc + Aply radiative BC's by group index in given direction + @enddesc + @calls + @calledby + @history + @endhistory + +@@*/ + +int BndRadiativeDirGI(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL speed, + int gi, + int gi_p) +{ + int numvars, first_vi, first_vi_p; + + first_vi = CCTK_FirstVarIndexI(gi); + first_vi_p = CCTK_FirstVarIndexI(gi_p); + numvars = CCTK_NumVarsInGroupI(gi); + + return ApplyBndRadiative(GH,stencil_size,NULL,var0, + speed,first_vi,first_vi_p,numvars,dir); + +} + +void CCTK_FCALL CCTK_FNAME(BndRadiativeDirGI) + (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, + CCTK_REAL *speed, int *gi, int *gi_p) +{ + *ierr = BndRadiativeDirGI(GH, *stencil_size, *dir, *var0, *speed, *gi, *gi_p); +} + /*@@ @routine BndRadiativeGI @date Tue Jul 18 18:04:07 2000 @@ -79,9 +120,9 @@ static int ApplyBndRadiative(cGH *GH, @@*/ int BndRadiativeGI(cGH *GH, - int *sw, + int *stencil_array, CCTK_REAL var0, - CCTK_REAL v0, + CCTK_REAL speed, int gi, int gi_p) { @@ -91,17 +132,72 @@ int BndRadiativeGI(cGH *GH, first_vi_p = CCTK_FirstVarIndexI(gi_p); numvars = CCTK_NumVarsInGroupI(gi); - return ApplyBndRadiative(GH,sw,var0,v0,first_vi,first_vi_p,numvars); + return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,first_vi,first_vi_p,numvars,0); } void CCTK_FCALL CCTK_FNAME(BndRadiativeGI) - (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, - CCTK_REAL *v0, int *gi, int *gi_p) + (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, + CCTK_REAL *speed, int *gi, int *gi_p) { - *ierr = BndRadiativeGI(GH, sw, *var0, *v0, *gi, *gi_p); + *ierr = BndRadiativeGI(GH, stencil_array, *var0, *speed, *gi, *gi_p); } + +/* ===================================================================== */ + + +/*@@ + @routine BndRadiativeDirGN + @date Sat Jan 20 2001 + @author Gabrielle Allen + @desc + Apply radiative BC's by group name in given direction + @enddesc + @calls + @calledby + @history + @endhistory + +@@*/ + +int BndRadiativeDirGN(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL speed, + const char *gn, + const char *gn_p) +{ + int retval; + int gi = CCTK_GroupIndex(gn); + int gi_p = CCTK_GroupIndex(gn_p); + + if ((gi>-1)&&(gi_p>-1)) + retval = BndRadiativeDirGI(GH, stencil_size, dir, var0, speed, gi, gi_p); + else + { + CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", + "BndRadiativeDirGN: Grid variable groups %s or %s not found", + gn, gn_p); + retval = -1; + } + + return retval; +} + +void CCTK_FCALL CCTK_FNAME(BndRadiativeDirGN) + (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, + CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS) +{ + TWO_FORTSTRINGS_CREATE(gn, gn_p) + *ierr = BndRadiativeDirGN(GH, *stencil_size, *dir, *var0, *speed, gn, gn_p); + free(gn); + free(gn_p); + return; +} + + /*@@ @routine BndRadiativeGI @date Tue Jul 18 18:04:07 2000 @@ -116,9 +212,9 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeGI) @@*/ int BndRadiativeGN(cGH *GH, - int *sw, + int *stencil_array, CCTK_REAL var0, - CCTK_REAL v0, + CCTK_REAL speed, const char *gn, const char *gn_p) { @@ -127,7 +223,7 @@ int BndRadiativeGN(cGH *GH, int gi_p = CCTK_GroupIndex(gn_p); if ((gi>-1)&&(gi_p>-1)) - retval = BndRadiativeGI(GH, sw, var0, v0, gi, gi_p); + retval = BndRadiativeGI(GH, stencil_array, var0, speed, gi, gi_p); else { CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", @@ -140,18 +236,58 @@ int BndRadiativeGN(cGH *GH, } void CCTK_FCALL CCTK_FNAME(BndRadiativeGN) - (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, - CCTK_REAL *v0, TWO_FORTSTRINGS_ARGS) + (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, + CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS) { TWO_FORTSTRINGS_CREATE(gn, gn_p) - *ierr = BndRadiativeGN(GH, sw, *var0, *v0, gn, gn_p); + *ierr = BndRadiativeGN(GH, stencil_array, *var0, *speed, gn, gn_p); free(gn); free(gn_p); return; } + +/* ===================================================================== */ + + /*@@ - @routine BndRadiativeGI + @routine BndRadiativeDirVI + @date Sat Jan 20 2001 + @author Gabrielle Allen + @desc + Apply radiative BC's by variable index in given direction + @enddesc + @calls + @calledby + @history + @endhistory + +@@*/ + +int BndRadiativeDirVI(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL speed, + int vi, + int vi_p) +{ + return ApplyBndRadiative(GH,stencil_size,NULL,var0,speed,vi,vi_p,1,dir); +} + +void CCTK_FCALL CCTK_FNAME(BndRadiativeDirVI) + (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, + CCTK_REAL *speed, int *vi, int *vi_p) +{ + *ierr = BndRadiativeDirVI(GH, *stencil_size, *dir, *var0, + *speed, *vi, *vi_p); + return; +} + + + +/*@@ + @routine BndRadiativeVI @date Tue Jul 18 18:04:07 2000 @author Gerd Lanfermann @desc @@ -165,23 +301,81 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeGN) @@*/ int BndRadiativeVI(cGH *GH, - int *sw, + int *stencil_array, CCTK_REAL var0, - CCTK_REAL v0, + CCTK_REAL speed, int vi, int vi_p) { - return ApplyBndRadiative(GH,sw,var0,v0,vi,vi_p,1); + return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,vi,vi_p,1,0); } void CCTK_FCALL CCTK_FNAME(BndRadiativeVI) - (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, - CCTK_REAL *v0, int *vi, int *vi_p) + (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, + CCTK_REAL *speed, int *vi, int *vi_p) { - *ierr = BndRadiativeVI(GH, sw, *var0, *v0, *vi, *vi_p); + *ierr = BndRadiativeVI(GH, stencil_array, *var0, *speed, *vi, *vi_p); return; } + +/* ======================================================================= */ + + +/*@@ + @routine BndRadiativeDirGI + @date Sat Jan 20 2001 + @author Gabrielle Allen + @desc + Apply radiative BC's by variable name in given direction + @enddesc + @calls + @calledby + @history + @endhistory + +@@*/ + +int BndRadiativeDirVN(cGH *GH, + int stencil_size, + int dir, + CCTK_REAL var0, + CCTK_REAL speed, + const char *vn, + const char *vn_p) +{ + int retval; + int vi = CCTK_VarIndex(vn); + int vi_p = CCTK_VarIndex(vn_p); + + if (vi>=0 && vi_p>=0) + { + retval = BndRadiativeDirVI(GH, stencil_size, dir, var0, speed, vi, vi_p); + } + else + { + CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", + "BndRadiativeDirVN: Grid variable %s or %s not found", + vn, vn_p); + retval = -1; + } + + return retval; + +} + +void CCTK_FCALL CCTK_FNAME(BndRadiativeDirVN) + (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, + CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS) +{ + TWO_FORTSTRINGS_CREATE(vn,vn_p) + *ierr = BndRadiativeDirVN(GH, *stencil_size, *dir, *var0, *speed, vn, vn_p); + free(vn); + free(vn_p); + return; +} + + /*@@ @routine BndRadiativeGI @date Tue Jul 18 18:04:07 2000 @@ -197,9 +391,9 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVI) @@*/ int BndRadiativeVN(cGH *GH, - int *sw, + int *stencil_array, CCTK_REAL var0, - CCTK_REAL v0, + CCTK_REAL speed, const char *vn, const char *vn_p) { @@ -209,7 +403,7 @@ int BndRadiativeVN(cGH *GH, if (vi>=0 && vi_p>=0) { - retval = BndRadiativeVI(GH, sw, var0, v0, vi, vi_p); + retval = BndRadiativeVI(GH, stencil_array, var0, speed, vi, vi_p); } else { @@ -224,11 +418,11 @@ int BndRadiativeVN(cGH *GH, } void CCTK_FCALL CCTK_FNAME(BndRadiativeVN) - (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, - CCTK_REAL *v0, TWO_FORTSTRINGS_ARGS) + (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, + CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS) { TWO_FORTSTRINGS_CREATE(vn,vn_p) - *ierr = BndRadiativeVN(GH, sw, *var0, *v0, vn, vn_p); + *ierr = BndRadiativeVN(GH, stencil_array, *var0, *speed, vn, vn_p); free(vn); free(vn_p); return; @@ -240,7 +434,7 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVN) ********************* Local Routines ************************* ********************************************************************/ /*@@ - @routine BndRadiativeGI + @routine ApplyBndRadiative @date Tue Jul 18 18:04:07 2000 @author Gerd Lanfermann @desc @@ -254,44 +448,117 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVN) @@*/ static int ApplyBndRadiative(cGH *GH, - int *sw, + int stencil, + int *stencil_array, CCTK_REAL var0, - CCTK_REAL v0, + CCTK_REAL speed, int var_current, int var_previous, - int num_vars) + int num_vars, + int dir) { - SymmetryGHex *sGHex; int xi=-1; int yi=-1; int zi=-1; int ri=-1; + int count; int gi; int vi1, vi2; int symmetry_handle; int idim, dim; int time, time_p; int berr,ierr; - int *doBC, *dstag, *lssh; + int dirloop1,dirloop2; + int doBC; + int retval; + int *sw; + int *dstag; + int *lssh; CCTK_REAL dx[3], dt; + SymmetryGHex *sGHex; + +#ifdef DEBUG_BOUNDARY + printf("Input arguments for ApplyBndRadiative:\n"); + printf("GH = %x\n",GH); + printf("stencil = %d\n",stencil); + printf("stencil_array = %x\n",stencil_array); + printf("var0 = %f\n",var0); + printf("speed = %f\n",speed); + printf("var_current = %d\n",var_current); + printf("var_previous = %d\n",var_previous); + printf("num_vars = %d\n",num_vars); + printf("dir = %d\n",dir); + printf("-----------------------------------\n"); +#endif + + /* Get group dimension */ + dim = CCTK_GroupDimFromVarI(var_current); - /* See if we have a symmetry array */ - symmetry_handle = CCTK_GHExtensionHandle("Symmetry"); - if (symmetry_handle < 0) + if (dim<3 || dim>3) { - sGHex = NULL; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndRadiative: Invalid grid function dimension %d", + dim); + retval=-1; } else { - sGHex = (SymmetryGHex*)GH->extensions[symmetry_handle]; - } + /* See if we have a symmetry array */ + symmetry_handle = CCTK_GHExtensionHandle("Symmetry"); + if (symmetry_handle < 0) + { + sGHex = NULL; + } + else + { + sGHex = (SymmetryGHex*)GH->extensions[symmetry_handle]; + } - /* Get group dimension */ - dim = CCTK_GroupDimFromVarI(var_current); + /* Set up stencil width array if needed */ + if (stencil_array && dir == 0) + { + sw = stencil_array; + } + else if (dir==0) + { + CCTK_WARN(0,"ApplyBndRadiative: Direction 0 invalid for directional BC"); + } + else if (dir>dim||-dir>dim) + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndRadiative: Direction %d greater than dimension %d", + dir,dim); + } + else + { - /* Radiative boundaries need the underlying Cartesian coordinates */ - switch (dim) - { + int i; + /* Should be single direction */ + sw = (int *)malloc(dim*sizeof(int)); + if (sw) + { + for (i=0;i 0) + { + sw[dir-1] = stencil; + } + else + { + sw[-dir-1] = stencil; + } + } + else + { + CCTK_WARN(0,"ApplyBndRadiative: Malloc sw failed"); + } + } + + /* Radiative boundaries need the underlying Cartesian coordinates */ + switch (dim) + { case 1: xi = CCTK_CoordIndex(-1,"x","cart1d"); break; @@ -306,130 +573,148 @@ static int ApplyBndRadiative(cGH *GH, ri = CCTK_CoordIndex(-1,"r","spher3d"); break; default: - CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "Radiative boundary only for dim=3: grid varible '%s'", + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndRadiative: Radiative boundary only for dim=3: grid variable '%s'", CCTK_VarName(var_current)); break; - } + } + /* Allocate temporary arrays */ + dstag = (int *)malloc(dim*sizeof(int)); + lssh = (int *)malloc(dim*sizeof(int)); - /* Allocate temporary arrays */ - doBC = (int *)malloc((2*dim)*sizeof(int)); - dstag = (int *)malloc(dim*sizeof(int)); - lssh = (int *)malloc(dim*sizeof(int)); + /* get the directional staggering of the group */ + gi = CCTK_GroupIndexFromVarI(var_current); + ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi); - /* get the directional staggering of the group */ - gi = CCTK_GroupIndexFromVarI(var_current); - ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi); + + /* Use next time level, if available */ + time = CCTK_NumTimeLevelsFromVarI(var_current) - 1; + if (time < 0) + { + time = 0; + } + /* Use current time level, if available */ + time_p = CCTK_NumTimeLevelsFromVarI(var_previous) - 2; + if (time_p < 0) + { + time_p = 0; + } + /* Treat all boundaries or just one? */ + if (dir>0) + { + dirloop1 = 2*(dir-1)+1; + dirloop2 = dirloop1+1; + } + else if (dir<0) + { + dirloop1 = 2*(-dir-1); + dirloop2 = dirloop1+1; + } + else + { + dirloop1 = 0; + dirloop2 = 2*dim; + } + + for (vi1=var_current; vi1GFSym[vi1][idim*2]==GFSYM_NOSYM)|| - (sGHex->GFSym[vi1][idim*2]==GFSYM_UNSET)) && - GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]); - doBC[idim*2+1] = (((sGHex->GFSym[vi1][idim*2+1]==GFSYM_NOSYM)|| - (sGHex->GFSym[vi1][idim*2+1]==GFSYM_UNSET)) && - GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]); /* FIXME: THIS LINE FAILS WAVE TESTSUITES */ /*lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];*/ lssh[idim] = GH->cctk_lsh[idim]; } - } - else - { - for (idim=0;idimcctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]); - doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]); - lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)]; + if (sGHex) + { + doBC = ( + ( (sGHex->GFSym[vi1][count]==GFSYM_NOSYM) || + (sGHex->GFSym[vi1][count]==GFSYM_UNSET) ) && + GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? + count : -1; + } + else + { + doBC = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? + count : -1 ; + } + + if (doBC >=0) + { + switch (dim) + { + case 1: + berr = -1; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndRadiative: Radiative boundary conditions" + "not implemented for dimension = 1, grid variable %s", + CCTK_VarName(vi1)); + break; + case 2: + berr = -1; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndRadiative: Radiative boundary conditions" + "not implemented for dimension = 2, grid variable %s", + CCTK_VarName(vi1)); + break; + case 3: + /* According to the Cactus spec, the true dx and dt values for a + grid are calculated as follows: */ + dx[0] = GH->cctk_delta_space[0] / GH->cctk_levfac[0]; + dx[1] = GH->cctk_delta_space[1] / GH->cctk_levfac[1]; + dx[2] = GH->cctk_delta_space[2] / GH->cctk_levfac[2]; + dt = GH->cctk_delta_time; + + berr = BndApplyRadiative3Di + (GH, + doBC, /* apply BC on which face */ + lssh, /* local shape, respects staggering */ + sw, /* stencil width array */ + dx, /* dx,dy,dz */ + dt, /* dt */ + GH->data[vi1][time ], /* start of data array GF[] */ + GH->data[vi2][time_p], /* start of prev data array */ + GH->data[xi][0], /* pointer to the X field */ + GH->data[yi][0], /* pointer to the Y field */ + GH->data[zi][0], /* pointer to the Z field */ + GH->data[ri][0], /* pointer to the R field */ + var0, /* asymptotic limit */ + speed); /* wave speed */ + break; + default : + berr=-1; + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndRadiative: No Radiative BC for dim>3, grid variable %s", + CCTK_VarName(vi1)); + } + } } + retval=(berr<0)?-1:0; } + free(dstag); + free(lssh); - switch (dim) + if (!stencil_array) { - case 1: - berr = -1; - CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "ApplyBndRadiative: Radiative boundary conditions" - "not implemented for dimension = 1, grid variable '%s'", - CCTK_VarName(vi1)); - break; - case 2: - berr = -1; - CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "ApplyBndRadiative: Radiative boundary conditions" - "not implemented for dimension = 2, grid variable '%s'", - CCTK_VarName(vi1)); - break; - case 3: - /* According to the Cactus spec, the true dx and dt values for a - grid are calculated as follows: */ - dx[0] = GH->cctk_delta_space[0] / GH->cctk_levfac[0]; - dx[1] = GH->cctk_delta_space[1] / GH->cctk_levfac[1]; - dx[2] = GH->cctk_delta_space[2] / GH->cctk_levfac[2]; - dt = GH->cctk_delta_time; - - berr = BndApplyRadiative3Di - (GH, - dim, - sw, /* stencil_width */ - doBC, /* flag applying low/up BCs */ - lssh, /* local shape, respects staggering */ - dx, /* dx,dy,dz */ - dt, /* dt */ - GH->data[vi1][time ], /* start of data array GF[] */ - GH->data[vi2][time_p], /* start of prev data array */ - GH->data[xi][0], /* pointer to the X field */ - GH->data[yi][0], /* pointer to the Y field */ - GH->data[zi][0], /* pointer to the Z field */ - GH->data[ri][0], /* pointer to the R field */ - var0, /* asymptotic limit */ - v0); /* wave speed */ - break; - default : - berr=-1; - CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "No Radiative BC for dim>3, grid variable '%s'", - CCTK_VarName(vi1)); + free(sw); } - berr=(berr<0)?-1:0; } - - free(lssh); - free(dstag); - free(doBC); - - return ierr; -} - + return retval; +} /*@@ @routine BndApplyRadiative3Di @@ -515,10 +800,9 @@ static int ApplyBndRadiative(cGH *GH, @@*/ static int BndApplyRadiative3Di(cGH *GH, - int gdim, - int *sw, - int *doBC, + int doBC, int *lssh, + int *sw, CCTK_REAL *dxyz, CCTK_REAL dt, CCTK_REAL *var_n, @@ -528,7 +812,7 @@ static int BndApplyRadiative3Di(cGH *GH, CCTK_REAL *z, CCTK_REAL *r, CCTK_REAL var0, - CCTK_REAL v0) + CCTK_REAL speed) { DECLARE_CCTK_PARAMETERS @@ -588,7 +872,7 @@ static int BndApplyRadiative3Di(cGH *GH, /* Find Courant parameters. */ - dtv = v0*dt; + dtv = speed*dt; dtvh = half*dtv; dtvvar0 = dtv*var0; @@ -610,11 +894,17 @@ static int BndApplyRadiative3Di(cGH *GH, /* Lower x-bound: x(2,:,:) --> x[xgp2] */ - if (doBC[0]==1) { - for (k=0;k=0;i--) { - + if (doBC == 0) + { +#ifdef DEBUG_BOUNDARY + printf("Applying radiative boundary (lower x)\n"); +#endif + for (k=0;k=0;i--) + { xgp0 = CCTK_GFINDEX3D(GH,i ,j,k); xgp1 = CCTK_GFINDEX3D(GH,i+1,j,k); xgp2 = CCTK_GFINDEX3D(GH,i+2,j,k); @@ -633,11 +923,14 @@ static int BndApplyRadiative3Di(cGH *GH, ri0 = 1.0/r0; ri1 = 1.0/r1; - if (fac==0) { + if (fac==0) + { H = 0; - } else { + } + else + { H = 0; @@ -670,10 +963,17 @@ static int BndApplyRadiative3Di(cGH *GH, /* Upper x-bound: x(nx,:,:) --> xgp[xgp0] */ - if (doBC[1]==1){ - for (k=0;k=0;k--) { + if (doBC==4) + { +#ifdef DEBUG_BOUNDARY + printf("Applying radiative boundary (lower z)\n"); +#endif + for (j=0;j=0;k--) + { zgp0 = CCTK_GFINDEX3D(GH,i,j,k ); zgp1 = CCTK_GFINDEX3D(GH,i,j,k+1); @@ -883,11 +1202,14 @@ static int BndApplyRadiative3Di(cGH *GH, ri0 = 1.0/r0; ri1 = 1.0/r1; - if (fac==0) { + if (fac==0) + { H = 0; - } else { + } + else + { H = 0; @@ -920,10 +1242,17 @@ static int BndApplyRadiative3Di(cGH *GH, /* Upper z-bound */ - if (doBC[5] == 1) { - for (j=0;j 0) { sw[dir-1] = stencil; - coord = dir-1; } else { sw[-dir-1] = stencil; - coord = -dir-1; } } else @@ -606,12 +604,12 @@ static int ApplyBndScalar(cGH *GH, doBC = ( ( (sGHex->GFSym[vi][count]==GFSYM_NOSYM) || (sGHex->GFSym[vi][count]==GFSYM_UNSET) ) && - GH->cctk_lsh[coord]>1 && GH->cctk_bbox[count]) ? + GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? count : -1; } else { - doBC = (GH->cctk_lsh[coord]>1 && GH->cctk_bbox[count]) ? + doBC = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ? count : -1 ; } @@ -645,8 +643,8 @@ static int ApplyBndScalar(cGH *GH, break; default : berr = -1; - CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", - "No scalar boundary for dim>3: grid variable '%s'", + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "ApplyBndScalar: No scalar boundary for dim>3: grid variable '%s'", CCTK_VarName(vi)); } } -- cgit v1.2.3