diff options
| author | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-04-27 11:17:59 -0500 |
|---|---|---|
| committer | Barry Wardell <barry.wardell@gmail.com> | 2011-12-14 18:17:06 +0000 |
| commit | 8d73536d782c9b1c419b89170767ee7ba6b5820f (patch) | |
| tree | 46b55f7403c05387f2e16c519e698c87de21607e /Carpet | |
| parent | 580147ad1465531a2951d2e224b64d9bd3dcc7d8 (diff) | |
| parent | 31aa0a30ac88e00db2e149c33e4e3c54237ce9d3 (diff) | |
CarpetInterp: Merge
Diffstat (limited to 'Carpet')
59 files changed, 1537 insertions, 1360 deletions
diff --git a/Carpet/Carpet/options/README b/Carpet/Carpet/options/README deleted file mode 100644 index 982309470..000000000 --- a/Carpet/Carpet/options/README +++ /dev/null @@ -1,7 +0,0 @@ -This directory contains sample options files for various architectures -and compilers. The syntax of the file names is - - carpet-MACHINE-COMPILER - -so that e. g. carpet-lilypong-lf95 is a configuration for the machine -Lilypond, using the LF95 compiler. diff --git a/Carpet/Carpet/options/carpet-because b/Carpet/Carpet/options/carpet-because deleted file mode 100644 index 2aa561626..000000000 --- a/Carpet/Carpet/options/carpet-because +++ /dev/null @@ -1,27 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -F77 = g77 -F90 = g77 - -CXXFLAGS = -DCARPET_REAL -ftemplate-depth-30 - -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -F90_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double -F77_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = -Wall -F77_WARN_FLAGS = -Wall - -SYS_INC_DIRS = /usr/include/hdf /home/eschnett/FlexIO -LIBDIRS = /home/eschnett/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt g2c m - -HDF5 = yes -MPI = MPICH - -DEBUG = yes - -# PTHREADS -PTHREADS = no diff --git a/Carpet/Carpet/options/carpet-harpo-gcc b/Carpet/Carpet/options/carpet-harpo-gcc deleted file mode 100644 index a760e61aa..000000000 --- a/Carpet/Carpet/options/carpet-harpo-gcc +++ /dev/null @@ -1,36 +0,0 @@ -CC = gcc -CXX = g++ -F77 = g77 -F90 = g77 -LD = g++ - -#CC = /u1/eschnett/bin/gcc -#CXX = /u1/eschnett/bin/g++ -#F77 = /u1/eschnett/bin/g77 -#F90 = /u1/eschnett/bin/g77 -#LD = /u1/eschnett/bin/g++ - -CFLAGS = -CXXFLAGS = -ftemplate-depth-30 -DMPIPP_H -FFLAGS = -F90FLAGS = - -C_OPTIMISE_FLAGS = -O3 -funroll-loops -CXX_OPTIMISE_FLAGS = -O3 -funroll-loops -F90_OPTIMISE_FLAGS = -O3 -funroll-loops -malign-double -F77_OPTIMISE_FLAGS = -O3 -funroll-loops -malign-double - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -Winline -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -Winline -F90_WARN_FLAGS = -Wall -F77_WARN_FLAGS = -Wall - -SYS_INC_DIRS = /u1/eschnett/include /u1/eschnett/scratch/FlexIO -LIBDIRS = /u1/eschnett/lib /u1/eschnett/scratch/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z - -HDF5 = yes -MPI = NATIVE - -OPTIMISE = no -DEBUG = yes diff --git a/Carpet/Carpet/options/carpet-harpo-sgi b/Carpet/Carpet/options/carpet-harpo-sgi deleted file mode 100644 index 7844edeb3..000000000 --- a/Carpet/Carpet/options/carpet-harpo-sgi +++ /dev/null @@ -1,29 +0,0 @@ -CPP = /lib/cpp -DMPIPP_H - -CFLAGS = -mips4 -r12000 -64 -LANG:restrict=ON -CXXFLAGS = -mips4 -r12000 -64 -LANG:restrict=ON -LANG:std -LANG:vla=ON -no_auto_include -ptused -DMPIPP_H -F90FLAGS = -mips4 -r12000 -64 -F77FLAGS = -mips4 -r12000 -64 -LDFLAGS = -LANG:std -LANG:vla=ON -64 -mips4 -r12000 -Wl,"-woff 84","-woff 85" - -C_OPTIMISE_FLAGS = -O3 -OPT:Olimit=0 -CXX_OPTIMISE_FLAGS = -O3 -OPT:Olimit=0 -F77_OPTIMISE_FLAGS = -O3 -OPT:Olimit=0 -F90_OPTIMISE_FLAGS = -O3 -OPT:Olimit=0 - -C_DEBUG_FLAGS = -g3 -CXX_DEBUG_FLAGS = -g3 -F77_DEBUG_FLAGS = -g3 -F90_DEBUG_FLAGS = -g3 - -SYS_INC_DIRS = /usr/local/hdf/include /u1/eschnett/FlexIO -LIBDIRS = /usr/local/hdf/lib /u1/eschnett/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio ffio fpe fortran ftn ftn90 hdf5 mfhdf df jpeg z - -HDF5 = yes -HDF5_DIR = /usr/local/hdf5_64 - -MPI = NATIVE - -OPTIMISE = yes -DEBUG = yes diff --git a/Carpet/Carpet/options/carpet-hawaii-sgi b/Carpet/Carpet/options/carpet-hawaii-sgi deleted file mode 100644 index deaf2a33b..000000000 --- a/Carpet/Carpet/options/carpet-hawaii-sgi +++ /dev/null @@ -1,25 +0,0 @@ -CPP = /lib/cpp -DMPIPP_H - -CXXFLAGS = -LANG:std -LANG:vla=ON -FE:eliminate_duplicate_inline_copies -FE:template_in_elf_section -64 -mips4 -r12000 -no_auto_include -ptused -DMPIPP_H -LDFLAGS = -LANG:std -LANG:vla=ON -64 -mips4 -r12000 -Wl,"-woff 84","-woff 85" - -C_OPTIMISE_FLAGS = -O3 -INLINE -LNO -OPT:Olimit=100000 -CXX_OPTIMISE_FLAGS = -O3 -INLINE -LNO -OPT:Olimit=100000 -F77_OPTIMISE_FLAGS = -O3 -INLINE -LNO -OPT:Olimit=100000 -F90_OPTIMISE_FLAGS = -O3 -INLINE -LNO -OPT:Olimit=100000 - -C_DEBUG_FLAGS = -g3 -CXX_DEBUG_FLAGS = -g3 -F77_DEBUG_FLAGS = -g3 -F90_DEBUG_FLAGS = -g3 - -SYS_INC_DIRS = /usr/center/maui/eriks/include /usr/center/maui/eriks/FlexIO -LIBDIRS = /usr/center/maui/eriks/lib /usr/center/maui/eriks/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z - -HDF5 = yes -HDF5_DIR = /usr/center/maui/eriks -MPI = NATIVE - -OPTIMISE = yes -DEBUG = yes diff --git a/Carpet/Carpet/options/carpet-lilypond b/Carpet/Carpet/options/carpet-lilypond deleted file mode 100644 index 36d708f39..000000000 --- a/Carpet/Carpet/options/carpet-lilypond +++ /dev/null @@ -1,27 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -F77 = g77 -F90 = g77 - -CXXFLAGS = -DCARPET_REAL -ftemplate-depth-30 - -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -F90_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double -F77_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = -Wall -F77_WARN_FLAGS = -Wall - -SYS_INC_DIRS = /home/eschnett/proj/FlexIO -LIBDIRS = /home/eschnett/proj/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt g2c m - -HDF5 = yes -MPI = MPICH - -DEBUG = yes - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-lilypond-debug b/Carpet/Carpet/options/carpet-lilypond-debug deleted file mode 100644 index 5d12914b6..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-debug +++ /dev/null @@ -1,29 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -F77 = g77 -F90 = g77 - -CXXFLAGS = -DCARPET_REAL -ftemplate-depth-30 - -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -F90_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double -F77_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = -Wall -F77_WARN_FLAGS = -Wall - -#FLEXIO_DIR = /home/eschnett/proj/FlexIO -#SYS_INC_DIRS = /home/eschnett/proj/FlexIO -#LIBDIRS = /home/eschnett/proj/FlexIO -#LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt g2c m - -HDF5 = yes -MPI = MPICH - -DEBUG = yes -OPTIMISE = no - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-lilypond-gcc32 b/Carpet/Carpet/options/carpet-lilypond-gcc32 deleted file mode 100644 index 94ac69a44..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-gcc32 +++ /dev/null @@ -1,33 +0,0 @@ -CPP = /home/eschnett/gcc32/bin/cpp -traditional -DCARPET_REAL -CC = /home/eschnett/gcc32/bin/gcc -CXX = /home/eschnett/gcc32/bin/g++ -F77 = /home/eschnett/gcc32/bin/g77 -F90 = /home/eschnett/gcc32/bin/g77 - -CFLAGS = -m128bit-long-double -fmessage-length=0 -CXXFLAGS = -m128bit-long-double -fmessage-length=0 -DCARPET_REAL -ftemplate-depth-30 -F77FLAGS = -malign-double -fmessage-length=0 -F90FLAGS = -malign-double -fmessage-length=0 -LDFLAGS = -Wl,-rpath,/home/eschnett/gcc32/lib - -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentium2 -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentium2 -F90_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentium2 -F77_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentium2 - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = -Wall -F77_WARN_FLAGS = -Wall - -SYS_INC_DIRS = /home/eschnett/proj/FlexIO -LIBDIRS = /home/eschnett/proj/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt g2c m - -HDF5 = yes -MPI = MPICH - -DEBUG = yes - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-lilypond-ic b/Carpet/Carpet/options/carpet-lilypond-ic deleted file mode 100644 index 53d4fbf08..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-ic +++ /dev/null @@ -1,48 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -CC = icc -CXX = icc -F77 = ifc -F90 = ifc -LD = icc - -LIBDIRS = -LIBS = intrins IEPCF90 F90 imf m irc cxa cprts cxa - -# basic -CFLAGS = -DCARPET_REAL -c99 -CXXFLAGS = -DCARPET_REAL -c99 -F90FLAGS = -DCARPET_REAL -w95 -F77FLAGS = -DCARPET_REAL -w95 -LDFLAGS = -DCARPET_REAL -g - -# debugging -DEBUG = yes -C_DEBUG_FLAGS = -g -fp -CXX_DEBUG_FLAGS = -g -fp -F77_DEBUG_FLAGS = -g -fp -F90_DEBUG_FLAGS = -g -fp - -# warnings -#WARN=yes -C_WARN_FLAGS = -w1 -CXX_WARN_FLAGS = -w1 -F77_WARN_FLAGS = -F90_WARN_FLAGS = - -# optimisation -OPTIMISE = yes -C_OPTIMISE_FLAGS = -O3 -tpp6 -xi -CXX_OPTIMISE_FLAGS = -O3 -tpp6 -xi -F77_OPTIMISE_FLAGS = -O3 -tpp6 -F90_OPTIMISE_FLAGS = -O3 -tpp6 - -# HDF -HDF5 = YES -HDF5_DIR = /usr/local - -# MPI -MPI = MPICH -MPICH_DIR = /usr/lib/mpich - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-lilypond-ic-debug b/Carpet/Carpet/options/carpet-lilypond-ic-debug deleted file mode 100644 index 88e0f8abb..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-ic-debug +++ /dev/null @@ -1,45 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -CC = icc -CXX = icc -F77 = ifc -F90 = ifc -LD = icc - -LIBDIRS = -LIBS = intrins IEPCF90 F90 imf m irc cxa cprts cxa - -# basic -CFLAGS = -DCARPET_REAL -c99 -CXXFLAGS = -DCARPET_REAL -c99 -F90FLAGS = -DCARPET_REAL -w95 -F77FLAGS = -DCARPET_REAL -w95 -LDFLAGS = -DCARPET_REAL -g - -# debugging -DEBUG = yes -C_DEBUG_FLAGS = -g -fp -CXX_DEBUG_FLAGS = -g -fp -F77_DEBUG_FLAGS = -g -fp -F90_DEBUG_FLAGS = -g -fp - -# warnings -#WARN=yes -C_WARN_FLAGS = -w1 -CXX_WARN_FLAGS = -w1 -F77_WARN_FLAGS = -F90_WARN_FLAGS = - -# optimisation -OPTIMISE = yes -C_OPTIMISE_FLAGS = -O0 -tpp6 -xi -CXX_OPTIMISE_FLAGS = -O0 -tpp6 -xi -F77_OPTIMISE_FLAGS = -O0 -tpp6 -F90_OPTIMISE_FLAGS = -O0 -tpp6 - -# HDF -HDF5 = YES -HDF5_DIR = /usr/local - -# MPI -MPI = MPICH -MPICH_DIR = /usr/lib/mpich diff --git a/Carpet/Carpet/options/carpet-lilypond-ic7 b/Carpet/Carpet/options/carpet-lilypond-ic7 deleted file mode 100644 index 287525288..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-ic7 +++ /dev/null @@ -1,48 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -CC = /opt/intel/compiler70/ia32/bin/icc -CXX = /opt/intel/compiler70/ia32/bin/icc -F77 = /opt/intel/compiler70/ia32/bin/ifc -F90 = /opt/intel/compiler70/ia32/bin/ifc -LD = /opt/intel/compiler70/ia32/bin/icc - -LIBDIRS = -LIBS = intrins IEPCF90 F90 imf m irc cxa cprts cxa - -# basic -CFLAGS = -DCARPET_REAL -c99 -CXXFLAGS = -DCARPET_REAL -c99 -F90FLAGS = -DCARPET_REAL -w95 -F77FLAGS = -DCARPET_REAL -w95 -LDFLAGS = -g -Qoption,ld,-rpath,/opt/intel/compiler70/ia32/lib - -# debugging -DEBUG = yes -C_DEBUG_FLAGS = -g -fp -CXX_DEBUG_FLAGS = -g -fp -F77_DEBUG_FLAGS = -g -fp -F90_DEBUG_FLAGS = -g -fp - -# warnings -#WARN=yes -C_WARN_FLAGS = -w1 -CXX_WARN_FLAGS = -w1 -F77_WARN_FLAGS = -F90_WARN_FLAGS = - -# optimisation -OPTIMISE = yes -C_OPTIMISE_FLAGS = -O3 -tpp6 -xi -CXX_OPTIMISE_FLAGS = -O3 -tpp6 -xi -F77_OPTIMISE_FLAGS = -O3 -tpp6 -F90_OPTIMISE_FLAGS = -O3 -tpp6 - -# HDF -HDF5 = YES -HDF5_DIR = /usr/local - -# MPI -MPI = MPICH -MPICH_DIR = /usr/lib/mpich - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-lilypond-lf95 b/Carpet/Carpet/options/carpet-lilypond-lf95 deleted file mode 100644 index a6c6abbd8..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-lf95 +++ /dev/null @@ -1,37 +0,0 @@ -LAHEY = /usr/local/lf9560 - -CPP = cpp -traditional -DCARPET_REAL -F77 = lf95 -F90 = lf95 - -CFLAGS = -g3 -CXXFLAGS = -g3 -DCARPET_REAL -ftemplate-depth-30 -F77FLAGS = --pca -F90FLAGS = --pca -LDFLAGS = -g3 - -DEBUG = no -C_DEBUG_FLAGS = -CXX_DEBUG_FLAGS = -F77_DEBUG_FLAGS = -g --chk -F90_DEBUG_FLAGS = -g --chk - -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -F90_OPTIMISE_FLAGS = -O --tpp -F77_OPTIMISE_FLAGS = -O --tpp - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = --warn -F77_WARN_FLAGS = --warn - -SYS_INC_DIRS = /home/eschnett/proj/FlexIO -LIBDIRS = /home/eschnett/proj/FlexIO /usr/local/lf9560/lib -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt fj9i6 fj9f6 fj9e6 fccx86_6a m - -HDF5 = yes -MPI = MPICH - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-lilypond-lf95-debug b/Carpet/Carpet/options/carpet-lilypond-lf95-debug deleted file mode 100644 index 454afbb4a..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-lf95-debug +++ /dev/null @@ -1,38 +0,0 @@ -LAHEY = /usr/local/lf9560 - -CPP = cpp -traditional -DCARPET_REAL -F77 = lf95 -F90 = lf95 - -CFLAGS = -g3 -CXXFLAGS = -g3 -DCARPET_REAL -ftemplate-depth-30 -F77FLAGS = --pca -F90FLAGS = --pca -LDFLAGS = -g3 - -DEBUG = yes -C_DEBUG_FLAGS = -CXX_DEBUG_FLAGS = -F77_DEBUG_FLAGS = -g --chk -F90_DEBUG_FLAGS = -g --chk - -OPTIMISE = no -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -F90_OPTIMISE_FLAGS = -O --tpp -F77_OPTIMISE_FLAGS = -O --tpp - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = --warn -F77_WARN_FLAGS = --warn - -SYS_INC_DIRS = /home/eschnett/proj/FlexIO -LIBDIRS = /home/eschnett/proj/FlexIO /usr/local/lf9560/lib -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt fj9i6 fj9f6 fj9e6 fccx86_6a m - -HDF5 = yes -MPI = MPICH - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-lilypond-prof b/Carpet/Carpet/options/carpet-lilypond-prof deleted file mode 100644 index 18d3a97f6..000000000 --- a/Carpet/Carpet/options/carpet-lilypond-prof +++ /dev/null @@ -1,31 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -F77 = g77 -F90 = g77 - -CFLAGS = -pg -CXXFLAGS = -pg -DCARPET_REAL -ftemplate-depth-30 -F77FLAGS = -pg -F90FLAGS = -pg -LDFLAGS = -pg - -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -F90_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double -F77_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = -Wall -F77_WARN_FLAGS = -Wall - -SYS_INC_DIRS = /home/eschnett/proj/FlexIO -LIBDIRS = /home/eschnett/proj/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt g2c m - -HDF5 = yes -MPI = MPICH - -DEBUG = yes - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-mintaka b/Carpet/Carpet/options/carpet-mintaka deleted file mode 100644 index a761fce3c..000000000 --- a/Carpet/Carpet/options/carpet-mintaka +++ /dev/null @@ -1,30 +0,0 @@ -CPP = cpp -traditional -DCARPET_REAL -F77 = g77 -F90 = g77 - -CXXFLAGS = -DCARPET_REAL -ftemplate-depth-30 -LDFLAGS = -Wl,-rpath,/home/schnette/lib - -C_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -CXX_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -F90_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double -F77_OPTIMISE_FLAGS = -O2 -funroll-loops -march=pentiumpro -malign-double - -C_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -CXX_WARN_FLAGS = -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -F90_WARN_FLAGS = -Wall -F77_WARN_FLAGS = -Wall - -SYS_INC_DIRS = /home/schnette/lilypond/proj/FlexIO -LIBDIRS = /home/schnette/lilypond/proj/FlexIO -LIBS = AMR h5io hdfio hlio ieeeio hdf5 mfhdf df jpeg z crypt g2c m - -HDF5 = yes -HDF5_DIR = /home/schnette -MPI = MPICH -MPICH_DIR = /usr/lib/mpich - -DEBUG = yes - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/options/carpet-tat-ic b/Carpet/Carpet/options/carpet-tat-ic deleted file mode 100644 index 09b9bd87a..000000000 --- a/Carpet/Carpet/options/carpet-tat-ic +++ /dev/null @@ -1,48 +0,0 @@ -CPP = /home/schnette/gcc/bin/cpp -traditional -DCARPET_REAL -CC = icc -CXX = icc -F77 = ifc -F90 = ifc -LD = icc - -LIBDIRS = -LIBS = intrins IEPCF90 F90 imf m irc cxa cprts cxa - -# basic -CFLAGS = -DCARPET_REAL -c99 -CXXFLAGS = -DCARPET_REAL -c99 -F90FLAGS = -DCARPET_REAL -w95 -F77FLAGS = -DCARPET_REAL -w95 -LDFLAGS = -DCARPET_REAL -g -Qoption,ld,-rpath,/home/schnette/opt/intel/compiler60/ia32/lib - -# debugging -DEBUG = yes -C_DEBUG_FLAGS = -g -fp -CXX_DEBUG_FLAGS = -g -fp -F77_DEBUG_FLAGS = -g -fp -F90_DEBUG_FLAGS = -g -fp - -# warnings -#WARN=yes -C_WARN_FLAGS = -w1 -CXX_WARN_FLAGS = -w1 -F77_WARN_FLAGS = -F90_WARN_FLAGS = - -# optimisation -OPTIMISE = yes -C_OPTIMISE_FLAGS = -O3 -tpp6 -CXX_OPTIMISE_FLAGS = -O3 -tpp6 -F77_OPTIMISE_FLAGS = -O3 -tpp6 -F90_OPTIMISE_FLAGS = -O3 -tpp6 - -# HDF -HDF5 = YES -HDF5_DIR = /usr/local - -# MPI -MPI = MPICH -MPICH_DIR = /usr/lib/mpich - -# PTHREADS -PTHREADS = yes diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc index 3ed5ae848..38efe537f 100644 --- a/Carpet/Carpet/src/Evolve.cc +++ b/Carpet/Carpet/src/Evolve.cc @@ -467,6 +467,8 @@ namespace Carpet { for (int ml=mglevels-1; ml>=0; --ml) { bool have_done_global_mode = false; + bool have_done_early_global_mode = false; + bool have_done_late_global_mode = false; bool have_done_anything = false; for (int rl=0; rl<reflevels; ++rl) { @@ -477,14 +479,20 @@ namespace Carpet { ENTER_LEVEL_MODE (cctkGH, rl) { BeginTimingLevel (cctkGH); - do_early_global_mode = not have_done_global_mode; + do_early_global_mode = not have_done_early_global_mode; do_late_global_mode = reflevel==reflevels-1; do_early_meta_mode = do_early_global_mode and mglevel==mglevels-1; do_late_meta_mode = do_late_global_mode and mglevel==0; do_global_mode = do_late_global_mode; do_meta_mode = do_global_mode and do_late_meta_mode; assert (not (have_done_global_mode and do_global_mode)); + assert (not (have_done_early_global_mode and + do_early_global_mode)); + assert (not (have_done_late_global_mode and + do_late_global_mode)); have_done_global_mode |= do_global_mode; + have_done_early_global_mode |= do_early_global_mode; + have_done_late_global_mode |= do_late_global_mode; have_done_anything = true; if (use_tapered_grids and reflevel > 0) { @@ -544,6 +552,8 @@ namespace Carpet { } // for rl if (have_done_anything) assert (have_done_global_mode); + if (have_done_anything) assert (have_done_early_global_mode); + if (have_done_anything) assert (have_done_late_global_mode); } // for ml diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc index bb5f3ba20..5ff48751c 100644 --- a/Carpet/Carpet/src/Recompose.cc +++ b/Carpet/Carpet/src/Recompose.cc @@ -556,8 +556,8 @@ namespace Carpet { for (int c=0; c<hh.components(rl); ++c) { const rvect origin = domainspecs.AT(m).exterior_min; const rvect delta = (domainspecs.AT(m).exterior_max - domainspecs.AT(m).exterior_min) / rvect (domainspecs.AT(m).npoints - 1); - const ivect lower = dd.boxes.AT(ml).AT(rl).AT(c).exterior.lower(); - const ivect upper = dd.boxes.AT(ml).AT(rl).AT(c).exterior.upper(); + const ivect lower = dd.light_boxes.AT(ml).AT(rl).AT(c).exterior.lower(); + const ivect upper = dd.light_boxes.AT(ml).AT(rl).AT(c).exterior.upper(); const int convfact = ipow(mgfact, ml); const ivect levfact = spacereffacts.AT(rl); cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" @@ -850,11 +850,7 @@ namespace Carpet { for (int ml=0; ml<(int)regsss.size(); ++ml) { file << m << " " << ml << " reflevels " << regsss.AT(ml).size() << eol; for (int rl=0; rl<(int)regsss.AT(ml).size(); ++rl) { - ibset extents; - for (int c=0; c<(int)regsss.AT(ml).AT(rl).size(); ++c) { - extents += regsss.AT(ml).AT(rl).AT(c).extent; - } - extents.normalize(); + ibset const extents (regsss.AT(ml).AT(rl), & region_t::extent); file << m << " " << ml << " " << rl << " regions " << extents.setsize() << eol; int c=0; for (ibset::const_iterator @@ -964,7 +960,7 @@ namespace Carpet { gh const * const hh = arrdata.AT(g).AT(0).hh; dh const * const dd = arrdata.AT(g).AT(0).dd; for (int c=0; c<hh->components(0); ++c) { - dh::dboxes const & b = dd->boxes.AT(0).AT(0).AT(c); + dh::light_dboxes const & b = dd->light_boxes.AT(0).AT(0).AT(c); num_active_array_points += num_tl * num_vars * b.active_size; num_total_array_points += num_tl * num_vars * b.exterior_size; size_total_array_points += num_tl * size_vars * b.exterior_size; @@ -1001,7 +997,7 @@ namespace Carpet { } for (int c=0; c<hh->components(rl); ++c) { ++ num_comps; - dh::dboxes const & b = dd->boxes.AT(ml).AT(rl).AT(c); + dh::light_dboxes const & b = dd->light_boxes.AT(ml).AT(rl).AT(c); num_active_mem_points += num_gfs * b.active_size; num_owned_mem_points += num_gfs * b.owned_size; num_total_mem_points += num_gfs * b.exterior_size; @@ -1032,7 +1028,7 @@ namespace Carpet { for (int ml=0; ml<mglevels; ++ml) { for (int c=0; c<hh->components(rl); ++c) { int const p = hh->processor(rl,c); - dh::dboxes const & b = dd->boxes.AT(ml).AT(rl).AT(c); + dh::light_dboxes const & b = dd->light_boxes.AT(ml).AT(rl).AT(c); num_active_per_proc.AT(p) += num_gfs * b.active_size; } } diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc index 4961cec3b..f8fd5a7fc 100644 --- a/Carpet/Carpet/src/SetupGH.cc +++ b/Carpet/Carpet/src/SetupGH.cc @@ -1583,7 +1583,8 @@ namespace Carpet { baseextents.AT(ml).AT(rl) = baseextent; } else { // Refine next coarser refinement level - if (refcentering == vertex_centered) { + switch (refcentering) { + case vertex_centered: { ibbox const & cbox = baseextents.AT(ml).AT(rl-1); assert (not any (any (is_staggered))); i2vect const bnd_shift = @@ -1595,7 +1596,9 @@ namespace Carpet { ibbox (cbox_phys.lower(), cbox_phys.upper(), fstride); ibbox const fbox = fbox_phys.expand (bnd_shift); baseextents.AT(ml).AT(rl) = fbox; - } else { + break; + } + case cell_centered: { ibbox const & cbox = baseextents.AT(ml).AT(rl-1); assert (all (all (is_staggered))); ivect const cstride = cbox.stride(); @@ -1621,6 +1624,10 @@ namespace Carpet { fbox_phys.upper() + bnd_shift_fstride[1], fbox_phys.stride()); baseextents.AT(ml).AT(rl) = fbox; + break; + } + default: + assert (0); } } } else { diff --git a/Carpet/Carpet/src/modes.cc b/Carpet/Carpet/src/modes.cc index 038815393..d101a509b 100644 --- a/Carpet/Carpet/src/modes.cc +++ b/Carpet/Carpet/src/modes.cc @@ -113,7 +113,7 @@ namespace Carpet { const int c = dist::rank(); const ibbox& baseext = arrdata.AT(group).AT(m).hh->baseextents.AT(ml).AT(rl); - const ibbox& ext = arrdata.AT(group).AT(m).dd->boxes.AT(ml).AT(rl).AT(c).exterior; + const ibbox& ext = arrdata.AT(group).AT(m).dd->light_boxes.AT(ml).AT(rl).AT(c).exterior; const b2vect& obnds = arrdata.AT(group).AT(m).hh->outer_boundaries(rl,c); cGroupDynamicData & info = groupdata.AT(group).info; @@ -545,7 +545,7 @@ namespace Carpet { // Set cGH fields const ibbox& baseext = vhh.AT(map)->baseextents.AT(mglevel).AT(reflevel); - const ibbox& ext = vdd.AT(map)->boxes.AT(mglevel).AT(reflevel).AT(component).exterior; + const ibbox& ext = vdd.AT(map)->light_boxes.AT(mglevel).AT(reflevel).AT(component).exterior; const b2vect& obnds = vhh.AT(map)->outer_boundaries(reflevel,component); ivect_ref(cctkGH->cctk_lbnd) = @@ -1352,6 +1352,30 @@ namespace Carpet { // + // Loop over a bounding box + // + + void ibbox2iminimax (ibbox const& ext, // component extent + ibbox const& box, // this bbox + ivect& imin, ivect& imax) + { + ivect const izero = ivect(0); + + assert (all ((box.lower() - ext.lower() ) >= 0)); + assert (all ((box.upper() - ext.lower() + ext.stride()) >= 0)); + assert (all ((box.lower() - ext.lower() ) % ext.stride() == 0)); + assert (all ((box.upper() - ext.lower() + ext.stride()) % ext.stride() == 0)); + + imin = (box.lower() - ext.lower() ) / ext.stride(); + imax = (box.upper() - ext.lower() + ext.stride()) / ext.stride(); + + assert (all (0 <= imin)); + assert (box.empty() xor all (imin <= imax)); + } + + + + // // Call a scheduling group // diff --git a/Carpet/Carpet/src/modes.hh b/Carpet/Carpet/src/modes.hh index 868c54172..37d3ec646 100644 --- a/Carpet/Carpet/src/modes.hh +++ b/Carpet/Carpet/src/modes.hh @@ -434,6 +434,48 @@ namespace Carpet { + // Loop over a bounding box + + void ibbox2iminimax (ibbox const& ext, // component extent + ibbox const& box, // this bbox + ivect& imin, ivect& imax); + +#define LOOP_OVER_BBOX(cctkGH, box_, box, imin, imax) \ + do { \ + bool bbox_loop_ = true; \ + assert (Carpet::is_local_mode()); \ + dh::light_mboxes const& light_boxes_ = \ + Carpet::vdd.AT(Carpet::map)->light_boxes; \ + ibbox const& ext_ = \ + light_boxes_.AT(mglevel).AT(reflevel).AT(component).exterior; \ + ibbox const box = (box_) & ext_; \ + if (not box.empty()) { \ + ivect imin; \ + ivect imax; \ + Carpet::ibbox2iminimax (ext_, box, imin, imax); \ + { +#define END_LOOP_OVER_BBOX \ + } \ + } /* if not empty */ \ + assert (bbox_loop_); \ + bbox_loop_ = false; \ + } while (false) + +#define LOOP_OVER_BSET(cctkGH, set_, box, imin, imax) \ + do { \ + bool bset_loop_ = true; \ + ibset const& set1_ (set_); \ + for (ibset::const_iterator bi = set1_.begin(); bi != set1_.end(); ++bi) { \ + LOOP_OVER_BBOX(cctkGH, *bi, box, imin, imax) { +#define END_LOOP_OVER_BSET \ + } END_LOOP_OVER_BBOX; \ + } /* for */ \ + assert (bset_loop_); \ + bset_loop_ = false; \ + } while (false) + + + } // namespace Carpet #endif // #ifndef MODES_HH diff --git a/Carpet/CarpetEvolutionMask/src/evolution_mask.cc b/Carpet/CarpetEvolutionMask/src/evolution_mask.cc index 527ff33f2..a6b3c844d 100644 --- a/Carpet/CarpetEvolutionMask/src/evolution_mask.cc +++ b/Carpet/CarpetEvolutionMask/src/evolution_mask.cc @@ -84,7 +84,6 @@ namespace CarpetEvolutionMask { refcomp = refcomp.expand(expand_left,expand_right); refined |= refcomp; } - refined.normalize(); // cout << "refined: " << refined << endl; @@ -101,15 +100,10 @@ namespace CarpetEvolutionMask { for (int d=0; d<dim;d++) { antishrinkby[d] = cctkGH->cctk_nghostzones[d] + buffer_width; } - ibset antishrunk; - for (ibset::const_iterator bi = antirefined.begin(); - bi != antirefined.end(); - ++bi) - { - antishrunk |= (*bi).expand(antishrinkby, antishrinkby).expanded_for(coarsebase); - } - antishrunk.normalize(); - + ibset antishrunk + (antirefined.expand(antishrinkby, antishrinkby). + expanded_for(coarsebase)); + // cout << "antishrunk1: " << antishrunk << endl; // now cut away dangling edges @@ -118,7 +112,7 @@ namespace CarpetEvolutionMask { // cout << "antishrunk2: " << antishrunk << endl; // cut holes into coarsebase - ibset shrunk = coarsebase - antishrunk; + ibset const shrunk = coarsebase - antishrunk; // cout << "shrunk: " << shrunk << endl; @@ -127,13 +121,11 @@ namespace CarpetEvolutionMask { for (int c=0; c<hh.components(reflevel-1); ++c) { parent |= hh.extent(mglevel,reflevel-1,c); } - parent.normalize(); // cout << "parent: " << parent << endl; // Subtract the refined region - ibset notrefined = parent - shrunk; - notrefined.normalize(); + ibset const notrefined = parent - shrunk; // cout << "notrefined: " << notrefined << endl; @@ -143,23 +135,14 @@ namespace CarpetEvolutionMask { enlargeby[d] = cctkGH->cctk_nghostzones[d] + buffer_width + stencil_size; } - ibset enlarged; - for (ibset::const_iterator bi = notrefined.begin(); - bi != notrefined.end(); - ++bi) - { - enlarged |= (*bi).expand(enlargeby, enlargeby); - } - enlarged.normalize(); + ibset const enlarged (notrefined.expand(enlargeby, enlargeby)); // cout << "enlarged: " << enlarged << endl; // Intersect with the original union - ibset evolveon = parent & enlarged; - evolveon.normalize(); + ibset const evolveon = parent & enlarged; - ibset notevolveon = parent - evolveon; - notevolveon.normalize(); + ibset const notevolveon = parent - evolveon; // cout << "notevolveon: " << notevolveon << endl; @@ -178,7 +161,7 @@ namespace CarpetEvolutionMask { DECLARE_CCTK_ARGUMENTS; ibbox const & ext - = dd.boxes.at(mglevel).at(reflevel).at(component).exterior; + = dd.light_boxes.at(mglevel).at(reflevel).at(component).exterior; for (ibset::const_iterator bi = notevolveon.begin(); bi != notevolveon.end(); diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc index 349c18e9b..216de8ae3 100644 --- a/Carpet/CarpetIOASCII/src/ioascii.cc +++ b/Carpet/CarpetIOASCII/src/ioascii.cc @@ -570,7 +570,7 @@ namespace CarpetIOASCII { int const proc = hh->processor(rl,c); if (dist::rank() == proc or dist::rank() == ioproc) { - const ibbox& data_ext = dd->boxes.at(ml).at(rl).at(c).exterior; + const ibbox& data_ext = dd->light_boxes.at(ml).at(rl).at(c).exterior; const ibbox ext = GetOutputBBox (cctkGH, group, rl, m, c, data_ext); CCTK_REAL coord_time; diff --git a/Carpet/CarpetIOHDF5/src/Input.cc b/Carpet/CarpetIOHDF5/src/Input.cc index 3fbf02192..10e71036f 100644 --- a/Carpet/CarpetIOHDF5/src/Input.cc +++ b/Carpet/CarpetIOHDF5/src/Input.cc @@ -431,9 +431,8 @@ int Recover (cGH* cctkGH, const char *basefilename, int called_from) for (int lc=0; lc<data.hh->local_components(reflevel); ++lc) { int const c = data.hh->get_component(reflevel,lc); group_bboxes[gindex][m] |= - data.dd->boxes.at(mglevel).at(reflevel).at(c).exterior; + data.dd->light_boxes.at(mglevel).at(reflevel).at(c).exterior; } - group_bboxes[gindex][m].normalize(); } } else { if (mglevel==0 and reflevel==0) { @@ -441,8 +440,7 @@ int Recover (cGH* cctkGH, const char *basefilename, int called_from) struct arrdesc& data = arrdata.at(gindex).at(m); int const c=dist::rank(); group_bboxes[gindex][m] |= - data.dd->boxes.at(mglevel).at(reflevel).at(c).exterior; - group_bboxes[gindex][m].normalize(); + data.dd->light_boxes.at(mglevel).at(reflevel).at(c).exterior; } } } @@ -658,7 +656,6 @@ int Recover (cGH* cctkGH, const char *basefilename, int called_from) "variable '%s' timelevel %d has been read only partially", fullname, tl); // for (int m = 0; m < maps; m++) { - // bboxes_read[vindex][tl][m].normalize(); // const int gindex = CCTK_GroupIndexFromVarI (vindex); // cout << "Need: " << group_bboxes[gindex][m] << endl; // cout << "Read: " << bboxes_read[vindex][tl][m] << endl; @@ -1168,7 +1165,7 @@ static int ReadVar (const cGH* const cctkGH, } #endif const ibbox& exterior_membox = - data.dd->boxes.at(mglevel).at(reflevel).at(component).exterior; + data.dd->light_boxes.at(mglevel).at(reflevel).at(component).exterior; // skip this dataset if it doesn't overlap with this component's // exterior, or if it only reads what has already been read @@ -1189,11 +1186,10 @@ static int ReadVar (const cGH* const cctkGH, // get the overlap with this component's exterior const ibbox& membox = - data.dd->boxes.at(mglevel).at(reflevel).at(component).exterior; + data.dd->light_boxes.at(mglevel).at(reflevel).at(component).exterior; const ibbox overlap = membox & filebox; // cout << "Overlap: " << overlap << endl; bboxes_read.at(Carpet::map) |= overlap; - bboxes_read.at(Carpet::map).normalize(); // cout << "New read: " << bboxes_read.at(Carpet::map) << endl; // calculate hyperslab selection parameters diff --git a/Carpet/CarpetIOHDF5/src/Output.cc b/Carpet/CarpetIOHDF5/src/Output.cc index 70aa35af9..28516b76b 100644 --- a/Carpet/CarpetIOHDF5/src/Output.cc +++ b/Carpet/CarpetIOHDF5/src/Output.cc @@ -77,15 +77,12 @@ int WriteVarUnchunked (const cGH* const cctkGH, // Using "interior" removes ghost zones and refinement boundaries. #if 0 bboxes += arrdata.at(gindex).at(Carpet::map).dd-> - boxes.at(mglevel).at(refinementlevel).at(c).interior; + light_boxes.at(mglevel).at(refinementlevel).at(c).interior; #endif bboxes += arrdata.at(gindex).at(Carpet::map).dd-> - boxes.at(mglevel).at(refinementlevel).at(c).exterior; + light_boxes.at(mglevel).at(refinementlevel).at(c).exterior; } - // Normalise the set, i.e., try to represent the set with fewer bboxes. - bboxes.normalize(); - #if 0 // According to Cactus conventions, DISTRIB=CONSTANT arrays // (including grid scalars) are assumed to be the same on all @@ -177,10 +174,10 @@ int WriteVarUnchunked (const cGH* const cctkGH, dh const * const dd = arrdata.at(gindex).at(Carpet::map).dd; #if 0 ibbox const overlap = *bbox & - dd->boxes.at(mglevel).at(refinementlevel).at(component).interior; + dd->light_boxes.at(mglevel).at(refinementlevel).at(component).interior; #endif ibbox const overlap = *bbox & - dd->boxes.at(mglevel).at(refinementlevel).at(component).exterior; + dd->light_boxes.at(mglevel).at(refinementlevel).at(component).exterior; // Continue if this component is not part of this combination if (overlap.empty()) continue; @@ -341,7 +338,7 @@ int WriteVarChunkedSequential (const cGH* const cctkGH, gh const * const hh = arrdata.at(gindex).at(Carpet::map).hh; dh const * const dd = arrdata.at(gindex).at(Carpet::map).dd; ibbox const& bbox = - dd->boxes.at(mglevel).at(refinementlevel).at(component).exterior; + dd->light_boxes.at(mglevel).at(refinementlevel).at(component).exterior; // Get the shape of the HDF5 dataset (in Fortran index order) hsize_t shape[dim]; @@ -405,7 +402,7 @@ int WriteVarChunkedSequential (const cGH* const cctkGH, datasetname << " rl=" << refinementlevel; } if (arrdata.at(gindex).at(Carpet::map).dd-> - boxes.at(mglevel).at(refinementlevel).size () > 1 and + light_boxes.at(mglevel).at(refinementlevel).size () > 1 and group.disttype != CCTK_DISTRIB_CONSTANT) { datasetname << " c=" << component; } @@ -517,7 +514,7 @@ int WriteVarChunkedParallel (const cGH* const cctkGH, // 0 : component, mglevel)->extent(); const dh* dd = arrdata.at(gindex).at(Carpet::map).dd; const ibbox& bbox = - dd->boxes.AT(mglevel).AT(refinementlevel) + dd->light_boxes.AT(mglevel).AT(refinementlevel) .AT(group.disttype == CCTK_DISTRIB_CONSTANT ? 0 : component).exterior; // Don't create zero-sized components @@ -566,7 +563,7 @@ int WriteVarChunkedParallel (const cGH* const cctkGH, datasetname << " rl=" << refinementlevel; } if (arrdata.at(gindex).at(Carpet::map).dd-> - boxes.at(mglevel).at(refinementlevel).size () > 1 and + light_boxes.at(mglevel).at(refinementlevel).size () > 1 and group.disttype != CCTK_DISTRIB_CONSTANT) { datasetname << " c=" << component; } diff --git a/Carpet/CarpetIOHDF5/src/OutputSlice.cc b/Carpet/CarpetIOHDF5/src/OutputSlice.cc index 49e75d7c7..d6459c094 100644 --- a/Carpet/CarpetIOHDF5/src/OutputSlice.cc +++ b/Carpet/CarpetIOHDF5/src/OutputSlice.cc @@ -556,7 +556,7 @@ namespace CarpetIOHDF5 { int const proc = hh->processor(rl,c); if (dist::rank() == proc or dist::rank() == ioproc) { - const ibbox& data_ext = dd->boxes.at(ml).at(rl).at(c).exterior; + const ibbox& data_ext = dd->light_boxes.at(ml).at(rl).at(c).exterior; const ibbox ext = GetOutputBBox (cctkGH, group, rl, m, c, data_ext); CCTK_REAL coord_time; diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc index b04d268dd..0f809926e 100644 --- a/Carpet/CarpetInterp/src/interp.cc +++ b/Carpet/CarpetInterp/src/interp.cc @@ -1342,8 +1342,7 @@ namespace CarpetInterp { ENTER_LEVEL_MODE (cctkGH, rl) { // Number of neccessary time levels - CCTK_REAL const level_time = - cctkGH->cctk_time; + CCTK_REAL const level_time = cctkGH->cctk_time; vector<int> num_tl (N_input_arrays, 0); vector<bool> need_time_interp (N_output_arrays); for (int m=0; m<N_output_arrays; ++m) { diff --git a/Carpet/CarpetInterp2/interface.ccl b/Carpet/CarpetInterp2/interface.ccl index 15b3919c7..0b1363a6e 100644 --- a/Carpet/CarpetInterp2/interface.ccl +++ b/Carpet/CarpetInterp2/interface.ccl @@ -1,4 +1,4 @@ -# Interface definition for thorn CarpetInterp +# Interface definition for thorn CarpetInterp2 IMPLEMENTS: interp2 @@ -34,6 +34,22 @@ REQUIRES FUNCTION GetCoordRange +CCTK_INT FUNCTION \ + MultiPatch_LocalToGlobal \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN ndims, \ + CCTK_INT IN npoints, \ + CCTK_INT ARRAY IN patch, \ + CCTK_POINTER_TO_CONST IN localcoords, \ + CCTK_POINTER IN globalcoords, \ + CCTK_POINTER IN dxda, \ + CCTK_POINTER IN det_dxda, \ + CCTK_POINTER IN dadx, \ + CCTK_POINTER IN ddxdada, \ + CCTK_POINTER IN ddadxdx, \ + CCTK_POINTER IN dddxdadada) +USES FUNCTION MultiPatch_LocalToGlobal + CCTK_INT FUNCTION \ MultiPatch_GlobalToLocal \ (CCTK_POINTER_TO_CONST IN cctkGH, \ diff --git a/Carpet/CarpetInterp2/param.ccl b/Carpet/CarpetInterp2/param.ccl index 569c4f1f2..f7b5fffef 100644 --- a/Carpet/CarpetInterp2/param.ccl +++ b/Carpet/CarpetInterp2/param.ccl @@ -1,4 +1,4 @@ -# Parameter definitions for thorn CarpetInterp +# Parameter definitions for thorn CarpetInterp2 BOOLEAN verbose "Produce info output" STEERABLE=always { diff --git a/Carpet/CarpetInterp2/schedule.ccl b/Carpet/CarpetInterp2/schedule.ccl index 96d18609b..2712c53bb 100644 --- a/Carpet/CarpetInterp2/schedule.ccl +++ b/Carpet/CarpetInterp2/schedule.ccl @@ -1 +1 @@ -# Schedule definitions for thorn CarpetInterp +# Schedule definitions for thorn CarpetInterp2 diff --git a/Carpet/CarpetInterp2/src/fasterp.cc b/Carpet/CarpetInterp2/src/fasterp.cc index 4b5126b78..c59330556 100644 --- a/Carpet/CarpetInterp2/src/fasterp.cc +++ b/Carpet/CarpetInterp2/src/fasterp.cc @@ -643,19 +643,33 @@ namespace CarpetInterp2 { ostringstream msg; msg << "Interpolation point " << n << " on map " << m << " " << "at " << pos << " is outside of the grid hierarchy"; + msg << "\n" + << "rl=" << rl << " c=" << c << "\n" + << "rpos=" << rpos << "\n" + << "ipos=" << ipos << "\n" + << "lower=" << lower << "\n" + << "upper=" << upper << "\n" + << "delta=" << delta << "\n" + << "idelta=" << idelta << "\n" + << "hh=" << *hh << "\n"; CCTK_WARN (CCTK_WARN_ABORT, msg.str().c_str()); } } assert (rl>=0 and c>=0); ibbox const & ext = - Carpet::vdd.AT(m)->boxes.AT(Carpet::mglevel).AT(rl).AT(c).exterior; + Carpet::vdd.AT(m)->light_boxes.AT(Carpet::mglevel).AT(rl).AT(c).exterior; rvect dpos = rpos - rvect(ipos); // Convert from Carpet indexing to grid point indexing assert (all (ipos % ext.stride() == ivect(0))); ipos /= ext.stride(); dpos /= rvect(ext.stride()); + if (not (all (abs(dpos) <= rvect(0.5)))) { + cout << "fasterp.cc:659\n" + << " dpos=" << dpos << "\n" + << " ext=" << ext << "\n"; + } assert (all (abs(dpos) <= rvect(0.5))); ivect const ind = ipos - ext.lower() / ext.stride(); @@ -905,7 +919,7 @@ namespace CarpetInterp2 { int const c = themrc.c; assert (Carpet::vhh.AT(m)->is_local(rl,c)); ibbox const & ext = - Carpet::vdd.AT(m)->boxes.AT(Carpet::mglevel).AT(rl).AT(c).exterior; + Carpet::vdd.AT(m)->light_boxes.AT(Carpet::mglevel).AT(rl).AT(c).exterior; send_comp.lsh = ext.shape() / ext.stride(); send_comp.offset = offset; diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl index 689553c93..b4421f097 100644 --- a/Carpet/CarpetLib/param.ccl +++ b/Carpet/CarpetLib/param.ccl @@ -99,9 +99,9 @@ BOOLEAN use_ipm_timing_regions "Call IPM (via MPI_Pcontrol) to define regions" S INT print_memstats_every "Report periodically how much memory is used per process" STEERABLE=always { -1 :: "don't report" - 0 :: "don't report" + 0 :: "report after setting up initial data" 1:* :: "report every so many iterations" -} 0 +} -1 INT max_allowed_memory_MB "Maximum allowed amount of memory per process that can be allocated for grid variables (in Megabytes)" STEERABLE=always { diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc index 164d63884..e0b72f1a6 100644 --- a/Carpet/CarpetLib/src/bbox.cc +++ b/Carpet/CarpetLib/src/bbox.cc @@ -19,7 +19,7 @@ using namespace std; // Consistency checks -template<class T, int D> +template<typename T, int D> void bbox<T,D>::assert_bbox_limits () const { ASSERT_BBOX (all(_stride>T(0))); @@ -42,7 +42,7 @@ void bbox<T,D>::assert_bbox_limits () const // Poison -template<class T, int D> +template<typename T, int D> bbox<T,D> bbox<T,D>::poison () { DECLARE_CCTK_PARAMETERS; @@ -54,7 +54,7 @@ bbox<T,D> bbox<T,D>::poison () // Accessors -template<class T, int D> +template<typename T, int D> typename bbox<T,D>::size_type bbox<T,D>::size () const { if (empty()) return 0; const vect<T,D> sh(shape()/stride()); @@ -78,12 +78,12 @@ typename bbox<T,D>::size_type bbox<T,D>::size () const { // Queries // Containment -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::contains (const vect<T,D>& x) const { return all(x>=lower() and x<=upper()); } -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::is_contained_in (const bbox& b) const { if (empty()) return true; // no alignment check @@ -91,7 +91,7 @@ bool bbox<T,D>::is_contained_in (const bbox& b) const { } // Intersection -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::intersects (const bbox& b) const { if (empty()) return false; if (b.empty()) return false; @@ -100,27 +100,27 @@ bool bbox<T,D>::intersects (const bbox& b) const { } // Alignment check -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::is_aligned_with (const bbox& b) const { return all(stride()==b.stride() and (lower()-b.lower()) % stride() == T(0)); } // Operators -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::operator== (const bbox& b) const { if (empty() and b.empty()) return true; ASSERT_BBOX (all(stride()==b.stride())); return all(lower()==b.lower() and upper()==b.upper()); } -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::operator!= (const bbox& b) const { return not (*this == b); } #if 0 // Introduce an ordering on bboxes -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::operator< (const bbox& b) const { // An arbitraty order: empty boxes come first, then sorted by lower // bound, then by upper bound, then by coarseness @@ -142,28 +142,28 @@ bool bbox<T,D>::operator< (const bbox& b) const { } #endif -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::operator<= (const bbox& b) const { return is_contained_in (b); } -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::operator>= (const bbox& b) const { return b <= *this; } -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::operator< (const bbox& b) const { return *this <= b and *this != b; } -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::operator> (const bbox& b) const { return b < *this; } // Expand the bbox a little by multiples of the stride -template<class T, int D> +template<typename T, int D> bbox<T,D> bbox<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const { // Allow expansion only into directions where the extent is not negative // ASSERT_BBOX (all(lower()<=upper() or (lo==T(0) and hi==T(0)))); @@ -175,7 +175,7 @@ bbox<T,D> bbox<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const { } // Find the smallest b-compatible box around *this -template<class T, int D> +template<typename T, int D> bbox<T,D> bbox<T,D>::expanded_for (const bbox& b) const { if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride()); const vect<T,D> str = b.stride(); @@ -187,7 +187,7 @@ bbox<T,D> bbox<T,D>::expanded_for (const bbox& b) const { } // Find the largest b-compatible box inside *this -template<class T, int D> +template<typename T, int D> bbox<T,D> bbox<T,D>::contracted_for (const bbox& b) const { if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride()); const vect<T,D> str = b.stride(); @@ -199,7 +199,7 @@ bbox<T,D> bbox<T,D>::contracted_for (const bbox& b) const { } // Smallest bbox containing both boxes -template<class T, int D> +template<typename T, int D> bbox<T,D> bbox<T,D>::expanded_containing (const bbox& b) const { if (empty()) return b; if (b.empty()) return *this; @@ -211,18 +211,18 @@ bbox<T,D> bbox<T,D>::expanded_containing (const bbox& b) const { } // Iterators -template<class T, int D> +template<typename T, int D> bbox<T,D>::iterator::iterator (const bbox& box_, const vect<T,D>& pos_) : box(box_), pos(pos_) { if (box.empty()) pos=box.upper(); } -template<class T, int D> +template<typename T, int D> bool bbox<T,D>::iterator::operator!= (const iterator& i) const { return any(pos!=i.pos); } -template<class T, int D> +template<typename T, int D> typename bbox<T,D>::iterator& bbox<T,D>::iterator::operator++ () { for (int d=0; d<D; ++d) { pos[d]+=box.stride()[d]; @@ -232,12 +232,12 @@ typename bbox<T,D>::iterator& bbox<T,D>::iterator::operator++ () { return *this; } -template<class T, int D> +template<typename T, int D> typename bbox<T,D>::iterator bbox<T,D>::begin () const { return iterator(*this, lower()); } -template<class T, int D> +template<typename T, int D> typename bbox<T,D>::iterator bbox<T,D>::end () const { return iterator(*this, lower()); } @@ -245,7 +245,7 @@ typename bbox<T,D>::iterator bbox<T,D>::end () const { // Input -template<class T,int D> +template<typename T,int D> void bbox<T,D>::input (istream& is) { try { skipws (is); @@ -300,7 +300,7 @@ void bbox<T,D>::input (istream& is) { // Output -template<class T,int D> +template<typename T,int D> void bbox<T,D>::output (ostream& os) const { os << "(" << lower() << ":" << upper() << ":" << stride() << "/" << lower() / stride() << ":" << upper() / stride() diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh index 5b5029bdc..d2dac83cc 100644 --- a/Carpet/CarpetLib/src/bbox.hh +++ b/Carpet/CarpetLib/src/bbox.hh @@ -22,12 +22,12 @@ using namespace std; // Forward declaration -template<class T, int D> class bbox; +template<typename T, int D> class bbox; // Input/Output -template<class T, int D> +template<typename T, int D> istream& operator>> (istream& is, bbox<T,D>& b); -template<class T, int D> +template<typename T, int D> ostream& operator<< (ostream& os, const bbox<T,D>& b); @@ -36,7 +36,7 @@ ostream& operator<< (ostream& os, const bbox<T,D>& b); * A bounding box, i.e., a rectangle with lower and upper bound and a * stride. */ -template<class T, int D> +template<typename T, int D> class bbox { // Fields @@ -215,9 +215,9 @@ public: // Memory usage -template<class T, int D> +template<typename T, int D> inline size_t memoryof (bbox<T,D> const & b) CCTK_ATTRIBUTE_CONST; -template<class T, int D> +template<typename T, int D> inline size_t memoryof (bbox<T,D> const & b) { return b.memory(); } @@ -225,7 +225,7 @@ inline size_t memoryof (bbox<T,D> const & b) { return b.memory(); } // Input /** Read a formatted bbox from a stream. */ -template<class T,int D> +template<typename T,int D> inline istream& operator>> (istream& is, bbox<T,D>& b) { b.input(is); return is; @@ -236,7 +236,7 @@ inline istream& operator>> (istream& is, bbox<T,D>& b) { // Output /** Write a bbox formatted to a stream. */ -template<class T,int D> +template<typename T,int D> inline ostream& operator<< (ostream& os, const bbox<T,D>& b) { b.output(os); return os; diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc index a5748a7c6..7ffd61699 100644 --- a/Carpet/CarpetLib/src/bboxset.cc +++ b/Carpet/CarpetLib/src/bboxset.cc @@ -14,44 +14,62 @@ using namespace std; // Constructors -template<class T, int D> +template<typename T, int D> bboxset<T,D>::bboxset () { assert (invariant()); } -template<class T, int D> +template<typename T, int D> bboxset<T,D>::bboxset (const box& b) { //S if (!b.empty()) bs.insert(b); if (!b.empty()) bs.push_back(b); assert (invariant()); } -template<class T, int D> +template<typename T, int D> bboxset<T,D>::bboxset (const bboxset& s): bs(s.bs) { assert (invariant()); } -template<class T, int D> +template<typename T, int D> bboxset<T,D>::bboxset (const list<box>& lb) { for (typename list<box>::const_iterator li = lb.begin(); li != lb.end(); ++ li) { *this |= *li; } - normalize(); } -template<class T, int D> +template<typename T, int D> bboxset<T,D>::bboxset (const vector<list<box> >& vlb) { for (typename vector<list<box> >::const_iterator vli = vlb.begin(); vli != vlb.end(); ++ vli) { *this |= bboxset (*vli); } - normalize(); } -template<class T, int D> +template<typename T, int D> +template<typename U> +bboxset<T,D>::bboxset (const vector<U>& vb, const bbox<T,D> U::* const v) { + for (typename vector<U>::const_iterator + vi = vb.begin(); vi != vb.end(); ++ vi) + { + *this |= (*vi).*v; + } +} + +template<typename T, int D> +template<typename U> +bboxset<T,D>::bboxset (const vector<U>& vb, const bboxset U::* const v) { + for (typename vector<U>::const_iterator + vi = vb.begin(); vi != vb.end(); ++ vi) + { + *this |= (*vi).*v; + } +} + +template<typename T, int D> bboxset<T,D> bboxset<T,D>::poison () { return bboxset (bbox<T,D>::poison()); } @@ -59,7 +77,7 @@ bboxset<T,D> bboxset<T,D>::poison () { // Invariant -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::invariant () const { // This is very slow when there are many bboxes #if 0 && defined(CARPET_DEBUG) @@ -78,7 +96,7 @@ bool bboxset<T,D>::invariant () const { // Normalisation -template<class T, int D> +template<typename T, int D> void bboxset<T,D>::normalize () { assert (invariant()); @@ -200,14 +218,15 @@ void bboxset<T,D>::normalize () size_type const newsize = this->size(); - assert (*this == oldbs); + // Can't use operators on *this since these would call normalize again + // assert (*this == oldbs); assert (newsize <= oldsize); } // Accessors -template<class T, int D> +template<typename T, int D> typename bboxset<T,D>::size_type bboxset<T,D>::size () const { size_type s=0; for (const_iterator bi=begin(); bi!=end(); ++bi) { @@ -224,7 +243,7 @@ typename bboxset<T,D>::size_type bboxset<T,D>::size () const { // Queries // Intersection -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::intersects (const box& b) const { for (const_iterator bi=begin(); bi!=end(); ++bi) { if ((*bi).intersects(b)) return true; @@ -235,7 +254,7 @@ bool bboxset<T,D>::intersects (const box& b) const { // Add (bboxes that don't overlap) -template<class T, int D> +template<typename T, int D> bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s) { for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) { *this += *bi; @@ -244,14 +263,15 @@ bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s) { return *this; } -template<class T, int D> +template<typename T, int D> bboxset<T,D>& bboxset<T,D>::add_transfer (bboxset& s) { bs.splice (bs.end(), s.bs); assert (invariant()); + normalize(); return *this; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator+ (const box& b) const { bboxset r(*this); r += b; @@ -259,7 +279,7 @@ bboxset<T,D> bboxset<T,D>::operator+ (const box& b) const { return r; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator+ (const bboxset& s) const { bboxset r(*this); r += s; @@ -270,7 +290,7 @@ bboxset<T,D> bboxset<T,D>::operator+ (const bboxset& s) const { // Union -template<class T, int D> +template<typename T, int D> bboxset<T,D>& bboxset<T,D>::operator|= (const box& b) { bboxset tmp = b - *this; add_transfer (tmp); @@ -278,7 +298,7 @@ bboxset<T,D>& bboxset<T,D>::operator|= (const box& b) { return *this; } -template<class T, int D> +template<typename T, int D> bboxset<T,D>& bboxset<T,D>::operator|= (const bboxset& s) { bboxset tmp = s - *this; add_transfer (tmp); @@ -286,7 +306,7 @@ bboxset<T,D>& bboxset<T,D>::operator|= (const bboxset& s) { return *this; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator| (const box& b) const { bboxset r(*this); r |= b; @@ -294,7 +314,7 @@ bboxset<T,D> bboxset<T,D>::operator| (const box& b) const { return r; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator| (const bboxset& s) const { bboxset r(*this); r |= s; @@ -305,7 +325,7 @@ bboxset<T,D> bboxset<T,D>::operator| (const bboxset& s) const { // Intersection -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator& (const box& b) const { // start with an empty set bboxset r; @@ -318,7 +338,7 @@ bboxset<T,D> bboxset<T,D>::operator& (const box& b) const { return r; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator& (const bboxset& s) const { // start with an empty set bboxset r; @@ -332,14 +352,14 @@ bboxset<T,D> bboxset<T,D>::operator& (const bboxset& s) const { return r; } -template<class T, int D> +template<typename T, int D> bboxset<T,D>& bboxset<T,D>::operator&= (const box& b) { *this = *this & b; assert (invariant()); return *this; } -template<class T, int D> +template<typename T, int D> bboxset<T,D>& bboxset<T,D>::operator&= (const bboxset& s) { *this = *this & s; assert (invariant()); @@ -349,7 +369,7 @@ bboxset<T,D>& bboxset<T,D>::operator&= (const bboxset& s) { // Difference -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2) { assert (b1.is_aligned_with(b2)); if (b1.empty()) return bboxset<T,D>(); @@ -387,7 +407,7 @@ bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2) { return r; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator- (const box& b) const { // start with an empty set bboxset r; @@ -401,7 +421,7 @@ bboxset<T,D> bboxset<T,D>::operator- (const box& b) const { return r; } -template<class T, int D> +template<typename T, int D> bboxset<T,D>& bboxset<T,D>::operator-= (const bboxset& s) { for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) { *this -= *bi; @@ -410,7 +430,7 @@ bboxset<T,D>& bboxset<T,D>::operator-= (const bboxset& s) { return *this; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::operator- (const bboxset& s) const { bboxset r(*this); r -= s; @@ -418,7 +438,7 @@ bboxset<T,D> bboxset<T,D>::operator- (const bboxset& s) const { return r; } -template<class T, int D> +template<typename T, int D> bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b, const bboxset<T,D>& s) { bboxset<T,D> r = bboxset<T,D>(b) - s; assert (r.invariant()); @@ -427,33 +447,80 @@ bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b, const bboxset<T,D>& s) { +template<typename T, int D> +typename bboxset<T,D>::box bboxset<T,D>::container () const { + box b; + for (const_iterator bi=begin(); bi!=end(); ++bi) { + b = b.expanded_containing(*bi); + } + return b; +} + +template<typename T, int D> +bboxset<T,D> bboxset<T,D>::pseudo_inverse (const int n) const { + assert (not empty()); + box const c = container().expand(n,n); + return c - *this; +} + +template<typename T, int D> +bboxset<T,D> bboxset<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const { + // We don't know (yet?) how to shrink a set + assert (all (lo>=0 and hi>=0)); + bboxset res; + for (const_iterator bi=begin(); bi!=end(); ++bi) { + res |= (*bi).expand(lo,hi); + } + return res; +} + +template<typename T, int D> +bboxset<T,D> bboxset<T,D>::expanded_for (const box& b) const { + bboxset res; + for (const_iterator bi=begin(); bi!=end(); ++bi) { + res |= (*bi).expanded_for(b); + } + return res; +} + +template<typename T, int D> +bboxset<T,D> bboxset<T,D>::contracted_for (const box& b) const { + bboxset res; + for (const_iterator bi=begin(); bi!=end(); ++bi) { + res |= (*bi).contracted_for(b); + } + return res; +} + + + // Equality -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::operator<= (const bboxset<T,D>& s) const { return (*this - s).empty(); } -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::operator< (const bboxset<T,D>& s) const { return (*this - s).empty() && ! (s - *this).empty(); } -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::operator>= (const bboxset<T,D>& s) const { return s <= *this; } -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::operator> (const bboxset<T,D>& s) const { return s < *this; } -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::operator== (const bboxset<T,D>& s) const { return (*this <= s) && (*this >= s); } -template<class T, int D> +template<typename T, int D> bool bboxset<T,D>::operator!= (const bboxset<T,D>& s) const { return ! (*this == s); } @@ -461,7 +528,7 @@ bool bboxset<T,D>::operator!= (const bboxset<T,D>& s) const { // Input -template<class T,int D> +template<typename T,int D> istream& bboxset<T,D>::input (istream& is) { T Tdummy; try { @@ -491,13 +558,14 @@ istream& bboxset<T,D>::input (istream& is) { cout << "Input error while reading a bboxset<" << typestring(Tdummy) << "," << D << ">" << endl; throw err; } + normalize(); return is; } // Output -template<class T,int D> +template<typename T,int D> ostream& bboxset<T,D>::output (ostream& os) const { T Tdummy; os << "bboxset<" << typestring(Tdummy) << "," << D << ">:{" @@ -511,3 +579,13 @@ ostream& bboxset<T,D>::output (ostream& os) const { template class bboxset<int,dim>; +template size_t memoryof (const bboxset<int,dim>& s); +template istream& operator>> (istream& is, bboxset<int,dim>& s); +template ostream& operator<< (ostream& os, const bboxset<int,dim>& s); + +#include "dh.hh" +#include "region.hh" + +template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bbox<int,dim> dh::full_dboxes::* const v); +template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bboxset<int,dim> dh::full_dboxes::* const v); +template bboxset<int,dim>::bboxset (const vector<region_t>& vb, const bbox<int,dim> region_t::* const v); diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh index 5d206da33..5624ef4b1 100644 --- a/Carpet/CarpetLib/src/bboxset.hh +++ b/Carpet/CarpetLib/src/bboxset.hh @@ -16,30 +16,30 @@ using namespace std; // Forward declaration -template<class T, int D> class bboxset; +template<typename T, int D> class bboxset; -// template<class T,int D> +// template<typename T,int D> // bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2); -// template<class T,int D> +// template<typename T,int D> // bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s); -// template<class T,int D> +// template<typename T,int D> // bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2); -// template<class T,int D> +// template<typename T,int D> // bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s); // Input -template<class T,int D> +template<typename T,int D> istream& operator>> (istream& is, bboxset<T,D>& s); // Output -template<class T,int D> +template<typename T,int D> ostream& operator<< (ostream& os, const bboxset<T,D>& s); // Bounding box class -template<class T, int D> +template<typename T, int D> class bboxset { // Types @@ -63,15 +63,21 @@ public: bboxset (const list<box>& lb); bboxset (const vector<list<box> >& vlb); + template<typename U> + bboxset (const vector<U>& vb, const bbox<T,D> U::* const v); + template<typename U> + bboxset (const vector<U>& vb, const bboxset U::* const v); static bboxset poison (); // Invariant bool invariant () const; +private: // Normalisation void normalize (); +public: // Accessors bool empty () const { return bs.empty(); } // T size () const; @@ -131,6 +137,23 @@ public: // friend bboxset operator- <T,D>(const box& b, const bboxset& s); static bboxset minus (const box& b, const bboxset& s); + /** Find a bbox containing the whole set. */ + box container () const; + /** Find the pseudo-inverse. */ + bboxset pseudo_inverse (const int n) const; + + /** Expand (enlarge) the bbox by multiples of the stride. */ + bboxset expand (const vect<T,D>& lo, const vect<T,D>& hi) const; + bboxset expand (const vect<vect<T,D>,2>& lohi) const + { return expand (lohi[0], lohi[1]); } + + /** Find the smallest b-compatible box around this bbox. + ("compatible" means having the same stride.) */ + bboxset expanded_for (const box& b) const; + + /** Find the largest b-compatible box inside this bbox. */ + bboxset contracted_for (const box& b) const; + // Equality bool operator== (const bboxset& s) const; bool operator!= (const bboxset& s) const; @@ -160,89 +183,89 @@ public: -template<class T,int D> +template<typename T,int D> inline bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2) { return bboxset<T,D>(b1) + bboxset<T,D>(b2); } -template<class T,int D> +template<typename T,int D> inline bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>(b) + s; } -template<class T,int D> +template<typename T,int D> inline bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) { return bboxset<T,D>::minus(b1,b2); } -template<class T,int D> +template<typename T,int D> inline bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>::minus(b,s); } -template<class T,int D> +template<typename T,int D> inline bboxset<T,D> operator| (const bbox<T,D>& b, const bboxset<T,D>& s) { return s | b; } -template<class T,int D> +template<typename T,int D> inline bboxset<T,D> operator& (const bbox<T,D>& b, const bboxset<T,D>& s) { return s & b; } -template<class T,int D> +template<typename T,int D> inline bool operator== (const bbox<T,D>& b, const bboxset<T,D>& s) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator== (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>(b) == s; } -template<class T,int D> +template<typename T,int D> inline bool operator!= (const bbox<T,D>& b, const bboxset<T,D>& s) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator!= (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>(b) != s; } -template<class T,int D> +template<typename T,int D> inline bool operator< (const bbox<T,D>& b, const bboxset<T,D>& s) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator< (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>(b) < s; } -template<class T,int D> +template<typename T,int D> inline bool operator<= (const bbox<T,D>& b, const bboxset<T,D>& s) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator<= (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>(b) <= s; } -template<class T,int D> +template<typename T,int D> inline bool operator> (const bbox<T,D>& b, const bboxset<T,D>& s) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator> (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>(b) > s; } -template<class T,int D> +template<typename T,int D> inline bool operator>= (const bbox<T,D>& b, const bboxset<T,D>& s) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator>= (const bbox<T,D>& b, const bboxset<T,D>& s) { return bboxset<T,D>(b) >= s; @@ -250,55 +273,55 @@ inline bool operator>= (const bbox<T,D>& b, const bboxset<T,D>& s) -template<class T,int D> +template<typename T,int D> inline bool operator== (const bboxset<T,D>& s, const bbox<T,D>& b) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator== (const bboxset<T,D>& s, const bbox<T,D>& b) { return s == bboxset<T,D>(b); } -template<class T,int D> +template<typename T,int D> inline bool operator!= (const bboxset<T,D>& s, const bbox<T,D>& b) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator!= (const bboxset<T,D>& s, const bbox<T,D>& b) { return s != bboxset<T,D>(b); } -template<class T,int D> +template<typename T,int D> inline bool operator< (const bboxset<T,D>& s, const bbox<T,D>& b) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator< (const bboxset<T,D>& s, const bbox<T,D>& b) { return s < bboxset<T,D>(b); } -template<class T,int D> +template<typename T,int D> inline bool operator<= (const bboxset<T,D>& s, const bbox<T,D>& b) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator<= (const bboxset<T,D>& s, const bbox<T,D>& b) { return s <= bboxset<T,D>(b); } -template<class T,int D> +template<typename T,int D> inline bool operator> (const bboxset<T,D>& s, const bbox<T,D>& b) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator> (const bboxset<T,D>& s, const bbox<T,D>& b) { return s > bboxset<T,D>(b); } -template<class T,int D> +template<typename T,int D> inline bool operator>= (const bboxset<T,D>& s, const bbox<T,D>& b) CCTK_ATTRIBUTE_PURE; -template<class T,int D> +template<typename T,int D> inline bool operator>= (const bboxset<T,D>& s, const bbox<T,D>& b) { return s >= bboxset<T,D>(b); @@ -307,17 +330,17 @@ inline bool operator>= (const bboxset<T,D>& s, const bbox<T,D>& b) // Memory usage -template<class T, int D> +template<typename T, int D> inline size_t memoryof (bboxset<T,D> const & s) CCTK_ATTRIBUTE_PURE; -template<class T, int D> +template<typename T, int D> inline size_t memoryof (bboxset<T,D> const & s) { return s.memory(); } // Input -template<class T,int D> +template<typename T,int D> inline istream& operator>> (istream& is, bboxset<T,D>& s) { return s.input(is); } @@ -325,7 +348,7 @@ inline istream& operator>> (istream& is, bboxset<T,D>& s) { // Output -template<class T,int D> +template<typename T,int D> inline ostream& operator<< (ostream& os, const bboxset<T,D>& s) { return s.output(os); } diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 7b6f4e0a8..484582073 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -94,7 +94,6 @@ call_operator (void } # if ! defined (NDEBUG) && ! defined (CARPET_OPTIMISE) if (not (allregbboxes == ibset (regbbox))) { - allregbboxes.normalize(); cout << "allregbboxes=" << allregbboxes << endl << "regbbox=" << regbbox << endl; } @@ -163,7 +162,6 @@ call_operator (void } # if ! defined (NDEBUG) && ! defined (CARPET_OPTIMISE) if (not (allregbboxes == ibset (regbbox))) { - allregbboxes.normalize(); cout << "allregbboxes=" << allregbboxes << endl << "regbbox=" << regbbox << endl; } diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc index a2b366ea2..f620a0b0d 100644 --- a/Carpet/CarpetLib/src/defs.cc +++ b/Carpet/CarpetLib/src/defs.cc @@ -330,6 +330,8 @@ template CCTK_REAL ipow (CCTK_REAL x, int y); template vect<int,dim> ipow (vect<int,dim> x, int y); template vect<CCTK_REAL,dim> ipow (vect<CCTK_REAL,dim> x, int y); +template bool equals (vector<ibset> const& v, vector<ibset> const& w); + template size_t memoryof (list<ibbox> const & l); template size_t memoryof (list<ivect> const & l); template size_t memoryof (list<dh*> const & l); @@ -342,6 +344,7 @@ template size_t memoryof (vector<bool> const & v); template size_t memoryof (vector<int> const & v); template size_t memoryof (vector<CCTK_REAL> const & v); template size_t memoryof (vector<ibbox> const & v); +template size_t memoryof (vector<ibset> const & v); template size_t memoryof (vector<ivect> const & v); template size_t memoryof (vector<i2vect> const & v); template size_t memoryof (vector<fulltree <int,dim,pseudoregion_t> *> const & f); @@ -351,12 +354,16 @@ template size_t memoryof (vector<sendrecv_pseudoregion_t> const & v); template size_t memoryof (vector<vector<int> > const & v); template size_t memoryof (vector<vector<CCTK_REAL> > const & v); template size_t memoryof (vector<vector<ibbox> > const & v); -template size_t memoryof (vector<vector<dh::dboxes> > const & v); template size_t memoryof (vector<vector<dh::fast_dboxes> > const & v); +template size_t memoryof (vector<vector<dh::full_dboxes> > const & v); +template size_t memoryof (vector<vector<dh::light_dboxes> > const & v); +template size_t memoryof (vector<vector<dh::local_dboxes> > const & v); template size_t memoryof (vector<vector<region_t> > const & v); template size_t memoryof (vector<vector<vector<CCTK_REAL> > > const & v); -template size_t memoryof (vector<vector<vector<dh::dboxes> > > const & v); template size_t memoryof (vector<vector<vector<dh::fast_dboxes> > > const & v); +template size_t memoryof (vector<vector<vector<dh::full_dboxes> > > const & v); +template size_t memoryof (vector<vector<vector<dh::light_dboxes> > > const & v); +template size_t memoryof (vector<vector<vector<dh::local_dboxes> > > const & v); template size_t memoryof (vector<vector<vector<region_t> > > const & v); template size_t memoryof (vector<vector<vector<gdata*> > > const & v); template size_t memoryof (vector<vector<vector<vector<gdata*> > > > const & v); @@ -367,6 +374,7 @@ template istream& input (istream& os, vector<int>& v); template istream& input (istream& os, vector<CCTK_REAL>& v); template istream& input (istream& os, vector<ibbox>& v); template istream& input (istream& os, vector<rbbox>& v); +template istream& input (istream& os, vector<ibset>& v); template istream& input (istream& os, vector<ivect>& v); template istream& input (istream& os, vector<bbvect>& v); template istream& input (istream& os, vector<i2vect>& v); @@ -394,12 +402,15 @@ template ostream& output (ostream& os, const vector<int>& v); template ostream& output (ostream& os, const vector<CCTK_REAL>& v); template ostream& output (ostream& os, const vector<ibbox>& v); template ostream& output (ostream& os, const vector<rbbox>& v); +template ostream& output (ostream& os, const vector<ibset>& v); template ostream& output (ostream& os, const vector<ivect>& v); template ostream& output (ostream& os, const vector<rvect>& v); template ostream& output (ostream& os, const vector<bbvect>& v); template ostream& output (ostream& os, const vector<i2vect>& v); -template ostream& output (ostream& os, const vector<dh::dboxes> & v); template ostream& output (ostream& os, const vector<dh::fast_dboxes> & v); +template ostream& output (ostream& os, const vector<dh::full_dboxes> & v); +template ostream& output (ostream& os, const vector<dh::light_dboxes> & v); +template ostream& output (ostream& os, const vector<dh::local_dboxes> & v); template ostream& output (ostream& os, const vector<region_t>& v); template ostream& output (ostream& os, const vector<pseudoregion_t>& v); template ostream& output (ostream& os, const vector<sendrecv_pseudoregion_t>& v); @@ -410,11 +421,15 @@ template ostream& output (ostream& os, const vector<vector<ibbox> >& v); template ostream& output (ostream& os, const vector<vector<rbbox> >& v); template ostream& output (ostream& os, const vector<vector<bbvect> >& v); template ostream& output (ostream& os, const vector<vector<i2vect> >& v); -template ostream& output (ostream& os, const vector<vector<dh::dboxes> > & b); template ostream& output (ostream& os, const vector<vector<dh::fast_dboxes> > & b); +template ostream& output (ostream& os, const vector<vector<dh::full_dboxes> > & b); +template ostream& output (ostream& os, const vector<vector<dh::light_dboxes> > & b); +template ostream& output (ostream& os, const vector<vector<dh::local_dboxes> > & b); template ostream& output (ostream& os, const vector<vector<region_t> >& v); template ostream& output (ostream& os, const vector<vector<vector<CCTK_REAL> > >& v); template ostream& output (ostream& os, const vector<vector<vector<ibbox> > >& v); -template ostream& output (ostream& os, const vector<vector<vector<dh::dboxes> > > & b); template ostream& output (ostream& os, const vector<vector<vector<dh::fast_dboxes> > > & b); +template ostream& output (ostream& os, const vector<vector<vector<dh::full_dboxes> > > & b); +template ostream& output (ostream& os, const vector<vector<vector<dh::light_dboxes> > > & b); +template ostream& output (ostream& os, const vector<vector<vector<dh::local_dboxes> > > & b); template ostream& output (ostream& os, const vector<vector<vector<region_t> > >& v); diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh index d4ff68505..65a1edce5 100644 --- a/Carpet/CarpetLib/src/defs.hh +++ b/Carpet/CarpetLib/src/defs.hh @@ -103,6 +103,8 @@ typedef vect<vect<CCTK_INT,2>,dim> jjvect; typedef vect<vect<bool,dim>,2> b2vect; typedef vect<vect<int,dim>,2> i2vect; +typedef vect<vect<CCTK_INT,dim>,2> j2vect; +typedef vect<vect<CCTK_REAL,dim>,2> r2vect; @@ -212,23 +214,23 @@ inline const char * typestring (const CCTK_COMPLEX32&) // Capture the system's fpclassify, isfinite, isinf, isnan, and // isnormal functions -#ifdef HAVE_CCTK_REAL4 -inline int myfpclassify (CCTK_REAL4 const & x) CCTK_ATTRIBUTE_CONST; -int myfpclassify (CCTK_REAL4 const & x) -{ return fpclassify (x); } -#endif -#ifdef HAVE_CCTK_REAL8 -inline int myfpclassify (CCTK_REAL8 const & x) CCTK_ATTRIBUTE_CONST; -inline int myfpclassify (CCTK_REAL8 const & x) -{ return fpclassify (x); } -#endif -#ifdef HAVE_CCTK_REAL16 -inline int myfpclassify (CCTK_REAL16 const & x) CCTK_ATTRIBUTE_CONST; -inline int myfpclassify (CCTK_REAL16 const & x) -{ return fpclassify (x); } -#endif - -#undef fpclassify +// #ifdef HAVE_CCTK_REAL4 +// inline int myfpclassify (CCTK_REAL4 const & x) CCTK_ATTRIBUTE_CONST; +// int myfpclassify (CCTK_REAL4 const & x) +// { return fpclassify (x); } +// #endif +// #ifdef HAVE_CCTK_REAL8 +// inline int myfpclassify (CCTK_REAL8 const & x) CCTK_ATTRIBUTE_CONST; +// inline int myfpclassify (CCTK_REAL8 const & x) +// { return fpclassify (x); } +// #endif +// #ifdef HAVE_CCTK_REAL16 +// inline int myfpclassify (CCTK_REAL16 const & x) CCTK_ATTRIBUTE_CONST; +// inline int myfpclassify (CCTK_REAL16 const & x) +// { return fpclassify (x); } +// #endif +// +// #undef fpclassify #ifdef HAVE_CCTK_REAL4 inline int myisfinite (CCTK_REAL4 const & x) CCTK_ATTRIBUTE_CONST; @@ -363,48 +365,48 @@ namespace CarpetLib { { return CCTK_Cmplx32Abs (x); } #endif - // - // fpclassify - // - - // Default implementation, only good for integers - template <typename T> - inline int fpclassify (T const & x) CCTK_ATTRIBUTE_CONST; - template <typename T> - inline int fpclassify (T const & x) - { return x ? FP_NORMAL : FP_ZERO; } - -#ifdef HAVE_CCTK_REAL4 - template<> inline int fpclassify (CCTK_REAL4 const & x) CCTK_ATTRIBUTE_CONST; - template<> inline int fpclassify (CCTK_REAL4 const & x) - { return myfpclassify (x); } -#endif -#ifdef HAVE_CCTK_REAL8 - template<> inline int fpclassify (CCTK_REAL8 const & x) CCTK_ATTRIBUTE_CONST; - template<> inline int fpclassify (CCTK_REAL8 const & x) - { return myfpclassify (x); } -#endif -#ifdef HAVE_CCTK_REAL16 - template<> inline int fpclassify (CCTK_REAL16 const & x) CCTK_ATTRIBUTE_CONST; - template<> inline int fpclassify (CCTK_REAL16 const & x) - { return myfpclassify (x); } -#endif - -#ifdef HAVE_CCTK_COMPLEX8 - template<> inline int fpclassify (CCTK_COMPLEX8 const & x) CCTK_ATTRIBUTE_CONST; - template<> inline int fpclassify (CCTK_COMPLEX8 const & x) - { assert(0); } -#endif -#ifdef HAVE_CCTK_COMPLEX16 - template<> inline int fpclassify (CCTK_COMPLEX16 const & x) CCTK_ATTRIBUTE_CONST; - template<> inline int fpclassify (CCTK_COMPLEX16 const & x) - { assert(0); } -#endif -#ifdef HAVE_CCTK_COMPLEX32 - template<> inline int fpclassify (CCTK_COMPLEX32 const & x) CCTK_ATTRIBUTE_CONST; - template<> inline int fpclassify (CCTK_COMPLEX32 const & x) - { assert(0); } -#endif +// // +// // fpclassify +// // +// +// // Default implementation, only good for integers +// template <typename T> +// inline int fpclassify (T const & x) CCTK_ATTRIBUTE_CONST; +// template <typename T> +// inline int fpclassify (T const & x) +// { return x ? FP_NORMAL : FP_ZERO; } +// +// #ifdef HAVE_CCTK_REAL4 +// template<> inline int fpclassify (CCTK_REAL4 const & x) CCTK_ATTRIBUTE_CONST; +// template<> inline int fpclassify (CCTK_REAL4 const & x) +// { return myfpclassify (x); } +// #endif +// #ifdef HAVE_CCTK_REAL8 +// template<> inline int fpclassify (CCTK_REAL8 const & x) CCTK_ATTRIBUTE_CONST; +// template<> inline int fpclassify (CCTK_REAL8 const & x) +// { return myfpclassify (x); } +// #endif +// #ifdef HAVE_CCTK_REAL16 +// template<> inline int fpclassify (CCTK_REAL16 const & x) CCTK_ATTRIBUTE_CONST; +// template<> inline int fpclassify (CCTK_REAL16 const & x) +// { return myfpclassify (x); } +// #endif +// +// #ifdef HAVE_CCTK_COMPLEX8 +// template<> inline int fpclassify (CCTK_COMPLEX8 const & x) CCTK_ATTRIBUTE_CONST; +// template<> inline int fpclassify (CCTK_COMPLEX8 const & x) +// { assert(0); } +// #endif +// #ifdef HAVE_CCTK_COMPLEX16 +// template<> inline int fpclassify (CCTK_COMPLEX16 const & x) CCTK_ATTRIBUTE_CONST; +// template<> inline int fpclassify (CCTK_COMPLEX16 const & x) +// { assert(0); } +// #endif +// #ifdef HAVE_CCTK_COMPLEX32 +// template<> inline int fpclassify (CCTK_COMPLEX32 const & x) CCTK_ATTRIBUTE_CONST; +// template<> inline int fpclassify (CCTK_COMPLEX32 const & x) +// { assert(0); } +// #endif // // isfinite @@ -583,6 +585,52 @@ namespace CarpetLib { +// Use #defines to access the CarpetLib::good functions from the +// global namespace. Because macros are not scoped, this makes the +// corresponding functions in CarpetLib::good inaccessible, so we +// rename them by adding an underscore. + +namespace CarpetLib { + namespace good { + +#define ALIAS_GOOD_FUNCTION(typ,func) \ + template <typename T> \ + inline typ func##_ (T const & x) CCTK_ATTRIBUTE_CONST; \ + template <typename T> \ + inline typ func##_ (T const & x) \ + { return func (x); } + + // ALIAS_GOOD_FUNCTION(int,fpclassify) + ALIAS_GOOD_FUNCTION(int,isfinite) + ALIAS_GOOD_FUNCTION(int,isinf) + ALIAS_GOOD_FUNCTION(int,isnan) + ALIAS_GOOD_FUNCTION(int,isnormal) + +#undef ALIAS_GOOD_FUNCTION + + } // namespace good +} // namespace CarpetLib + +#define isfinite (::CarpetLib::good::isfinite_) +#define isinf (::CarpetLib::good::isinf_) +#define isnan (::CarpetLib::good::isnan_) +#define isnormal (::CarpetLib::good::isnormal_) + + + +// Container equality +template <typename T> +bool equals (vector<T> const& v, vector<T> const& w) +{ + if (v.size() != w.size()) return false; + for (size_t i=0; i<v.size(); ++i) { + if (v.AT(i) != w.AT(i)) return false; + } + return true; +} + + + // Container memory usage inline size_t memoryof (char const e) CCTK_ATTRIBUTE_CONST; inline size_t memoryof (short const e) CCTK_ATTRIBUTE_CONST; diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index d16fde8eb..9679b3fd7 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -1,8 +1,9 @@ #include <cassert> #include <cstddef> +#include <sstream> -#include "cctk.h" -#include "cctk_Parameters.h" +#include <cctk.h> +#include <cctk_Parameters.h> #include "CarpetTimers.hh" @@ -120,36 +121,39 @@ bool there_was_an_error = false; static void assert_error (char const * restrict const checkstring, + char const * restrict const file, int const line, int const ml, int const rl, char const * restrict const message) { CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, - "[ml=%d rl=%d] The following grid structure consistency check failed:\n %s\n %s", - ml, rl, message, checkstring); + "\n%s:%d:\n [ml=%d rl=%d] The following grid structure consistency check failed:\n %s\n %s", + file, line, ml, rl, message, checkstring); there_was_an_error = true; } static void assert_error (char const * restrict const checkstring, + char const * restrict const file, int const line, int const ml, int const rl, int const c, char const * restrict const message) { CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, - "[ml=%d rl=%d c=%d] The following grid structure consistency check failed:\n %s\n %s", - ml, rl, c, message, checkstring); + "\n%s:%d:\n [ml=%d rl=%d c=%d] The following grid structure consistency check failed:\n %s\n %s", + file, line, ml, rl, c, message, checkstring); there_was_an_error = true; } static void assert_error (char const * restrict const checkstring, + char const * restrict const file, int const line, int const ml, int const rl, int const c, int const cc, char const * restrict const message) { CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, - "[ml=%d rl=%d c=%d cc=%d] The following grid structure consistency check failed:\n %s\n %s", - ml, rl, c, cc, message, checkstring); + "\n%s:%d:\n [ml=%d rl=%d c=%d cc=%d] The following grid structure consistency check failed:\n %s\n %s", + file, line, ml, rl, c, cc, message, checkstring); there_was_an_error = true; } @@ -162,25 +166,25 @@ assert_error (char const * restrict const checkstring, #else -#define ASSERT_rl(check, message) \ - do { \ - if (not (check)) { \ - assert_error (#check, ml, rl, message); \ - } \ +#define ASSERT_rl(check, message) \ + do { \ + if (not (check)) { \ + assert_error (#check, __FILE__, __LINE__, ml, rl, message); \ + } \ } while (false) -#define ASSERT_c(check, message) \ - do { \ - if (not (check)) { \ - assert_error (#check, ml, rl, c, message); \ - } \ +#define ASSERT_c(check, message) \ + do { \ + if (not (check)) { \ + assert_error (#check, __FILE__, __LINE__, ml, rl, c, message); \ + } \ } while (false) -#define ASSERT_cc(check, message) \ - do { \ - if (not (check)) { \ - assert_error (#check, ml, rl, c, cc, message); \ - } \ +#define ASSERT_cc(check, message) \ + do { \ + if (not (check)) { \ + assert_error (#check, __FILE__, __LINE__, ml, rl, c, cc, message); \ + } \ } while (false) #endif @@ -201,8 +205,8 @@ regrid (bool const do_init) static Timer total ("CarpetLib::dh::regrid"); total.start (); - mboxes oldboxes; - swap (boxes, oldboxes); + light_mboxes old_light_boxes; + swap (light_boxes, old_light_boxes); full_mboxes full_boxes; @@ -210,21 +214,22 @@ regrid (bool const do_init) - // cerr << "QQQ: regrid[1]" << endl; - boxes.resize (h.mglevels()); + light_boxes.resize (h.mglevels()); + local_boxes.resize (h.mglevels()); full_boxes.resize (h.mglevels()); fast_boxes.resize (h.mglevels()); for (int ml = 0; ml < h.mglevels(); ++ ml) { - // cerr << "QQQ: regrid[2] ml=" << ml << endl; - boxes.AT(ml).resize (h.reflevels()); + light_boxes.AT(ml).resize (h.reflevels()); + local_boxes.AT(ml).resize (h.reflevels()); full_boxes.AT(ml).resize (h.reflevels()); fast_boxes.AT(ml).resize (h.reflevels()); for (int rl = 0; rl < h.reflevels(); ++ rl) { - // cerr << "QQQ: regrid[3] rl=" << rl << endl; - boxes.AT(ml).AT(rl).resize (h.components(rl)); + light_boxes.AT(ml).AT(rl).resize (h.components(rl)); + local_boxes.AT(ml).AT(rl).resize (h.local_components(rl)); full_boxes.AT(ml).AT(rl).resize (h.components(rl)); - cboxes & level = boxes.AT(ml).AT(rl); + light_cboxes & light_level = light_boxes.AT(ml).AT(rl); + local_cboxes & local_level = local_boxes.AT(ml).AT(rl); full_cboxes & full_level = full_boxes.AT(ml).AT(rl); fast_dboxes & fast_level = fast_boxes.AT(ml).AT(rl); @@ -238,7 +243,6 @@ regrid (bool const do_init) // Domain: - // cerr << "QQQ: regrid[a]" << endl; static Carpet::Timer timer_domain ("CarpetLib::dh::regrid::domain"); timer_domain.start(); @@ -259,8 +263,7 @@ regrid (bool const do_init) ASSERT_rl (domain_active <= domain_exterior, "The active part of the domain must be contained in the exterior part of the domain"); - ibset domain_boundary = domain_exterior - domain_active; - domain_boundary.normalize(); + ibset const domain_boundary = domain_exterior - domain_active; timer_domain.stop(); @@ -269,7 +272,6 @@ regrid (bool const do_init) static Carpet::Timer timer_region ("CarpetLib::dh::regrid::region"); timer_region.start(); - // cerr << "QQQ: regrid[b]" << endl; for (int c = 0; c < h.components(rl); ++ c) { full_dboxes & box = full_level.AT(c); @@ -330,7 +332,7 @@ regrid (bool const do_init) // (The exterior must not be empty) // Variables may have size zero // ASSERT_c (not extr.empty(), - // "The experior must not be empty"); + // "The exterior must not be empty"); // The exterior must be contained in the domain ASSERT_c (extr <= domain_exterior, @@ -344,7 +346,6 @@ regrid (bool const do_init) ghosts = ibset::poison(); ghosts = extr - intr; - ghosts.normalize(); // The ghosts must be contained in the domain. Different from // the boundaries, the ghost can include part of the outer @@ -379,7 +380,6 @@ regrid (bool const do_init) outer_boundaries = ibset::poison(); outer_boundaries = extr - comm; - outer_boundaries.normalize(); // The outer boundary must be contained in the outer boundary // of the domain @@ -422,7 +422,6 @@ regrid (bool const do_init) boundaries = ibset::poison(); boundaries = comm - owned; - boundaries.normalize(); // The boundary must be contained in the active part of the // domain. This prevents that a region is too close to the @@ -438,7 +437,6 @@ regrid (bool const do_init) // Conjunction of all buffer zones: - // cerr << "QQQ: regrid[c]" << endl; static Carpet::Timer timer_buffers ("CarpetLib::dh::regrid::buffers"); timer_buffers.start(); @@ -448,37 +446,67 @@ regrid (bool const do_init) ibbox const domain_enlarged = domain_active.expand (safedist); // All owned regions - ibset allowned; - for (int c = 0; c < h.components(rl); ++ c) { - full_dboxes const & box = full_level.AT(c); - allowned += box.owned; - } - allowned.normalize(); + ibset const allowned (full_level, & full_dboxes::owned); ASSERT_rl (allowned <= domain_active, "The owned regions must be contained in the active part of the domain"); // All not-owned regions - ibset notowned = domain_enlarged - allowned; - notowned.normalize(); + ibset const notowned = domain_enlarged - allowned; // All not-active points - ibset notactive; - for (ibset::const_iterator - ri = notowned.begin(); ri != notowned.end(); ++ ri) - { - ibbox const & r = * ri; - ibbox const r_enlarged = r.expand (buffer_width); - notactive |= r_enlarged; + ibset const notactive = notowned.expand (buffer_width); + + // All not-active points, in stages + int const num_substeps = + any (any (ghost_width == 0)) ? + 0 : + minval (minval (buffer_width / ghost_width)); + if (not all (all (buffer_width == num_substeps * ghost_width))) { + ostringstream buf; + buf << "The buffer width " << buffer_width << " is not a multiple of the ghost width " << ghost_width << " on level " << rl; + CCTK_WARN (CCTK_WARN_COMPLAIN, buf.str().c_str()); + } + + vector<ibset> notactive_stepped (num_substeps+1); + notactive_stepped.AT(0) = notowned; + for (int substep = 1; substep <= num_substeps; ++ substep) { + notactive_stepped.AT(substep) = + notactive_stepped.AT(substep-1).expand (ghost_width); + } + if (all (all (buffer_width == num_substeps * ghost_width))) { + ASSERT_rl (notactive_stepped.AT(num_substeps) == notactive, + "The stepped not-active region must be equal to the not-active region"); } - notactive.normalize(); // All buffer zones - ibset allbuffers = allowned & notactive; - allbuffers.normalize(); + ibset const allbuffers = allowned & notactive; + + // All active points + ibset const allactive = allowned - notactive; + + // All stepped buffer zones + vector<ibset> allbuffers_stepped (num_substeps); + ibset allbuffers_stepped_combined; + for (int substep = 0; substep < num_substeps; ++ substep) { + allbuffers_stepped.AT(substep) = + allowned & + (notactive_stepped.AT(substep+1) - notactive_stepped.AT(substep)); + allbuffers_stepped_combined += allbuffers_stepped.AT(substep); + } + if (all (all (buffer_width == num_substeps * ghost_width))) { + ASSERT_rl (allbuffers_stepped_combined == allbuffers, + "The stepped buffer zones must be equal to the buffer zones"); + } // Buffer zones must be in the active part of the domain + ASSERT_rl (allactive <= domain_active, + "The active region must be in the active part of the domain"); ASSERT_rl (allbuffers <= domain_active, "The buffer zones must be in the active part of the domain"); + ASSERT_rl ((allactive & allbuffers).empty(), + "The active points and the buffer zones cannot overlap"); + ASSERT_rl (allactive + allbuffers == allowned, + "The active points and the buffer points together must be exactly the owned region"); @@ -487,25 +515,34 @@ regrid (bool const do_init) // Buffer zones: box.buffers = box.owned & allbuffers; - box.buffers.normalize(); // Active region: - box.active = box.owned - box.buffers; - box.active.normalize(); - + box.active = box.owned & allactive; + ASSERT_c (box.active == box.owned - box.buffers, + "The active region must equal the owned region minus the buffer zones"); } // for c + for (int lc = 0; lc < h.local_components(rl); ++ lc) { + int const c = h.get_component (rl, lc); + local_dboxes & local_box = local_level.AT(lc); + full_dboxes const& box = full_level.AT(c); + + // Stepped buffer zones: + local_box.buffers = box.buffers; + + local_box.buffers_stepped.resize (num_substeps); + for (int substep = 0; substep < num_substeps; ++ substep) { + local_box.buffers_stepped.AT(substep) = + box.owned & allbuffers_stepped.AT(substep); + } + + local_box.active = box.active; + } // for lc + // The conjunction of all buffer zones must equal allbuffers - // cerr << "QQQ: regrid[d]" << endl; - - ibset allbuffers1; - for (int c = 0; c < h.components(rl); ++ c) { - full_dboxes const & box = full_level.AT(c); - allbuffers1 += box.buffers; - } - allbuffers1.normalize(); + ibset const allbuffers1 (full_level, & full_dboxes::buffers); ASSERT_rl (allbuffers1 == allbuffers, "Buffer zone consistency check"); @@ -514,7 +551,6 @@ regrid (bool const do_init) // Test constituency relations: - // cerr << "QQQ: regrid[e]" << endl; static Carpet::Timer timer_test ("CarpetLib::dh::regrid::test"); timer_test.start(); @@ -548,18 +584,16 @@ regrid (bool const do_init) } // for c timer_test.stop(); - + // Communication schedule: - // cerr << "QQQ: regrid[4]" << endl; static Carpet::Timer timer_comm ("CarpetLib::dh::regrid::comm"); timer_comm.start(); for (int lc = 0; lc < h.local_components(rl); ++ lc) { int const c = h.get_component (rl, lc); - // cerr << "QQQ: regrid[4a] lc=" << lc << " c=" << c << endl; full_dboxes & box = full_level.AT(c); @@ -580,17 +614,9 @@ regrid (bool const do_init) ibset needrecv = box.active; - ibset contracted_oactive; - for (ibset::const_iterator - ai = obox.active.begin(); ai != obox.active.end(); ++ ai) - { - ibbox const & oactive = * ai; - contracted_oactive += oactive.contracted_for (box.interior); - } - contracted_oactive.normalize(); - - ibset ovlp = needrecv & contracted_oactive; - ovlp.normalize(); + ibset const contracted_oactive + (obox.active.contracted_for (box.interior)); + ibset const ovlp = needrecv & contracted_oactive; for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) @@ -604,7 +630,6 @@ regrid (bool const do_init) } needrecv -= ovlp; - needrecv.normalize(); // All points must have been received ASSERT_c (needrecv.empty(), @@ -617,7 +642,6 @@ regrid (bool const do_init) // Multigrid prolongation: - // cerr << "QQQ: regrid[f]" << endl; static Carpet::Timer timer_comm_mgprol ("CarpetLib::dh::regrid::comm::mprol"); @@ -635,17 +659,8 @@ regrid (bool const do_init) i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); - ibset expanded_active; - for (ibset::const_iterator - ai = box.active.begin(); ai != box.active.end(); ++ ai) - { - ibbox const & active = * ai; - expanded_active += active.expanded_for (obox.interior); - } - expanded_active.normalize(); - - ibset ovlp = oneedrecv & expanded_active; - ovlp.normalize(); + ibset const expanded_active (box.active.expanded_for (obox.interior)); + ibset const ovlp = oneedrecv & expanded_active; for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) @@ -660,7 +675,6 @@ regrid (bool const do_init) } oneedrecv -= ovlp; - oneedrecv.normalize(); // All points must have been received ASSERT_c (oneedrecv.empty(), @@ -673,7 +687,6 @@ regrid (bool const do_init) // Refinement prolongation: - // cerr << "QQQ: regrid[g]" << endl; static Carpet::Timer timer_comm_refprol ("CarpetLib::dh::regrid::comm::refprol"); @@ -696,19 +709,15 @@ regrid (bool const do_init) for (int cc = 0; cc < h.components(orl); ++ cc) { full_dboxes const & obox = full_boxes.AT(ml).AT(orl).AT(cc); - ibset contracted_oactive; - for (ibset::const_iterator - ai = obox.active.begin(); ai != obox.active.end(); ++ ai) - { - ibbox const & oactive = * ai; - // untested for cell centering - contracted_oactive += - oactive.contracted_for (box.interior).expand (reffact); - } - contracted_oactive.normalize(); - - ibset ovlp = needrecv & contracted_oactive; - ovlp.normalize(); +#if 0 + // untested for cell centering + ibset const expanded_oactive + (obox.active.contracted_for (box.interior).expand (reffact)); +#else + ibset const expanded_oactive + (obox.active.expanded_for (box.interior).expand (reffact)); +#endif + ibset const ovlp = needrecv & expanded_oactive; for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) @@ -733,7 +742,6 @@ regrid (bool const do_init) } // for cc // All points must have been received - needrecv.normalize(); ASSERT_c (needrecv.empty(), "Refinement prolongation: All points must have been received"); @@ -744,7 +752,6 @@ regrid (bool const do_init) // Synchronisation: - // cerr << "QQQ: regrid[h]" << endl; static Carpet::Timer timer_comm_sync ("CarpetLib::dh::regrid::comm::sync"); @@ -765,7 +772,6 @@ regrid (bool const do_init) // compatibility. ibset needrecv = box.ghosts; #endif - ibset const needrecv_orig = needrecv; ibset & sync = box.sync; @@ -773,11 +779,10 @@ regrid (bool const do_init) full_dboxes const & obox = full_level.AT(cc); #if 0 - ibset ovlp = needrecv & obox.owned; + ibset const ovlp = needrecv & obox.owned; #else - ibset ovlp = needrecv & obox.interior; + ibset const ovlp = needrecv & obox.interior; #endif - ovlp.normalize(); if (cc == c) { ASSERT_cc (ovlp.empty(), @@ -804,7 +809,6 @@ regrid (bool const do_init) } // for cc - sync.normalize(); } @@ -813,7 +817,6 @@ regrid (bool const do_init) // Boundary prolongation: - // cerr << "QQQ: regrid[i]" << endl; static Carpet::Timer timer_comm_refbndprol ("CarpetLib::dh::regrid::comm::refbndprol"); @@ -844,9 +847,6 @@ regrid (bool const do_init) // also all buffer zones needrecv += box.buffers; - needrecv.normalize(); - ibset const needrecv_orig = needrecv; - ibset & bndref = box.bndref; i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); @@ -860,19 +860,15 @@ regrid (bool const do_init) for (int cc = 0; cc < h.components(orl); ++ cc) { full_dboxes const & obox = full_boxes.AT(ml).AT(orl).AT(cc); - ibset contracted_oactive; - for (ibset::const_iterator - ai = obox.active.begin(); ai != obox.active.end(); ++ ai) - { - ibbox const & oactive = * ai; - // untested for cell centering - contracted_oactive += - oactive.contracted_for (box.interior).expand (reffact); - } - contracted_oactive.normalize(); - - ibset ovlp = needrecv & contracted_oactive; - ovlp.normalize(); +#if 0 + // untested for cell centering + ibset const expanded_oactive + (obox.active.contracted_for (box.interior).expand (reffact)); +#else + ibset const expanded_oactive + (obox.active.expanded_for (box.interior).expand (reffact)); +#endif + ibset const ovlp = needrecv & expanded_oactive; for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) @@ -897,11 +893,8 @@ regrid (bool const do_init) } // for cc - bndref.normalize(); - // All points must now have been received, either through // synchronisation or through boundary prolongation - needrecv.normalize(); ASSERT_c (needrecv.empty(), "Synchronisation and boundary prolongation: All points must have been received"); @@ -914,7 +907,6 @@ regrid (bool const do_init) // Refinement restriction: - // cerr << "QQQ: regrid[j]" << endl; static Carpet::Timer timer_comm_refrest ("CarpetLib::dh::regrid::comm::refrest"); @@ -934,29 +926,14 @@ regrid (bool const do_init) // Refinement restriction may fill all active points, and // must use all active points - ibset needrecv; - for (ibset::const_iterator - ai = box.active.begin(); ai != box.active.end(); ++ ai) - { - ibbox const & active = * ai; - needrecv += active.contracted_for (obox0.interior); - } - needrecv.normalize(); + ibset needrecv (box.active.contracted_for (obox0.interior)); for (int cc = 0; cc < h.components(orl); ++ cc) { full_dboxes & obox = full_boxes.AT(ml).AT(orl).AT(cc); - ibset contracted_active; - for (ibset::const_iterator - ai = box.active.begin(); ai != box.active.end(); ++ ai) - { - ibbox const & active = * ai; - contracted_active += active.contracted_for (obox0.interior); - } - contracted_active.normalize(); - - ibset ovlp = obox.active & contracted_active; - ovlp.normalize(); + ibset const contracted_active + (box.active.contracted_for (obox0.interior)); + ibset const ovlp = obox.active & contracted_active; for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) @@ -980,7 +957,6 @@ regrid (bool const do_init) } // for cc // All points must have been received - needrecv.normalize(); ASSERT_rl (needrecv.empty(), "Refinement restriction: All points must have been received"); @@ -995,8 +971,260 @@ regrid (bool const do_init) + // Mask: + static Carpet::Timer timer_mask ("CarpetLib::dh::regrid::mask"); + timer_mask.start(); + + if (rl > 0) { + int const orl = rl - 1; + full_cboxes const& full_olevel = full_boxes.AT(ml).AT(orl); + // Local boxes are not communicated or post-processed, and + // thus can be modified even for coarser levels + local_cboxes& local_olevel = local_boxes.AT(ml).AT(orl); + + ivect const reffact = h.reffacts.AT(rl) / h.reffacts.AT(orl); + + // This works only when the refinement factor is 2 + if (all (reffact == 2)) { + + ibbox const& base = domain_exterior; + ibbox const& obase = h.baseextent(ml,orl); + + // Calculate the union of all coarse regions + ibset const allointr (full_olevel, & full_dboxes::interior); + + // Project to current level + ivect const rf(reffact); + // TODO: Why expand by rf-1? This expansion shouldn't + // matter, since the coarse grid is larger than the fine + // grid anyway, except at the outer boundary + // ibset const parent (allointr.expand(rf-1,rf-1).contracted_for(base)); + ibset const parent (allointr.expanded_for(base)); + + // Subtract the active region + ibset const notrefined = parent - allactive; + + // Enlarge this set + vect<ibset,dim> enlarged; + for (int d=0; d<dim; ++d) { + switch (h.refcent) { + case vertex_centered: { + ivect const dir = ivect::dir(d); + enlarged[d] = ibset (notrefined.expand(dir, dir)); + break; + } + case cell_centered: { + enlarged[d] = notrefined; +#warning "TODO: restriction boundaries are wrong (they are empty, but should not be) with cell centring when fine cell cut coarse cells" + bool const old_there_was_an_error = there_was_an_error; + ASSERT_rl (notrefined.contracted_for(obase).expanded_for(base) == + notrefined, + "Refinement mask: Fine grid boundaries must be aligned with coarse grid cells"); + // Do not abort because of this problem + there_was_an_error = old_there_was_an_error; + break; + } + default: + assert (0); + } + } + + // Intersect with the active region + vect<ibset,dim> all_boundaries; + for (int d=0; d<dim; ++d) { + all_boundaries[d] = allactive & enlarged[d]; + } + +#warning "TODO: Ensure that the prolongation boundaries all_boundaries are contained in the boundary prolongated region" + +#warning "TODO: Ensure that the restriction boundaries and the restricted region are contained in the restricted region" + + // Subtract the boundaries from the refined region + ibset all_refined = allactive; + for (int d=0; d<dim; ++d) { + all_refined -= all_boundaries[d]; + } + + for (int lc = 0; lc < h.local_components(rl); ++ lc) { + int const c = h.get_component(rl, lc); + full_dboxes & box = full_level.AT(c); + full_dboxes const& obox = full_olevel.AT(c); + local_dboxes & local_box = local_level.AT(lc); + // Local boxes are not communicated or post-processed, and + // thus can be modified even for coarser levels + local_dboxes & local_obox = local_olevel.AT(lc); + + // Set restriction information for next coarser level + local_obox.restricted_region = + all_refined.contracted_for(obox.exterior) & obox.owned; + + // Set prolongation information for current level + for (int d=0; d<dim; ++d) { + local_obox.restriction_boundaries[d] = + all_boundaries[d].contracted_for(obox.exterior) & obox.owned; + } + + // Set prolongation information for current level + for (int d=0; d<dim; ++d) { + local_box.prolongation_boundaries[d] = + all_boundaries[d] & box.owned; + } + + } // for lc + + } // if reffact != 2 + + } // if not coarsest level + + timer_mask.stop(); + + + + // Refluxing: + static Carpet::Timer timer_reflux ("CarpetLib::dh::regrid::reflux"); + timer_reflux.start(); + + // If there is no coarser level, do nothing + if (rl > 0) { + int const orl = rl - 1; + full_cboxes const& full_olevel = full_boxes.AT(ml).AT(orl); + local_cboxes & local_olevel = local_boxes.AT(ml).AT(orl); + + // This works only with cell centred grids + if (h.refcent == cell_centered) { + + // Fine grids + ibset const& fine_level = allactive; + + ivect const izero = ivect (0); + ivect const ione = ivect (1); + + vect<vect<ibset,2>,dim> all_fine_boundary; + + for (int dir=0; dir<dim; ++dir) { + // Unit vector + ivect const idir = ivect::dir(dir); + + // Fine grids with boundary + ibset fine_level_plus; + for (ibset::const_iterator fine_level_i = fine_level.begin(); + fine_level_i != fine_level.end(); + ++ fine_level_i) + { + // Expand region to the left (but not to the right), + // since the fluxes are staggered by one half to the + // right + fine_level_plus |= fine_level_i->expand (idir, izero); + } + + // Fine grids without boundary + ibset fine_level_minus; + for (ibset::const_iterator fine_level_i = fine_level.begin(); + fine_level_i != fine_level.end(); + ++ fine_level_i) + { + // Shrink region at the left boundary (but not at the + // right boundary), since the fluxes are staggered by + // one half to the right + fine_level_minus |= fine_level_i->expand (izero, -idir); + } + + // Fine boundary only + all_fine_boundary[dir][0] = fine_level_plus - fine_level ; + all_fine_boundary[dir][1] = fine_level - fine_level_minus; + } // for dir + + + + // Coarse grids + ibbox const& coarse_domain_exterior = h.baseextent(ml, orl); + ibbox const& coarse_ext = coarse_domain_exterior; + + // Coarse equivalent of fine boundary + vect<vect<ibset,2>,dim> all_coarse_boundary; + for (int dir=0; dir<3; ++dir) { + // Unit vector + ivect const idir = ivect::dir(dir); + for (int face=0; face<2; ++face) { + for (ibset::const_iterator + all_fine_boundary_i = all_fine_boundary[dir][face].begin(); + all_fine_boundary_i != all_fine_boundary[dir][face].end(); + ++ all_fine_boundary_i) + { + ibbox const& fb = *all_fine_boundary_i; + + ivect const cstr = coarse_ext.stride(); + + ivect const flo = fb.lower(); + ivect const fup = fb.upper(); + ivect const fstr = fb.stride(); + + assert (all (h.reffacts.AT(rl) % h.reffacts.AT(orl) == 0)); + ivect const reffact = h.reffacts.AT(rl) / h.reffacts.AT(orl); + assert (all (reffact % 2 == 0)); + assert (all (fstr % 2 == 0)); + assert (all (cstr % 2 == 0)); + + ibbox const ftmp (flo + idir*(fstr/2-cstr/2), + fup + idir*(fstr/2-cstr/2), + fstr); + ibbox const ctmp = ftmp.contracted_for (coarse_ext); + + ivect const clo = ctmp.lower(); + ivect const cup = ctmp.upper(); + ibbox const cb (clo, cup, cstr); + + all_coarse_boundary[dir][face] += cb; + } + } // for face + } // for dir + + + + for (int lc = 0; lc < h.local_components(rl); ++ lc) { + int const c = h.get_component(rl, lc); + full_dboxes & box = full_level.AT(c); + local_dboxes & local_box = local_level.AT(lc); + + for (int dir = 0; dir < dim; ++ dir) { + for (int face = 0; face < 2; ++ face) { + + // This is not really used (only for debugging) + local_box.fine_boundary[dir][face] = + box.exterior & all_fine_boundary[dir][face]; + + } // for face + } // for dir + + } // for lc + + for (int lc = 0; lc < h.local_components(orl); ++ lc) { + int const c = h.get_component(orl, lc); + full_dboxes const & obox = full_olevel.AT(c); + local_dboxes & local_obox = local_olevel.AT(lc); + + for (int dir = 0; dir < dim; ++ dir) { + for (int face = 0; face < 2; ++ face) { + +#warning "TODO: Set up communication schedule for refluxing" + + local_obox.coarse_boundary[dir][face] = + obox.exterior & all_coarse_boundary[dir][face]; + + } // for face + } // for dir + + } // for lc + + } // if cell centered + + } // if rl > 0 + + timer_reflux.stop(); + + + // Regridding schedule: - // cerr << "QQQ: regrid[5]" << endl; fast_level.do_init = do_init; if (do_init) { @@ -1006,7 +1234,6 @@ regrid (bool const do_init) for (int lc = 0; lc < h.local_components(rl); ++ lc) { int const c = h.get_component (rl, lc); - // cerr << "QQQ: regrid[5a] lc=" << lc << " c=" << c << endl; full_dboxes & box = full_level.AT(c); @@ -1014,27 +1241,26 @@ regrid (bool const do_init) - // Synchronisation: - // cerr << "QQQ: regrid[k]" << endl; + // Regridding synchronisation: static Carpet::Timer timer_regrid_sync ("CarpetLib::dh::regrid::regrid::sync"); timer_regrid_sync.start(); - if (int (oldboxes.size()) > ml and int (oldboxes.AT(ml).size()) > rl) + if (int (old_light_boxes.size()) > ml and + int (old_light_boxes.AT(ml).size()) > rl) { - int const oldcomponents = oldboxes.AT(ml).AT(rl).size(); + int const oldcomponents = old_light_boxes.AT(ml).AT(rl).size(); // Synchronisation copies from the same level of the old // grid structure. It should fill as many active points // as possible. for (int cc = 0; cc < oldcomponents; ++ cc) { - dboxes const & obox = oldboxes.AT(ml).AT(rl).AT(cc); + light_dboxes const & obox = old_light_boxes.AT(ml).AT(rl).AT(cc); - ibset ovlp = needrecv & obox.owned; - ovlp.normalize(); + ibset const ovlp = needrecv & obox.owned; for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) @@ -1055,16 +1281,13 @@ regrid (bool const do_init) } // for cc - needrecv.normalize(); - - } // if not oldboxes.empty + } // if not old_light_boxes.empty timer_regrid_sync.stop(); - // Prolongation: - // cerr << "QQQ: regrid[l]" << endl; + // Regridding prolongation: static Carpet::Timer timer_regrid_prolongate ("CarpetLib::dh::regrid::regrid::prolongate"); @@ -1087,19 +1310,15 @@ regrid (bool const do_init) for (int cc = 0; cc < h.components(orl); ++ cc) { full_dboxes const & obox = full_boxes.AT(ml).AT(orl).AT(cc); - ibset contracted_oactive; - for (ibset::const_iterator - ai = obox.active.begin(); ai != obox.active.end(); ++ ai) - { - ibbox const & oactive = * ai; - // untested for cell centering - contracted_oactive += - oactive.contracted_for (box.interior).expand (reffact); - } - contracted_oactive.normalize(); - - ibset ovlp = needrecv & contracted_oactive; - ovlp.normalize(); +#if 0 + // untested for cell centering + ibset const expanded_oactive + (obox.active.contracted_for (box.interior).expand (reffact)); +#else + ibset const expanded_oactive + (obox.active.expanded_for (box.interior).expand (reffact)); +#endif + ibset const ovlp = needrecv & expanded_oactive; for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) @@ -1123,11 +1342,11 @@ regrid (bool const do_init) } // for cc - needrecv.normalize(); - } // if rl > 0 - if (int (oldboxes.size()) > ml and int (oldboxes.AT(ml).size()) > 0) { + if (int (old_light_boxes.size()) > ml and + int (old_light_boxes.AT(ml).size()) > 0) + { // All points must now have been received, either through // synchronisation or through prolongation ASSERT_c (needrecv.empty(), @@ -1144,26 +1363,27 @@ regrid (bool const do_init) - // cerr << "QQQ: regrid[6]" << endl; for (int lc = 0; lc < h.local_components(rl); ++ lc) { int const c = h.get_component (rl, lc); - level.AT(c).exterior = full_level.AT(c).exterior; - level.AT(c).owned = full_level.AT(c).owned; - level.AT(c).interior = full_level.AT(c).interior; + light_level.AT(c).exterior = full_level.AT(c).exterior; + light_level.AT(c).owned = full_level.AT(c).owned; + light_level.AT(c).interior = full_level.AT(c).interior; +#if 0 dh::dboxes::ibset2ibboxs (full_level.AT(c).active, - level.AT(c).active, level.AT(c).numactive); + light_level.AT(c).active, + light_level.AT(c).numactive); +#endif - level.AT(c).exterior_size = full_level.AT(c).exterior.size(); - level.AT(c).owned_size = full_level.AT(c).owned.size(); - level.AT(c).active_size = full_level.AT(c).active.size(); + light_level.AT(c).exterior_size = full_level.AT(c).exterior.size(); + light_level.AT(c).owned_size = full_level.AT(c).owned.size(); + light_level.AT(c).active_size = full_level.AT(c).active.size(); } // for lc // Broadcast grid structure and communication schedule - // cerr << "QQQ: regrid[7]" << endl; { @@ -1172,20 +1392,18 @@ regrid (bool const do_init) timer_bcast_boxes.start(); int const count_send = h.local_components(rl); - vector<dboxes> level_send (count_send); + vector<light_dboxes> level_send (count_send); for (int lc = 0; lc < h.local_components(rl); ++ lc) { int const c = h.get_component (rl, lc); - level_send.AT(lc) = level.AT(c); + level_send.AT(lc) = light_level.AT(c); } - // cerr << "QQQ: regrid[7a]" << endl; - vector<vector<dboxes> > const level_recv = + vector<vector<light_dboxes> > const level_recv = allgatherv (dist::comm(), level_send); - // cerr << "QQQ: regrid[7b]" << endl; vector<int> count_recv (dist::size(), 0); for (int c = 0; c < h.components(rl); ++ c) { int const p = this_proc (rl, c); if (p != dist::rank()) { - level.AT(c) = level_recv.AT(p).AT(count_recv.AT(p)); + light_level.AT(c) = level_recv.AT(p).AT(count_recv.AT(p)); ++ count_recv.AT(p); } } @@ -1194,7 +1412,6 @@ regrid (bool const do_init) assert (count_recv.AT(p) == int(level_recv.AT(p).size())); } } - // cerr << "QQQ: regrid[7c]" << endl; timer_bcast_boxes.stop(); @@ -1263,17 +1480,41 @@ regrid (bool const do_init) // Output: if (output_bboxes or there_was_an_error) { + cout << eol; + cout << "ml=" << ml << " rl=" << rl << eol; + cout << "baseextent=" << h.baseextent(ml,rl) << eol; + + for (int c = 0; c < h.components(rl); ++ c) { + cout << eol; + cout << "ml=" << ml << " rl=" << rl << " c=" << c << eol; + cout << "extent=" << h.extent(ml,rl,c) << eol; + cout << "outer_boundaries=" << h.outer_boundaries(rl,c) << eol; + cout << "processor=" << h.outer_boundaries(rl,c) << eol; + } // for c + for (int c = 0; c < h.components(rl); ++ c) { full_dboxes const & box = full_boxes.AT(ml).AT(rl).AT(c); - cout << eol; cout << "ml=" << ml << " rl=" << rl << " c=" << c << eol; cout << box; - } // for c - fast_dboxes const & fast_box = fast_boxes.AT(ml).AT(rl); + for (int c = 0; c < h.components(rl); ++ c) { + light_dboxes const & box = light_boxes.AT(ml).AT(rl).AT(c); + cout << eol; + cout << "ml=" << ml << " rl=" << rl << " c=" << c << eol; + cout << box; + } // for c + for (int lc = 0; lc < h.local_components(rl); ++ lc) { + int const c = h.get_component (rl, lc); + local_dboxes const & box = local_boxes.AT(ml).AT(rl).AT(lc); + cout << eol; + cout << "ml=" << ml << " rl=" << rl << " lc=" << lc << " c=" << c << eol; + cout << box; + } // for lc + + fast_dboxes const & fast_box = fast_boxes.AT(ml).AT(rl); cout << eol; cout << "ml=" << ml << " rl=" << rl << eol; cout << fast_box; @@ -1283,8 +1524,10 @@ regrid (bool const do_init) // Free memory early to save space - if (int (oldboxes.size()) > ml and int (oldboxes.AT(ml).size()) > rl) { - oldboxes.AT(ml).AT(rl).clear(); + if (int (old_light_boxes.size()) > ml and + int (old_light_boxes.AT(ml).size()) > rl) + { + old_light_boxes.AT(ml).AT(rl).clear(); } if (ml > 0) { @@ -1323,7 +1566,8 @@ regrid (bool const do_init) cout << eol; cout << "memoryof(gh)=" << memoryof(h) << eol; cout << "memoryof(dh)=" << memoryof(*this) << eol; - cout << "memoryof(dh.boxes)=" << memoryof(boxes) << eol; + cout << "memoryof(dh.light_boxes)=" << memoryof(light_boxes) << eol; + cout << "memoryof(dh.local_boxes)=" << memoryof(local_boxes) << eol; cout << "memoryof(dh.fast_boxes)=" << memoryof(fast_boxes) << eol; int gfcount = 0; size_t gfmemory = 0; @@ -1357,7 +1601,6 @@ broadcast_schedule (vector<fast_dboxes> & fast_level_otherprocs, fast_dboxes & fast_level, srpvect fast_dboxes::* const schedule_item) { - // cerr << "QQQ: broadcast_schedule[1]" << endl; static Carpet::Timer timer_bs1 ("CarpetLib::dh::bs1"); timer_bs1.start(); vector <srpvect> send (dist::size()); @@ -1376,7 +1619,6 @@ broadcast_schedule (vector<fast_dboxes> & fast_level_otherprocs, (fast_level.*schedule_item).insert ((fast_level.*schedule_item).end(), recv.begin(), recv.end()); timer_bs3.stop(); - // cerr << "QQQ: broadcast_schedule[2]" << endl; } @@ -1512,12 +1754,12 @@ operator== (full_dboxes const & b) const // MPI datatypes MPI_Datatype -mpi_datatype (dh::dboxes const &) +mpi_datatype (dh::light_dboxes const &) { static bool initialised = false; static MPI_Datatype newtype; if (not initialised) { - static dh::dboxes s; + static dh::light_dboxes s; #define ENTRY(type, name) \ { \ sizeof s.name / sizeof(type), /* count elements */ \ @@ -1530,17 +1772,19 @@ mpi_datatype (dh::dboxes const &) ENTRY(int, exterior), ENTRY(int, owned), ENTRY(int, interior), +#if 0 ENTRY(int, numactive), ENTRY(int, active), - ENTRY(dh::dboxes::size_type, exterior_size), - ENTRY(dh::dboxes::size_type, owned_size), - ENTRY(dh::dboxes::size_type, active_size), +#endif + ENTRY(dh::light_dboxes::size_type, exterior_size), + ENTRY(dh::light_dboxes::size_type, owned_size), + ENTRY(dh::light_dboxes::size_type, active_size), {1, sizeof s, MPI_UB, "MPI_UB", "MPI_UB"} }; #undef ENTRY newtype = dist::create_mpi_datatype (sizeof descr / sizeof descr[0], descr, - "dh::dboxes", sizeof s); + "dh::light::dboxes", sizeof s); #if 0 int type_size; MPI_Type_size (newtype, & type_size); @@ -1606,7 +1850,7 @@ memory () memoryof (ghost_widths) + memoryof (buffer_widths) + memoryof (prolongation_orders_space) + - memoryof (boxes) + + memoryof (light_boxes) + memoryof (fast_boxes) + memoryof (gfs); } @@ -1624,10 +1868,11 @@ allmemory () return mem; } -int const dh::dboxes::maxactive; +#if 0 +int const dh::light::dboxes::maxactive; void -dh::dboxes:: +dh::light_dboxes:: ibset2ibboxs (ibset const& s, ibbox* const bs, int& nbs) { assert (s.setsize() <= maxactive); @@ -1638,7 +1883,7 @@ ibset2ibboxs (ibset const& s, ibbox* const bs, int& nbs) } void -dh::dboxes:: +dh::light_dboxes:: ibboxs2ibset (ibbox const* const bs, int const& nbs, ibset& s) { s = ibset(); @@ -1647,9 +1892,10 @@ ibboxs2ibset (ibbox const* const bs, int const& nbs, ibset& s) s += bs[n]; } } +#endif size_t -dh::dboxes:: +dh::light_dboxes:: memory () const { @@ -1657,7 +1903,29 @@ memory () memoryof (exterior) + memoryof (owned) + memoryof (interior) + - memoryof (active); +#if 0 + memoryof (numactive) + + memoryof (active) + +#endif + memoryof (exterior_size) + + memoryof (owned_size) + + memoryof (active_size); +} + +size_t +dh::local_dboxes:: +memory () + const +{ + return + memoryof (buffers) + + memoryof (buffers_stepped) + + memoryof (active) + + memoryof (restricted_region) + + memoryof (restriction_boundaries) + + memoryof (prolongation_boundaries) + + memoryof (coarse_boundary) + + memoryof (fine_boundary); } size_t @@ -1701,13 +1969,13 @@ memory () // Input istream & -dh::dboxes:: +dh::light_dboxes:: input (istream & is) { // Regions: try { skipws (is); - consume (is, "dh::dboxes:{"); + consume (is, "dh::light_dboxes:{"); skipws (is); consume (is, "exterior:"); is >> exterior; @@ -1725,7 +1993,48 @@ input (istream & is) skipws (is); consume (is, "}"); } catch (input_error & err) { - cout << "Input error while reading a dh::full_dboxes" << endl; + cout << "Input error while reading a dh::light_dboxes" << endl; + throw err; + } + return is; +} + +istream & +dh::local_dboxes:: +input (istream & is) +{ + // Regions: + try { + skipws (is); + consume (is, "dh::local_dboxes:{"); + skipws (is); + consume (is, "buffers:"); + is >> buffers; + skipws (is); + consume (is, "buffers_stepped:"); + is >> buffers_stepped; + skipws (is); + consume (is, "active:"); + is >> active; + skipws (is); + consume (is, "restricted_region:"); + is >> restricted_region; + skipws (is); + consume (is, "restriction_boundaries:"); + is >> restriction_boundaries; + skipws (is); + consume (is, "prolongation_boundaries:"); + is >> prolongation_boundaries; + skipws (is); + consume (is, "coarse_boundary:"); + is >> coarse_boundary; + skipws (is); + consume (is, "fine_boundary:"); + is >> fine_boundary; + skipws (is); + consume (is, "}"); + } catch (input_error & err) { + cout << "Input error while reading a dh::local_dboxes" << endl; throw err; } return is; @@ -1838,7 +2147,7 @@ output (ostream & os) << "ghost_widths=" << ghost_widths << "," << "buffer_widths=" << buffer_widths << "," << "prolongation_orders_space=" << prolongation_orders_space << "," - << "boxes=" << boxes << "," + << "light_boxes=" << light_boxes << "," << "fast_boxes=" << fast_boxes << "," << "gfs={"; { @@ -1855,12 +2164,12 @@ output (ostream & os) } ostream & -dh::dboxes:: +dh::light_dboxes:: output (ostream & os) const { // Regions: - os << "dh::dboxes:{" << eol + os << "dh::light_dboxes:{" << eol << " exterior: " << exterior << eol << " owned: " << owned << eol << " interior: " << interior << eol @@ -1870,6 +2179,25 @@ output (ostream & os) } ostream & +dh::local_dboxes:: +output (ostream & os) + const +{ + // Regions: + os << "dh::local_dboxes:{" << eol + << " buffers: " << buffers << eol + << " buffers_stepped: " << buffers_stepped << eol + << " active: " << active << eol + << " restricted_region: " << restricted_region << eol + << " restriction_boundaries: " << restriction_boundaries << eol + << " prolongation_boundaries: " << prolongation_boundaries << eol + << " coarse_boundary: " << coarse_boundary << eol + << " fine_boundary: " << fine_boundary << eol + << "}" << eol; + return os; +} + +ostream & dh::full_dboxes:: output (ostream & os) const diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index d509f1bda..83e44e0f0 100644 --- a/Carpet/CarpetLib/src/dh.hh +++ b/Carpet/CarpetLib/src/dh.hh @@ -41,7 +41,7 @@ public: - struct dboxes { + struct light_dboxes { // Region description: @@ -49,18 +49,44 @@ public: ibbox owned; // evolved in time ibbox interior; // interior (without ghost zones) +#if 0 // TODO: Create a new datatype bboxarr for this? Or get rid of // it? int numactive; static int const maxactive = 4; ibbox active[maxactive]; // owned minus buffers +#endif // Region statistics: typedef ibbox::size_type size_type; size_type exterior_size, owned_size, active_size; +#if 0 static void ibset2ibboxs (ibset const& s, ibbox* bs, int& nbs); static void ibboxs2ibset (ibbox const* bs, int const& nbs, ibset& s); +#endif + + size_t memory () const CCTK_ATTRIBUTE_PURE; + istream & input (istream & is); + ostream & output (ostream & os) const; + }; + + struct local_dboxes { + + // Information about the processor-local region: + + ibset buffers; // buffer zones + vector<ibset> buffers_stepped; // buffer zones [substep] + ibset active; // owned minus buffers + + // Mask + ibset restricted_region; // filled by restriction + vect<ibset,dim> restriction_boundaries; // partly filled by restriction + vect<ibset,dim> prolongation_boundaries; // partly used by prolongation + + // Refluxing + vect<vect<ibset,2>,dim> coarse_boundary; + vect<vect<ibset,2>,dim> fine_boundary; size_t memory () const CCTK_ATTRIBUTE_PURE; istream & input (istream & is); @@ -131,9 +157,13 @@ public: ostream & output (ostream & os) const; }; - typedef vector<dboxes> cboxes; // ... for each component - typedef vector<cboxes> rboxes; // ... for each refinement level - typedef vector<rboxes> mboxes; // ... for each multigrid level + typedef vector<light_dboxes> light_cboxes; // ... for each component + typedef vector<light_cboxes> light_rboxes; // ... for each refinement level + typedef vector<light_rboxes> light_mboxes; // ... for each multigrid level + + typedef vector<local_dboxes> local_cboxes; // ... for each component + typedef vector<local_cboxes> local_rboxes; // ... for each refinement level + typedef vector<local_rboxes> local_mboxes; // ... for each multigrid level typedef vector<full_dboxes> full_cboxes; // ... for each component typedef vector<full_cboxes> full_rboxes; // ... for each refinement level @@ -163,8 +193,9 @@ public: // should be readonly vector<int> prolongation_orders_space; // order of spatial // prolongation operator [rl] - mboxes boxes; // grid hierarchy - fast_mboxes fast_boxes; // grid hierarchy + light_mboxes light_boxes; // grid hierarchy [ml][rl][c] + local_mboxes local_boxes; // grid hierarchy [ml][rl][lc] + fast_mboxes fast_boxes; // grid hierarchy [ml][rl][p] typedef list<ggf*>::iterator ggf_handle; list<ggf*> gfs; // list of all grid functions @@ -213,21 +244,27 @@ public: -MPI_Datatype mpi_datatype (dh::dboxes const &) CCTK_ATTRIBUTE_CONST; +MPI_Datatype mpi_datatype (dh::light_dboxes const &) CCTK_ATTRIBUTE_CONST; MPI_Datatype mpi_datatype (dh::fast_dboxes const &); namespace dist { - template<> inline MPI_Datatype mpi_datatype<dh::dboxes> () + template<> inline MPI_Datatype mpi_datatype<dh::light_dboxes> () CCTK_ATTRIBUTE_CONST; - template<> inline MPI_Datatype mpi_datatype<dh::dboxes> () - { dh::dboxes dummy; return mpi_datatype(dummy); } + template<> inline MPI_Datatype mpi_datatype<dh::light_dboxes> () + { dh::light_dboxes dummy; return mpi_datatype(dummy); } template<> inline MPI_Datatype mpi_datatype<dh::fast_dboxes> () CCTK_ATTRIBUTE_CONST; template<> inline MPI_Datatype mpi_datatype<dh::fast_dboxes> () { dh::fast_dboxes dummy; return mpi_datatype(dummy); } } -inline size_t memoryof (dh::dboxes const & b) CCTK_ATTRIBUTE_PURE; -inline size_t memoryof (dh::dboxes const & b) +inline size_t memoryof (dh::light_dboxes const & b) CCTK_ATTRIBUTE_PURE; +inline size_t memoryof (dh::light_dboxes const & b) +{ + return b.memory (); +} + +inline size_t memoryof (dh::local_dboxes const & b) CCTK_ATTRIBUTE_PURE; +inline size_t memoryof (dh::local_dboxes const & b) { return b.memory (); } @@ -250,7 +287,12 @@ inline size_t memoryof (dh const & d) return d.memory (); } -inline istream & operator>> (istream & is, dh::dboxes & b) +inline istream & operator>> (istream & is, dh::light_dboxes & b) +{ + return b.input (is); +} + +inline istream & operator>> (istream & is, dh::local_dboxes & b) { return b.input (is); } @@ -265,7 +307,12 @@ inline istream & operator>> (istream & is, dh::fast_dboxes & b) return b.input (is); } -inline ostream & operator<< (ostream & os, dh::dboxes const & b) +inline ostream & operator<< (ostream & os, dh::light_dboxes const & b) +{ + return b.output (os); +} + +inline ostream & operator<< (ostream & os, dh::local_dboxes const & b) { return b.output (os); } diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index 8a7c643da..02b90bc38 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -93,7 +93,7 @@ void ggf::set_timelevels (const int ml, const int rl, const int new_timelevels) for (int tl=timelevels(ml,rl); tl<new_timelevels; ++tl) { storage.AT(ml).AT(rl).AT(lc).AT(tl) = typed_data(tl,rl,lc,ml); storage.AT(ml).AT(rl).AT(lc).AT(tl)->allocate - (d.boxes.AT(ml).AT(rl).AT(c).exterior, dist::rank()); + (d.light_boxes.AT(ml).AT(rl).AT(c).exterior, dist::rank()); } // for tl } // for lc @@ -157,7 +157,7 @@ void ggf::recompose_allocate (const int rl) for (int tl=0; tl<timelevels(ml,rl); ++tl) { storage.AT(ml).AT(rl).AT(lc).AT(tl) = typed_data(tl,rl,lc,ml); storage.AT(ml).AT(rl).AT(lc).AT(tl)->allocate - (d.boxes.AT(ml).AT(rl).AT(c).exterior, dist::rank()); + (d.light_boxes.AT(ml).AT(rl).AT(c).exterior, dist::rank()); } // for tl } // for lc } // for ml diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc index 5ecc90a5c..22ac4ae06 100644 --- a/Carpet/CarpetLib/src/gh.cc +++ b/Carpet/CarpetLib/src/gh.cc @@ -161,7 +161,6 @@ regrid (rregs const & superregs, mregs const & regs, bool const do_init) for (int c=0; c<components(rl-1); ++c) { coarse += extent(ml,rl-1,c); } - coarse.normalize(); // then check all fine grids for (int c=0; c<components(rl); ++c) { ibbox const & fine = diff --git a/Carpet/CarpetLib/src/mpi_string.cc b/Carpet/CarpetLib/src/mpi_string.cc index bf794881a..e0022e3e0 100644 --- a/Carpet/CarpetLib/src/mpi_string.cc +++ b/Carpet/CarpetLib/src/mpi_string.cc @@ -471,9 +471,9 @@ namespace CarpetLib template - vector <vector <dh::dboxes> > + vector <vector <dh::light_dboxes> > allgatherv (MPI_Comm comm, - vector <dh::dboxes> const & data); + vector <dh::light_dboxes> const & data); template vector <sendrecv_pseudoregion_t> diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc index 860d24c6b..00dabc410 100644 --- a/Carpet/CarpetLib/src/region.cc +++ b/Carpet/CarpetLib/src/region.cc @@ -130,12 +130,6 @@ combine_regions (vector<region_t> const & oldregs, } } } - comps.normalize(); - for (int f = 0; f < 2; ++ f) { - for (int d = 0; d < dim; ++ d) { - cobnds[f][d].normalize(); - } - } // Reserve (generous) memory for the result size_t const needsize = newregs.size() + comps.setsize(); diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc index 02b857557..b02a46979 100644 --- a/Carpet/CarpetLib/src/vect.cc +++ b/Carpet/CarpetLib/src/vect.cc @@ -4,6 +4,7 @@ #include "cctk.h" #include "defs.hh" +#include "bboxset.hh" #include "vect.hh" @@ -12,7 +13,7 @@ using namespace std; // Input -template<class T,int D> +template<typename T,int D> void vect<T,D>::input (istream& is) { skipws (is); consume (is, '['); @@ -31,7 +32,7 @@ void vect<T,D>::input (istream& is) { // Output -template<class T,int D> +template<typename T,int D> void vect<T,D>::output (ostream& os) const { os << "["; for (int d=0; d<D; ++d) { @@ -83,3 +84,37 @@ template void vect<vect<bool,2>,dim>::output (ostream& os) const; template void vect<vect<int,2>,dim>::output (ostream& os) const; template void vect<vect<bool,dim>,2>::output (ostream& os) const; template void vect<vect<int,dim>,2>::output (ostream& os) const; +template void vect<vect<CCTK_REAL,dim>,2>::output (ostream& os) const; + + + +// Instantiate for bboxset class + +#define DEFINE_FAKE_VECT_OPERATIONS(T,D) \ +template<> vect<T,D> vect<T,D>::dir (const int d) { assert(0); } \ +template<> vect<T,D> vect<T,D>::seq () { assert(0); } \ +template<> vect<T,D> vect<T,D>::seq (const int n) { assert(0); } \ +template<> vect<T,D> vect<T,D>::seq (const int n, const int s) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator*= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator*= (const T&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator/= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator/= (const T&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator%= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator%= (const T&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator^= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator^= (const T&) { assert(0); } \ +template<> vect<T,D> vect<T,D>::operator+ () const { assert(0); } \ +template<> vect<T,D> vect<T,D>::operator- () const { assert(0); } \ +template<> vect<T,D> vect<T,D>::operator~ () const { assert(0); } \ +template class vect<T,D>; \ +template size_t memoryof (const vect<T,D>&); \ +template istream& operator>> (istream& is, vect<T,D>&); \ +template ostream& operator<< (ostream& os, const vect<T,D>&); + +typedef bboxset<int,dim> T1; +typedef vect<bboxset<int,dim>,2> T2; + +DEFINE_FAKE_VECT_OPERATIONS(T1,dim) +DEFINE_FAKE_VECT_OPERATIONS(T2,dim) + +#undef DEFINE_FAKE_VECT_OPERATIONS diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh index d7d648b7a..27b4ffad0 100644 --- a/Carpet/CarpetLib/src/vect.hh +++ b/Carpet/CarpetLib/src/vect.hh @@ -394,11 +394,11 @@ DECLARE_FUNCTION_1 (floor) DECLARE_FUNCTION_1 (sqrt) namespace CarpetLib { namespace good { - DECLARE_FUNCTION_1_RET (fpclassify, int) - DECLARE_FUNCTION_1_RET (isfinite, int) - DECLARE_FUNCTION_1_RET (isinf, int) - DECLARE_FUNCTION_1_RET (isnan, int) - DECLARE_FUNCTION_1_RET (isnormal, int) + // DECLARE_FUNCTION_1_RET (fpclassify, int) + DECLARE_FUNCTION_1_RET (isfinite_, int) + DECLARE_FUNCTION_1_RET (isinf_, int) + DECLARE_FUNCTION_1_RET (isnan_, int) + DECLARE_FUNCTION_1_RET (isnormal_, int) } } @@ -474,6 +474,18 @@ inline int maxloc (const vect<T,D>& a) return r; } +/** Return the index of the last maximum element. */ +template<typename T,int D> +inline int maxloc1 (const vect<T,D>& a) CCTK_ATTRIBUTE_PURE; +template<typename T,int D> +inline int maxloc1 (const vect<T,D>& a) +{ + ASSERT_VECT (D>0); + int r(D-1); + for (int d=D-2; d>=0; --d) if (a[d]>a[r]) r=d; + return r; +} + /** Return the index of the first minimum element. */ template<typename T,int D> inline int minloc (const vect<T,D>& a) CCTK_ATTRIBUTE_PURE; @@ -486,6 +498,18 @@ inline int minloc (const vect<T,D>& a) return r; } +/** Return the index of the last minimum element. */ +template<typename T,int D> +inline int minloc1 (const vect<T,D>& a) CCTK_ATTRIBUTE_PURE; +template<typename T,int D> +inline int minloc1 (const vect<T,D>& a) +{ + ASSERT_VECT (D>0); + int r(D-1); + for (int d=D-2; d>=0; --d) if (a[d]<a[r]) r=d; + return r; +} + /** Return the n-dimensional linear array index. */ template<typename T,int D> inline T index (const vect<T,D>& lsh, const vect<T,D>& ind) CCTK_ATTRIBUTE_PURE; diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc index 260afb23d..e469d105b 100644 --- a/Carpet/CarpetReduce/src/mask_carpet.cc +++ b/Carpet/CarpetReduce/src/mask_carpet.cc @@ -22,28 +22,6 @@ namespace CarpetMask { - void ibbox2iminimax (ibbox const& ext, // component extent - ibbox const& box, // this bbox - ivect const& lsh, // component's lsh - ivect& imin, ivect& imax) - { - ivect const izero = ivect(0); - - assert (all ((box.lower() - ext.lower() ) >= 0)); - assert (all ((box.upper() - ext.lower() + ext.stride()) >= 0)); - assert (all ((box.lower() - ext.lower() ) % ext.stride() == 0)); - assert (all ((box.upper() - ext.lower() + ext.stride()) % ext.stride() == 0)); - - imin = (box.lower() - ext.lower() ) / ext.stride(); - imax = (box.upper() - ext.lower() + ext.stride()) / ext.stride(); - - assert (all (izero <= imin)); - assert (box.empty() xor all (imin <= imax)); - assert (all (imax <= lsh)); - } - - - /** * Reduce the weight on the current and the next coarser level to * make things consistent. Set the weight to 0 inside the @@ -76,86 +54,27 @@ namespace CarpetMask { - // Calculate the union of all refined regions - ibset active; - for (int c=0; c<hh.components(reflevel); ++c) { - // refined |= hh.extent(mglevel,reflevel,c); - ibset this_active; - dh::dboxes const& box = dd.boxes.AT(mglevel).AT(reflevel).AT(c); - dh::dboxes::ibboxs2ibset (box.active, box.numactive, this_active); - active += this_active; - } - active.normalize(); - - // Calculate the union of all coarse regions - ibset parent; - for (int c=0; c<hh.components(reflevel-1); ++c) { - parent |= hh.extent(mglevel,reflevel-1,c) - .expand(ivect(reffact-1),ivect(reffact-1)).contracted_for(base); - } - parent.normalize(); - - // Subtract the refined region - ibset notrefined = parent - active; - notrefined.normalize(); - - // Enlarge this set - ibset enlarged[dim]; - for (int d=0; d<dim; ++d) { - for (ibset::const_iterator - bi = notrefined.begin(); bi != notrefined.end(); ++bi) - { - if (hh.refcent == vertex_centered) { - enlarged[d] |= (*bi).expand(ivect::dir(d), ivect::dir(d)); - } else { - enlarged[d] |= *bi; - } - } - enlarged[d].normalize(); - } - - // Intersect with the union of refined regions - ibset boundaries[dim]; - for (int d=0; d<dim; ++d) { - boundaries[d] = active & enlarged[d]; - boundaries[d].normalize(); - } - - // Subtract the boundaries from the refined region - ibset refined = active; - for (int d=0; d<dim; ++d) { - refined -= boundaries[d]; - } - refined.normalize(); - - - // Set prolongation boundaries of this level { BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { DECLARE_CCTK_ARGUMENTS; - ibbox const & ext - = dd.boxes.at(mglevel).at(reflevel).at(component).exterior; + ibbox const & ext = + dd.light_boxes.AT(mglevel).AT(reflevel).AT(component).exterior; + ibset const & active = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).active; + + vect<ibset,dim> const & boundaries = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).prolongation_boundaries; - ibset notactive = ext - active; - notactive.normalize(); + ibset const notactive = ext - active; for (int d=0; d<dim; ++d) { assert ((notactive & boundaries[d]).empty()); } - for (ibset::const_iterator bi = notactive.begin(); - bi != notactive.end(); - ++bi) - { - ibbox const& box = *bi; - assert (box <= ext); - assert (not box.empty()); - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); + LOOP_OVER_BSET (cctkGH, notactive, box, imin, imax) { if (verbose) { ostringstream buf; @@ -175,41 +94,31 @@ namespace CarpetMask { weight[ind] = 0.0; } LC_ENDLOOP3(CarpetMaskSetup_buffers); - } // for box + } END_LOOP_OVER_BSET; for (int d=0; d<dim; ++d) { - for (ibset::const_iterator bi = boundaries[d].begin(); - bi != boundaries[d].end(); - ++bi) - { + LOOP_OVER_BSET (cctkGH, boundaries[d], box, imin, imax) { - ibbox const & box = (*bi) & ext; - if (not box.empty()) { - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); - - if (verbose) { - ostringstream buf; - buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; - CCTK_INFO (buf.str().c_str()); - } - - // Set weight on the boundary to 1/2 - assert (dim == 3); + if (verbose) { + ostringstream buf; + buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; + CCTK_INFO (buf.str().c_str()); + } + + // Set weight on the boundary to 1/2 + assert (dim == 3); #pragma omp parallel - LC_LOOP3(CarpetMaskSetup_prolongation, - i,j,k, - imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) - { - int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - weight[ind] *= 0.5; - } LC_ENDLOOP3(CarpetMaskSetup_prolongation); - - } // if box not empty + LC_LOOP3(CarpetMaskSetup_prolongation, + i,j,k, + imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + weight[ind] *= 0.5; + } LC_ENDLOOP3(CarpetMaskSetup_prolongation); + - } // for box + } END_LOOP_OVER_BSET; } // for d } END_LOCAL_COMPONENT_LOOP; @@ -231,41 +140,32 @@ namespace CarpetMask { DECLARE_CCTK_ARGUMENTS; - ibbox const & ext - = dd.boxes.at(mglevel).at(reflevel).at(component).exterior; + ibset const & refined = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).restricted_region; + vect<ibset,dim> const & boundaries = + dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).restriction_boundaries; - for (ibset::const_iterator bi = refined.begin(); - bi != refined.end(); - ++bi) - { + LOOP_OVER_BSET (cctkGH, refined, box, imin, imax) { - ibbox const & box = (*bi).contracted_for(ext) & ext; - if (not box.empty()) { - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); - - if (verbose) { - ostringstream buf; - buf << "Setting restricted region on level " << reflevel << " to weight 0: " << imin << ":" << imax-ione; - CCTK_INFO (buf.str().c_str()); - } - - // Set weight in the restricted region to 0 - assert (dim == 3); + if (verbose) { + ostringstream buf; + buf << "Setting restricted region on level " << reflevel << " to weight 0: " << imin << ":" << imax-ione; + CCTK_INFO (buf.str().c_str()); + } + + // Set weight in the restricted region to 0 + assert (dim == 3); #pragma omp parallel - LC_LOOP3(CarpetMaskSetup_restriction, - i,j,k, - imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) - { - int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - weight[ind] = 0; - } LC_ENDLOOP3(CarpetMaskSetup_restriction); - - } // if box not empty - - } // for box + LC_LOOP3(CarpetMaskSetup_restriction, + i,j,k, + imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + weight[ind] = 0; + } LC_ENDLOOP3(CarpetMaskSetup_restriction); + + } END_LOOP_OVER_BSET; assert (dim == 3); vector<int> mask (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]); @@ -282,41 +182,30 @@ namespace CarpetMask { } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_init); for (int d=0; d<dim; ++d) { - for (ibset::const_iterator bi = boundaries[d].begin(); - bi != boundaries[d].end(); - ++bi) - { + LOOP_OVER_BSET (cctkGH, boundaries[d], box, imin, imax) { - ibbox const & box = (*bi).contracted_for(ext) & ext; - if (not box.empty()) { - - ivect imin, imax; - ibbox2iminimax (ext, box, ivect::ref(cctk_lsh), imin, imax); - - if (verbose) { - ostringstream buf; - buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; - CCTK_INFO (buf.str().c_str()); - } - - // Set weight on the boundary to 1/2 - assert (dim == 3); + if (verbose) { + ostringstream buf; + buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione; + CCTK_INFO (buf.str().c_str()); + } + + // Set weight on the boundary to 1/2 + assert (dim == 3); #pragma omp parallel - LC_LOOP3(CarpetMaskSetup_restriction_boundary_partial, - i,j,k, - imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) - { - int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - if (mask[ind] == 0) { - mask[ind] = 1; - } - mask[ind] *= 2; - } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_partial); - - } // if box not empty + LC_LOOP3(CarpetMaskSetup_restriction_boundary_partial, + i,j,k, + imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + if (mask[ind] == 0) { + mask[ind] = 1; + } + mask[ind] *= 2; + } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_partial); - } // for box + } END_LOOP_OVER_BSET; } // for d assert (dim == 3); diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc index 577fe2478..7e28a9ddc 100644 --- a/Carpet/CarpetReduce/src/reduce.cc +++ b/Carpet/CarpetReduce/src/reduce.cc @@ -250,10 +250,10 @@ namespace CarpetReduce { #ifdef HAVE_CCTK_REAL4 - template<> inline int myisnan (complex<CCTK_REAL4> const x) - { - return isnan (x.real()) or isnan (x.imag()); - } + // template<> inline int myisnan (complex<CCTK_REAL4> const x) + // { + // return isnan (x.real()) or isnan (x.imag()); + // } template<> inline complex<CCTK_REAL4> mymin (const complex<CCTK_REAL4> x, const complex<CCTK_REAL4> y) @@ -299,10 +299,10 @@ namespace CarpetReduce { #ifdef HAVE_CCTK_REAL8 - template<> inline int myisnan (complex<CCTK_REAL8> const x) - { - return isnan (x.real()) or isnan (x.imag()); - } + // template<> inline int myisnan (complex<CCTK_REAL8> const x) + // { + // return isnan (x.real()) or isnan (x.imag()); + // } template<> inline complex<CCTK_REAL8> mymin (const complex<CCTK_REAL8> x, const complex<CCTK_REAL8> y) @@ -348,10 +348,10 @@ namespace CarpetReduce { #ifdef HAVE_CCTK_REAL16 - template<> inline int myisnan (complex<CCTK_REAL16> const x) - { - return isnan (x.real()) or isnan (x.imag()); - } + // template<> inline int myisnan (complex<CCTK_REAL16> const x) + // { + // return isnan (x.real()) or isnan (x.imag()); + // } template<> inline complex<CCTK_REAL16> mymin (const complex<CCTK_REAL16> x, const complex<CCTK_REAL16> y) diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl index c87da7a8f..46bd9f6c2 100644 --- a/Carpet/CarpetRegrid2/param.ccl +++ b/Carpet/CarpetRegrid2/param.ccl @@ -4,6 +4,10 @@ BOOLEAN verbose "Display regridding information on the terminal" STEERABLE=alway { } "no" +BOOLEAN veryverbose "Display much regridding information on the terminal" STEERABLE=always +{ +} "no" + CCTK_INT min_distance "Minimum distance (in grid points) between coarse and fine grid boundaries" STEERABLE=always diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc index 76e4563ab..99c1832ff 100644 --- a/Carpet/CarpetRegrid2/src/regrid.cc +++ b/Carpet/CarpetRegrid2/src/regrid.cc @@ -126,22 +126,17 @@ namespace CarpetRegrid2 { rvect const & origin, rvect const & scale, gh const & hh, int const rl) { - DECLARE_CCTK_PARAMETERS; - - ivect const istride = hh.baseextents.at(0).at(rl).stride(); - ivect const cistride = - snap_to_coarse and rl > 0 - ? hh.baseextents.at(0).at(rl-1).stride() - : istride; + ivect const istride = hh.baseextents.at(0).at(rl).stride(); + ivect const bistride = hh.baseextents.at(0).at(0).stride(); if (hh.refcent == cell_centered) { assert (all (istride % 2 == 0)); } - + ivect const ipos = hh.refcent == vertex_centered - ? ivect (floor ((rpos - origin) * scale / rvect(cistride) + rvect(0.5) )) * cistride - : ivect (floor ((rpos - origin) * scale / rvect(cistride) + rvect(0.5) - rvect(istride / 2) / rvect(cistride))) * cistride + istride / 2; + ? ivect (floor (((rpos - origin) * scale ) / rvect(istride) + rvect(0.5))) * istride + : ivect (floor (((rpos - origin) * scale - rvect(bistride/2)) / rvect(istride) )) * istride + istride/2 + bistride/2; return ipos; } @@ -154,13 +149,8 @@ namespace CarpetRegrid2 { rvect const & origin, rvect const & scale, gh const & hh, int const rl) { - DECLARE_CCTK_PARAMETERS; - - ivect const istride = hh.baseextents.at(0).at(rl).stride(); - ivect const cistride = - snap_to_coarse and rl > 0 - ? hh.baseextents.at(0).at(rl-1).stride() - : istride; + ivect const istride = hh.baseextents.at(0).at(rl).stride(); + ivect const bistride = hh.baseextents.at(0).at(0).stride(); if (hh.refcent == cell_centered) { assert (all (istride % 2 == 0)); @@ -168,8 +158,8 @@ namespace CarpetRegrid2 { ivect const ipos = hh.refcent == vertex_centered - ? ivect (ceil ((rpos - origin) * scale / rvect(cistride) - rvect(0.5) )) * cistride - : ivect (ceil ((rpos - origin) * scale / rvect(cistride) - rvect(0.5) + rvect(istride / 2) / rvect(cistride))) * cistride - istride / 2; + ? ivect (ceil (((rpos - origin) * scale ) / rvect(istride) - rvect(0.5))) * istride + : ivect (ceil (((rpos - origin) * scale - rvect(bistride/2)) / rvect(istride) )) * istride - istride/2 + bistride/2; return ipos; } @@ -195,6 +185,31 @@ namespace CarpetRegrid2 { ipos2rpos (ib.stride(), zero , scale, hh, rl)); } + // Snap (enlarge) a bbox to the next coarser level, if desired + ibbox + snap_ibbox (ibbox const & ib, + gh const & hh, int const rl) + { + DECLARE_CCTK_PARAMETERS; + + assert (not ib.empty()); + assert (rl > 0); + + if (not snap_to_coarse) return ib; + + ibbox const & base = hh.baseextents.at(0).at(rl); + ibbox const & cbase = hh.baseextents.at(0).at(rl-1); + assert (all (cbase.stride() % base.stride() == 0)); + ivect const reffact = cbase.stride() / base.stride(); + + ivect const lo = ib.lower(); + ivect const up = ib.upper(); + ivect const str = ib.stride(); + assert (all (str == base.stride())); + + return ib.expand(reffact-1).contracted_for(cbase).expanded_for(base); + } + void @@ -309,7 +324,7 @@ namespace CarpetRegrid2 { { DECLARE_CCTK_PARAMETERS; - if (verbose) CCTK_INFO ("Regridding"); + if (verbose or veryverbose) CCTK_INFO ("Regridding"); assert (is_singlemap_mode()); gh const & hh = * vhh.at (Carpet::map); @@ -381,7 +396,7 @@ namespace CarpetRegrid2 { if (is_active and in_sync) break; ++ min_rl; } - if (verbose) { + if (verbose or veryverbose) { CCTK_VInfo (CCTK_THORNSTRING, "Regridding levels %d and up", min_rl); } } @@ -415,12 +430,17 @@ namespace CarpetRegrid2 { // Set up coarse levels for (int rl = 0; rl < min_rl; ++ rl) { + if (veryverbose) { + cout << "Refinement level " << rl << ": will not be changed" << endl; + } for (size_t c = 0; c < regss.at(rl).size(); ++ c) { regions.at(rl) += regss.at(rl).at(c).extent; } - regions.at(rl).normalize(); + if (veryverbose) { + cout << "Refinement level " << rl << ": regions are " << regions.at(rl) << endl; + } } - + // Refine only patch 0 if (Carpet::map == 0) { // Loop over all centres @@ -431,6 +451,10 @@ namespace CarpetRegrid2 { // Loop over all levels for this centre for (int rl = min_rl; rl < centre.num_levels; ++ rl) { + if (veryverbose) { + cout << "Refinement level " << rl << ": importing refined region settings..." << endl; + } + // Calculate a bbox for this region rvect const rmin = centre.position - centre.radius.at(rl); rvect const rmax = centre.position + centre.radius.at(rl); @@ -444,14 +468,18 @@ namespace CarpetRegrid2 { rpos2ipos1 (rmax, origin, scale, hh, rl) + boundary_shiftout * istride; - ibbox const region (imin, imax, istride); + ibbox const region = + snap_ibbox (ibbox (imin, imax, istride), hh, rl); // Add this region to the list of regions if (static_cast <int> (regions.size()) < rl+1) { regions.resize (rl+1); } regions.at(rl) |= region; - regions.at(rl).normalize(); + + if (veryverbose) { + cout << "Refinement level " << rl << ": preliminary regions are " << regions.at(rl) << endl; + } } // for rl @@ -477,9 +505,21 @@ namespace CarpetRegrid2 { if (ensure_proper_nesting) { if (rl < int(regions.size()) - 1) { + if (veryverbose) { + cout << "Refinement level " << rl << ": ensuring proper nesting..." << endl; + } + assert (not regions.at(rl).empty()); ibbox const coarse0 = * regions.at(rl).begin(); + // This is the location of the outermost grid points. For + // cell centring, these are 1/2 grid spacing inside of the + // boundary. + ivect const level_physical_ilower = + rpos2ipos (physical_lower, origin, scale, hh, rl); + ivect const level_physical_iupper = + rpos2ipos1 (physical_upper, origin, scale, hh, rl); + i2vect const fdistance = dd.ghost_widths.at(rl); i2vect const cdistance = i2vect (min_distance + dd.prolongation_stencil_size(rl)); @@ -490,8 +530,8 @@ namespace CarpetRegrid2 { { ibbox const & fbb = * ibb; - bvect const lower_is_outer = fbb.lower() <= physical_ilower; - bvect const upper_is_outer = fbb.upper() >= physical_iupper; + bvect const lower_is_outer = fbb.lower() <= level_physical_ilower; + bvect const upper_is_outer = fbb.upper() >= level_physical_iupper; b2vect const ob (lower_is_outer, upper_is_outer); ibbox const ebb = fbb.expand (i2vect (not ob) * fdistance); @@ -499,10 +539,12 @@ namespace CarpetRegrid2 { ibbox const ecbb = cbb.expand (i2vect (not ob) * cdistance); // Enlarge this level - regions.at(rl) |= ecbb; + regions.at(rl) |= snap_ibbox (ecbb, hh, rl); } - regions.at(rl).normalize(); + if (veryverbose) { + cout << "Refinement level " << rl << ": enlarged regions to " << regions.at(rl) << endl; + } } } // if ensure proper nesting @@ -513,6 +555,10 @@ namespace CarpetRegrid2 { // Add buffer zones // { + if (veryverbose) { + cout << "Refinement level " << rl << ": adding buffer zones..." << endl; + } + ibboxset buffered; for (ibboxset::const_iterator ibb = regions.at(rl).begin(); ibb != regions.at(rl).end(); @@ -520,19 +566,15 @@ namespace CarpetRegrid2 { { ibbox const & bb = * ibb; - // bvect const lower_is_outer = bb.lower() <= physical_ilower; - // bvect const upper_is_outer = bb.upper() >= physical_iupper; - // - // b2vect const ob (lower_is_outer, upper_is_outer); - // - // ibbox const bbb = bb.expand (i2vect (not ob) * dd.buffer_width); - ibbox const bbb = bb.expand (dd.buffer_widths.at(rl)); buffered |= bbb; } regions.at(rl) = buffered; - regions.at(rl).normalize(); + + if (veryverbose) { + cout << "Refinement level " << rl << ": enlarged regions to " << regions.at(rl) << endl; + } } @@ -545,6 +587,10 @@ namespace CarpetRegrid2 { // // TODO: Check after clipping { + if (veryverbose) { + cout << "Refinement level " << rl << ": checking whether regions should be combined..." << endl; + } + ibbox single; for (ibboxset::const_iterator ibb = regions.at(rl).begin(); ibb != regions.at(rl).end(); @@ -563,15 +609,30 @@ namespace CarpetRegrid2 { if (regions_size >= min_fraction * single_size) { // Combine the boxes regions.at(rl) = ibboxset (single); + if (veryverbose) { + cout << "Refinement level " << rl << ": combining regions to " << regions.at(rl) << endl; + } + } else { + if (veryverbose) { + cout << "Refinement level " << rl << ": not combining" << endl; + } } } // Find the location of the outer boundary + if (veryverbose) { + cout << "Refinement level " << rl << ": determining outer boundary..." << endl; + } + rvect const level_physical_lower = physical_lower; rvect const level_physical_upper = physical_upper; rvect const level_spacing = spacing / rvect (hh.reffacts.at(rl)); + if (veryverbose) { + cout << "Refinement level " << rl << ": physical coordinate boundary is at " << r2vect(level_physical_lower, level_physical_upper) << endl; + cout << "Refinement level " << rl << ": spacing is " << level_spacing << endl; + } #if 0 rvect level_interior_lower, level_interior_upper; rvect level_exterior_lower, level_exterior_upper; @@ -590,9 +651,28 @@ namespace CarpetRegrid2 { calculate_exterior_boundary (level_physical_lower, level_physical_upper, level_exterior_lower, level_exterior_upper, level_spacing); + if (veryverbose) { + cout << "Refinement level " << rl << ": exterior coordinate boundary is at " << r2vect(level_exterior_lower, level_exterior_upper) << endl; + } ibbox const & baseextent = hh.baseextents.at(0).at(rl); ivect const istride = baseextent.stride(); + if (veryverbose) { + cout << "Refinement level " << rl << ": stride is " << istride << endl; + } + + // This is the location of the outermost grid points. For cell + // centring, these are 1/2 grid spacing inside of the boundary. + ivect const level_physical_ilower = + rpos2ipos (physical_lower, origin, scale, hh, rl); + ivect const level_physical_iupper = + rpos2ipos1 (physical_upper, origin, scale, hh, rl); + if (veryverbose) { + cout << "Refinement level " << rl << ": physical boundary is at " << i2vect(level_physical_ilower, level_physical_iupper) << endl; + cout << "Refinement level " << rl << ": reconstructed physical coordinate boundary is at " + << r2vect(ipos2rpos(level_physical_ilower - (hh.refcent == cell_centered ? istride/2 : 0), origin, scale, hh, rl), + ipos2rpos(level_physical_iupper + (hh.refcent == cell_centered ? istride/2 : 0), origin, scale, hh, rl)) << endl; + } ivect const level_exterior_ilower = rpos2ipos (level_exterior_lower, origin, scale, hh, rl); @@ -600,6 +680,12 @@ namespace CarpetRegrid2 { rpos2ipos1 (level_exterior_upper, origin, scale, hh, rl); assert (all (level_exterior_ilower >= baseextent.lower())); assert (all (level_exterior_iupper <= baseextent.upper())); + if (veryverbose) { + cout << "Refinement level " << rl << ": exterior boundary is at " << i2vect(level_exterior_ilower, level_exterior_iupper) << endl; + cout << "Refinement level " << rl << ": reconstructed exterior coordinate boundary is at " + << r2vect(ipos2rpos(level_exterior_ilower, origin, scale, hh, rl), + ipos2rpos(level_exterior_iupper, origin, scale, hh, rl)) << endl; + } // Find the minimum necessary distance away from the outer // boundary due to buffer and ghost zones. This is e.g. the @@ -619,18 +705,28 @@ namespace CarpetRegrid2 { // Make the boxes rotating-90 symmetric // if (symmetry_rotating90) { + if (veryverbose) { + cout << "Refinement level " << rl << ": making regions rotating-90 symmetric" << endl; + } + ibboxset added; for (ibboxset::const_iterator ibb = regions.at(rl).begin(); ibb != regions.at(rl).end(); ++ ibb) { ibbox const & bb = * ibb; + if (veryverbose) { + cout << " considering box " << bb << "..." << endl; + } bvect const lower_is_outside_lower = - bb.lower() - min_bnd_dist_away[0] * bb.stride() <= physical_ilower; + bb.lower() - min_bnd_dist_away[0] * bb.stride() <= level_physical_ilower; // Treat both x and y directions for (int dir=0; dir<=1; ++dir) { + if (veryverbose) { + cout << " symmetrising in " << "xy"[dir] << " direction..." << endl; + } if (lower_is_outside_lower[dir]) { ivect const ilo = bb.lower(); ivect const iup = bb.upper(); @@ -640,9 +736,11 @@ namespace CarpetRegrid2 { rvect const axis (physical_lower[0], physical_lower[1], physical_lower[2]); // z component is unused - ivect const iaxis0 = rpos2ipos (axis, origin, scale, hh, rl); + ivect const iaxis0 = + rpos2ipos (axis, origin, scale, hh, rl); assert (all (iaxis0 % istr == 0)); - ivect const iaxis1 = rpos2ipos1 (axis, origin, scale, hh, rl); + ivect const iaxis1 = + rpos2ipos1 (axis, origin, scale, hh, rl); assert (all (iaxis1 % istr == 0)); ivect const offset = iaxis1 - iaxis0; assert (all (offset % istr == 0)); @@ -674,13 +772,18 @@ namespace CarpetRegrid2 { //assert (new_bb.is_contained_in (baseextent)); added |= new_bb; + if (veryverbose) { + cout << " adding box " << new_bb << endl; + } } } } regions.at(rl) |= added; - regions.at(rl).normalize(); + if (veryverbose) { + cout << "Refinement level " << rl << ": symmetrised regions are " << regions.at(rl) << endl; + } } @@ -689,6 +792,9 @@ namespace CarpetRegrid2 { // Make the boxes rotating-180 symmetric // if (symmetry_rotating180) { + if (veryverbose) { + cout << "Refinement level " << rl << ": making regions rotating-180 symmetric" << endl; + } ibboxset added; for (ibboxset::const_iterator ibb = regions.at(rl).begin(); @@ -696,12 +802,18 @@ namespace CarpetRegrid2 { ++ ibb) { ibbox const & bb = * ibb; + if (veryverbose) { + cout << " considering box " << bb << "..." << endl; + } bvect const lower_is_outside_lower = - bb.lower() - min_bnd_dist_away[0] * bb.stride() <= physical_ilower; + bb.lower() - min_bnd_dist_away[0] * bb.stride() <= level_physical_ilower; // Treat x direction int const dir = 0; + if (veryverbose) { + cout << " symmetrising in " << "x"[dir] << " direction..." << endl; + } if (lower_is_outside_lower[dir]) { ivect const ilo = bb.lower(); ivect const iup = bb.upper(); @@ -713,9 +825,11 @@ namespace CarpetRegrid2 { rvect const axis (physical_lower[0], (physical_lower[1] + physical_upper[1]) / 2, physical_lower[2]); // z component is unused - ivect const iaxis0 = rpos2ipos (axis, origin, scale, hh, rl); + ivect const iaxis0 = + rpos2ipos (axis, origin, scale, hh, rl); assert (all ((iaxis0 - baseextent.lower()) % istr == 0)); - ivect const iaxis1 = rpos2ipos1 (axis, origin, scale, hh, rl); + ivect const iaxis1 = + rpos2ipos1 (axis, origin, scale, hh, rl); assert (all ((iaxis1 - baseextent.lower()) % istr == 0)); ivect const offset = iaxis1 - iaxis0; assert (all (offset % istr == 0)); @@ -744,12 +858,17 @@ namespace CarpetRegrid2 { //assert (new_bb.is_contained_in (baseextent)); added |= new_bb; + if (veryverbose) { + cout << " adding box " << new_bb << endl; + } } } regions.at(rl) |= added; - regions.at(rl).normalize(); + if (veryverbose) { + cout << "Refinement level " << rl << ": symmetrised regions are " << regions.at(rl) << endl; + } } // if symmetry_rotating180 @@ -758,32 +877,39 @@ namespace CarpetRegrid2 { // Clip at the outer boundary // { + if (veryverbose) { + cout << "Refinement level " << rl << ": clipping at outer boundary..." << endl; + } + ibboxset clipped; for (ibboxset::const_iterator ibb = regions.at(rl).begin(); ibb != regions.at(rl).end(); ++ ibb) { ibbox const & bb = * ibb; + if (veryverbose) { + cout << " considering box " << bb << "..." << endl; + } // Clip boxes that extend outside the boundary. Enlarge // boxes that are inside but too close to the outer // boundary. bvect const lower_is_outside_lower = - bb.lower() - min_bnd_dist_away[0] * bb.stride() <= physical_ilower; + bb.lower() - min_bnd_dist_away[0] * bb.stride() <= level_physical_ilower; // Remove bboxes that are completely outside. bvect const upper_is_outside_lower = - bb.upper() < physical_ilower; + bb.upper() < level_physical_ilower; // Enlarge bboxes that extend not far enough inwards. bvect const upper_is_almost_outside_lower = - bb.upper() < physical_ilower + min_bnd_dist_incl[0] * bb.stride(); + bb.upper() < level_physical_ilower + min_bnd_dist_incl[0] * bb.stride(); // Ditto for the upper boundary. bvect const upper_is_outside_upper = - bb.upper() + min_bnd_dist_away[1] * bb.stride() >= physical_iupper; + bb.upper() + min_bnd_dist_away[1] * bb.stride() >= level_physical_iupper; bvect const lower_is_outside_upper = - bb.lower() > physical_iupper; + bb.lower() > level_physical_iupper; bvect const lower_is_almost_outside_upper = - bb.lower() > physical_iupper - min_bnd_dist_incl[1] * bb.stride(); + bb.lower() > level_physical_iupper - min_bnd_dist_incl[1] * bb.stride(); assert (not any (lower_is_almost_outside_upper and lower_is_outside_lower)); @@ -792,6 +918,9 @@ namespace CarpetRegrid2 { if (any (upper_is_outside_lower or lower_is_outside_upper)) { // The box is completely outside. Ignore it. + if (veryverbose) { + cout << " box is completely outside; dropping it" << endl; + } continue; } @@ -806,8 +935,8 @@ namespace CarpetRegrid2 { << " lower_is_outside_lower=" << lower_is_outside_lower << " upper_is_outside_upper=" << upper_is_outside_upper << " boundary_staggering_mismatch=" << boundary_staggering_mismatch - << " physical_ilower=" << physical_ilower - << " physical_iupper=" << physical_iupper + << " level_physical_ilower=" << level_physical_ilower + << " level_physical_iupper=" << level_physical_iupper << " baseextent=" << baseextent; CCTK_WARN (CCTK_WARN_ABORT, msg.str().c_str()); } @@ -816,13 +945,13 @@ namespace CarpetRegrid2 { (either (lower_is_outside_lower, level_exterior_ilower, either (lower_is_almost_outside_upper, - (physical_iupper - + (level_physical_iupper - min_bnd_dist_incl[1] * bb.stride()), bb.lower())), either (upper_is_outside_upper, level_exterior_iupper, either (upper_is_almost_outside_lower, - (physical_ilower + + (level_physical_ilower + min_bnd_dist_incl[0] * bb.stride()), bb.upper())), bb.stride()); @@ -837,19 +966,24 @@ namespace CarpetRegrid2 { << " upper_is_almost_outside_lower=" << upper_is_almost_outside_lower << " level_exterior_ilower=" << level_exterior_ilower << " level_exterior_iupper=" << level_exterior_iupper - << " physical_ilower=" << physical_ilower - << " physical_iupper=" << physical_iupper + << " level_physical_ilower=" << level_physical_ilower + << " level_physical_iupper=" << level_physical_iupper << " baseextent=" << baseextent; CCTK_WARN (CCTK_WARN_ABORT, msg.str().c_str()); } assert (clipped_bb.is_contained_in (baseextent)); clipped |= clipped_bb; + if (veryverbose) { + cout << " Clipped box is " << clipped_bb << endl; + } } // for ibb regions.at(rl) = clipped; - regions.at(rl).normalize(); + if (veryverbose) { + cout << "Refinement level " << rl << ": clipped regions are " << regions.at(rl) << endl; + } } @@ -874,8 +1008,8 @@ namespace CarpetRegrid2 { ibbox const & bb = * ibb; assert (bb.is_contained_in (hh.baseextents.at(0).at(rl))); - bvect const lower_is_outer = bb.lower() <= physical_ilower; - bvect const upper_is_outer = bb.upper() >= physical_iupper; + bvect const lower_is_outer = bb.lower() <= level_physical_ilower; + bvect const upper_is_outer = bb.upper() >= level_physical_iupper; b2vect const ob (lower_is_outer, upper_is_outer); @@ -899,6 +1033,13 @@ namespace CarpetRegrid2 { assert (not regions.at(rl-1).empty()); ibbox const coarse0 = * regions.at(rl-1).begin(); + // This is the location of the outermost grid points. For cell + // centring, these are 1/2 grid spacing inside of the boundary. + ivect const level_physical_ilower = + rpos2ipos (physical_lower, origin, scale, hh, rl); + ivect const level_physical_iupper = + rpos2ipos1 (physical_upper, origin, scale, hh, rl); + i2vect const fdistance = dd.ghost_widths.at(rl); i2vect const cdistance = i2vect (min_distance + dd.prolongation_stencil_size(rl)); @@ -911,8 +1052,8 @@ namespace CarpetRegrid2 { { ibbox const & fbb = * ibb; - bvect const lower_is_outer = fbb.lower() <= physical_ilower; - bvect const upper_is_outer = fbb.upper() >= physical_iupper; + bvect const lower_is_outer = fbb.lower() <= level_physical_ilower; + bvect const upper_is_outer = fbb.upper() >= level_physical_iupper; b2vect const ob (lower_is_outer, upper_is_outer); ibbox const cbb = fbb.expanded_for (coarse0); @@ -969,7 +1110,7 @@ namespace CarpetRegrid2 { } } - if (verbose) { + if (verbose or veryverbose) { if (do_recompose) { for (int n = 0; n < num_centres; ++ n) { if (active[n]) { @@ -1051,7 +1192,7 @@ namespace CarpetRegrid2 { if (do_recompose) break; } // for n - if (verbose) { + if (verbose or veryverbose) { if (not do_recompose) { CCTK_INFO ("Refined regions have not changed sufficiently; skipping regridding"); @@ -1141,7 +1282,7 @@ namespace CarpetRegrid2 { } } - if (verbose) { + if (verbose or veryverbose) { if (do_recompose) { for (int n = 0; n < num_centres; ++ n) { if (active[n]) { @@ -1223,7 +1364,7 @@ namespace CarpetRegrid2 { if (do_recompose) break; } // for n - if (verbose) { + if (verbose or veryverbose) { if (not do_recompose) { CCTK_INFO ("Refined regions have not changed sufficiently; skipping regridding"); diff --git a/Carpet/CarpetSlab/src/GetHyperslab.cc b/Carpet/CarpetSlab/src/GetHyperslab.cc index 8f86d8568..9b1e08f77 100644 --- a/Carpet/CarpetSlab/src/GetHyperslab.cc +++ b/Carpet/CarpetSlab/src/GetHyperslab.cc @@ -186,9 +186,8 @@ namespace CarpetSlab { mydata = (*myff)(tl, rl, component, mglevel); // Calculate overlapping extents - bboxset<int,dim> myextents = - mydd->boxes.at(mglevel).at(rl).at(component).interior & hextent; - myextents.normalize(); + bboxset<int,dim> const myextents = + mydd->light_boxes.at(mglevel).at(rl).at(component).interior & hextent; // Loop over overlapping extents for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin(); diff --git a/Carpet/CarpetSlab/src/slab.cc b/Carpet/CarpetSlab/src/slab.cc index b20566d47..d9d587a98 100644 --- a/Carpet/CarpetSlab/src/slab.cc +++ b/Carpet/CarpetSlab/src/slab.cc @@ -220,7 +220,7 @@ namespace CarpetSlab { // Calculate overlapping extents const bboxset<int,dim> myextents = - mydd->boxes.at(mglevel).at(rl).at(component).interior & hextent; + mydd->light_boxes.at(mglevel).at(rl).at(component).interior & hextent; // Loop over overlapping extents for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin(); diff --git a/Carpet/CarpetWeb/Makefile b/Carpet/CarpetWeb/Makefile index 7a7d2ad67..ea359dfd7 100644 --- a/Carpet/CarpetWeb/Makefile +++ b/Carpet/CarpetWeb/Makefile @@ -1,7 +1,6 @@ all: sync: - rsync -a -v -z --exclude .DS_Store --exclude "*~" --exclude doxygen --exclude Makefile --delete --delete-excluded ./ carpet@carpetcode.dyndns.org:www.carpetcode.org/htdocs - rsync -a -v -z --exclude .DS_Store --exclude "*~" --exclude doxygen --exclude Makefile --delete --delete-excluded -e ssh --rsync-path=/home/schnette/bin/rsync ./ carpet@www.carpetcode.org:www.carpetcode.org/htdocs + rsync -a -v -z --exclude .DS_Store --exclude "*~" --exclude doxygen --exclude Makefile --delete --delete-excluded ./ carpet@www.carpetcode.org:www.carpetcode.org/htdocs .PHONY: all sync diff --git a/Carpet/CarpetWeb/get-carpet.html b/Carpet/CarpetWeb/get-carpet.html index 319d8291f..3783c958d 100644 --- a/Carpet/CarpetWeb/get-carpet.html +++ b/Carpet/CarpetWeb/get-carpet.html @@ -113,7 +113,7 @@ <p>The development version of Carpet is available via <a href="http://git.or.cz/">git</a>:</p> <pre> cd Cactus - git clone -o carpet git://carpetcode.dyndns.org/carpet.git + git clone -o carpet git://carpetcode.org/carpet.git cd arrangements ln -s ../carpet/Carpet* .</pre> <p>(Don't miss the dot after the <code>Carpet*</code> in the last @@ -121,7 +121,7 @@ given <a href="#git">below</a>.</p> <!-- This doesn't work yet You can also have a look at - the <a href="http://carpetcode.dyndns.org/~carpet/git/">development + the <a href="http://www.carpetcode.org/~carpet/git/">development source tree</a> in your web browser.</p> --> @@ -163,7 +163,7 @@ ssh. Once you have an account set up, you obtain e.g. the development version with</p> <pre> cd Cactus - git clone carpetgit@carpetcode.dyndns.org:carpet.git + git clone carpetgit@carpetcode.org:carpet.git cd arrangements ln -s ../carpet/Carpet* .</pre> <p>(Don't miss the dot after the <code>Carpet*</code> in the last diff --git a/Carpet/CarpetWeb/publications/carpet-publications.bib b/Carpet/CarpetWeb/publications/carpet-publications.bib index 81cbcafd9..30b7f0f40 100644 --- a/Carpet/CarpetWeb/publications/carpet-publications.bib +++ b/Carpet/CarpetWeb/publications/carpet-publications.bib @@ -1232,7 +1232,7 @@ Article{CarpetResult-Grandclement2007a, doi = {10.1103/PhysRevD.79.084010}, } -@Article{Carpet-Baiotti2008b, +@Article{Carpet-Baiotti2008b, status = {refereed}, author = {Luca Baiotti and Sebastiano Bernuzzi and Giovanni Corvino and Roberto De Pietri and Alessandro Nagar}, @@ -1250,6 +1250,22 @@ Article{CarpetResult-Grandclement2007a, doi = {10.1103/PhysRevD.79.024002}, } +@Article{Carpet-Ott2008a, + status = {refereed}, + author = {Christian D. Ott}, + title = {The gravitational-wave signature of core-collapse + supernovae}, + journal = {Class. Quantum Grav.}, + year = 2009, + volume = 26, + pages = 063001, + annote = {Class. Quantum Grav. highlight of 2008/2009}, + eprint = {arXiv:0809.0695 [astro-ph]}, + url = {http://arxiv.org/abs/0809.0695}, + receiveddate = {2008-09-03}, + doi = {10.1088/0264-9381/26/6/063001} +} + @Article{Carpet-Brown2007b, status = {refereed}, author = {David Brown and Peter Diener and Olivier Sarbach and @@ -1722,6 +1738,19 @@ Article{Carpet-Cadonati2009b, doi = {10.1088/0264-9381/27/7/075014} } +@Article{Carpet-Zink2010a, + status = {refereed}, + author = {Burkhard Zink and Oleg Korobkin and Erik Schnetter + and Nikolaos Stergioulas}, + title = {On the frequency band of the $f$-mode {CFS} + instability}, + journal = {Phys. Rev. D (accepted)}, + year = 2010, + eprint = {arXiv:1003.0779 [astro-ph.SR]}, + url = {http://arxiv.org/abs/1003.0779}, + receiveddate = {2010-03-03}, +} + ******************************************************************************** @@ -2231,6 +2260,19 @@ REPORT: receiveddate = {2009-12-17}, } +@Unpublished{Carpet-Kelly2009a, + status = {report}, + author = {Bernard J. Kelly and Wolfgang Tichy and Yosef + Zlochower and Manuela Campanelli and Bernard + Whiting}, + title = {Post-Newtonian Initial Data with Waves: Progress in + Evolution}, + note = {arXiv:0912.5311 [gr-qc]}, + year = 2009, + url = {http://arxiv.org/abs/0912.5311}, + receiveddate = {2009-12-29}, +} + @Unpublished{Carpet-Zilhao2010a, status = {report}, author = {Miguel Zilh{\~a}o and Helvi Witek and Ulrich diff --git a/Carpet/LoopControl/param.ccl b/Carpet/LoopControl/param.ccl index 229d1f489..f529229e2 100644 --- a/Carpet/LoopControl/param.ccl +++ b/Carpet/LoopControl/param.ccl @@ -15,30 +15,6 @@ BOOLEAN debug "Output debug information" STEERABLE=always { } "no" -###################################### -# LoopControl demonstration parameters - -BOOLEAN run_demo "Run the embedded wavetoy as a LoopControl demo application" -{ -} "no" - -CCTK_INT nx "Number of grid points in X dimension (for the embedded wavetoy example)" -{ - 1:* :: "a positive integer" -} 100 -CCTK_INT ny "Number of grid points in Y dimension (for the embedded wavetoy example)" -{ - 1:* :: "a positive integer" -} 100 -CCTK_INT nz "Number of grid points in Z dimension (for the embedded wavetoy example)" -{ - 1:* :: "a positive integer" -} 100 -CCTK_INT nsteps "Number of time steps (for the embedded wavetoy example)" -{ - 1:* :: "a positive integer" -} 100 - ################# diff --git a/Carpet/LoopControl/schedule.ccl b/Carpet/LoopControl/schedule.ccl index 1a7c28d10..fc14c7b0d 100644 --- a/Carpet/LoopControl/schedule.ccl +++ b/Carpet/LoopControl/schedule.ccl @@ -6,11 +6,3 @@ if (printstats) { LANG: C } "Output loop control statistics" } - -if (run_demo) -{ - SCHEDULE lc_demo AT startup - { - LANG:C - } "Run embedded wavetoy as a LoopControl demo application" -} diff --git a/Carpet/LoopControl/src/make.code.defn b/Carpet/LoopControl/src/make.code.defn index 961ce4edd..426f7ee52 100644 --- a/Carpet/LoopControl/src/make.code.defn +++ b/Carpet/LoopControl/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn LoopControl # Source files in this directory -SRCS = loopcontrol.c loopcontrol.F90 loopcontrol_types.F90 lc_auto.c lc_siman.c lc_hill.c wavetoy-loopcontrol.c +SRCS = loopcontrol.c loopcontrol.F90 loopcontrol_types.F90 lc_auto.c lc_siman.c lc_hill.c # Subdirectories containing source files SUBDIRS = |