Commit 93340c63 authored by Tuelin Kaman's avatar Tuelin Kaman

MG and AMR Solver Capabilities in OPAL

parent 6e7e1ad3
......@@ -36,6 +36,8 @@ FIND_PACKAGE (GSL REQUIRED)
FIND_PACKAGE (MPI REQUIRED)
# Handle options
OPTION (ENABLE_AMR_SOLVER "Enable BoxLib based AMR solver" OFF)
OPTION (ENABLE_ML_SOLVER "Enable iteartive SA-AMG-PCG self field solver" OFF)
OPTION (DBG_SCALARFIELD "Enable dump of scalar field rho_m" OFF)
......@@ -47,6 +49,62 @@ OPTION (NOCPLUSPLUS11_NULLPTR "Disable C++11 nullptr support" OFF)
OPTION (NO_FIELD_ASSIGN_OPTIMIZATION "Disable compiler optimization of IPPL field assignment" OFF)
IF (ENABLE_AMR_SOLVER)
MESSAGE (STATUS "----> Enable AMR")
#
FIND_PACKAGE (CCSE REQUIRED HINTS $ENV{CCSE_INCLUDE_DIRS} $ENV{CCSE_LIBRARY_DIR})
#
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")
SET(BL_PRECISION "DOUBLE" CACHE INT "Precision of BoxLib build")
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 1 CACHE INT "Include backtrace information in BoxLib build")
INCLUDE($ENV{BOXLIB_HOME}/Tools/CMake/CCSEOptions.cmake)
SET (CMAKE_CXX_FLAGS "-Wno-sign-compare -DHAVE_AMR_SOLVER -DTULIN ${CMAKE_CXX_FLAGS}")
SET (CCSE_INCLUDE_DIRS $ENV{CCSE_INCLUDE_DIRS})
SET (CCSE_LIBRARY_DIR $ENV{CCSE_LIBRARY_DIR})
#
SET(CCSE_LIBRARIES box_camrdata;cboxlib;cfboxlib;fboxlib)
FOREACH (L ${CCSE_LIBRARIES})
find_library(CCSE_LIBRARY
NAMES ${L}
HINTS ${CCSE_LIBRARY_DIR}
NO_DEFAULT_PATH)
IF ( NOT CCSE_LIBRARY )
message(SEND_ERROR "Cannot locate CCSE library: ${L}")
ENDIF()
ENDFOREACH()
MESSAGE (STATUS "----> Enable AMR done")
# Get Trilinos as one entity
FIND_PACKAGE(Trilinos REQUIRED HINTS $ENV{TRILINOS_PREFIX} $ENV{TRILINOS_DIR} $ENV{TRILINOS})
# Echo trilinos build info just for fun
MESSAGE (STATUS "Found Trilinos: ${Trilinos_DIR}")
MESSAGE (STATUS " Trilinos version: ${Trilinos_VERSION}")
MESSAGE (STATUS " Trilinos package list: ${Trilinos_PACKAGE_LIST}")
MESSAGE (STATUS " Trilinos libraries: ${Trilinos_LIBRARIES}")
MESSAGE (STATUS " Trilinos TPL libraries: ${Trilinos_TPL_LIBRARIES}")
# Make sure to use same compilers and flags as Trilinos
IF (NOT ${CMAKE_CXX_COMPILER} STREQUAL ${Trilinos_CXX_COMPILER} )
MESSAGE (STATUS "Compiler mismatch:")
MESSAGE (STATUS " Trilinos was compiled with:")
MESSAGE (STATUS " ${Trilinos_C_COMPILER}")
MESSAGE (STATUS " ${Trilinos_CXX_COMPILER}")
MESSAGE (STATUS " ${Trilinos_Fortran_COMPILER}")
MESSAGE (STATUS " You are using:")
MESSAGE (STATUS " ${CMAKE_C_COMPILER}")
MESSAGE (STATUS " ${CMAKE_CXX_COMPILER}")
MESSAGE (STATUS " ${CMAKE_Fortran_COMPILER}")
MESSAGE (FATAL_ERROR "Oops ...")
ENDIF ()
ENDIF (ENABLE_AMR_SOLVER)
IF (DBG_SCALARFIELD)
MESSAGE (STATUS "\nWrite scalar rho_m field is enabled ")
SET (CMAKE_CXX_FLAGS "-DDBG_SCALARFIELD ${CMAKE_CXX_FLAGS}")
......
......@@ -82,12 +82,12 @@ void AlignWrapper::setElement(ElementBase *elem) {
}
Geometry &AlignWrapper::getGeometry() {
BGeometryBase &AlignWrapper::getGeometry() {
return itsElement->getGeometry();
}
const Geometry &AlignWrapper::getGeometry() const {
const BGeometryBase &AlignWrapper::getGeometry() const {
return itsElement->getGeometry();
}
......
......@@ -80,12 +80,12 @@ public:
/// Get geometry.
// Return the element geometry.
// Version for non-constant object.
virtual Geometry &getGeometry();
virtual BGeometryBase &getGeometry();
/// Get geometry.
// Return the element geometry.
// Version for constant object.
virtual const Geometry &getGeometry() const;
virtual const BGeometryBase &getGeometry() const;
/// Return the offset.
// This method can be used to get or set the offset. The offset is
......
......@@ -144,12 +144,12 @@ public:
/// Get geometry.
// Return the element geometry.
// Version for non-constant object.
virtual Geometry &getGeometry() = 0;
virtual BGeometryBase &getGeometry() = 0;
/// Get geometry.
// Return the element geometry
// Version for constant object.
virtual const Geometry &getGeometry() const = 0;
virtual const BGeometryBase &getGeometry() const = 0;
/// Get arc length.
// Return the entire arc length measured along the design orbit
......
......@@ -102,11 +102,11 @@ bool SBend3D::bends() const {
return true;
}
Geometry& SBend3D::getGeometry() {
BGeometryBase& SBend3D::getGeometry() {
return planarArcGeometry_m;
}
const Geometry& SBend3D::getGeometry() const {
const BGeometryBase& SBend3D::getGeometry() const {
return planarArcGeometry_m;
}
......
......@@ -118,10 +118,10 @@ class SBend3D : public Component {
void getDimensions(double &zBegin, double &zEnd) const {}
/** Return the cell geometry */
Geometry& getGeometry();
BGeometryBase& getGeometry();
/** Return the cell geometry */
const Geometry& getGeometry() const;
const BGeometryBase& getGeometry() const;
/** Return a dummy (0.) field value (what is this for?) */
EMField &getField();
......
......@@ -818,11 +818,13 @@ void PartBunch::computeSelfFields(int binNumber) {
}
void PartBunch::resizeMesh() {
double scaleFactor = 1.0;
//scaleFactor = (Physics::c * dt_m);
//get x, y range and scale to unit-less
double xmin = fs_m->solver_m->getXRangeMin() / (Physics::c * dt_m);
double xmax = fs_m->solver_m->getXRangeMax() / (Physics::c * dt_m);
double ymin = fs_m->solver_m->getYRangeMin() / (Physics::c * dt_m);
double ymax = fs_m->solver_m->getYRangeMax() / (Physics::c * dt_m);
double xmin = fs_m->solver_m->getXRangeMin() * scaleFactor ;
double xmax = fs_m->solver_m->getXRangeMax() * scaleFactor;
double ymin = fs_m->solver_m->getYRangeMin() * scaleFactor;
double ymax = fs_m->solver_m->getYRangeMax() * scaleFactor;
// Check if the new domain is larger than bunch max, mins
......@@ -848,8 +850,8 @@ void PartBunch::resizeMesh() {
get_bounds(rmin_m, rmax_m);
}
hr_m[0] = (xmax - xmin) / (nr_m[0] - 1);
hr_m[1] = (ymax - ymin) / (nr_m[1] - 1);
hr_m[0] = (xmax - xmin) / (nr_m[0]-1);
hr_m[1] = (ymax - ymin) / (nr_m[1]-1);
//hr_m[2] = (rmax_m[2] - rmin_m[2]) / (nr_m[2] - 1);
// we cannot increase the number of mesh points
......@@ -884,6 +886,7 @@ void PartBunch::computeSelfFields() {
if(fs_m->getFieldSolverType() == "MG") // || fs_m->getFieldSolverType() == "FFTBOX") {
resizeMesh();
std::cout << "after resizeMesh" << hr_m << std::endl;
//scatter charges onto grid
this->Q *= this->dt;
......@@ -906,6 +909,30 @@ void PartBunch::computeSelfFields() {
//divide charge by a 'grid-cube' volume to get [C/m^3]
rho_m *= tmp2;
#define DBG_SCALARFIELD
#ifdef DBG_SCALARFIELD
INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
ofstream fstr1;
fstr1.precision(9);
std::ostringstream istr;
istr << fieldDBGStep_m;
string SfileName = OpalData::getInstance()->getInputBasename();
string rho_fn = string("data/") + SfileName + string("-rho_scalar-") + string(istr.str());
fstr1.open(rho_fn.c_str(), ios::out);
NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
for(int x = myidx1[0].first(); x <= myidx1[0].last(); x++) {
for(int y = myidx1[1].first(); y <= myidx1[1].last(); y++) {
for(int z = myidx1[2].first(); z <= myidx1[2].last(); z++) {
fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << endl;
}
}
}
fstr1.close();
INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif
// charge density is in rho_m
fs_m->solver_m->computePotential(rho_m, hr_scaled);
......@@ -924,37 +951,25 @@ void PartBunch::computeSelfFields() {
//write out rho
// #define DBG_SCALARFIELD
#ifdef DBG_SCALARFIELD
INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
ostringstream oss;
// MPI_File file;
//MPI_Status status;
//MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL, &file);
ofstream fstr2;
fstr2.precision(9);
std::ostringstream istr;
istr << fieldDBGStep_m;
string SfileName = OpalData::getInstance()->getInputBasename();
string rho_fn = string("data/") + SfileName + string("-rho_scalar-") + string(istr.str());
fstr2.open(rho_fn.c_str(), ios::out);
string phi_fn = string("data/") + SfileName + string("-phi_scalar-") + string(istr.str());
fstr2.open(phi_fn.c_str(), ios::out);
NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
for(int x = myidx[0].first(); x <= myidx[0].last(); x++) {
for(int y = myidx[1].first(); y <= myidx[1].last(); y++) {
for(int z = myidx[2].first(); z <= myidx[2].last(); z++) {
fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << endl;
//oss << x+1 << " " << y+1 << " " << z+1 << " " << rho_m[x][y][z].get() << endl;
}
}
}
fstr2.close();
//MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
//MPI_File_close(&file);
INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif
......@@ -972,8 +987,6 @@ void PartBunch::computeSelfFields() {
//MPI_Status status;
//MPI_Info fileinfo;
//MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, fileinfo, &file);
ofstream fstr;
fstr.precision(9);
......
// ------------------------------------------------------------------------
// $RCSfile: Geometry.cpp,v $
// $RCSfile: BGeometryBase.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Geometry
// Class: BGeometryBase
// Pure virtual base class for all Beamline Geometries
//
// ------------------------------------------------------------------------
// Class category: BeamlineGeometry
// Class category: BeamlineBGeometryBase
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:34 $
......@@ -22,57 +22,57 @@
#include "BeamlineGeometry/Euclid3D.h"
// Class Geometry.
// Class BGeometryBase.
// ------------------------------------------------------------------------
Geometry::~Geometry()
BGeometryBase::~BGeometryBase()
{}
void Geometry::setElementLength(double)
void BGeometryBase::setElementLength(double)
{}
double Geometry::getOrigin() const {
double BGeometryBase::getOrigin() const {
return getArcLength() / 2.0;
}
double Geometry::getEntrance() const {
double BGeometryBase::getEntrance() const {
return - getOrigin();
}
double Geometry::getExit() const {
double BGeometryBase::getExit() const {
return getArcLength() - getOrigin();
}
Euclid3D Geometry::getTotalTransform() const {
Euclid3D BGeometryBase::getTotalTransform() const {
return getTransform(getExit(), getEntrance());
}
Euclid3D Geometry::getTransform(double s) const {
Euclid3D BGeometryBase::getTransform(double s) const {
return getTransform(0.0, s);
}
Euclid3D Geometry::getEntranceFrame() const {
Euclid3D BGeometryBase::getEntranceFrame() const {
return getTransform(0.0, getEntrance());
}
Euclid3D Geometry::getExitFrame() const {
Euclid3D BGeometryBase::getExitFrame() const {
return getTransform(0.0, getExit());
}
Euclid3D Geometry::getEntrancePatch() const {
Euclid3D BGeometryBase::getEntrancePatch() const {
return Euclid3D::identity();
}
Euclid3D Geometry::getExitPatch() const {
Euclid3D BGeometryBase::getExitPatch() const {
return Euclid3D::identity();
}
#ifndef CLASSIC_Geometry_HH
#define CLASSIC_Geometry_HH
#ifndef CLASSIC_BGeometryBase_HH
#define CLASSIC_BGeometryBase_HH
// ------------------------------------------------------------------------
// $RCSfile: Geometry.h,v $
// $RCSfile: BGeometryBase.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Geometry
// Class: BGeometryBase
//
// ------------------------------------------------------------------------
// Class category: BeamlineGeometry
// Class category: BeamlineBGeometryBase
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:34 $
......@@ -24,29 +24,29 @@
class Euclid3D;
// Class Geometry
// Class BGeometryBase
// ------------------------------------------------------------------------
/// Abstract base class for accelerator geometry classes.
// A Geometry can be considered a 3-dimensional space line parameterised by
// A BGeometryBase can be considered a 3-dimensional space line parameterised by
// the distance along the line (arc length) s. All Geometries have an exit
// and entrance plane and an origin. At any position s, a Geometry can define
// and entrance plane and an origin. At any position s, a BGeometryBase can define
// a unique 3-d rectilinear coordinate frame whose origin is on the geometry
// at s, and whose local z-axis is tangential to the geometry at s. The
// orientation of the local x- and y-axes are arbitrarilly specified by
// the Geometry. A special frame, referred to as the Geometry Local Frame
// the BGeometryBase. A special frame, referred to as the BGeometryBase Local Frame
// (or Local Frame when it is unambiguous) is specified at s = origin. The
// Local Frame is is used to define that frame about which translations and
// rotations can be applied to the Geometry. The entrance and exit planes
// rotations can be applied to the BGeometryBase. The entrance and exit planes
// are defined as those x-y planes (z=0, s=constant) in the frames defined
// at s=entrance and s=exit.
class Geometry {
class BGeometryBase {
public:
Geometry();
Geometry(const Geometry &right);
virtual ~Geometry();
const Geometry &operator=(const Geometry &right);
BGeometryBase();
BGeometryBase(const BGeometryBase &right);
virtual ~BGeometryBase();
const BGeometryBase &operator=(const BGeometryBase &right);
/// Get arc length.
// Return the length of the geometry, measured along the design arc.
......@@ -126,13 +126,13 @@ public:
};
// inlined (trivial) member functions
inline Geometry::Geometry()
inline BGeometryBase::BGeometryBase()
{ }
inline Geometry::Geometry(const Geometry &)
inline BGeometryBase::BGeometryBase(const BGeometryBase &)
{ }
inline const Geometry &Geometry::operator=(const Geometry &)
inline const BGeometryBase &BGeometryBase::operator=(const BGeometryBase &)
{ return *this; }
#endif // CLASSIC_Geometry_HH
#endif // CLASSIC_BGeometryBase_HH
......@@ -29,7 +29,7 @@ class Euclid3D;
// ------------------------------------------------------------------------
/// Geometry representing an identity transform.
class NullGeometry: public Geometry {
class NullGeometry: public BGeometryBase {
public:
NullGeometry();
......
......@@ -25,13 +25,13 @@
// Class OffsetGeometry.
// ------------------------------------------------------------------------
OffsetGeometry::OffsetGeometry(const Geometry &g, const Geometry &l,
OffsetGeometry::OffsetGeometry(const BGeometryBase &g, const BGeometryBase &l,
const Euclid3D &t):
global(g), local(l), g2l(t)
{}
OffsetGeometry::OffsetGeometry(const Geometry &g, const Euclid3D &t):
OffsetGeometry::OffsetGeometry(const BGeometryBase &g, const Euclid3D &t):
global(g), local(g), g2l(t)
{}
......
......@@ -39,19 +39,19 @@
// an alignment error, the local and global geometries can both refer to
// the same geometry.
class OffsetGeometry : public Geometry {
class OffsetGeometry : public BGeometryBase {
public:
/// Constructor.
// Assign the [b]global[/b] and [b]local[/b] geometries,
// and the displacement [b]euclid[/b] between their origins.
OffsetGeometry
(const Geometry &global, const Geometry &local, const Euclid3D &euclid);
(const BGeometryBase &global, const BGeometryBase &local, const Euclid3D &euclid);
/// Constructor.
// Both global and local geometries are set to [b]geom[/b],
// and the displacement is [b]euclid[/b].
OffsetGeometry(const Geometry &geom, const Euclid3D &euclid);
OffsetGeometry(const BGeometryBase &geom, const Euclid3D &euclid);
OffsetGeometry(const OffsetGeometry &);
virtual ~OffsetGeometry();
......@@ -138,8 +138,8 @@ public:
private:
const Geometry &global;
const Geometry &local;
const BGeometryBase &global;
const BGeometryBase &local;
Euclid3D g2l;
};
......
......@@ -53,7 +53,7 @@
// [/ul]
// [/ul]
class PlanarArcGeometry : public Geometry {
class PlanarArcGeometry : public BGeometryBase {
public:
/// Constructor.
......
......@@ -25,13 +25,13 @@
// Class SRotatedGeometry.
// ------------------------------------------------------------------------
SRotatedGeometry::SRotatedGeometry(const Geometry &g,
SRotatedGeometry::SRotatedGeometry(const BGeometryBase &g,
double sin, double sout):
srotIn(sin), srotOut(sout), geom(g)
{}
SRotatedGeometry::SRotatedGeometry(const Geometry &g,
SRotatedGeometry::SRotatedGeometry(const BGeometryBase &g,
double srot, BalanceMode mode):
srotIn(srot), srotOut(0), geom(g) {
balanceSrots(mode);
......
......@@ -38,7 +38,7 @@
// within the geometry (i.e. from s1 to s2, where s1 and/or s2 are not
// the entrance or exit planes) do not contain the s-rotations.
class SRotatedGeometry : public Geometry {
class SRotatedGeometry : public BGeometryBase {
public:
/// Balance mode.
......@@ -58,13 +58,13 @@ public:
/// Constructor.
// Use the wrapped geometry [b]geom[/b],
// and the two angles [b]srotIn[/b] and [b]srotOut[/b].
SRotatedGeometry(const Geometry &geom, double srotIn, double srotOut);
SRotatedGeometry(const BGeometryBase &geom, double srotIn, double srotOut);
/// Constructor.
// Use the wrapped geometry [b]geom[/b],
// the entrance rotation [b]srotIn[/b],
// and the balanc mode [b]mode[/b].
SRotatedGeometry(const Geometry &geom, double srotIn, BalanceMode mode = tilt);
SRotatedGeometry(const BGeometryBase &geom, double srotIn, BalanceMode mode = tilt);
SRotatedGeometry(const SRotatedGeometry &);
virtual ~SRotatedGeometry();
......@@ -162,7 +162,7 @@ private:
double srotIn;
double srotOut;
const Geometry &geom;
const BGeometryBase &geom;
};
#endif // CLASSIC_SRotatedGeometry_HH
......
......@@ -31,7 +31,7 @@
// transformations are correspondingly only simple translations along
// the z-axis.
class StraightGeometry : public Geometry {
class StraightGeometry : public BGeometryBase {
public:
/// Constructor.
......
......@@ -29,7 +29,7 @@ class Beamline;
// ------------------------------------------------------------------------
/// Implements the composite geometry of a beam line.
class BeamlineGeometry: public Geometry {
class BeamlineGeometry: public BGeometryBase {
public:
......
......@@ -54,12 +54,12 @@ MPSplitIntegrator *MPSplitIntegrator::clone() const {
}
Geometry &MPSplitIntegrator::getGeometry() {
BGeometryBase &MPSplitIntegrator::getGeometry() {
return itsElement->getGeometry();
}
const Geometry &MPSplitIntegrator::getGeometry() const {
const BGeometryBase &MPSplitIntegrator::getGeometry() const {
return itsElement->getGeometry();
}
......
......@@ -81,12 +81,12 @@ public:
/// Get geometry.
// Return the element geometry.
// Version for non-constant object.
virtual Geometry &getGeometry();
virtual BGeometryBase &getGeometry();
/// Get geometry.
// Return the element geometry
// Version for constant object.
virtual const Geometry &getGeometry() const;
virtual const BGeometryBase &getGeometry() const;
/// Get element type string.
virtual const string &getType() const;
......
......@@ -121,7 +121,6 @@ SeyNum_m(0.0),
sphys_m(NULL){
}
ParallelTTracker::ParallelTTracker(const Beamline &beamline,
PartBunch &bunch,
DataSink &ds,
......@@ -182,6 +181,66 @@ SeyNum_m(0.0) {
}
#ifdef HAVE_AMR_SOLVER
ParallelTTracker::ParallelTTracker(const Beamline &beamline,
PartBunch &bunch,
DataSink &ds,
const PartData &reference,
bool revBeam,
bool revTrack,
int maxSTEPS,
double zstop,
int timeIntegrator,
size_t N,
Amr* amrptr_in):
Tracker(beamline, reference, revBeam, revTrack),
itsBunch(&bunch),
specifiedNPart_m(N),
itsDataSink_m(&ds),
bgf_m(NULL),
itsOpalBeamline_m(),
lineDensity_m(),
RefPartR_zxy_m(0.0),
RefPartP_zxy_m(0.0),
RefPartR_suv_m(0.0),
RefPartP_suv_m(0.0),
globalEOL_m(false),
wakeStatus_m(false),
surfaceStatus_m(false),
secondaryFlg_m(false),
mpacflg_m(true),
nEmissionMode_m(false),
zStop_m(zstop),
scaleFactor_m(itsBunch->getdT() * Physics::c),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
dtTrack_m(0.0),
surfaceEmissionStop_m(-1),
minStepforReBin_m(-1),
minBinEmitted_m(std::numeric_limits<size_t>::max()),
repartFreq_m(-1),
lastVisited_m(-1),
numRefs_m(-1),
gunSubTimeSteps_m(-1),
emissionSteps_m(numeric_limits<unsigned int>::max()),
localTrackSteps_m(maxSTEPS),
maxNparts_m(0),
numberOfFieldEmittedParticles_m(numeric_limits<size_t>::max()),
bends_m(0),
numParticlesInSimulation_m(0),
space_orientation_m(0.0),