Commit 3876559f authored by adelmann's avatar adelmann 🎗
Browse files

Merge branch 'master' of gitlab.psi.ch:OPAL/src

parents 612b636b c0ad5e53
CMAKE_MINIMUM_REQUIRED (VERSION 2.8.10)
PROJECT (OPAL)
SET (OPAL_VERSION_MAJOR 1)
SET (OPAL_VERSION_MINOR 5.00.1)
SET (OPAL_VERSION_MINOR 5.00.2)
set (PACKAGE \"opal\")
set (PACKAGE_BUGREPORT \"opal@lists.psi.ch\")
set (PACKAGE_NAME \"OPAL\")
......@@ -185,11 +185,6 @@ if (NOT ${HOSTNAME_BASE} MATCHES "edison" AND NOT ${HOSTNAME_BASE} MATCHES "cori
endif ()
endif ()
IF (AMR_DUMMY_SOLVE)
MESSAGE (STATUS "\nSolve Poisson equation with rhs equal 1")
SET (CMAKE_CXX_FLAGS "-DAMR_DUMMY_SOLVE ${CMAKE_CXX_FLAGS}")
ENDIF (AMR_DUMMY_SOLVE)
IF (DBG_SCALARFIELD)
MESSAGE (STATUS "\nWrite scalar rho_m field is enabled ")
SET (CMAKE_CXX_FLAGS "-DDBG_SCALARFIELD ${CMAKE_CXX_FLAGS}")
......
......@@ -771,7 +771,7 @@ WARN_LOGFILE =
# spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
# Note: If this tag is empty the current directory is searched.
INPUT = .
INPUT = . amr_plot_tools
# 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
......
......@@ -12,13 +12,19 @@ AmrOpal::AmrOpal(const RealBox* rb, int max_level_in, const Array<int>& n_cell_i
initBaseLevel();
nPartPerCell_m.resize(max_level_in + 1);//, PArrayManage);
// nPartPerCell_m.set(0, new MultiFab(this->boxArray(0), 1, 1));
#ifdef UNIQUE_PTR
nPartPerCell_m[0] = std::unique_ptr<MultiFab>(new MultiFab(this->boxArray(0), 1, 1, this->DistributionMap(0)));
nPartPerCell_m[0]->setVal(0.0);
chargeOnGrid_m.resize(max_level_in + 1);
chargeOnGrid_m[0] = std::unique_ptr<MultiFab>(new MultiFab(this->boxArray(0), 1, 0, this->DistributionMap(0)));
#else
nPartPerCell_m.set(0, new MultiFab(this->boxArray(0), 1, 1, this->DistributionMap(0)));
nPartPerCell_m[0].setVal(0.0);
chargeOnGrid_m.resize(max_level_in + 1);
chargeOnGrid_m.set(0, new MultiFab(this->boxArray(0), 1, 0, this->DistributionMap(0)));
#endif
// bunch_m->AssignDensitySingleLevel(0, *nPartPerCell_m[0], 0);
// bunch_m->myUpdate();
}
......@@ -42,23 +48,7 @@ void AmrOpal::initBaseLevel() {
// void AmrOpal::initFineLevels() { }
void AmrOpal::writePlotFile(std::string filename, int step) {
// Array<std::string> varnames(1, "rho");
//
// Array<const MultiFab*> tmp(finest_level + 1/*nPartPerCell_m.size()*/);
// for (/*unsigned*/ int i = 0; i < finest_level + 1/*nPartPerCell_m.size()*/; ++i) {
// tmp[i] = chargeOnGrid_m[i].get();
// }
//
// const auto& mf = tmp;
//
// Array<int> istep(finest_level+1, step);
//
// BoxLib::WriteMultiLevelPlotfile(filename, finest_level + 1, mf, varnames,
// Geom(), 0.0, istep, refRatio());
//
// return;
void AmrOpal::writePlotFileYt(std::string filename, int step) {
std::string dir = filename;
int nLevels = chargeOnGrid_m.size();
......@@ -85,7 +75,11 @@ void AmrOpal::writePlotFile(std::string filename, int step) {
HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
#ifdef UNIQUE_PTR
int nData = chargeOnGrid_m[0]->nComp();
#else
int nData = chargeOnGrid_m[0].nComp();
#endif
if (ParallelDescriptor::IOProcessor())
{
......@@ -100,7 +94,11 @@ void AmrOpal::writePlotFile(std::string filename, int step) {
HeaderFile << nData << '\n';
// variable names
#ifdef UNIQUE_PTR
for (int ivar = 1; ivar <= chargeOnGrid_m[0]->nComp(); ivar++) {
#else
for (int ivar = 1; ivar <= chargeOnGrid_m[0].nComp(); ivar++) {
#endif
HeaderFile << "rho\n";
}
......@@ -174,12 +172,21 @@ void AmrOpal::writePlotFile(std::string filename, int step) {
if (ParallelDescriptor::IOProcessor())
{
#ifdef UNIQUE_PTR
HeaderFile << lev << ' ' << chargeOnGrid_m[lev]->boxArray().size() << ' ' << 0 << '\n';
HeaderFile << 0 << '\n'; // # time steps at this level
for (int i = 0; i < chargeOnGrid_m[lev]->boxArray().size(); ++i)
{
RealBox loc = RealBox(chargeOnGrid_m[lev]->boxArray()[i],geom[lev].CellSize(),geom[lev].ProbLo());
#else
HeaderFile << lev << ' ' << chargeOnGrid_m[lev].boxArray().size() << ' ' << 0 << '\n';
HeaderFile << 0 << '\n'; // # time steps at this level
for (int i = 0; i < chargeOnGrid_m[lev].boxArray().size(); ++i)
{
RealBox loc = RealBox(chargeOnGrid_m[lev].boxArray()[i],geom[lev].CellSize(),geom[lev].ProbLo());
#endif
for (int n = 0; n < BL_SPACEDIM; n++)
HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n';
}
......@@ -189,8 +196,11 @@ void AmrOpal::writePlotFile(std::string filename, int step) {
HeaderFile << PathNameInHeader << '\n';
}
#ifdef UNIQUE_PTR
MultiFab data(chargeOnGrid_m[lev]->boxArray(), nData, 0);
#else
MultiFab data(chargeOnGrid_m[lev].boxArray(), nData, 0);
#endif
// dst, src, srccomp, dstcomp, numcomp, nghost
/*
* srccomp: the component to copy
......@@ -198,8 +208,11 @@ void AmrOpal::writePlotFile(std::string filename, int step) {
* numcomp: how many components to copy
*/
// data.copy(*nPartPerCell_m[lev],0,0,0);
#ifdef UNIQUE_PTR
MultiFab::Copy(data, *chargeOnGrid_m[lev], 0, 0, 1, 0);
#else
MultiFab::Copy(data, chargeOnGrid_m[lev], 0, 0, 1, 0);
#endif
//
// Use the Full pathname when naming the MultiFab.
//
......@@ -211,20 +224,54 @@ void AmrOpal::writePlotFile(std::string filename, int step) {
}
void AmrOpal::writePlotFile(std::string filename, int step) {
Array<std::string> varnames(1, "rho");
Array<const MultiFab*> tmp(finest_level + 1/*nPartPerCell_m.size()*/);
for (/*unsigned*/ int i = 0; i < finest_level + 1/*nPartPerCell_m.size()*/; ++i) {
#ifdef UNIQUE_PTR
tmp[i] = chargeOnGrid_m[i].get();
#else
tmp[i] = &chargeOnGrid_m[i];
#endif
}
const auto& mf = tmp;
Array<int> istep(finest_level+1, step);
BoxLib::WriteMultiLevelPlotfile(filename, finest_level + 1, mf, varnames,
Geom(), 0.0, istep, refRatio());
}
void AmrOpal::ErrorEst(int lev, TagBoxArray& tags, Real time, int /*ngrow*/) {
for (int i = lev; i <= finest_level; ++i) {
#ifdef UNIQUE_PTR
nPartPerCell_m[i]->setVal(0.0);
bunch_m->AssignDensitySingleLevel(0, *nPartPerCell_m[i], i);
#else
nPartPerCell_m[i].setVal(0.0);
bunch_m->AssignDensitySingleLevel(0, nPartPerCell_m[i], i);
#endif
}
#ifdef UNIQUE_PTR
for (int i = finest_level-1; i >= lev; --i) {
MultiFab tmp(nPartPerCell_m[i]->boxArray(), 1, 0, nPartPerCell_m[i]->DistributionMap());
tmp.setVal(0.0);
BoxLib::average_down(*nPartPerCell_m[i+1], tmp, 0, 1, refRatio(i));
MultiFab::Add(*nPartPerCell_m[i], tmp, 0, 0, 1, 0);
}
#else
for (int i = finest_level-1; i >= lev; --i) {
MultiFab tmp(nPartPerCell_m[i].boxArray(), 1, 0, nPartPerCell_m[i].DistributionMap());
tmp.setVal(0.0);
BoxLib::average_down(nPartPerCell_m[i+1], tmp, 0, 1, refRatio(i));
MultiFab::Add(nPartPerCell_m[i], tmp, 0, 0, 1, 0);
}
#endif
// std::cout << lev << " " << nPartPerCell_m[lev]->min(0) << " " << nPartPerCell_m[lev]->max(0) << std::endl;
......@@ -242,8 +289,11 @@ void AmrOpal::ErrorEst(int lev, TagBoxArray& tags, Real time, int /*ngrow*/) {
#endif
{
Array<int> itags;
#ifdef UNIQUE_PTR
for (MFIter mfi(*nPartPerCell_m[lev],false/*true*/); mfi.isValid(); ++mfi) {
#else
for (MFIter mfi(nPartPerCell_m[lev],false/*true*/); mfi.isValid(); ++mfi) {
#endif
const Box& tilebx = mfi.validbox();//mfi.tilebox();
TagBox& tagfab = tags[mfi];
......@@ -258,7 +308,11 @@ void AmrOpal::ErrorEst(int lev, TagBoxArray& tags, Real time, int /*ngrow*/) {
const int* thi = tilebx.hiVect();
state_error(tptr, ARLIM_3D(tlo), ARLIM_3D(thi),
#ifdef UNIQUE_PTR
BL_TO_FORTRAN_3D((*nPartPerCell_m[lev])[mfi]),
#else
BL_TO_FORTRAN_3D((nPartPerCell_m[lev])[mfi]),
#endif
&tagval, &clearval,
ARLIM_3D(tilebx.loVect()), ARLIM_3D(tilebx.hiVect()),
ZFILL(dx), ZFILL(prob_lo), &time, &nPart);
......@@ -346,9 +400,17 @@ AmrOpal::RemakeLevel (int lev, Real time,
SetBoxArray(lev, new_grids);
SetDistributionMap(lev, new_dmap);
#ifdef UNIQUE_PTR
nPartPerCell_m[lev].reset(new MultiFab(new_grids, 1, 1, new_dmap));
chargeOnGrid_m[lev].reset(new MultiFab(new_grids, 1, 0, new_dmap));
#else
nPartPerCell_m.clear(lev);
nPartPerCell_m.set(lev, new MultiFab(new_grids, 1, 1, new_dmap));
chargeOnGrid_m.clear(lev);
chargeOnGrid_m.set(lev, new MultiFab(new_grids, 1, 0, new_dmap));
#endif
}
void
......@@ -358,15 +420,26 @@ AmrOpal::MakeNewLevel (int lev, Real time,
SetBoxArray(lev, new_grids);
SetDistributionMap(lev, new_dmap);
#ifdef UNIQUE_PTR
nPartPerCell_m[lev] = std::unique_ptr<MultiFab>(new MultiFab(new_grids, 1, 1, dmap[lev]));
chargeOnGrid_m[lev] = std::unique_ptr<MultiFab>(new MultiFab(new_grids, 1, 0, dmap[lev]));
#else
nPartPerCell_m.set(lev, new MultiFab(new_grids, 1, 1, dmap[lev]));
chargeOnGrid_m.set(lev, new MultiFab(new_grids, 1, 0, dmap[lev]));
#endif
}
void AmrOpal::ClearLevel(int lev) {
#ifdef UNIQUE_PTR
nPartPerCell_m[lev].reset(nullptr);
chargeOnGrid_m[lev].reset(nullptr);
#else
nPartPerCell_m.clear(lev);
chargeOnGrid_m.clear(lev);
#endif
ClearBoxArray(lev);
ClearDistributionMap(lev);
}
......@@ -24,7 +24,11 @@
class AmrOpal : public AmrCore {
private:
typedef Array<std::unique_ptr<MultiFab> > mfs_mt; ///< instead of using PArray<MultiFab>
// typedef Array<std::unique_ptr<MultiFab> > mfs_mt; ///< instead of using PArray<MultiFab>
typedef PArray<MultiFab > mfs_mt;
// typedef PArray<MultiFab> mp_mt;
public:
......@@ -76,13 +80,25 @@ public:
void info() {
for (int i = 0; i < finest_level; ++i)
std::cout << "density level " << i << ": "
#ifdef UNIQUE_PTR
<< nPartPerCell_m[i]->min(0) << " "
<< nPartPerCell_m[i]->max(0) << std::endl;
#else
<< nPartPerCell_m[i].min(0) << " "
<< nPartPerCell_m[i].max(0) << std::endl;
#endif
}
/*!
* Write a timestamp file for displaying with yt.
*/
void writePlotFileYt(std::string filename, int step);
/*!
* Write a timestamp file for displaying with AmrVis.
*/
void writePlotFile(std::string filename, int step);
mfs_mt* getPartPerCell() {
......@@ -92,7 +108,11 @@ public:
void assignDensity() {
for (int i = 0; i < finest_level; ++i)
#ifdef UNIQUE_PTR
chargeOnGrid_m[i]->setVal(0.0);
#else
chargeOnGrid_m[i].setVal(0.0);
#endif
bunch_m->AssignDensity(0, false, chargeOnGrid_m, 0, 1, finest_level);
......
......@@ -105,7 +105,7 @@ public:
void destroyAll() {
for (unsigned int i = 0; i < m_particles.size(); ++i)
for (std::size_t i = 0; i < m_particles.size(); ++i)
this->RemoveParticlesAtLevel(i);
nLocalParticles_m = 0;
}
......@@ -153,6 +153,40 @@ public:
}
}
Vector_t interpolate(int i, MultiFab& quantity) {
int lev, grid, dq;
std::tie(lev, grid, dq) = idxMap_m[i];
const Real strttime = ParallelDescriptor::second();
MultiFab* ac_pointer;
if (OnSameGrids(lev,quantity)) {
ac_pointer = &quantity;
} else {
ac_pointer = new MultiFab(m_gdb->ParticleBoxArray(lev),quantity.nComp(),quantity.nGrow(),
m_gdb->ParticleDistributionMap(lev),Fab_allocate);
for (MFIter mfi(*ac_pointer); mfi.isValid(); ++mfi)
ac_pointer->setVal(0.);
ac_pointer->copy(quantity,0,0,quantity.nComp());
ac_pointer->FillBoundary(); // DO WE NEED GHOST CELLS FILLED ???
}
const FArrayBox& gfab = (*ac_pointer)[grid];
Real grav[1] = { 0.0 };
int idx[1] = { 0 };
ParticleBase::Interp(m_particles[lev][grid][dq],
m_gdb->Geom(lev),
gfab,
idx,
grav,
1);
return Vector_t(grav[0], 0, 0);
}
private:
/// Create the index mapping in order to have random access
void buildIndexMapping_m();
......
......@@ -14,7 +14,7 @@ IF (ENABLE_AMR_SOLVER)
)
set (CMAKE_CXX_FLAGS
"${IPPL_CMAKE_CXX_FLAGS} -DUSEH5FEDV2 -DPARALLEL_IO ${CMAKE_CXX_FLAGS}"
"${IPPL_CMAKE_CXX_FLAGS} -std=c++11 -DUSEH5FEDV2 -DPARALLEL_IO ${CMAKE_CXX_FLAGS}"
)
include_directories (
......@@ -25,6 +25,7 @@ IF (ENABLE_AMR_SOLVER)
${H5Hut_INCLUDE_DIR}
${HDF5_INCLUDE_DIR}
${CCSE_INCLUDE_DIR}
${HYPRE_INCLUDE_DIR}
)
link_directories (
......@@ -32,6 +33,7 @@ IF (ENABLE_AMR_SOLVER)
${CMAKE_SOURCE_DIR}/src
${IPPL_LIBRARY_DIR}
${CCSE_LIBRARY_DIR}
${HYPRE_LIBRARY_DIR}
)
# remove digits from hostname: edison03 -> edison
......@@ -62,20 +64,42 @@ IF (ENABLE_AMR_SOLVER)
set(F90_source_files Tagging_nd.f90)
add_executable (testAmrPartBunch testAmrPartBunch.cpp AmrOpal.cpp AmrPartBunch.cpp H5Reader.cpp Distribution.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
add_executable (testSolver testSolver.cpp H5Reader.cpp Distribution.cpp Solver.cpp)
# add_executable (testAmrPartBunch testAmrPartBunch.cpp AmrOpal.cpp AmrPartBunch.cpp H5Reader.cpp Distribution.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
# add_executable (testSolver testSolver.cpp H5Reader.cpp Distribution.cpp Solver.cpp)
add_executable (iterative iterative.cpp)
add_executable (testError error.cpp AmrPartBunch.cpp H5Reader.cpp Distribution.cpp PoissonProblems.cpp Solver.cpp AmrOpal.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp)
# add_executable (testError error.cpp AmrPartBunch.cpp H5Reader.cpp Distribution.cpp PoissonProblems.cpp Solver.cpp AmrOpal.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp)
add_executable (testGaussian testGaussian.cpp AmrOpal.cpp AmrPartBunch.cpp H5Reader.cpp Distribution.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
add_executable (testH5Read testH5Read.cpp H5Reader.cpp Distribution.cpp)
add_executable (testUnifSphere testUnifSphere.cpp AmrOpal.cpp AmrPartBunch.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)# HypreABecLap.H HypreABecLap.cpp)
add_executable (testGridSolve testGridSolve.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
add_executable (testUnifSphereGrid testUnifSphereGrid.cpp AmrOpal.cpp AmrPartBunch.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
add_executable(testReal testReal.cpp H5Reader.cpp Distribution.cpp AmrOpal.cpp AmrPartBunch.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
add_executable(testMultiBunch testMultiBunch.cpp H5Reader.cpp Distribution.cpp AmrOpal.cpp AmrPartBunch.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
add_executable (testH5Read testH5Read.cpp H5Reader.cpp)
add_executable(toFile toFile.cpp H5Reader.cpp Distribution.cpp)
# Boxlib has a circular dependency between -lcboxlib and -lfboxlib --> resolve by: -lcboxlib -lfboxlib -lcboxlib
target_link_libraries (testAmrPartBunch ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testSolver ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testError ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
# target_link_libraries (testAmrPartBunch ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
# target_link_libraries (testSolver ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
# target_link_libraries (testError ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testH5Read ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testGaussian ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
ENDIF(ENABLE_AMR_SOLVER)
\ No newline at end of file
target_link_libraries (testUnifSphere ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testGridSolve ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testUnifSphereGrid ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testReal ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (testMultiBunch ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
target_link_libraries (toFile ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
ENDIF(ENABLE_AMR_SOLVER)
#include "Distribution.h"
#include "H5Reader.h"
#include <fstream>
// ----------------------------------------------------------------------------
// PUBLIC MEMBER FUNCTIONS
// ----------------------------------------------------------------------------
Distribution::Distribution():
x_m(0),
y_m(0),
z_m(0),
px_m(0),
py_m(0),
pz_m(0),
nloc_m(0),
ntot_m(0)
{ }
void Distribution::uniform(double lower, double upper, size_t nloc, int seed) {
nloc_m = nloc;
......@@ -78,6 +91,8 @@ void Distribution::readH5(const std::string& filename, int step) {
nloc_m = h5.getNumParticles();
ntot_m = nloc_m;
size_t numParticlesPerNode = nloc_m / Ippl::getNodes();
size_t firstParticle = numParticlesPerNode * Ippl::myNode();
......@@ -108,19 +123,22 @@ void Distribution::readH5(const std::string& filename, int step) {
}
void Distribution::injectBeam(PartBunchBase& bunch) {
void Distribution::injectBeam(PartBunchBase& bunch, bool doDelete, std::array<double, 3> shift) {
// destroy all partcles
if ( bunch.getLocalNum() )
if ( doDelete && bunch.getLocalNum() )
bunch.destroyAll();
// previous number of particles
int prevnum = bunch.getLocalNum();
// create memory space
bunch.create(nloc_m);
for (int i = bunch.getLocalNum() - 1; i >= 0; --i) {
bunch.setR(Vector_t(x_m[i], y_m[i], z_m[i]), i);
bunch.setP(Vector_t(px_m[i], py_m[i], pz_m[i]), i);
bunch.setQM(q_m[i], i);
for (int i = nloc_m - 1; i >= 0; --i) {
bunch.setR(Vector_t(x_m[i] + shift[0], y_m[i] + shift[1], z_m[i] + shift[2]), i + prevnum);
bunch.setP(Vector_t(px_m[i], py_m[i], pz_m[i]), i + prevnum);
bunch.setQM(q_m[i], i + prevnum);
x_m.pop_back();
y_m.pop_back();
......@@ -141,3 +159,32 @@ void Distribution::setDistribution(PartBunchBase& bunch, const std::string& file
for (unsigned int i = 0; i < bunch.getLocalNum(); ++i)
bunch.setR(Vector_t(x_m[i], y_m[i], z_m[i]), i);
}
void Distribution::print2file(std::string pathname) {
for (int n = 0; n < Ippl::getNodes(); ++n) {
if ( n == Ippl::myNode() ) {
std::ofstream out;
switch (n) {
case 0:
out.open(pathname);
out << ntot_m << std::endl;
break;
default:
out.open(pathname, std::ios::app);
break;
}
for (std::size_t i = 0; i < x_m.size(); ++i)
out << x_m[i] << " " << px_m[i] << " "
<< y_m[i] << " " << py_m[i] << " "
<< z_m[i] << " " << pz_m[i] << std::endl;
out.close();
}
Ippl::Comm->barrier();
}
}
......@@ -15,6 +15,7 @@
#include <random>
#include <iostream>
#include <array>
#include "PartBunchBase.h"
......@@ -28,6 +29,8 @@ public:
public:
Distribution();
/// Generate an uniform particle distribution
/*!
* @param lower boundary
......@@ -55,8 +58,11 @@ public:
/// Transfer distribution to particle bunch object.
/*! @param bunch is either an AmrPartBunch or an PartBunch object
* @param doDelete removes all particles already in bunch before
* injection.
* @param shift all particles, each direction independently
*/
void injectBeam(PartBunchBase& bunch);
void injectBeam(PartBunchBase& bunch, bool doDelete = true, std::array<double, 3> shift = {{0.0, 0.0, 0.0}});
/// Update a distribution (only single-core)
/*! @param bunch is either an AmrPartBunch or an PartBunch object
......@@ -65,6 +71,12 @@ public:
*/
void setDistribution(PartBunchBase& bunch, const std::string& filename, int step);
/// Write the particles to a text file that can be read by OPAL. (sec. 11.3 in OPAL manual)
/*!
* @param pathname where to store.
*/
void print2file(std::string pathname);
private:
container_t x_m; ///< Horizontal particle positions [m]
container_t y_m; ///< Vertical particle positions [m]
......@@ -76,6 +88,7 @@ private:
container_t q_m; ///< Particle charge (always set to 1.0, except for Distribution::readH5)
size_t nloc_m; ///< Local number of particles
size_t ntot_m; ///< Total number of particles
};
#endif
\ No newline at end of file
#endif
#include "H5Reader.h"
#include <vector>
#include <cassert>
#include "Ippl.h"
......@@ -8,9 +9,28 @@ H5Reader::H5Reader(const std::string& filename)
: filename_m(filename), file_m(0)
{ }
H5Reader::H5Reader()
: filename_m(""), file_m(0)
{ }
void H5Reader::open(int step) {
close();
file_m = H5OpenFile(filename_m.c_str(), H5_O_RDONLY, Ippl::<