Commit 595b4b83 authored by Christof Metzger-Kraus's avatar Christof Metzger-Kraus
Browse files

All changes from project 'OPAL3D'. Couldn't keep DKS stuff in

ParallelTTracker. Particle-Matter interaction doesn't work either at the moment.

I hope that I didn't create too much of a mess when merging!
parent dfc106b6
CMAKE_MINIMUM_REQUIRED (VERSION 2.8.10)
PROJECT (OPAL)
SET (OPAL_VERSION_MAJOR 1)
SET (OPAL_VERSION_MINOR 5.00.2)
SET (OPAL_VERSION_MINOR 9.00)
set (PACKAGE \"opal\")
set (PACKAGE_BUGREPORT \"opal@lists.psi.ch\")
set (PACKAGE_NAME \"OPAL\")
set (PACKAGE_TARNAME \"opal\")
set (PACKAGE_VERSION "\"${OPAL_VERSION_MAJOR}.${OPAL_VERSION_MINOR}\"")
set (OPAL_VERSION "${OPAL_VERSION_MAJOR}.${OPAL_VERSION_MINOR}")
STRING (REGEX REPLACE "\\.([0-9])\\." ".0\\1." PACKAGE_VERSION ${OPAL_VERSION})
set (PACKAGE_VERSION_STR "\"${PACKAGE_VERSION}\"")
STRING (REGEX REPLACE "\\." "" PACKAGE_VERSION ${PACKAGE_VERSION})
IF (NOT CMAKE_CONFIGURATION_TYPES AND NOT CMAKE_BUILD_TYPE)
SET (CMAKE_BUILD_TYPE RelWithDebInfo CACHE STRING
......@@ -61,7 +64,7 @@ set (BOOSTROOT $ENV{BOOST_DIR})
set(Boost_USE_STATIC_LIBS ON)
set(Boost_USE_MULTITHREADED OFF)
set(Boost_USE_STATIC_RUNTIME OFF)
find_package (Boost 1.55.0 REQUIRED COMPONENTS regex filesystem system)
find_package (Boost 1.54.0 REQUIRED COMPONENTS regex filesystem system iostreams)
if (Boost_INCLUDE_DIRS)
message (STATUS "Found boost include dir: ${Boost_INCLUDE_DIR}")
message (STATUS "Found boost library dir: ${Boost_LIBRARY_DIR}")
......@@ -85,7 +88,7 @@ IF (ENABLE_DKS)
### CUDA compiler flags ###
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -lcudart -lcufft -lcublas -lnvToolsExt -DDKS_CUDA")
### if any accelerator enabled set flag to use DKS ###
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DIPPL_DKS -DIPPL_DKS_CUDA -DOPAL_DKS")
ENDIF (ENABLE_DKS)
......@@ -105,6 +108,11 @@ OPTION (NOCPLUSPLUS11_NULLPTR "Disable C++11 nullptr support" OFF)
OPTION (NO_FIELD_ASSIGN_OPTIMIZATION "Disable compiler optimization of IPPL field assignment" OFF)
IF (BUILD_OPAL_UNIT_TESTS)
FIND_PACKAGE (GTest REQUIRED)
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DWITH_UNIT_TESTS")
ENDIF (BUILD_OPAL_UNIT_TESTS)
IF (USE_H5HUT2)
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DUSE_H5HUT2")
ENDIF (USE_H5HUT2)
......@@ -161,8 +169,8 @@ STRING(REGEX MATCH "[^0-9]*" HOSTNAME_BASE "${HOSTNAME}")
# cray-tpsl/16.06.1 libraries -> Linker error.
# We can avoid this issue by not going into the if-statement
if (NOT ${HOSTNAME_BASE} MATCHES "edison" AND NOT ${HOSTNAME_BASE} MATCHES "cori" AND (ENABLE_SAAMG_SOLVER OR ENABLE_AMR_SOLVER) )
MESSAGE ("Enable SAAMG_SOLVER " ${ENABLE_SAAMG_SOLVER})
MESSAGE ("Enable AMR_SOLVER " ${ENABLE_AMR_SOLVER})
MESSAGE (STATUS "Enable SAAMG_SOLVER " ${ENABLE_SAAMG_SOLVER})
MESSAGE (STATUS "Enable AMR_SOLVER " ${ENABLE_AMR_SOLVER})
find_package (Trilinos REQUIRED HINTS $ENV{TRILINOS_PREFIX} $ENV{TRILINOS_DIR} $ENV{TRILINOS})
message (STATUS "Found Trilinos: ${Trilinos_DIR}")
......@@ -237,6 +245,8 @@ ENDIF()
SET (IPPL_USED_FROM_OPAL "TRUE")
ADD_SUBDIRECTORY (ippl)
ADD_SUBDIRECTORY (src)
ADD_SUBDIRECTORY (tools)
ADD_SUBDIRECTORY (doc/user_guide)
IF (BUILD_OPAL_UNIT_TESTS)
ADD_SUBDIRECTORY(tests)
......@@ -249,4 +259,4 @@ set(OPAL_CXX_FLAGS ${OPAL_CXX_FLAGS}
set(OPAL_LIBS ${OPAL_LIBS}
CACHE INTERNAL "" FORCE
)
)
\ No newline at end of file
......@@ -21,7 +21,7 @@ elseif(FortranCInterface_GLOBAL_SUFFIX STREQUAL "_" AND FortranCInterface_GLOBAL
# message(STATUS " Fortran name mangling scheme to UNDERSCORE (lower case, append underscore)")
set(BL_FORTLINK "UNDERSCORE")
#else()
# message(AUTHOR_WARNING "Fortran to C mangling not backward compatible with older style BoxLib code")
# message(AUTHOR_WARNING "Fortran to C mangling not backward compatible with older style BoxLib code")
endif()
set(BL_MACHINE ${CMAKE_SYSTEM_NAME})
......
......@@ -6,7 +6,7 @@
# Control the search through CCSE_DIR or setting environment variable
# CCSE_ROOT to the ccse installation prefix.
#
# This module does not search default paths!
# This module does not search default paths!
#
# Following variables are set:
# CCSE_FOUND (BOOL) Flag indicating if CCSE was found
......@@ -46,12 +46,12 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
if(CCSE_LIBRARY_DIR)
set(CCSE_LIBRARY_DIR "${CCSE_LIBRARY_DIR}" CACHE PATH "Path to search for CCSE library files")
endif()
if(CCSE_PERL_DIR)
set(CCSE_PERL_DIR "${CCSE_PERL_DIR}" CACHE PATH "Path to search for CCSE perl scripts")
endif()
# Search for include files
# Search order preference:
# (1) CCSE_INCLUDE_DIR - check existence of path AND if the include files exist
......@@ -77,7 +77,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
set(CCSE_INCLUDE_DIR "CCSE_INCLUDE_DIR-NOTFOUND")
endif()
else()
else()
set(ccse_inc_suffixes "include")
if(CCSE_DIR)
......@@ -93,7 +93,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
else()
message(SEND_ERROR "CCSE_DIR=${CCSE_DIR} does not exist")
set(CCSE_INCLUDE_DIR "CCSE_INCLUDE_DIR-NOTFOUND")
endif()
endif()
else()
......@@ -112,7 +112,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
set(CCSE_INCLUDE_DIRS ${CCSE_INCLUDE_DIR})
endif()
# Search for libraries
# Search for libraries
# Search order preference:
# (1) CCSE_LIBRARY_DIR - check existence of path AND if the include files exist
# (2) CCSE_DIR/<lib,Lib>
......@@ -132,7 +132,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
set(CCSE_LIBRARY "CCSE_LIBRARY-NOTFOUND")
endif()
else()
else()
if(CCSE_DIR)
......@@ -148,7 +148,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
else()
message(SEND_ERROR "CCSE libs not found in CCSE_DIR/lib${CCSE_LIBDIR_MPI_SUFFIX}${CCSE_LIBDIR_OMP_SUFFIX}, (CCSE_DIR=${CCSE_DIR})")
set(CCSE_LIBRARY "CCSE_LIBRARY-NOTFOUND")
endif()
endif()
else()
......@@ -156,7 +156,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
NAMES ${ccse_lib_name})
get_filename_component(CCSE_LIBRARY_DIR "${CCSE_LIBRARY}" PATH)
endif()
endif()
......@@ -199,7 +199,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
set(CCSE_PERL "CCSE_PERL-NOTFOUND")
endif()
else()
else()
if(CCSE_DIR)
......@@ -214,7 +214,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
message(SEND_ERROR "CCSE Perl scripts not in CCSE_DIR/perl=${CCSE_DIR}/perl")
else()
set(CCSE_PERL_DIR ${CCSE_PERL})
endif()
endif()
else()
......@@ -230,7 +230,7 @@ else(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
set(CCSE_EXT_LIBRARIES "")
endif(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
endif(CCSE_LIBRARIES AND CCSE_INCLUDE_DIRS AND CCSE_PERL_DIR)
# Send useful message if everything is found
find_package_handle_standard_args(CCSE DEFAULT_MSG
......
......@@ -765,7 +765,8 @@ WARN_LOGFILE =
# Note: If this tag is empty the current directory is searched.
INPUT = ./src \
./ippl
./src/Classic \
./ippl/src
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
......@@ -994,7 +995,6 @@ USE_HTAGS = NO
# The default value is: YES.
VERBATIM_HEADERS = YES
#---------------------------------------------------------------------------
# Configuration options related to the alphabetical class index
#---------------------------------------------------------------------------
......@@ -1772,7 +1772,6 @@ RTF_EXTENSIONS_FILE =
# classes and files.
# The default value is: NO.
GENERATE_MAN = NO
# The MAN_OUTPUT tag is used to specify where the man pages will be put. If a
# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
......@@ -1836,7 +1835,7 @@ XML_OUTPUT = xml
XML_PROGRAMLISTING = YES
#---------------------------------------------------------------------------
# Configuration options related to the DOCBOOK output
# Configuration options for the AutoGen Definitions output
#---------------------------------------------------------------------------
# If the GENERATE_DOCBOOK tag is set to YES doxygen will generate Docbook files
......@@ -1983,7 +1982,7 @@ EXPAND_AS_DEFINED =
SKIP_FUNCTION_MACROS = YES
#---------------------------------------------------------------------------
# Configuration options related to external references
# Configuration additions related to external references
#---------------------------------------------------------------------------
# The TAGFILES tag can be used to specify one or more tag files. For each tag
......@@ -2326,4 +2325,4 @@ GENERATE_LEGEND = YES
# The default value is: YES.
# This tag requires that the tag HAVE_DOT is set to YES.
DOT_CLEANUP = YES
DOT_CLEANUP = YES
\ No newline at end of file
OPAL
\ No newline at end of file
OPAL
\ No newline at end of file
......@@ -6,8 +6,7 @@ echo "#define GIT_VERSION \"$REV\"" > revision.tmp
diff revision.tmp src/revision.h 2>&1 > /dev/null
if [ $? -eq 1 ];then
cp revision.tmp src/revision.h
mv revision.tmp classic/5.0/src/revision.h
mv revision.tmp src/revision.h
else
rm revision.tmp
fi
\ No newline at end of file
......@@ -6,8 +6,7 @@ echo "#define GIT_VERSION \"$REV\"" > revision.tmp
diff revision.tmp src/revision.h 2>&1 > /dev/null
if [ $? -eq 1 ];then
cp revision.tmp src/revision.h
mv revision.tmp classic/5.0/src/revision.h
mv revision.tmp src/revision.h
else
rm revision.tmp
fi
\ No newline at end of file
......@@ -6,8 +6,7 @@ echo "#define GIT_VERSION \"$REV\"" > revision.tmp
diff revision.tmp src/revision.h 2>&1 > /dev/null
if [ $? -eq 1 ];then
cp revision.tmp src/revision.h
mv revision.tmp classic/5.0/src/revision.h
mv revision.tmp src/revision.h
else
rm revision.tmp
fi
\ No newline at end of file
......@@ -6,8 +6,7 @@ echo "#define GIT_VERSION \"$REV\"" > revision.tmp
diff revision.tmp src/revision.h 2>&1 > /dev/null
if [ $? -eq 1 ];then
cp revision.tmp src/revision.h
mv revision.tmp classic/5.0/src/revision.h
mv revision.tmp src/revision.h
else
rm revision.tmp
fi
\ No newline at end of file
......@@ -76,7 +76,7 @@ IF (ENABLE_DKS)
### CUDA compiler flags ###
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -lcudart -lcufft -lcublas -lnvToolsExt -DDKS_CUDA")
### if any accelerator enabled set flag to use DKS ###
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DIPPL_DKS -DIPPL_DKS_CUDA -DOPAL_DKS")
ENDIF (ENABLE_DKS)
......@@ -99,19 +99,19 @@ OPTION (NO_FIELD_ASSIGN_OPTIMIZATION "Disable compiler optimization of IPPL fiel
IF (ENABLE_AMR_SOLVER)
# set specific parameters for BoxLib used in OPAL
ENABLE_LANGUAGE (Fortran)
SET ( CCSE_DIR $ENV{BOXLIB_PREFIX} )
SET ( CCSE_LIBRARY_DIR $ENV{BOXLIB_LIBRARY_DIR} )
FIND_PACKAGE (CCSE REQUIRED)
IF (CCSE_FOUND)
MESSAGE (STATUS "Found BoxLib include dir: ${CCSE_INCLUDE_DIR}")
MESSAGE (STATUS "Found BoxLib library dir: ${CCSE_LIBRARY_DIR}")
MESSAGE (STATUS "Found BoxLib perl dir: ${CCSE_PERL_DIR}")
INCLUDE_DIRECTORIES (${CCSE_INCLUDE_DIR})
ENDIF (CCSE_FOUND)
SET (BL_SPACEDIM 3 CACHE INT "Dimension of BoxLib build")
SET (ENABLE_MPI 1 CACHE INT "Enable build with MPI")
SET (ENABLE_OpenMP 0 CACHE INT "Enable build with OpenMP")
......@@ -119,10 +119,10 @@ IF (ENABLE_AMR_SOLVER)
SET (BL_USE_PARTICLES 1 CACHE INT "Include Particles classes in BoxLib build")
SET (ENABLE_PROFILING 0 CACHE INT "Include profiling information in BoxLib build")
SET (ENABLE_BACKTRACE 0 CACHE INT "Include backtrace information in BoxLib build")
# disable due to BoxLib
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unused-variable")
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DBL_PRECISION=${BL_PRECISION}")
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DENABLE_MPI=${ENABLE_MPI}")
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DENABLE_OpenMP=${ENABLE_OpenMP}")
......@@ -130,7 +130,7 @@ IF (ENABLE_AMR_SOLVER)
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DBL_DEBUG=${BL_DEBUG}")
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DENABLE_BACKTRACE=${ENABLE_BACKTRACE}")
SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DENABLE_PROFILING=${ENABLE_PROFILING}")
MESSAGE (STATUS "Use following BoxLib settings:")
MESSAGE (STATUS " BL_SPACEDIM = ${BL_SPACEDIM} (INT: 1, 2, 3)")
MESSAGE (STATUS " BL_MACHINE = ${BL_MACHINE} (STRING: <ARCH>)")
......@@ -217,8 +217,8 @@ ENDIF()
SET (IPPL_USED_FROM_OPAL "TRUE")
ADD_SUBDIRECTORY (ippl)
ADD_SUBDIRECTORY (classic/5.0/src)
ADD_SUBDIRECTORY (src)
#ADD_SUBDIRECTORY (classic/5.0/src)
#ADD_SUBDIRECTORY (src)
IF (BUILD_OPAL_UNIT_TESTS)
ADD_SUBDIRECTORY(tests)
......@@ -279,45 +279,6 @@ SET (IPPL_CXX_FLAGS
"-DIPPL_MPI -DMPICH_SKIP_MPICXX -DIPPL_DONT_POOL -DIPPL_USE_XDIV_RNG -DPETE_BITWISE_COPY -DIPPL_HAS_TEMPLATED_COMPLEX -DIPPL_USE_PARTIAL_SPECIALIZATION -Drestrict=__restrict__ -DNOCTAssert ${IPPL_CXX_FLAGS}"
)
### CXX_FLAGS common to all compilers and platforms which ###
### may *not* be written to the configure file ###
SET (OTHER_CXX_FLAGS
......@@ -443,4 +404,4 @@ IF (LINUXGCC)
RENAME ${PROJECT_NAME}Config.cmake
)
ENDIF (NOT IPPL_USED_FROM_OPAL)
ENDIF (NOT IPPL_USED_FROM_OPAL)
\ No newline at end of file
\chapter{Introduction}
\label{sec:Introduction}
One of \ippl 's most attractive features is its high performance on both single-processor and distributed-memory multiprocessor machines. As future releases of the library will also support shared-memory machines.
One of \ippl 's most attractive features is its high performance on both single-processor and distributed-memory multiprocessor machines. As future releases of the library will also support shared-memory machines.
The heart of the problem \ippl 's authors face is that while data-parallel programming is a natural way to express many scientific and numerical algorithms, straightforward implementations of it do exactly the wrong thing on modern architectures, whose performance depends critically on the re-use of data loaded into cache. If a program evaluates A+B+C for three arrays A, B, and C by adding A to B, then adding C to that calculation's result, performance suffers both because of the overhead of executing two loops instead of one, but also (and more importantly) because every value in the temporary array that stores the result of A+B has to be accessed twice: once to write it, and once to read it back in. As soon as this array is too large to fit into cache, the program's performance will drop dramatically.
......@@ -27,8 +27,8 @@ int main(int argc, char *argv[])
const unsigned N=8;
const unsigned Dim=2;
Index IGLOBAL(N); // Specify the global domain
Index JGLOBAL(N);
Index IGLOBAL(N); // Specify the global domain
Index JGLOBAL(N);
Index I(1, N-1); // Specify the interior domain
Index J(1, N-1);
......@@ -41,7 +41,7 @@ int main(int argc, char *argv[])
assign(A,0.0); // Assign initial conditions
assign(b,0.0);
b[N/2][N/2] = -1.0; // put a spike on the RHS
b[N/2][N/2] = -1.0; // put a spike on the RHS
double fact = 0.25;
// Iterate 200 times
......@@ -50,50 +50,50 @@ int main(int argc, char *argv[])
A[I-1][J] +
A[I][J+1] +
A[I][J-1] - b[I][J]));
}
}
msg << A << endl;
return 0;
}
\end{code}
The syntax is very similar to that of Fortran 90: a single assignment fills an entire array with a scalar value, subscripts express ranges as well as single points, and so on. In fact, the combination of C++ and \ippl provides so many of the features of Fortran 90 that one might well ask whether it wouldn't better to just use the latter language straight up. One answer comes down to economics. While the various flavors of Fortran are still used in scientific computing, Fortran's user base is shrinking, particularly in comparison to C++. Networking, graphics, database access, and operating system interfaces are available in C++ programmers long before they're available in Fortran (if they become available at all). What's more, support tools such as debuggers and memory inspectors are primarily targeted at C++ developers, as are hundreds of books, journal articles, and web sites.
Another answer is that the abstraction facilities of C++ are much more powerful that those in Fortran. While Fortran 90 supports an attractive array syntax for floating point arrays one could not, for example, efficiently extend this high level syntax to arrays of vectors or tensors. Until recently, Fortran has had two powerful arguments in its favor: legacy applications, and performance. However, the importance of the former is diminishing as the invention of new algorithms force programmers to rewrite old codes, while the invention of techniques such as {\it expression templates} has made it possible for C++ programs to match, or exceed, the performance of highly-optimized Fortran 77.
Another answer is that the abstraction facilities of C++ are much more powerful that those in Fortran. While Fortran 90 supports an attractive array syntax for floating point arrays one could not, for example, efficiently extend this high level syntax to arrays of vectors or tensors. Until recently, Fortran has had two powerful arguments in its favor: legacy applications, and performance. However, the importance of the former is diminishing as the invention of new algorithms force programmers to rewrite old codes, while the invention of techniques such as {\it expression templates} has made it possible for C++ programs to match, or exceed, the performance of highly-optimized Fortran 77.
%!TEX encoding = UTF-8 Unicode
\section{Example 2 Power Spectrum}
A sinussoidal field $\rho(i,j,k) = a_1sin(k_1 \frac{2\pi}{n_x} i) + a_5 sin(k_5 \frac{2\pi}{n_x} i)$,
A sinussoidal field $\rho(i,j,k) = a_1sin(k_1 \frac{2\pi}{n_x} i) + a_5 sin(k_5 \frac{2\pi}{n_x} i)$,
$i= 1 \dots n_x, ~ j= 1 \dots n_y, ~ k= 1 \dots n_z $ with $n_x,n_y$ and $n_z$ denoting the grid size is generated and
the power spectrum calculated. This examples shows how to initialise fields, compute discrete complex-complex FFT and
the power spectrum calculated. This examples shows how to initialise fields, compute discrete complex-complex FFT and
compute the resulting powerspectrum.
Assume a real density field is defined like
Assume a real density field is defined like
\begin{smallcode}
typedef Field<double,Dim,Mesh_t,Center_t> Field_t;
Field_t rho;
\end{smallcode}
we then can immediately initialize the field according to the above formula
we then can immediately initialize the field according to the above formula
\begin{smallcode}
assign(rho[I][J][K], a1*sin(2.0*pi/nr_m[0]*k1*I) +
assign(rho[I][J][K], a1*sin(2.0*pi/nr_m[0]*k1*I) +
a5*sin(2.0*pi/nr_m[0]*k5*I));
\end{smallcode}
Normalizing to $\max(\rho) \le 1.0 $ with
\begin{smallcode}
rho /= max(rho)
\end{smallcode}
we then assume to have defined a complex field "fC" and a complex-complex FFT.
\end{smallcode}
we then assume to have defined a complex field "fC" and a complex-complex FFT.
\begin{smallcode}
fC = rho;
fft->transform("forward" , fC);
\end{smallcode}
Here we used the in place version of the FFT to obtain $\rho $ in Fourier space. Now
we can compute the power spectrum:
\begin{smallcode}
we can compute the power spectrum:
\begin{smallcode}
pwrSpec = real(fC*conj(fC));
\end{smallcode}
%\pagebreak
and calculate the 1D pwr-spectrum (in x direction) by integrating over y and z: \\
\begin{code}
NDIndex<3> elem;
\begin{code}
NDIndex<3> elem;
for (int i=lDomain[0].min(); i<=(lDomain[0].max()-1)/2; ++i) {
elem[0]=Index(i,i);
for (int j=lDomain[1].min(); j<=(lDomain[1].max()-1)/2; ++j) {
......@@ -107,24 +107,24 @@ and calculate the 1D pwr-spectrum (in x direction) by integrating over y and z:
\end{code}
The power spectra of the local domain is stored in $f1$. We have to update all other node
so that each node has the full power spectrum by:
\begin{smallcode}
\begin{smallcode}
reduce(&(f1[0]),&(f1[0])+f1_lenght,OpAddAssign());
\end{smallcode}
\end{smallcode}
assuming the non local part of $f1$ is initialized with zero.
%The full code pwrspec-1.cpp is located at $\ippl_ROOT/test/simple.
\section{Example 3 Particle in Cell Code (PIC)}
This example discusses how to write a 3D Particle in Cell Code (PIC). The
complete source file can be found at {\em \$IPPL\_ROOT/test/particles}. The
This example discusses how to write a 3D Particle in Cell Code (PIC). The
complete source file can be found at {\em \$IPPL\_ROOT/test/particles}. The
this presentation details are omitted, only the structure and important issues are
highlighted.
highlighted.
\subsection{The {\em ChargedParticles} Class}
The base class {\tt IpplParticleBase} is augmented with attributes such as charge to mass ration
{\tt qm}, the vector momenta {\tt P} and the vector holding the electric field {\tt E}.
{\tt qm}, the vector momenta {\tt P} and the vector holding the electric field {\tt E}.
\begin{code}
ChargedParticles(PL* pl, Vector_t hr, Vector_t rmin,
ChargedParticles(PL* pl, Vector_t hr, Vector_t rmin,
Vector_t rmax, e_dim_tag decomp[Dim]) :
IpplParticleBase<PL>(pl),
hr_m(hr),
......@@ -173,7 +173,7 @@ int main(int argc, char *argv[]) {
decomp[d] = (d == serialDim) ? SERIAL : PARALLEL;
}
\end{code}
In the fist part of main, the discrete computational domain ({\tt domain}) and the
In the fist part of main, the discrete computational domain ({\tt domain}) and the
domain decomposition ({\tt decomp}) is constructed. We have choose a 2D domain decomposition
with $z$ serial i.e. not parallelized. \\
\begin{code}
......@@ -186,7 +186,7 @@ with $z$ serial i.e. not parallelized. \\
Vector_t rmax(nr);
partBunch=new ChargedParticles<playout_t>(PL,hr,rmin,rmax,decomp);
\end{code}
\end{code}
Here we construct the {\tt mesh} the field layout ({\tt FL}), describing how the fields are distributed
and finally the particle layout {\tt PL}. The latter is the used as a template argument to construct the
particle container. For this example the mesh size is set to unity and the computational domain is
......@@ -201,10 +201,10 @@ with $z$ serial i.e. not parallelized. \\
}
partBunch->qm = 1.0/totalP;
partBunch->myUpdate();
partBunch->myUpdate();
partBunch->initFields();
\end{code}
Now each node created {\tt nloc} particles and initialized the coordinates randomly in the
\end{code}
Now each node created {\tt nloc} particles and initialized the coordinates randomly in the
computational domain. A fixed charge to mass ration is assigned. The \texttt{myUpdate()} moves
all particles to their node defined by the domain decomposition and initialized the fields. In the last
call the fields gets initialized with the sinusoidal electric field.
......@@ -222,11 +222,11 @@ with $z$ serial i.e. not parallelized. \\
\end{code}
The last part of main consists of a simple integration scheme to advance the particles. The call
{\tt gather} interpolates the electric field at the particle position form the nearby grid points by a second
order {\it cloud in cell} (CIC) interpolation scheme.
order {\it cloud in cell} (CIC) interpolation scheme.
\subsection{{\em initFields}}
\begin{code}
void initFields() {
......@@ -237,13 +237,13 @@ void initFields() {
int nx = nr_m[0]; int ny = nr_m[1]; int nz = nr_m[2];
double phi0 = 0.1*nx;
double phi0 = 0.1*nx;
Index I(nx), J(ny), K(nz);
assign(EFD_m[I][J][K](0),
-2.0*pi*phi0/nx *
cos(2.0*pi*(I+0.5)/nx) *
assign(EFD_m[I][J][K](0),
-2.0*pi*phi0/nx *
cos(2.0*pi*(I+0.5)/nx) *
cos(4.0*pi*(J+0.5)/ny) * cos(pi*(K+0.5)/nz));
assign(EFD_m[I][J][K](1), ..... ;
......@@ -258,7 +258,7 @@ void initFields() {
\end{code}
\subsection{{\em myUpdate}}
\begin{code}
void myUpdate() {
......@@ -269,24 +269,24 @@ void myUpdate() {
EFD_m.initialize(getMesh(), getFieldLayout(), GuardCellSizes<Dim>(1), vbc_m);
EFDMag_m.initialize(getMesh(), getFieldLayout(), GuardCellSizes<Dim>(1), bc_m);
}
this->update();
this->update();
}
\end{code}
\subsection{{\em gather}}
\begin{code}
void gather() {
void gather() {
IntCIC myinterp;
E.gather(EFD_m, this->R, myinterp);
}
\end{code}
\section{Example 4 A Particle Particle - Particle Mesh (P$^3$M) Solver}
Particle Particle - Particle Mesh solvers take the close range interaction of particles
into account by combining a mesh solver as seen in the previous section with a quadratic
Particle Particle computation for particles that are closer than a given interaction
radius $r_i$ (see \url{http://en.wikipedia.org/wiki/P3M},
radius $r_i$ (see \url{http://en.wikipedia.org/wiki/P3M},
\url{http://arxiv.org/abs/astro-ph/9805096} and \url{http://arxiv.org/abs/astro-ph/0512030}).
To be able to combine these two solutions the Greens functions of the PIC solver and the
particle-particle interaction have to be modified such that they add up to the desired
......@@ -334,7 +334,7 @@ struct SpecializedGreensFunction<3> {
grn = where(lt(R*R, grn), 1./sqrt(grn),
((grn*sqrt(grn))/R-2*grn)/(R*R*R) + 2/R);
grn[0][0][0] = grn[0][0][1];
}
}
};
\end{code}
\vspace{5mm}
......@@ -350,7 +350,7 @@ function:
{
HashPairBuilder< ChargedParticles<playout_t> > HPB(*this);
//apply the field to each pair, the -1 is the field constant
HPB.for_each(RadiusCondition<double, Dim>(interaction_radius),
HPB.for_each(RadiusCondition<double, Dim>(interaction_radius),
ApplyField<double>(-1,interaction_radius));
}
\end{code}
......@@ -369,11 +369,11 @@ struct ApplyField {
double sqr = 0;
for(int d = 0;d<Dim;++d)
sqr += diff[d]*diff[d];
if(sqr!=0)
{
double r = std::sqrt(sqr);