diff options
Diffstat (limited to 'src/Compute_Norm.c')
| -rw-r--r-- | src/Compute_Norm.c | 184 |
1 files changed, 184 insertions, 0 deletions
diff --git a/src/Compute_Norm.c b/src/Compute_Norm.c new file mode 100644 index 0000000..083f498 --- /dev/null +++ b/src/Compute_Norm.c @@ -0,0 +1,184 @@ +/* XXX cleanup the include files ! */ + +#include <assert.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + + +static const char * rcsid = "$Header$"; + +CCTK_FILEVERSION(Norms_Compute_Norm_c); + + +void Norms_Compute_Norms (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int vind; /* variable index */ + CCTK_REAL *data; + int index,indexp, i,j,k; + int ierr; + int istart, iend, jstart, jend, kstart, kend; + int npoints; + CCTK_REAL hd, dx,dy,dz; + CCTK_REAL deriv; + int target_proc,reduction_handle; + CCTK_REAL result; + + if (verbose>1) + CCTK_INFO("computing norms"); + + *norm=0.; + result=0.; + + dx = CCTK_DELTA_SPACE(0); + dy = CCTK_DELTA_SPACE(1); + dz = CCTK_DELTA_SPACE(2); + + /* XXX ghost size */ + istart = 1; + jstart = 1; + kstart = 1; + + iend = cctk_lsh[0]-1; + jend = cctk_lsh[1]-1; + kend = cctk_lsh[2]-1; + + for (k=kstart;k<kend;k++) { + for (j=jstart;j<jend;j++) { + for (i=istart;i<iend;i++) { + index=CCTK_GFINDEX3D(cctkGH,i,j,k); + diff_term[index]=0.; + } + } + } + + target_proc = 0; + reduction_handle = CCTK_ReductionArrayHandle("norm2"); + if (reduction_handle < 0) + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "can't get sum-reduction handle!"); + + if (verbose>2) + fprintf(stderr," nr1stvars %d nr2ndvars %d\n", + *nr1stvars,*nr2ndvars); + + /* L2-norm of 1st order variables */ + for (vind=0; vind< *nr1stvars; vind++) { + ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1, + CCTK_VARIABLE_REAL,&result,1,varindices_1st[vind]); + if (ierr<0) + CCTK_WARN(1,"reduction failed"); + *norm=*norm+result*result; + } + + if (verbose>2) + fprintf(stderr," after 1st order vars norm=%16.8g\n",*norm); + + /* L2-norm of 2nd order variables */ + for (vind=0; vind< *nr2ndvars; vind++) { + ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1, + CCTK_VARIABLE_REAL,&result,1,varindices_2nd[vind]); + if (ierr<0) + CCTK_WARN(1,"reduction failed"); + *norm=*norm+result*result; + } + + if (verbose>2) + fprintf(stderr," after 2nd order vars norm=%16.8g\n",*norm); + + /* L2-norm of D_{+i} terms */ + /* d_x */ + for (vind=0; vind< *nr2ndvars; vind++) { + data=CCTK_VarDataPtrI(cctkGH,0,varindices_2nd[vind]); + for (k=kstart;k<kend;k++) { + for (j=jstart;j<jend;j++) { + for (i=istart;i<iend;i++) { + index =CCTK_GFINDEX3D(cctkGH,i,j,k); + indexp=CCTK_GFINDEX3D(cctkGH,i+1,j,k); + diff_term[index]=(data[indexp]-data[index])/dx; + } + } + } + ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1, + CCTK_VARIABLE_REAL,&result,1, + CCTK_VarIndex("norms::diff_term")); + if (ierr<0) + CCTK_WARN(1,"reduction failed"); + *norm=*norm+result*result; + } + + if (verbose>2) + fprintf(stderr," after dx 2nd order terms norm=%16.8g\n",*norm); + + /* d_y */ + for (vind=0; vind< *nr2ndvars; vind++) { + data=CCTK_VarDataPtrI(cctkGH,0,varindices_2nd[vind]); + for (k=kstart;k<kend;k++) { + for (j=jstart;j<jend;j++) { + for (i=istart;i<iend;i++) { + index =CCTK_GFINDEX3D(cctkGH,i,j,k); + indexp=CCTK_GFINDEX3D(cctkGH,i,j+1,k); + diff_term[index]=(data[indexp]-data[index])/dy; + } + } + } + ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1, + CCTK_VARIABLE_REAL,&result,1, + CCTK_VarIndex("norms::diff_term")); + if (ierr<0) + CCTK_WARN(1,"reduction failed"); + *norm=*norm+result*result; + } + + if (verbose>2) + fprintf(stderr," after dy 2nd order terms norm=%16.8g\n",*norm); + + /* d_z */ + for (vind=0; vind< *nr2ndvars; vind++) { + data=CCTK_VarDataPtrI(cctkGH,0,varindices_2nd[vind]); + for (k=kstart;k<kend;k++) { + for (j=jstart;j<jend;j++) { + for (i=istart;i<iend;i++) { + index =CCTK_GFINDEX3D(cctkGH,i,j,k); + indexp=CCTK_GFINDEX3D(cctkGH,i,j,k+1); + diff_term[index]=(data[indexp]-data[index])/dz; + } + } + } + ierr=CCTK_Reduce(cctkGH,target_proc,reduction_handle,1, + CCTK_VARIABLE_REAL,&result,1, + CCTK_VarIndex("norms::diff_term")); + if (ierr<0) + CCTK_WARN(1,"reduction failed"); + *norm=*norm+result*result; + } + + if (verbose>2) + fprintf(stderr," after dz 2nd order terms norm=%16.8g\n",*norm); + + if (verbose>2) + fprintf(stderr," after 2nd order vars (incl. D_{+i}): norm=%16.8g\n", + *norm); + if (verbose>2) + fprintf(stderr,"final norm^2=%16.8g on this CPU\n",*norm); + + *norm=sqrt(*norm); + + /* multiply with h^d - see discussion above Eq. (12) in + * gr-qc/XXX */ + hd=dx*dy*dz; + *norm=*norm*hd; + + CCTK_VInfo(CCTK_THORNSTRING,"norm= %16.8g",*norm); + +} |