diff --git a/src/Amr/AbstractAmrWriter.h b/src/Amr/AbstractAmrWriter.h
index ae2e3a38e8b1d392adeefbd71ad6309958892885..1462a91be3b3da40cf63c1d75f8bdae9dffb4065 100644
--- a/src/Amr/AbstractAmrWriter.h
+++ b/src/Amr/AbstractAmrWriter.h
@@ -26,7 +26,7 @@
#include "Algorithms/AmrPartBunch.h"
class AbstractAmrWriter {
-
+
public:
/*!
* @param rho is the charge density on all levels
@@ -46,7 +46,7 @@ public:
const int& nLevel,
const double& time,
const double& scale = 1.0) = 0;
-
+
/*!
* @param bunch_p
* @param time
@@ -55,9 +55,9 @@ public:
virtual void writeBunch(const AmrPartBunch* bunch_p,
const double& time,
const double& scale = 1.0) = 0;
-
+
virtual ~AbstractAmrWriter() { }
-
+
};
#endif
diff --git a/src/Amr/AmrBoxLib.cpp b/src/Amr/AmrBoxLib.cpp
index 18e678e7dd2aa8d52b2734ac2952c7ec9e589230..597b70dab5015c8b95e2be8ce79bd93dc43fbfa5 100644
--- a/src/Amr/AmrBoxLib.cpp
+++ b/src/Amr/AmrBoxLib.cpp
@@ -22,23 +22,24 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-#include "AmrBoxLib.h"
+#include "Amr/AmrBoxLib.h"
#include "Algorithms/AmrPartBunch.h"
+#include "Amr/AmrYtWriter.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"
-#include "Structure/FieldSolver.h"
#include "Solvers/PoissonSolver.h"
+#include "Structure/FieldSolver.h"
+#include "Utilities/OpalException.h"
+#include "Utilities/Options.h"
#include "Utility/PAssert.h"
-#include "Amr/AmrYtWriter.h"
-
+#include <AMReX_BCUtil.H>
#include <AMReX_MultiFabUtil.H>
-
#include <AMReX_ParmParse.H> // used in initialize function
-#include <AMReX_BCUtil.H>
-#include <map>
+#include <algorithm>
+#include <cmath>
extern Inform* gmsg;
@@ -71,8 +72,7 @@ AmrBoxLib::AmrBoxLib(const AmrDomain_t& domain,
std::unique_ptr<AmrBoxLib> AmrBoxLib::create(const AmrInfo& info,
- AmrPartBunch* bunch_p)
-{
+ AmrPartBunch* bunch_p) {
/* The bunch is initialized first with a Geometry,
* BoxArray and DistributionMapping on
* the base level (level = 0). Thus, we take the domain specified there in
@@ -135,8 +135,8 @@ void AmrBoxLib::regrid(double time) {
void AmrBoxLib::getGridStatistics(std::map<int, long>& gridPtsPerCore,
- std::vector<int>& gridsPerLevel) const
-{
+ std::vector<int>& gridsPerLevel) const {
+
typedef std::vector<int> container_t;
gridPtsPerCore.clear();
@@ -590,8 +590,8 @@ void AmrBoxLib::redistributeGrids(int /*how*/) {
void AmrBoxLib::RemakeLevel (int lev, AmrReal_t /*time*/,
const AmrGrid_t& new_grids,
- const AmrProcMap_t& new_dmap)
-{
+ const AmrProcMap_t& new_dmap) {
+
SetBoxArray(lev, new_grids);
SetDistributionMap(lev, new_dmap);
@@ -623,8 +623,7 @@ void AmrBoxLib::RemakeLevel (int lev, AmrReal_t /*time*/,
void AmrBoxLib::MakeNewLevel (int lev, AmrReal_t /*time*/,
const AmrGrid_t& new_grids,
- const AmrProcMap_t& new_dmap)
-{
+ const AmrProcMap_t& new_dmap) {
SetBoxArray(lev, new_grids);
SetDistributionMap(lev, new_dmap);
@@ -643,8 +642,6 @@ void AmrBoxLib::MakeNewLevel (int lev, AmrReal_t /*time*/,
efield_m[lev][j]->setVal(0.0, 1);
}
-
-
/*
* particles need to know the BoxArray
* and DistributionMapping
@@ -675,22 +672,22 @@ void AmrBoxLib::ErrorEst(int lev, TagBoxArray_t& tags,
*gmsg << level2 << "* Start tagging of level " << lev << endl;
switch ( tagging_m ) {
- case CHARGE_DENSITY:
+ case TaggingCriteria::CHARGE_DENSITY:
tagForChargeDensity_m(lev, tags, time, ngrow);
break;
- case POTENTIAL:
+ case TaggingCriteria::POTENTIAL:
tagForPotentialStrength_m(lev, tags, time, ngrow);
break;
- case EFIELD:
+ case TaggingCriteria::EFIELD:
tagForEfield_m(lev, tags, time, ngrow);
break;
- case MOMENTA:
+ case TaggingCriteria::MOMENTA:
tagForMomenta_m(lev, tags, time, ngrow);
break;
- case MAX_NUM_PARTICLES:
+ case TaggingCriteria::MAX_NUM_PARTICLES:
tagForMaxNumParticles_m(lev, tags, time, ngrow);
break;
- case MIN_NUM_PARTICLES:
+ case TaggingCriteria::MIN_NUM_PARTICLES:
tagForMinNumParticles_m(lev, tags, time, ngrow);
break;
default:
@@ -764,7 +761,7 @@ void AmrBoxLib::preRegrid_m() {
* So, we need to solve the Poisson problem first assuming
* a single bin only
*/
- if ( tagging_m == POTENTIAL || tagging_m == EFIELD ) {
+ if ( tagging_m == TaggingCriteria::POTENTIAL || tagging_m == TaggingCriteria::EFIELD ) {
this->solvePoisson_m();
}
}
@@ -1178,11 +1175,11 @@ void AmrBoxLib::tagForMinNumParticles_m(int lev, TagBoxArray_t& tags,
for (int k = tlo[2]; k <= thi[2]; ++k) {
AmrIntVect_t iv(i, j, k);
if ( cells.find(iv) != cells.end() &&
- cells[iv] >= minNumPart_m )
- {
+ cells[iv] >= minNumPart_m ) {
tagfab(iv) = tagval;
- } else
+ } else {
tagfab(iv) = clearval;
+ }
}
}
}
@@ -1529,21 +1526,18 @@ void AmrBoxLib::fillPhysbc_m(AmrField_t& mf, int lev) {
* Author: Weiqun Zhang <weiqunzhang@lbl.gov>
* Date: Mon Jul 2 10:40:21 2018 -0700
*/
- if (AmrGeometry_t::isAllPeriodic())
+ if (AmrGeometry_t::isAllPeriodic()) {
return;
+ }
+
// Set up BC; see Src/Base/AMReX_BC_TYPES.H for supported types
amrex::Vector<amrex::BCRec> bc(mf.nComp());
- for (int n = 0; n < mf.nComp(); ++n)
- {
- for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
- {
- if (AmrGeometry_t::isPeriodic(idim))
- {
+ for (int n = 0; n < mf.nComp(); ++n) {
+ for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+ if (AmrGeometry_t::isPeriodic(idim)) {
bc[n].setLo(idim, amrex::BCType::int_dir); // interior
bc[n].setHi(idim, amrex::BCType::int_dir);
- }
- else
- {
+ } else {
bc[n].setLo(idim, amrex::BCType::foextrap); // first-order extrapolation.
bc[n].setHi(idim, amrex::BCType::foextrap);
}
diff --git a/src/Amr/AmrBoxLib.h b/src/Amr/AmrBoxLib.h
index 8742ad72a2d707b48b67ad33927d7ff0251de255..699488ed02d67ecc3229dac1f8ca0e03a23900fb 100644
--- a/src/Amr/AmrBoxLib.h
+++ b/src/Amr/AmrBoxLib.h
@@ -27,16 +27,18 @@
#include "Amr/AmrObject.h"
-class AmrPartBunch;
-
// AMReX headers
#include <AMReX_AmrMesh.H>
#include <AMReX.H>
-class AmrBoxLib : public AmrObject,
- public amrex::AmrMesh
-{
-
+#include <map>
+#include <vector>
+
+class AmrPartBunch;
+
+class AmrBoxLib: public AmrObject,
+ public amrex::AmrMesh {
+
public:
typedef amr::AmrField_t AmrField_t;
typedef amr::AmrScalarFieldContainer_t AmrScalarFieldContainer_t;
@@ -51,14 +53,13 @@ public:
typedef amr::AmrProcMap_t AmrProcMap_t;
typedef amr::AmrGeometry_t AmrGeometry_t;
typedef amr::AmrIntVect_t AmrIntVect_t;
-
+
typedef amrex::FArrayBox FArrayBox_t;
typedef amrex::Box Box_t;
typedef amrex::TagBox TagBox_t;
typedef amrex::TagBoxArray TagBoxArray_t;
typedef amrex::MFIter MFIter_t;
-
-
+
/*!
* Used for the redistribution of grids
*/
@@ -68,9 +69,8 @@ public:
RANDOM = 2,
KNAPSACK = 3
};
-
+
public:
-
/*!
* See other constructors documentation for further info.
* @param domain is the physical domain of the problem
@@ -82,7 +82,7 @@ public:
const AmrIntArray_t& nGridPts,
int maxLevel,
AmrPartBunch* bunch_p);
-
+
/*!
* Create a new object
* @param info are the initial informations to construct an AmrBoxLib object.
@@ -92,7 +92,7 @@ public:
*/
static std::unique_ptr<AmrBoxLib> create(const AmrInfo& info,
AmrPartBunch* bunch_p);
-
+
/*!
* Inherited from AmrObject
*/
@@ -100,50 +100,50 @@ public:
void getGridStatistics(std::map<int, long>& gridPtsPerCore,
std::vector<int>& gridsPerLevel) const;
-
+
/*!
* Initial gridding. Sets up all levels.
*/
void initFineLevels();
VectorPair_t getEExtrema();
-
+
double getRho(int x, int y, int z);
-
+
void computeSelfFields();
-
+
void computeSelfFields(int bin);
-
+
void computeSelfFields_cycl(double gamma);
-
+
void computeSelfFields_cycl(int bin);
-
+
void updateMesh();
-
+
/*!
* Mesh scaling for solver (gamma factor)
* (in particle rest frame, the longitudinal length enlarged)
*/
const Vector_t& getMeshScaling() const;
-
+
Vektor<int, 3> getBaseLevelGridPoints() const;
-
+
const int& maxLevel() const;
const int& finestLevel() const;
-
+
/*!
* @returns the time of the bunch [s]
*/
double getT() const;
-
+
void redistributeGrids(int how);
-
-
+
+
protected:
/*
* AmrMesh functions
*/
-
+
/*!
* Update the grids and the distributionmapping for a specific level
* (inherited from AmrMesh)
@@ -154,7 +154,7 @@ protected:
*/
void RemakeLevel (int lev, AmrReal_t time,
const AmrGrid_t& new_grids, const AmrProcMap_t& new_dmap);
-
+
/*!
* Create completeley new grids for a level (inherited from AmrMesh)
* @param lev is the current level
@@ -164,20 +164,19 @@ protected:
*/
void MakeNewLevel (int lev, AmrReal_t time,
const AmrGrid_t& new_grids, const AmrProcMap_t& new_dmap);
-
+
/*!
* Clean up a level.
* @param lev to free allocated memory.
*/
void ClearLevel(int lev);
-
+
/*!
* Is called in the AmrMesh function for performing tagging. (inherited from AmrMesh)
*/
virtual void ErrorEst(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow) override;
-
-
+
/*!
* Make a new level from scratch using provided BoxArray and
* DistributionMapping.
@@ -190,9 +189,9 @@ protected:
* @param ba the boxes
* @param dm the grid distribution among cores
*/
- void MakeNewLevelFromScratch (int lev, AmrReal_t time,
- const AmrGrid_t& ba,
- const AmrProcMap_t& dm);
+ void MakeNewLevelFromScratch(int lev, AmrReal_t time,
+ const AmrGrid_t& ba,
+ const AmrProcMap_t& dm);
/*!
* Make a new level using provided BoxArray and
@@ -205,9 +204,10 @@ protected:
* @param ba the boxes
* @param dm the grid distribution among cores
*/
- void MakeNewLevelFromCoarse (int lev, AmrReal_t time,
- const AmrGrid_t& ba,
- const AmrProcMap_t& dm);
+ void MakeNewLevelFromCoarse(int lev, AmrReal_t time,
+ const AmrGrid_t& ba,
+ const AmrProcMap_t& dm);
+
private:
/*!
@@ -217,7 +217,6 @@ private:
* @param time of regrid
*/
void doRegrid_m(int lbase, double time);
-
/*!
* Called within doRegrid_m(). Especially used
@@ -233,15 +232,13 @@ private:
*/
void postRegrid_m(int old_finest);
-
double solvePoisson_m();
-
/* ATTENTION
* The tagging routines assume the particles to be in the
* AMR domain, i.e. [-1, 1]^3
*/
-
+
/*!
* Mark a cell for refinement if the value is greater equal
* than some amount of charge (AmrObject::chargedensity_m).
@@ -253,7 +250,7 @@ private:
*/
void tagForChargeDensity_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow);
-
+
/*!
* Mark a cell for refinement if the potential value is greater
* equal than the maximum value of the potential on the grid
@@ -266,7 +263,7 @@ private:
*/
void tagForPotentialStrength_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow);
-
+
/*!
* Mark a cell for refinement if one of the electric field components
* is greater equal the maximum electric field value per direction scaled
@@ -279,7 +276,7 @@ private:
*/
void tagForEfield_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow);
-
+
/*!
* Mark a cell for refinement if at least one particle has
* a high momentum. The lower bound is specified by the
@@ -293,7 +290,7 @@ private:
*/
void tagForMomenta_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow);
-
+
/*!
* Mark a cell for refinement if it contains at most
* AmrObject::maxNumPart_m particles.
@@ -305,7 +302,7 @@ private:
*/
void tagForMaxNumParticles_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow);
-
+
/*!
* Mark a cell for refinement if it contains at least
* AmrObject::minNumPart_m particles.
@@ -317,7 +314,7 @@ private:
*/
void tagForMinNumParticles_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow);
-
+
/*!
* Use particle BoxArray and DistributionMapping for AmrObject and
* reset geometry for bunch
@@ -325,7 +322,7 @@ private:
* @param nGridPts per dimension (nx, ny, nz / nt)
*/
void initBaseLevel_m(const AmrIntArray_t& nGridPts);
-
+
/*!
* AMReX uses the ParmParse object to initialize
* parameters like the maximum level etc.
@@ -336,15 +333,14 @@ private:
* @param layout_p of bunch
*/
static void initParmParse_m(const AmrInfo& info, AmrLayout_t* layout_p);
-
+
/*!
* Fill the physical / mesh boundary values
*
* @param mf field to fill boundary values
*/
void fillPhysbc_m(AmrField_t& mf, int lev = 0);
-
-
+
// void gradient(int lev) {
//
// phi_m[lev]->FillBoundary(geom[lev].periodicity());
@@ -382,23 +378,24 @@ private:
// }
// }
// }
-
+
+
private:
/// bunch used for tagging strategies
- AmrPartBunch *bunch_mp;
-
+ AmrPartBunch* bunch_mp;
+
// the layout of the bunch
- AmrLayout_t *layout_mp;
-
+ AmrLayout_t* layout_mp;
+
/// charge density on the grid for all levels
AmrScalarFieldContainer_t rho_m;
-
+
/// scalar potential on the grid for all levels
AmrScalarFieldContainer_t phi_m;
-
+
/// vector field on the grid for all levels
AmrVectorFieldContainer_t efield_m;
-
+
/// in particle rest frame, the longitudinal length enlarged
Vector_t meshScaling_m;
diff --git a/src/Amr/AmrObject.cpp b/src/Amr/AmrObject.cpp
index d6e15be3ebc5065e21be742d36658c3d8fb18570..a00e5fb911cc09e7688724f7868332e27517b103 100644
--- a/src/Amr/AmrObject.cpp
+++ b/src/Amr/AmrObject.cpp
@@ -21,10 +21,12 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-#include "AmrObject.h"
+#include "Amr/AmrObject.h"
+
+#include "Utilities/OpalException.h"
AmrObject::AmrObject()
- : AmrObject(CHARGE_DENSITY, 0.75, 1.0e-15)
+ : AmrObject(TaggingCriteria::CHARGE_DENSITY, 0.75, 1.0e-15)
{ }
@@ -52,23 +54,24 @@ void AmrObject::setTagging(TaggingCriteria tagging) {
void AmrObject::setTagging(const std::string& tagging) {
- if ( tagging == "POTENTIAL" )
- tagging_m = TaggingCriteria::POTENTIAL;
- else if (tagging == "EFIELD" )
- tagging_m = TaggingCriteria::EFIELD;
- else if ( tagging == "MOMENTA" )
- tagging_m = TaggingCriteria::MOMENTA;
- else if ( tagging == "MAX_NUM_PARTICLES" )
- tagging_m = TaggingCriteria::MAX_NUM_PARTICLES;
- else if ( tagging == "MIN_NUM_PARTICLES" )
- tagging_m = TaggingCriteria::MIN_NUM_PARTICLES;
- else if ( tagging == "CHARGE_DENSITY" )
- tagging_m = TaggingCriteria::CHARGE_DENSITY;
- else
- throw OpalException("AmrObject::setTagging(std::string)",
- "Not supported refinement criteria "
+ static const std::map<std::string, TaggingCriteria> stringTaggingCriteria_s = {
+ {"CHARGE_DENSITY", TaggingCriteria::CHARGE_DENSITY},
+ {"POTENTIAL", TaggingCriteria::POTENTIAL},
+ {"EFIELD", TaggingCriteria::EFIELD},
+ {"MOMENTA", TaggingCriteria::MOMENTA},
+ {"MIN_NUM_PARTICLES", TaggingCriteria::MIN_NUM_PARTICLES},
+ {"MAX_NUM_PARTICLES", TaggingCriteria::MAX_NUM_PARTICLES }
+ };
+
+ if (stringTaggingCriteria_s.count(tagging)) {
+ tagging_m = stringTaggingCriteria_s.at(tagging);
+ } else {
+ throw OpalException("AmrObject::setTagging",
+ "Not supported refinement criteria.\n"
+ "Check the accepted values: "
"[CHARGE_DENSITY | POTENTIAL | EFIELD | "
- "MOMENTA | MAX_NUM_PARTICLES | MIN_NUM_PARTICLES].");
+ "MOMENTA | MIN_NUM_PARTICLES | MAX_NUM_PARTICLES].");
+ }
}
@@ -97,29 +100,29 @@ const bool& AmrObject::isRefined() const {
}
-std::string AmrObject::enum2string(int number) {
+std::string AmrObject::getTaggingString(int number) {
std::string tagging;
switch ( number ) {
- case TaggingCriteria::CHARGE_DENSITY:
+ case static_cast<std::underlying_type_t<TaggingCriteria>>(TaggingCriteria::CHARGE_DENSITY):
tagging = "CHARGE_DENSITY";
break;
- case TaggingCriteria::POTENTIAL:
+ case static_cast<std::underlying_type_t<TaggingCriteria>>(TaggingCriteria::POTENTIAL):
tagging = "POTENTIAL";
break;
- case TaggingCriteria::EFIELD:
+ case static_cast<std::underlying_type_t<TaggingCriteria>>(TaggingCriteria::EFIELD):
tagging = "EFIELD";
break;
- case TaggingCriteria::MOMENTA:
+ case static_cast<std::underlying_type_t<TaggingCriteria>>(TaggingCriteria::MOMENTA):
tagging = "MOMENTA";
break;
- case TaggingCriteria::MIN_NUM_PARTICLES:
+ case static_cast<std::underlying_type_t<TaggingCriteria>>(TaggingCriteria::MIN_NUM_PARTICLES):
tagging = "MIN_NUM_PARTICLES";
break;
- case TaggingCriteria::MAX_NUM_PARTICLES:
+ case static_cast<std::underlying_type_t<TaggingCriteria>>(TaggingCriteria::MAX_NUM_PARTICLES):
tagging = "MAX_NUM_PARTICLES";
break;
default:
- throw OpalException("AmrObject::enum2string(int)",
+ throw OpalException("AmrObject::getTaggingString",
"Only numbers between 0 and 5 allowed.");
}
return tagging;
diff --git a/src/Amr/AmrObject.h b/src/Amr/AmrObject.h
index 1fa29341349d724a1074d96a530d390760fb94fa..21f44bf193a7faec0ef3a8624b136c985d2e6d1a 100644
--- a/src/Amr/AmrObject.h
+++ b/src/Amr/AmrObject.h
@@ -29,23 +29,22 @@
#include "Algorithms/PBunchDefs.h"
class AmrObject {
-
+
public:
// FIXME Why not using typedef of PartBunchBase
typedef std::pair<Vector_t, Vector_t> VectorPair_t;
// using VectorPair_t = typename AmrPartBunch::VectorPair_t;
-
+
/// Methods for tagging cells for refinement
- enum TaggingCriteria {
- CHARGE_DENSITY = 0, // default
+ enum class TaggingCriteria: unsigned int {
+ CHARGE_DENSITY = 0,
POTENTIAL,
EFIELD,
MOMENTA,
- MIN_NUM_PARTICLES, ///< min. #particles per cell
- MAX_NUM_PARTICLES ///< max. #particles per cell
+ MIN_NUM_PARTICLES, ///< min. #particles per cell
+ MAX_NUM_PARTICLES ///< max. #particles per cell
};
-
-
+
/*!
* This data structure is only used for creating an object
* via the static member function AmrBoxLib::create()
@@ -58,17 +57,16 @@ public:
int maxlevel; ///< Maximum level for AMR (0: single-level)
int refratio[3]; ///< Mesh refinement ratio in x-, y- and z-direction
};
-
+
public:
-
AmrObject();
-
+
AmrObject(TaggingCriteria tagging,
double scaling,
double chargedensity);
-
+
virtual ~AmrObject();
-
+
/*!
* Collect information about grid load balancing.
* @param gridPtsPerCore is filled.
@@ -76,121 +74,120 @@ public:
*/
virtual void getGridStatistics(std::map<int, long>& gridPtsPerCore,
std::vector<int>& gridsPerLevel) const = 0;
-
+
/*!
* Setup all fine levels after object creation.
*/
virtual void initFineLevels() = 0;
-
+
/*!
* Update of mesh according to chosen refinement strategy.
* @param time of regrid
*/
virtual void regrid(double time) = 0;
-
+
/*!
* Choose a new tagging strategy.
* Is used in src/Structure/FieldSolver.cpp
* @param tagging strategy
*/
void setTagging(TaggingCriteria tagging);
-
+
/*!
* Choose a new tagging strategy (string version).
* Is used in src/Structure/FieldSolver.cpp
* @param tagging strategy
*/
void setTagging(const std::string& tagging);
-
+
/*!
* Scaling factor for tagging.
* It is used with POTENTIAL and EFIELD
* @param scaling factor in [0, 1]
*/
void setScalingFactor(double scaling);
-
+
/*!
* Charge density for tagging with CHARGE_DENSITY
* @param chargedensity >= 0.0 (e.g. 1e-14)
*/
void setChargeDensity(double chargedensity);
-
+
/*!
* Maximum number of particles per cell for tagging
* @param maxNumPart is upper bound for a cell to be marked
* for refinement
*/
void setMaxNumParticles(size_t maxNumPart);
-
+
/*!
* Minimum number of particles per cell for tagging
* @param minNumPart is lower bound for a cell to be marked
* for refinement
*/
void setMinNumParticles(size_t minNumPart);
-
+
/* Methods that are needed by the
* bunch
*/
virtual VectorPair_t getEExtrema() = 0;
-
+
virtual double getRho(int x, int y, int z) = 0;
-
+
virtual void computeSelfFields() = 0;
-
+
virtual void computeSelfFields(int b) = 0;
-
+
virtual void computeSelfFields_cycl(double gamma) = 0;
-
+
virtual void computeSelfFields_cycl(int b) = 0;
-
+
virtual void updateMesh() = 0;
-
+
virtual Vektor<int, 3> getBaseLevelGridPoints() const = 0;
-
+
virtual const int& maxLevel() const = 0;
virtual const int& finestLevel() const = 0;
-
+
/*!
* @returns the time of the simulation
*/
virtual double getT() const = 0;
-
+
/*!
* Rebalance the grids among the
* cores
*/
virtual void redistributeGrids(int /*how*/) { }
-
+
/*!
* Used in AmrPartBunch to check if we need to refine
* first.
* @returns true fine grids are initialized
*/
const bool& isRefined() const;
-
+
/*!
* Used in Fieldsolver in order to convert a number that
* specifies the tagging to the corresponding string. Check
* enum TaggingCriteria for ordering.
* @param number of tagging
*/
- static std::string enum2string(int number);
+ static std::string getTaggingString(int number);
protected:
-
TaggingCriteria tagging_m; ///< Tagging strategy
-
+
double scaling_m; ///< Scaling factor for tagging [0, 1]
// (POTENTIAL, EFIELD)
double chargedensity_m; ///< Tagging value for CHARGE_DENSITY
-
+
size_t maxNumPart_m; ///< Tagging value for MAX_NUM_PARTICLES
-
+
size_t minNumPart_m; ///< Tagging value for MIN_NUM_PARTICLES
-
+
bool refined_m; ///< Only set to true in AmrObject::initFineLevels()
-
+
/// timer for selfField calculation (used in concrete AmrObject classes)
IpplTimings::TimerRef amrSolveTimer_m;
IpplTimings::TimerRef amrRegridTimer_m;
diff --git a/src/Amr/AmrYtWriter.cpp b/src/Amr/AmrYtWriter.cpp
index a4e7846a915dbe3fdcb2077269d96c7ec7a876de..0a5bed7de31d85bfaca15d3aabf71e38b37764b4 100644
--- a/src/Amr/AmrYtWriter.cpp
+++ b/src/Amr/AmrYtWriter.cpp
@@ -22,7 +22,7 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-#include "AmrYtWriter.h"
+#include "Amr/AmrYtWriter.h"
#include <AMReX_Utility.H>
#include <AMReX_IntVect.H>
@@ -36,20 +36,19 @@
#include "Utilities/OpalException.h"
-AmrYtWriter::AmrYtWriter(int step, int bin)
-{
+AmrYtWriter::AmrYtWriter(int step, int bin) {
intData_m.resize(1);
// intData_m[0] = "id";
// intData_m[1] = "cpu";
intData_m[0] = "energy_bin";
-
+
realData_m.resize(//3 + // coordinates
3 + // momenta
3 + // charge + mass + timestep
// 1 + // potential at particle location
3 + // electric field at particle location
3); // magnetic field at particle location
-
+
int idx = 0;
// realData_m[idx++] ="position_x";
// realData_m[idx++] ="position_y";
@@ -67,16 +66,16 @@ AmrYtWriter::AmrYtWriter(int step, int bin)
realData_m[idx++] ="magnetic_field_x";
realData_m[idx++] ="magnetic_field_y";
realData_m[idx++] ="magnetic_field_z";
-
+
namespace fs = boost::filesystem;
-
+
fs::path dir = OpalData::getInstance()->getInputBasename();
std::string dataDir = OpalData::getInstance()->getAuxiliaryOutputDirectory();
boost::filesystem::path path = dir.parent_path() / dataDir / "amr" / "yt";
dir_m = amrex::Concatenate((path / "plt").string(), step, 10);
dir_m += "-";
dir_m = amrex::Concatenate(dir_m, bin, 3);
-
+
if ( Ippl::myNode() == 0 && !fs::exists(path) ) {
try {
fs::create_directories( path );
@@ -85,7 +84,7 @@ AmrYtWriter::AmrYtWriter(int step, int bin)
ex.what());
}
}
-
+
Ippl::Comm->barrier();
}
@@ -97,12 +96,11 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
const amr::AmrGeomContainer_t& geom,
const int& nLevel,
const double& time,
- const double& scale)
-{
+ const double& scale) {
/* We need to scale the geometry and cell sizes according to the
* particle mapping
*/
-
+
//
// Only let 64 CPUs be writing at any one time.
//
@@ -131,9 +129,8 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
+ efield[0][0]->nComp()
+ efield[0][1]->nComp()
+ efield[0][2]->nComp();
-
- if ( Ippl::myNode() == 0 )
- {
+
+ if ( Ippl::myNode() == 0 ) {
//
// Only the IOProcessor() writes to the header file.
//
@@ -147,19 +144,19 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
// variable names
for (int ivar = 1; ivar <= rho[0]->nComp(); ivar++)
HeaderFile << "rho\n";
-
+
for (int ivar = 1; ivar <= phi[0]->nComp(); ivar++)
HeaderFile << "phi\n";
-
+
HeaderFile << "Ex\nEy\nEz\n";
-
+
// dimensionality
HeaderFile << AMREX_SPACEDIM << '\n';
-
+
// time
HeaderFile << time << '\n';
HeaderFile << nLevel - 1 << '\n'; // maximum level number (0=single level)
-
+
// physical domain
for (int i = 0; i < AMREX_SPACEDIM; i++)
HeaderFile << geom[0].ProbLo(i) / scale << ' ';
@@ -167,34 +164,34 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
for (int i = 0; i < AMREX_SPACEDIM; i++)
HeaderFile << geom[0].ProbHi(i) / scale << ' ';
HeaderFile << '\n';
-
+
// reference ratio
for (int i = 0; i < refRatio.size(); ++i)
HeaderFile << refRatio[i] << ' ';
HeaderFile << '\n';
-
+
// geometry domain for all levels
for (int i = 0; i < nLevel; ++i)
HeaderFile << geom[i].Domain() << ' ';
HeaderFile << '\n';
-
+
// number of time steps
for (int i = 0; i < nLevel; ++i)
HeaderFile << 0 << ' ';
HeaderFile << '\n';
-
+
// cell sizes for all level
for (int i = 0; i < nLevel; ++i) {
for (int k = 0; k < AMREX_SPACEDIM; k++)
HeaderFile << geom[i].CellSize()[k] / scale << ' ';
HeaderFile << '\n';
}
-
+
// coordinate system
HeaderFile << geom[0].Coord() << '\n';
HeaderFile << "0\n"; // write boundary data
}
-
+
for (int lev = 0; lev < nLevel; ++lev) {
// Build the directory to hold the MultiFab at this level.
// The name is relative to the directory containing the Header file.
@@ -203,32 +200,33 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
char buf[64];
snprintf(buf, sizeof(buf), "Level_%d", lev);
std::string sLevel = buf;
-
+
//
// Now for the full pathname of that directory.
//
std::string FullPath = dir_m;
- if (!FullPath.empty() && FullPath[FullPath.length()-1] != '/')
+ if (!FullPath.empty() && FullPath[FullPath.length()-1] != '/') {
FullPath += '/';
+ }
FullPath += sLevel;
//
// Only the I/O processor makes the directory if it doesn't already exist.
//
- if ( Ippl::myNode() == 0 )
- if (!amrex::UtilCreateDirectory(FullPath, 0755))
+ if ( Ippl::myNode() == 0 ) {
+ if (!amrex::UtilCreateDirectory(FullPath, 0755)) {
amrex::CreateDirectoryFailed(FullPath);
+ }
+ }
//
// Force other processors to wait till directory is built.
//
Ippl::Comm->barrier();
-
- if ( Ippl::myNode() == 0 )
- {
+
+ if ( Ippl::myNode() == 0 ) {
HeaderFile << lev << ' ' << rho[lev]->boxArray().size() << ' ' << 0 /*time*/ << '\n';
HeaderFile << 0 /* level steps */ << '\n';
-
- for (int i = 0; i < rho[lev]->boxArray().size(); ++i)
- {
+
+ for (int i = 0; i < rho[lev]->boxArray().size(); ++i) {
amrex::Real dx[3] = {
geom[lev].CellSize(0),
geom[lev].CellSize(1),
@@ -252,12 +250,12 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
for (int n = 0; n < AMREX_SPACEDIM; n++)
HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n';
}
-
+
std::string PathNameInHeader = sLevel;
PathNameInHeader += BaseName;
HeaderFile << PathNameInHeader << '\n';
}
-
+
//
// We combine all of the multifabs
//
@@ -265,7 +263,7 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
rho[lev]->DistributionMap(),
nData, 0,
amrex::MFInfo());
-
+
//
// Cull data -- use no ghost cells.
//
@@ -280,16 +278,16 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
amr::AmrField_t::Copy(data, *efield[lev][0], 0, 2, 1, 0); // Ex
amr::AmrField_t::Copy(data, *efield[lev][1], 0, 3, 1, 0); // Ey
amr::AmrField_t::Copy(data, *efield[lev][2], 0, 4, 1, 0); // Ez
-
+
//
// Use the Full pathname when naming the MultiFab.
//
std::string TheFullPath = FullPath;
TheFullPath += BaseName;
-
+
amrex::VisMF::Write(data, TheFullPath, amrex::VisMF::NFiles, true);
}
-
+
if ( Ippl::myNode() == 0 )
HeaderFile.close();
}
@@ -297,8 +295,7 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho,
void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p,
const double& /*time*/,
- const double& gamma)
-{
+ const double& gamma) {
/* According to
*
* template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
@@ -314,24 +311,24 @@ void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p,
*
* ATTENTION: We need to Lorentz transform the particles.
*/
-
+
const AmrLayout_t* layout_p = static_cast<const AmrLayout_t*>(&bunch_p->getLayout());
const AmrPartBunch::pbase_t* amrpbase_p = bunch_p->getAmrParticleBase();
-
+
const int NProcs = amrex::ParallelDescriptor::NProcs();
const int IOProcNumber = amrex::ParallelDescriptor::IOProcessorNumber();
-
+
bool doUnlink = true;
-
+
//
// We store the particles in a subdirectory of "dir".
//
std::string pdir = dir_m;
-
+
if ( ! pdir.empty() && pdir[pdir.size()-1] != '/') {
pdir += '/';
}
-
+
pdir += "opal";
//
// Make the particle directories if they don't already exist.
@@ -341,49 +338,46 @@ void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p,
amrex::CreateDirectoryFailed(pdir);
}
}
-
+
// Force other processors to wait until directory is built.
Ippl::Comm->barrier();
-
+
//
// The header contains the info we need to read back in the particles.
//
// Only the I/O processor writes to the header file.
//
std::ofstream HdrFile;
-
+
long nParticles = bunch_p->getTotalNum();
-
+
// do not modify LocalNumPerLevel in here!!!
auto LocalNumPerLevel = amrpbase_p->getLocalNumPerLevel();
-
+
int nLevel = (&amrpbase_p->getAmrLayout())->maxLevel() + 1;
-
+
std::unique_ptr<size_t[]> partPerLevel( new size_t[nLevel] );
std::unique_ptr<size_t[]> globalPartPerLevel( new size_t[nLevel] );
-
+
for (size_t i = 0; i < LocalNumPerLevel.size(); ++i)
partPerLevel[i] = LocalNumPerLevel[i];
-
+
allreduce(*partPerLevel.get(),
*globalPartPerLevel.get(),
nLevel, std::plus<size_t>());
-
+
int finest_level = layout_p->finestLevel();
- if ( Ippl::myNode() == 0 )
- {
+ if ( Ippl::myNode() == 0 ) {
+
std::string HdrFileName = pdir;
-
if ( ! HdrFileName.empty() && HdrFileName[HdrFileName.size()-1] != '/') {
HdrFileName += '/';
-
}
-
+
HdrFileName += "Header";
-
HdrFile.open(HdrFileName.c_str(), std::ios::out|std::ios::trunc);
-
+
if ( ! HdrFile.good()) {
amrex::FileOpenFailed(HdrFileName);
}
@@ -399,20 +393,20 @@ void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p,
// AMREX_SPACEDIM and N for sanity checking.
//
HdrFile << AMREX_SPACEDIM << '\n';
-
+
// The number of extra real parameters
HdrFile << realData_m.size() << '\n';
for (std::size_t j = 0; j < realData_m.size(); ++j)
HdrFile << realData_m[j] << '\n';
-
+
// The number of extra int parameters
HdrFile << intData_m.size() << '\n';
for (std::size_t j = 0; j < intData_m.size(); ++j)
HdrFile << intData_m[j] << '\n';
-
+
// is_checkpoint
HdrFile << true << '\n';
-
+
//
// The total number of particles.
//
@@ -441,21 +435,21 @@ void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p,
int nOutFiles(256);
nOutFiles = std::max(1, std::min(nOutFiles, NProcs));
-
+
for (int lev = 0; lev <= finest_level; ++lev) {
bool gotsome = (globalPartPerLevel[lev] > 0);
//
// We store the particles at each level in their own subdirectory.
//
std::string LevelDir = pdir;
-
+
if (gotsome) {
if ( ! LevelDir.empty() && LevelDir[LevelDir.size()-1] != '/') {
LevelDir += '/';
}
-
+
LevelDir = amrex::Concatenate(LevelDir + "Level_", lev, 1);
-
+
if ( Ippl::myNode() == 0 ) {
if ( ! amrex::UtilCreateDirectory(LevelDir, 0755)) {
amrex::CreateDirectoryFailed(LevelDir);
@@ -466,20 +460,20 @@ void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p,
//
Ippl::Comm->barrier();
}
-
+
// Write out the header for each particle
if ( gotsome && Ippl::myNode() == 0 ) {
std::string HeaderFileName = LevelDir;
HeaderFileName += "/Particle_H";
std::ofstream ParticleHeader(HeaderFileName);
-
+
layout_p->ParticleBoxArray(lev).writeOn(ParticleHeader);
ParticleHeader << '\n';
ParticleHeader.flush();
ParticleHeader.close();
}
-
+
amrex::MFInfo info;
info.SetAlloc(false);
amr::AmrField_t state(layout_p->ParticleBoxArray(lev),
@@ -492,7 +486,7 @@ void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p,
amrex::Vector<int> which(state.size(),0);
amrex::Vector<int > count(state.size(),0);
amrex::Vector<long> where(state.size(),0);
-
+
std::string filePrefix(LevelDir);
filePrefix += '/';
filePrefix += "DATA_";
@@ -582,20 +576,20 @@ void AmrYtWriter::writeParticles_m(int level,
*
* in AMReX_ParticleContainerI.H
*/
-
+
const AmrLayout_t* layout_p = static_cast<const AmrLayout_t*>(&bunch_p->getLayout());
const AmrPartBunch::pbase_t* amrpbase_p = bunch_p->getAmrParticleBase();
-
+
const auto& LocalNumPerLevel = amrpbase_p->getLocalNumPerLevel();
size_t lBegin = LocalNumPerLevel.begin(level);
size_t lEnd = LocalNumPerLevel.end(level);
-
+
for (size_t ip = lBegin; ip < lEnd; ++ip) {
const int grid = amrpbase_p->Grid[ip];
-
+
count[grid] += 1;
}
-
+
amrex::MFInfo info;
info.SetAlloc(false);
amr::AmrField_t state(layout_p->ParticleBoxArray(level),
@@ -604,26 +598,25 @@ void AmrYtWriter::writeParticles_m(int level,
for (amrex::MFIter mfi(state); mfi.isValid(); ++mfi) {
const int grid = mfi.index();
-
+
which[grid] = fnum;
where[grid] = amrex::VisMF::FileOffset(ofs);
-
+
if (count[grid] == 0) {
continue;
}
-
-
+
// First write out the integer data in binary.
const int iChunkSize = 2 + intData_m.size();
amrex::Vector<int> istuff(count[grid]*iChunkSize);
int* iptr = istuff.dataPtr();
-
+
// inefficient
for (size_t ip = lBegin; ip < lEnd; ++ip) {
const int pGrid = amrpbase_p->Grid[ip];
-
+
if ( grid == pGrid ) {
-
+
iptr[0] = bunch_p->ID[ip];
iptr[1] = Ippl::myNode();
iptr[2] = bunch_p->Bin[ip];
@@ -632,55 +625,55 @@ void AmrYtWriter::writeParticles_m(int level,
}
amrex::writeIntData(istuff.dataPtr(), istuff.size(), ofs);
-
+
// Write the Real data in binary.
const int rChunkSize = AMREX_SPACEDIM + realData_m.size();
amrex::Vector<double> rstuff(count[grid]*rChunkSize);
double* rptr = rstuff.dataPtr();
-
-
+
+
// inefficient
for (size_t ip = lBegin; ip < lEnd; ++ip) {
const int pGrid = amrpbase_p->Grid[ip];
-
+
if ( grid == pGrid ) {
-
+
int idx = 0;
-
+
// coordinates
for (int j = 0; j < AMREX_SPACEDIM; ++j)
rptr[idx++] = bunch_p->R[ip](j);
-
+
// Lorentz transform
if ( OpalData::getInstance()->isInOPALCyclMode() )
rptr[1] *= gamma; // y is longitudinal direction
else
rptr[2] *= gamma; // z is longitudinal direction
-
+
// momenta
for (int j = 0; j < AMREX_SPACEDIM; ++j)
rptr[idx++] = bunch_p->P[ip](j);
-
+
// charge
rptr[idx++] = bunch_p->Q[ip];
-
+
// mass
rptr[idx++] = bunch_p->M[ip];
-
+
// timestep
rptr[idx++] = bunch_p->dt[ip];
-
+
// // potential
// rptr[idx++] = bunch_p->Phi[ip];
-
+
// electric field
for (int j = 0; j < AMREX_SPACEDIM; ++j)
rptr[idx++] = bunch_p->Ef[ip](j);
-
+
// magnetic field
for (int j = 0; j < AMREX_SPACEDIM; ++j)
rptr[idx++] = bunch_p->Bf[ip](j);
-
+
rptr += idx; //realData_m.size();
}
}
diff --git a/src/Amr/AmrYtWriter.h b/src/Amr/AmrYtWriter.h
index 581985071f34bcec56bb0f1fd172e92105d62664..65cd04bd6c874674734518c45cbb1d900ea2f5d9 100644
--- a/src/Amr/AmrYtWriter.h
+++ b/src/Amr/AmrYtWriter.h
@@ -31,16 +31,15 @@
#include <vector>
-class AmrYtWriter : public AbstractAmrWriter {
-
+class AmrYtWriter: public AbstractAmrWriter {
+
public:
-
/*!
* @param step we write
* @param bin energy bin we write (multi-bunch simulation)
*/
explicit AmrYtWriter(int step, int bin = 0);
-
+
/*!
* Write yt files to the simulation subdirectory
* data/amr/yt.
@@ -56,12 +55,11 @@ public:
const int& nLevel,
const double& time,
const double& scale);
-
-
+
void writeBunch(const AmrPartBunch* bunch_p,
const double& time,
const double& gamma);
-
+
private:
/* Copied and slightely modified version of
* AMReX_ParticleContainerI.H
@@ -83,10 +81,9 @@ private:
amrex::Vector<long>& where,
const AmrPartBunch* bunch_p,
const double gamma) const;
-
+
private:
std::string dir_m; ///< directory where to write files
-
std::vector<std::string> intData_m; ///< integer bunch data
std::vector<std::string> realData_m; ///< real bunch data
};
diff --git a/src/Amr/BoxLibLayout.h b/src/Amr/BoxLibLayout.h
index 29bab5306f822154189931000809f606f474e95e..bef7bbe038bcd0514c6fceff432650ab5c03a081 100644
--- a/src/Amr/BoxLibLayout.h
+++ b/src/Amr/BoxLibLayout.h
@@ -44,23 +44,23 @@
#include <AMReX_ParGDB.H>
template<class T, unsigned Dim>
-class BoxLibLayout : public ParticleAmrLayout<T, Dim>,
- public amrex::ParGDB
+class BoxLibLayout: public ParticleAmrLayout<T, Dim>,
+ public amrex::ParGDB
{
-
+
public:
typedef typename ParticleAmrLayout<T, Dim>::pair_t pair_t;
typedef typename ParticleAmrLayout<T, Dim>::pair_iterator pair_iterator;
typedef typename ParticleAmrLayout<T, Dim>::SingleParticlePos_t SingleParticlePos_t;
typedef typename ParticleAmrLayout<T, Dim>::Index_t Index_t;
-
+
typedef amr::AmrField_t AmrField_t;
typedef amr::AmrVectorField_t AmrVectorField_t;
typedef amr::AmrScalarFieldContainer_t AmrScalarFieldContainer_t;
typedef amr::AmrVectorFieldContainer_t AmrVectorFieldContainer_t;
typedef typename ParticleAmrLayout<T, Dim>::ParticlePos_t ParticlePos_t;
typedef ParticleAttrib<Index_t> ParticleIndex_t;
-
+
typedef amr::AmrProcMap_t AmrProcMap_t;
typedef amr::AmrGrid_t AmrGrid_t;
typedef amr::AmrGeometry_t AmrGeometry_t;
@@ -73,10 +73,10 @@ public:
typedef amr::AmrDomain_t AmrDomain_t;
typedef amr::AmrBox_t AmrBox_t;
typedef amr::AmrReal_t AmrReal_t;
-
+
typedef amrex::BaseFab<int> basefab_t;
typedef amrex::FabArray<basefab_t> mask_t;
-
+
/*!
* Lower physical domain boundary (each dimension). It has to be
* smaller than -1 since all particles are within \f$[-1, 1]^3\f$.
@@ -84,7 +84,7 @@ public:
* enlargement factor (in [%]) in BoxLibLayout::initBaseBox_m().
*/
static Vector_t lowerBound;
-
+
/*! Upper physical domain boundary (each dimension). It has to be
* greater than 1 since all particles are within \f$[-1, 1]^3\f$.
* The real computational domain is multiplied with the mesh
@@ -92,13 +92,13 @@ public:
*/
static Vector_t upperBound;
+
public:
-
/*!
* Initializes default Geometry, DistributionMapping and BoxArray.
*/
BoxLibLayout();
-
+
/*!
* Given a layout it copies that.
*/
@@ -109,7 +109,7 @@ public:
* @param maxGridSize for all levels.
*/
BoxLibLayout(int nGridPoints, int maxGridSize);
-
+
/*!
* Single-level constructor.
*
@@ -117,10 +117,10 @@ public:
* @param dmap is the distribution map for grids
* @param ba is the array of boxes for a level
*/
- BoxLibLayout(const AmrGeometry_t &geom,
- const AmrProcMap_t &dmap,
- const AmrGrid_t &ba);
-
+ BoxLibLayout(const AmrGeometry_t& geom,
+ const AmrProcMap_t& dmap,
+ const AmrGrid_t& ba);
+
/*!
* Multi-level constructor.
*
@@ -131,16 +131,15 @@ public:
* @param rr is the refinement ratio among the levels
* (always the ratio from l to l+1)
*/
- BoxLibLayout(const AmrGeomContainer_t &geom,
- const AmrProcMapContainer_t &dmap,
- const AmrGridContainer_t &ba,
- const AmrIntArray_t &rr);
-
-
+ BoxLibLayout(const AmrGeomContainer_t& geom,
+ const AmrProcMapContainer_t& dmap,
+ const AmrGridContainer_t& ba,
+ const AmrIntArray_t& rr);
+
/*
* Overloaded functions of ParticleAmrLayout
*/
-
+
/*!
* This method is used when creating the AMR object. OPAL
* takes the input argument BBOXINCR that is specified in
@@ -150,7 +149,7 @@ public:
* @param dh is the mesh enlargement factor
*/
void setBoundingBox(double dh);
-
+
/*!
* The Poisson computation domain is per default [-1,1]^3.
* With this method this can be changed in order to account
@@ -160,18 +159,17 @@ public:
*/
void setDomainRatio(const std::vector<double>& ratio);
-
/*
* Functions of IpplParticleBase
*/
-
+
/*!
* This method shouldn't be called. Otherwise
* it throws an exception.
*/
void update(IpplParticleBase< BoxLibLayout<T,Dim> >& PData,
const ParticleAttrib<char>* canSwap = 0);
-
+
/*!
* The proper update method for AMR.
*
@@ -185,12 +183,11 @@ public:
*/
void update(AmrParticleBase< BoxLibLayout<T,Dim> >& PData,
int lev_min = 0, int lev_max = -1, bool isRegrid = false);
-
-
+
/*
* Functions from AMReX that are adjusted to work with Ippl AmrParticleBase class
*/
-
+
/*!
* Get the cell of a particle
*
@@ -200,7 +197,7 @@ public:
*/
AmrIntVect_t Index (AmrParticleBase< BoxLibLayout<T,Dim> >& p,
const unsigned int ip, int level) const;
-
+
/*!
* Get the cell of a particle
*
@@ -208,13 +205,11 @@ public:
* @param lev is the level
*/
AmrIntVect_t Index (SingleParticlePos_t &R, int lev) const;
-
-
+
/*
* Additional methods
*/
-
-
+
/*!
* Build mask for a level used for interpolation from
* grid to particles to reduce spurious self field
@@ -225,8 +220,7 @@ public:
void clearLevelMask(int lev);
const std::unique_ptr<mask_t>& getLevelMask(int lev) const;
-
-
+
/*!
* The particles live initially on the coarsest level.
* Furthermore, the order the OPAL input file is parsed
@@ -248,7 +242,6 @@ public:
// this->m_rr.resize(maxLevel);
this->maxLevel_m = maxLevel;
}
-
/*!
* Set the geometry of the problem. It is called in
@@ -260,8 +253,7 @@ public:
for (unsigned int i = 0; i < geom.size(); ++i)
this->m_geom[i] = geom[i];
}
-
-
+
/*!
* Set the refinement ratios. It is called in
* AmrBoxLib::initBaseLevel_m().
@@ -273,44 +265,44 @@ public:
refRatio_m[i] = refRatio[i];
}
}
-
+
/*
* ParGDB overwritten functions
*/
-
+
/*!
* Check if an AMR level is well defined
*
* @param level to check
*/
inline bool LevelDefined (int level) const;
-
+
/*!
* @returns the current finest level
*/
inline int finestLevel () const;
-
+
/*!
* @returns the maximum level of simulation
*/
inline int maxLevel () const;
-
+
/*!
* @param level
* @returns the refinement ratio of this level to the next
* higher one
*/
inline AmrIntVect_t refRatio (int level) const;
-
+
/*!
* @param level
* @returns the maximum refinement ratio among all directions
* for the given level.
*/
inline int MaxRefRatio (int level) const;
-
+
+
private:
-
/*!
* Set up the box for the whole computation.
* The AMR object owning the bunch is not yet initialized.
@@ -320,12 +312,12 @@ private:
* @param dh is the mesh enlargement factor
*/
void initBaseBox_m(int nGridPoints, int maxGridSize, double dh = 0.04);
-
-
+
+
/*
* Functions from AMReX that are adjusted to work with Ippl AmrParticleBase class
*/
-
+
/*!
* Function from AMReX adjusted to work with Ippl AmrParticleBase class
* Checks/sets a particles location on levels lev_min and higher.
@@ -365,7 +357,7 @@ private:
* @returns true if the particle was shifted.
*/
bool PeriodicShift (SingleParticlePos_t R) const;
-
+
/*!
* Function from AMReX adjusted to work with Ippl AmrParticleBase class
*
@@ -378,12 +370,12 @@ private:
void locateParticle(AmrParticleBase< BoxLibLayout<T,Dim> >& p,
const unsigned int ip,
int lev_min, int lev_max, int nGrow) const;
-
+
+
private:
-
// don't use m_rr from ParGDB since it is the same refinement in all directions
AmrIntVectContainer_t refRatio_m; /// Refinement ratios [0:finest_level-1]
-
+
/* mask to reduce spurious self-field forces at
* coarse-fine interfaces
*/
diff --git a/src/Amr/BoxLibLayout.hpp b/src/Amr/BoxLibLayout.hpp
index e067ca277c98548107a4f6cbe027c4899f302fef..ac0a810dae4933e38bbf2b8fd819e41017b4db99 100644
--- a/src/Amr/BoxLibLayout.hpp
+++ b/src/Amr/BoxLibLayout.hpp
@@ -36,14 +36,18 @@
#ifndef BoxLibLayout_HPP
#define BoxLibLayout_HPP
-#include "BoxLibLayout.h"
+#include "Amr/BoxLibLayout.h"
+#include "Algorithms/Vektor.h"
#include "Message/Format.h"
#include "Message/MsgBuffer.h"
#include "Utility/PAssert.h"
#include "Utilities/OpalException.h"
#include <cmath>
+#include <utility>
+#include <vector>
+
template <class T, unsigned Dim>
Vector_t BoxLibLayout<T, Dim>::lowerBound = - Vector_t(1.0, 1.0, 1.0);
@@ -280,14 +284,12 @@ void BoxLibLayout<T, Dim>::update(AmrParticleBase< BoxLibLayout<T,Dim> >& PData,
std::vector<MsgBuffer*> buffers;
//create a message and send particles to nodes they belong to
- while (i!=p2n.end())
- {
+ while (i!=p2n.end()) {
unsigned cur_destination = i->first;
MsgBuffer *msgbuf = new MsgBuffer(format, p2n.count(i->first));
- for (; i!=p2n.end() && i->first == cur_destination; ++i)
- {
+ for (; i!=p2n.end() && i->first == cur_destination; ++i) {
Message msg;
PData.putMessage(msg, i->second);
PData.destroy(1, i->second);
@@ -305,54 +307,52 @@ void BoxLibLayout<T, Dim>::update(AmrParticleBase< BoxLibLayout<T,Dim> >& PData,
}
//destroy the particles that are sent to other domains
- if ( LocalNum < PData.getDestroyNum() )
+ if ( LocalNum < PData.getDestroyNum() ) {
throw OpalException("BoxLibLayout::update()",
"Rank " + std::to_string(myN) +
" can't destroy more particles than possessed.");
- else {
+ } else {
LocalNum -= PData.getDestroyNum(); // update local num
PData.performDestroy();
}
for (int lev = lev_min; lev <= lev_max; ++lev) {
- if ( LocalNumPerLevel[lev] < 0 )
+ if ( LocalNumPerLevel[lev] < 0 ) {
throw OpalException("BoxLibLayout::update()",
"Negative particle level count.");
+ }
}
//receive new particles
- for (int k = 0; k<msgrecv[myN]; ++k)
- {
+ for (int k = 0; k<msgrecv[myN]; ++k) {
int node = Communicate::COMM_ANY_NODE;
char *buffer = 0;
int bufsize = Ippl::Comm->raw_probe_receive(buffer, node, tag);
MsgBuffer recvbuf(format, buffer, bufsize);
Message *msg = recvbuf.get();
- while (msg != 0)
- {
- /* pBeginIdx is the start index of the new particle data
- * pEndIdx is the last index of the new particle data
- */
- size_t pBeginIdx = LocalNum;
+ while (msg != 0) {
+ /* pBeginIdx is the start index of the new particle data
+ * pEndIdx is the last index of the new particle data
+ */
+ size_t pBeginIdx = LocalNum;
- LocalNum += PData.getSingleMessage(*msg);
+ LocalNum += PData.getSingleMessage(*msg);
- size_t pEndIdx = LocalNum;
+ size_t pEndIdx = LocalNum;
- for (size_t idx = pBeginIdx; idx < pEndIdx; ++idx)
- ++LocalNumPerLevel[ PData.Level[idx] ];
+ for (size_t idx = pBeginIdx; idx < pEndIdx; ++idx)
+ ++LocalNumPerLevel[ PData.Level[idx] ];
- delete msg;
- msg = recvbuf.get();
- }
+ delete msg;
+ msg = recvbuf.get();
+ }
}
//wait for communication to finish and clean up buffers
MPI_Request* requests_ptr = requests.empty()? static_cast<MPI_Request*>(0): &(requests[0]);
MPI_Waitall(requests.size(), requests_ptr, MPI_STATUSES_IGNORE);
- for (unsigned int j = 0; j<buffers.size(); ++j)
- {
+ for (unsigned int j = 0; j<buffers.size(); ++j) {
delete buffers[j];
}
@@ -564,12 +564,10 @@ bool BoxLibLayout<T, Dim>::Where(AmrParticleBase< BoxLibLayout<T,Dim> >& p,
if (lev == (int)p.Level[ip]) {
// The fact that we are here means this particle does not belong to any finer grids.
- if (0 <= p.Grid[ip] && p.Grid[ip] < ba.size())
- {
+ if (0 <= p.Grid[ip] && p.Grid[ip] < ba.size()) {
const AmrBox_t& bx = ba.getCellCenteredBox(p.Grid[ip]);
const AmrBox_t& gbx = amrex::grow(bx,nGrow);
- if (gbx.contains(iv))
- {
+ if (gbx.contains(iv)) {
// if (bx != pld.m_gridbox || !pld.m_tilebox.contains(iv)) {
// pld.m_tile = getTileIndex(iv, bx, pld.m_tilebox);
// pld.m_gridbox = bx;
@@ -581,8 +579,7 @@ bool BoxLibLayout<T, Dim>::Where(AmrParticleBase< BoxLibLayout<T,Dim> >& p,
ba.intersections(AmrBox_t(iv, iv), isects, true, nGrow);
- if (!isects.empty())
- {
+ if (!isects.empty()) {
p.Level[ip] = lev;
p.Grid[ip] = isects[0].first;
@@ -613,19 +610,16 @@ bool BoxLibLayout<T, Dim>::EnforcePeriodicWhere (AmrParticleBase< BoxLibLayout<T
//
SingleParticlePos_t R = p.R[ip];
- if (PeriodicShift(R))
- {
+ if (PeriodicShift(R)) {
std::vector< std::pair<int, AmrBox_t> > isects;
- for (int lev = lev_max; lev >= lev_min; lev--)
- {
+ for (int lev = lev_max; lev >= lev_min; lev--) {
const AmrIntVect_t& iv = Index(R, lev);
const AmrGrid_t& ba = ParticleBoxArray(lev);
ba.intersections(AmrBox_t(iv,iv),isects,true,0);
- if (!isects.empty())
- {
+ if (!isects.empty()) {
D_TERM(p.R[ip][0] = R[0];,
p.R[ip][1] = R[1];,
p.R[ip][2] = R[2];);
@@ -658,20 +652,18 @@ bool BoxLibLayout<T, Dim>::PeriodicShift (SingleParticlePos_t R) const
const AmrIntVect_t& iv = Index(R, 0);
bool shifted = false;
- for (int i = 0; i < AMREX_SPACEDIM; i++)
- {
+ for (int i = 0; i < AMREX_SPACEDIM; i++) {
if (!geom.isPeriodic(i)) continue;
- if (iv[i] > dmn.bigEnd(i))
- {
- if (R[i] == geom.ProbHi(i))
+ if (iv[i] > dmn.bigEnd(i)) {
+ if (R[i] == geom.ProbHi(i)) {
//
// Don't let particles lie exactly on the domain face.
// Force the particle to be outside the domain so the
// periodic shift will bring it back inside.
//
R[i] += .125*geom.CellSize(i);
-
+ }
R[i] -= geom.ProbLength(i);
if (R[i] <= geom.ProbLo(i))
@@ -683,25 +675,24 @@ bool BoxLibLayout<T, Dim>::PeriodicShift (SingleParticlePos_t R) const
PAssert(R[i] >= geom.ProbLo(i));
shifted = true;
- }
- else if (iv[i] < dmn.smallEnd(i))
- {
- if (R[i] == geom.ProbLo(i))
+
+ } else if (iv[i] < dmn.smallEnd(i)) {
+ if (R[i] == geom.ProbLo(i)) {
//
// Don't let particles lie exactly on the domain face.
// Force the particle to be outside the domain so the
// periodic shift will bring it back inside.
//
R[i] -= .125*geom.CellSize(i);
-
+ }
R[i] += geom.ProbLength(i);
- if (R[i] >= geom.ProbHi(i))
+ if (R[i] >= geom.ProbHi(i)) {
//
// This can happen due to precision issues.
//
R[i] -= .125*geom.CellSize(i);
-
+ }
PAssert(R[i] <= geom.ProbHi(i));
shifted = true;
@@ -730,12 +721,10 @@ void BoxLibLayout<T, Dim>::locateParticle(
bool success = false;
- if (outside)
- {
+ if (outside) {
// Note that EnforcePeriodicWhere may shift the particle if it is successful.
success = EnforcePeriodicWhere(p, ip, lev_min, lev_max);
- if (!success && lev_min == 0)
- {
+ if (!success && lev_min == 0) {
// The particle has left the domain; invalidate it.
p.destroy(1, ip);
success = true;
@@ -747,19 +736,15 @@ void BoxLibLayout<T, Dim>::locateParticle(
"We're losing particles although we shouldn't");
}
- }
- else
- {
+ } else {
success = Where(p, ip, lev_min, lev_max);
}
- if (!success)
- {
+ if (!success) {
success = (nGrow > 0) && Where(p, ip, lev_min, lev_min, nGrow);
}
- if (!success)
- {
+ if (!success) {
std::stringstream ss;
ss << "Invalid particle with ID " << ip << " at position " << p.R[ip] << ".";
throw OpalException("BoxLibLayout::locateParticle()", ss.str());
@@ -788,8 +773,7 @@ int BoxLibLayout<T, Dim>::maxLevel () const {
template <class T, unsigned Dim>
typename BoxLibLayout<T, Dim>::AmrIntVect_t
- BoxLibLayout<T, Dim>::refRatio (int level) const
-{
+ BoxLibLayout<T, Dim>::refRatio (int level) const {
return refRatio_m[level];
}
@@ -806,8 +790,7 @@ int BoxLibLayout<T, Dim>::MaxRefRatio (int level) const {
template <class T, unsigned Dim>
void BoxLibLayout<T, Dim>::initBaseBox_m(int nGridPoints,
int maxGridSize,
- double dh)
-{
+ double dh) {
// physical box (in meters)
AmrDomain_t real_box;
for (int d = 0; d < AMREX_SPACEDIM; ++d) {
diff --git a/src/Amr/BoxLibParticle.h b/src/Amr/BoxLibParticle.h
index 51b69e8b219092db1c382535e47cae92e8c6f95d..83276a94e81e81fef2e9d16c09fe63681b59710b 100644
--- a/src/Amr/BoxLibParticle.h
+++ b/src/Amr/BoxLibParticle.h
@@ -35,8 +35,8 @@
template<class PLayout>
-class BoxLibParticle : public virtual AmrParticleBase<PLayout>
-{
+class BoxLibParticle: public virtual AmrParticleBase<PLayout> {
+
public:
typedef typename AmrParticleBase<PLayout>::ParticlePos_t ParticlePos_t;
typedef typename AmrParticleBase<PLayout>::ParticleIndex_t ParticleIndex_t;
@@ -47,24 +47,24 @@ public:
typedef typename AmrParticleBase<PLayout>::AmrVectorFieldContainer_t AmrVectorFieldContainer_t;
typedef typename AmrParticleBase<PLayout>::AmrVectorField_t AmrVectorField_t;
typedef typename AmrParticleBase<PLayout>::ParticleLevelCounter_t ParticleLevelCounter_t;
-
+
typedef typename PLayout::AmrProcMap_t AmrProcMap_t;
typedef typename PLayout::AmrGrid_t AmrGrid_t;
typedef typename PLayout::AmrGeometry_t AmrGeometry_t;
typedef typename PLayout::AmrIntVect_t AmrIntVect_t;
typedef typename PLayout::AmrBox_t AmrBox_t;
typedef typename PLayout::AmrReal_t AmrReal_t;
-
+
typedef amrex::FArrayBox FArrayBox_t;
-
+
public:
BoxLibParticle();
-
+
/*!
* @param layout that does the particle-to-core management
*/
BoxLibParticle(PLayout *layout);
-
+
/*!
* Multi-level scatter.
* Scatter the data from the given attribute onto the given field, using
@@ -83,8 +83,7 @@ public:
ParticleAttrib<Vektor<PT, Dim> >& pp,
int lbase, int lfine,
const ParticleAttrib<int>& pbin, int bin = -1);
-
-
+
/*!
* Single-level scatter.
* Scatter the data from the given attribute onto the given field, using
@@ -102,7 +101,7 @@ public:
ParticleAttrib<Vektor<PT, Dim> >& pp,
const ParticleAttrib<int>& pbin, int bin = -1,
int level = 0);
-
+
/*!
* Multi-level gather.
* Gather the data from the given Field into the given attribute, using
@@ -118,12 +117,12 @@ public:
void gather(ParticleAttrib<FT>& attrib, AmrVectorFieldContainer_t& f,
ParticleAttrib<Vektor<PT, Dim> >& pp,
int lbase, int lfine);
-
+
private:
/*
* AMReX functions adjusted to work with Ippl
*/
-
+
/*!
* Multi-level scatter (adjusted from AMReX).
*
@@ -135,11 +134,11 @@ private:
* @param finest_level level we want to end
*/
template <class AType>
- void AssignDensityFort(ParticleAttrib<AType> &pa,
+ void AssignDensityFort(ParticleAttrib<AType>& pa,
AmrScalarFieldContainer_t& mf_to_be_filled,
int lev_min, int ncomp, int finest_level,
const ParticleAttrib<int>& pbin, int bin = -1) const;
-
+
/*!
* Multi-level gather (adjusted from AMReX).
*
@@ -152,7 +151,7 @@ private:
void InterpolateFort(ParticleAttrib<AType> &pa,
AmrVectorFieldContainer_t& mesh_data,
int lev_min, int lev_max);
-
+
/*!
* Single-level gather (adjusted from AMReX).
*
@@ -162,8 +161,7 @@ private:
*/
template <class AType>
void InterpolateSingleLevelFort(ParticleAttrib<AType> &pa, AmrVectorField_t& mesh_data, int lev);
-
-
+
/*!
* Multi-level gather.
*
@@ -172,11 +170,10 @@ private:
* @param lev for which we get the mesh data
*/
template <class AType>
- void InterpolateMultiLevelFort(ParticleAttrib<AType> &pa,
+ void InterpolateMultiLevelFort(ParticleAttrib<AType>& pa,
AmrVectorFieldContainer_t& mesh_data,
int lev);
-
-
+
/*!
* Single-level scatter (adjusted from AMReX).
*
@@ -187,15 +184,15 @@ private:
* @param particle_lvl_offset is zero
*/
template <class AType>
- void AssignCellDensitySingleLevelFort(ParticleAttrib<AType> &pa, AmrField_t& mf, int level,
+ void AssignCellDensitySingleLevelFort(ParticleAttrib<AType>& pa, AmrField_t& mf, int level,
const ParticleAttrib<int>& pbin, int bin = -1,
int ncomp=1, int particle_lvl_offset = 0) const;
-
+
private:
IpplTimings::TimerRef AssignDensityTimer_m;
};
-#include "BoxLibParticle.hpp"
+#include "Amr/BoxLibParticle.hpp"
#endif
diff --git a/src/Amr/BoxLibParticle.hpp b/src/Amr/BoxLibParticle.hpp
index 7b33099a64cd715737f36d3b028c82cfa94be571..1bbe3a6109b214835d5a09bf4c4ce326e82c15f2 100644
--- a/src/Amr/BoxLibParticle.hpp
+++ b/src/Amr/BoxLibParticle.hpp
@@ -25,6 +25,7 @@
#define BOXLIB_PARTICLE_HPP
#include "Utilities/OpalException.h"
+#include "Utility/IpplTimings.h"
#include <AMReX_BLFort.H>
#include <AMReX_MultiFabUtil.H>
@@ -57,23 +58,23 @@ void BoxLibParticle<PLayout>::scatter(ParticleAttrib<FT>& attrib, AmrScalarField
this->scatter(attrib, *(f[lbase].get()), pp, pbin, bin, lbase);
return;
}
-
+
const PLayout *layout_p = &this->getLayout();
int nGrow = layout_p->refRatio(lbase)[0];
-
+
AmrScalarFieldContainer_t tmp(lfine+1);
for (int lev = lbase; lev <= lfine; ++lev) {
-
+
f[lev]->setVal(0.0, f[lev]->nGrow());
-
+
const AmrGrid_t& ba = f[lev]->boxArray();
const AmrProcMap_t& dm = f[lev]->DistributionMap();
tmp[lev].reset(new AmrField_t(ba, dm, 1, nGrow));
tmp[lev]->setVal(0.0, nGrow);
}
-
+
this->AssignDensityFort(attrib, tmp, lbase, 1, lfine, pbin, bin);
-
+
for (int lev = lbase; lev <= lfine; ++lev)
AmrField_t::Copy(*f[lev], *tmp[lev], 0, 0, f[lev]->nComp(), f[lev]->nGrow());
}
@@ -91,11 +92,11 @@ void BoxLibParticle<PLayout>::scatter(ParticleAttrib<FT>& attrib, AmrField_t& f,
AmrField_t tmp(ba, dmap, f.nComp(), 1);
tmp.setVal(0.0, 1);
-
+
this->AssignCellDensitySingleLevelFort(attrib, tmp, level, pbin, bin);
-
+
f.setVal(0.0, f.nGrow());
-
+
AmrField_t::Copy(f, tmp, 0, 0, f.nComp(), f.nGrow());
}
@@ -114,7 +115,7 @@ void BoxLibParticle<PLayout>::gather(ParticleAttrib<FT>& attrib, AmrVectorFieldC
// Scatter the particle attribute pa on the grid
template<class PLayout>
template <class AType>
-void BoxLibParticle<PLayout>::AssignDensityFort(ParticleAttrib<AType> &pa,
+void BoxLibParticle<PLayout>::AssignDensityFort(ParticleAttrib<AType>& pa,
AmrScalarFieldContainer_t& mf_to_be_filled,
int lev_min, int ncomp, int finest_level,
const ParticleAttrib<int>& pbin, int bin) const
@@ -122,16 +123,16 @@ void BoxLibParticle<PLayout>::AssignDensityFort(ParticleAttrib<AType> &pa,
// BL_PROFILE("AssignDensityFort()");
IpplTimings::startTimer(AssignDensityTimer_m);
const PLayout *layout_p = &this->getLayout();
-
+
// not done in amrex
int rho_index = 0;
-
+
amrex::PhysBCFunct cphysbc, fphysbc;
int lo_bc[] = {INT_DIR, INT_DIR, INT_DIR}; // periodic boundaries
int hi_bc[] = {INT_DIR, INT_DIR, INT_DIR};
amrex::Vector<amrex::BCRec> bcs(1, amrex::BCRec(lo_bc, hi_bc));
amrex::PCInterp mapper;
-
+
AmrScalarFieldContainer_t tmp(finest_level+1);
for (int lev = lev_min; lev <= finest_level; ++lev) {
const AmrGrid_t& ba = mf_to_be_filled[lev]->boxArray();
@@ -139,7 +140,7 @@ void BoxLibParticle<PLayout>::AssignDensityFort(ParticleAttrib<AType> &pa,
tmp[lev].reset(new AmrField_t(ba, dm, 1, 0));
tmp[lev]->setVal(0.0);
}
-
+
for (int lev = lev_min; lev <= finest_level; ++lev) {
AssignCellDensitySingleLevelFort(pa, *mf_to_be_filled[lev], lev, pbin, bin, 1, 0);
@@ -160,15 +161,15 @@ void BoxLibParticle<PLayout>::AssignDensityFort(ParticleAttrib<AType> &pa,
rho_index, 1, layout_p->refRatio(lev-1),
layout_p->Geom(lev-1), layout_p->Geom(lev));
}
-
+
mf_to_be_filled[lev]->plus(*tmp[lev], rho_index, ncomp, 0);
}
-
+
for (int lev = finest_level - 1; lev >= lev_min; --lev) {
amrex::average_down(*mf_to_be_filled[lev+1],
*mf_to_be_filled[lev], rho_index, ncomp, layout_p->refRatio(lev));
}
-
+
IpplTimings::stopTimer(AssignDensityTimer_m);
}
@@ -176,7 +177,7 @@ void BoxLibParticle<PLayout>::AssignDensityFort(ParticleAttrib<AType> &pa,
// This is the single-level version for cell-centered density
template<class PLayout>
template <class AType>
-void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AType> &pa,
+void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AType>& pa,
AmrField_t& mf_to_be_filled,
int level,
const ParticleAttrib<int>& pbin,
@@ -185,9 +186,9 @@ void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AT
int /*particle_lvl_offset*/) const
{
// BL_PROFILE("ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt>::AssignCellDensitySingleLevelFort()");
-
+
const PLayout *layout_p = &this->getLayout();
-
+
AmrField_t* mf_pointer;
if (layout_p->OnSameGrids(level, mf_to_be_filled)) {
@@ -220,30 +221,30 @@ void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AT
throw OpalException("BoxLibParticle::AssignCellDensitySingleLevelFort()",
"Problem must be periodic in no or all directions");
}
-
+
for (amrex::MFIter mfi(*mf_pointer); mfi.isValid(); ++mfi) {
(*mf_pointer)[mfi].setVal(0);
}
-
+
//loop trough particles and distribute values on the grid
const ParticleLevelCounter_t& LocalNumPerLevel = this->getLocalNumPerLevel();
size_t lBegin = LocalNumPerLevel.begin(level);
size_t lEnd = LocalNumPerLevel.end(level);
-
+
AmrReal_t inv_dx[3] = { 1.0 / dx[0], 1.0 / dx[1], 1.0 / dx[2] };
double lxyz[3] = { 0.0, 0.0, 0.0 };
double wxyz_hi[3] = { 0.0, 0.0, 0.0 };
double wxyz_lo[3] = { 0.0, 0.0, 0.0 };
int ijk[3] = {0, 0, 0};
-
+
for (size_t ip = lBegin; ip < lEnd; ++ip) {
-
+
if ( bin > -1 && pbin[ip] != bin )
continue;
-
+
const int grid = this->Grid[ip];
FArrayBox_t& fab = (*mf_pointer)[grid];
-
+
// not callable:
// begin amrex_deposit_cic(pbx.data(), nstride, N, fab.dataPtr(), box.loVect(), box.hiVect(), plo, dx);
for (int i = 0; i < 3; ++i) {
@@ -252,11 +253,11 @@ void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AT
wxyz_hi[i] = lxyz[i] - ijk[i];
wxyz_lo[i] = 1.0 - wxyz_hi[i];
}
-
+
int& i = ijk[0];
int& j = ijk[1];
int& k = ijk[2];
-
+
AmrIntVect_t i1(i-1, j-1, k-1);
AmrIntVect_t i2(i-1, j-1, k);
AmrIntVect_t i3(i-1, j, k-1);
@@ -265,7 +266,7 @@ void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AT
AmrIntVect_t i6(i, j-1, k);
AmrIntVect_t i7(i, j, k-1);
AmrIntVect_t i8(i, j, k);
-
+
fab(i1, 0) += wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*pa[ip];
fab(i2, 0) += wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*pa[ip];
fab(i3, 0) += wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*pa[ip];
@@ -276,7 +277,7 @@ void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AT
fab(i8, 0) += wxyz_hi[0]*wxyz_hi[1]*wxyz_hi[2]*pa[ip];
// end of amrex_deposit_cic
}
-
+
mf_pointer->SumBoundary(gm.periodicity());
// Only multiply the first component by (1/vol) because this converts mass
@@ -313,7 +314,7 @@ void BoxLibParticle<PLayout>::InterpolateFort(ParticleAttrib<AType> &pa,
template<class PLayout>
template <class AType>
-void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &pa,
+void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType>& pa,
AmrVectorField_t& mesh_data, int lev)
{
for (std::size_t i = 0; i < mesh_data.size(); ++i) {
@@ -321,23 +322,23 @@ void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &
throw OpalException("BoxLibParticle::InterpolateSingleLevelFort()",
"Must have at least one ghost cell when in InterpolateSingleLevelFort");
}
-
+
PLayout *layout_p = &this->getLayout();
-
+
const AmrGeometry_t& gm = layout_p->Geom(lev);
const AmrReal_t* plo = gm.ProbLo();
const AmrReal_t* dx = gm.CellSize();
-
+
//loop trough particles and distribute values on the grid
const ParticleLevelCounter_t& LocalNumPerLevel = this->getLocalNumPerLevel();
size_t lBegin = LocalNumPerLevel.begin(lev);
- size_t lEnd = LocalNumPerLevel.end(lev);
-
+ size_t lEnd = LocalNumPerLevel.end(lev);
+
// make sure that boundaries are filled!
for (std::size_t i = 0; i < mesh_data.size(); ++i) {
mesh_data[i]->FillBoundary(gm.periodicity());
}
-
+
AmrReal_t inv_dx[3] = { 1.0 / dx[0], 1.0 / dx[1], 1.0 / dx[2] };
double lxyz[3] = { 0.0, 0.0, 0.0 };
double wxyz_hi[3] = { 0.0, 0.0, 0.0 };
@@ -349,8 +350,7 @@ void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &
FArrayBox_t& exfab = (*mesh_data[0])[grid];
FArrayBox_t& eyfab = (*mesh_data[1])[grid];
FArrayBox_t& ezfab = (*mesh_data[2])[grid];
-
-
+
// not callable
// begin amrex_interpolate_cic(pbx.data(), nstride, N, fab.dataPtr(), box.loVect(), box.hiVect(), nComp, plo, dx);
for (int i = 0; i < 3; ++i) {
@@ -359,11 +359,11 @@ void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &
wxyz_hi[i] = lxyz[i] - ijk[i];
wxyz_lo[i] = 1.0 - wxyz_hi[i];
}
-
+
int& i = ijk[0];
int& j = ijk[1];
int& k = ijk[2];
-
+
AmrIntVect_t i1(i-1, j-1, k-1);
AmrIntVect_t i2(i-1, j-1, k);
AmrIntVect_t i3(i-1, j, k-1);
@@ -372,7 +372,7 @@ void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &
AmrIntVect_t i6(i, j-1, k);
AmrIntVect_t i7(i, j, k-1);
AmrIntVect_t i8(i, j, k);
-
+
pa[ip](0) = wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*exfab(i1) +
wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*exfab(i2) +
wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*exfab(i3) +
@@ -381,7 +381,7 @@ void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &
wxyz_hi[0]*wxyz_lo[1]*wxyz_hi[2]*exfab(i6) +
wxyz_hi[0]*wxyz_hi[1]*wxyz_lo[2]*exfab(i7) +
wxyz_hi[0]*wxyz_hi[1]*wxyz_hi[2]*exfab(i8);
-
+
pa[ip](1) = wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*eyfab(i1) +
wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*eyfab(i2) +
wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*eyfab(i3) +
@@ -390,7 +390,7 @@ void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &
wxyz_hi[0]*wxyz_lo[1]*wxyz_hi[2]*eyfab(i6) +
wxyz_hi[0]*wxyz_hi[1]*wxyz_lo[2]*eyfab(i7) +
wxyz_hi[0]*wxyz_hi[1]*wxyz_hi[2]*eyfab(i8);
-
+
pa[ip](2) = wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*ezfab(i1) +
wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*ezfab(i2) +
wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*ezfab(i3) +
@@ -406,7 +406,7 @@ void BoxLibParticle<PLayout>::InterpolateSingleLevelFort(ParticleAttrib<AType> &
template<class PLayout>
template <class AType>
-void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &pa,
+void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType>& pa,
AmrVectorFieldContainer_t& mesh_data,
int lev)
{
@@ -415,32 +415,32 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
throw OpalException("BoxLibParticle::InterpolateMultiLevelFort()",
"Must have at least one ghost cell when in InterpolateMultiLevelFort");
}
-
+
PLayout *layout_p = &this->getLayout();
-
+
const AmrGeometry_t& gm = layout_p->Geom(lev);
const AmrReal_t* plo = gm.ProbLo();
const AmrReal_t* fdx = gm.CellSize();
const AmrReal_t* cdx = layout_p->Geom(lev-1).CellSize();
const AmrReal_t* cplo = layout_p->Geom(lev-1).ProbLo();
-
+
//loop trough particles and distribute values on the grid
const ParticleLevelCounter_t& LocalNumPerLevel = this->getLocalNumPerLevel();
size_t lBegin = LocalNumPerLevel.begin(lev);
- size_t lEnd = LocalNumPerLevel.end(lev);
-
+ size_t lEnd = LocalNumPerLevel.end(lev);
+
// make sure that boundaries are filled!
for (std::size_t i = 0; i < mesh_data[lev].size(); ++i) {
mesh_data[lev][i]->FillBoundary(gm.periodicity());
}
-
+
AmrReal_t inv_fdx[3] = { 1.0 / fdx[0], 1.0 / fdx[1], 1.0 / fdx[2] };
AmrReal_t inv_cdx[3] = { 1.0 / cdx[0], 1.0 / cdx[1], 1.0 / cdx[2] };
double lxyz[3] = { 0.0, 0.0, 0.0 };
double wxyz_hi[3] = { 0.0, 0.0, 0.0 };
double wxyz_lo[3] = { 0.0, 0.0, 0.0 };
int ijk[3] = { 0, 0, 0 };
-
+
const AmrGrid_t& fba = mesh_data[lev][0]->boxArray();
const AmrProcMap_t& fdmap = mesh_data[lev][0]->DistributionMap();
AmrGrid_t cba = fba;
@@ -449,18 +449,18 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
AmrField_t cmesh_exdata(cba, fdmap,
mesh_data[lev][0]->nComp(),
mesh_data[lev][0]->nGrow());
-
+
cmesh_exdata.setVal(0.0, 0, 1, mesh_data[lev][0]->nGrow());
cmesh_exdata.copy(*mesh_data[lev-1][0], 0, 0, 1, 1, 1);
cmesh_exdata.FillBoundary(gm.periodicity());
-
+
AmrField_t cmesh_eydata(cba, fdmap,
mesh_data[lev][1]->nComp(),
mesh_data[lev][1]->nGrow());
cmesh_eydata.setVal(0.0, 0, 1, mesh_data[lev][1]->nGrow());
cmesh_eydata.copy(*mesh_data[lev-1][1], 0, 0, 1, 1, 1);
cmesh_eydata.FillBoundary(gm.periodicity());
-
+
AmrField_t cmesh_ezdata(cba, fdmap,
mesh_data[lev][2]->nComp(),
mesh_data[lev][2]->nGrow());
@@ -469,19 +469,19 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
cmesh_ezdata.FillBoundary(gm.periodicity());
for (size_t ip = lBegin; ip < lEnd; ++ip) {
-
+
const int grid = this->Grid[ip];
-
+
FArrayBox_t& exfab = (*(mesh_data[lev][0]))[grid];
FArrayBox_t& eyfab = (*(mesh_data[lev][1]))[grid];
FArrayBox_t& ezfab = (*(mesh_data[lev][2]))[grid];
-
+
FArrayBox_t& cexfab = cmesh_exdata[grid];
FArrayBox_t& ceyfab = cmesh_eydata[grid];
FArrayBox_t& cezfab = cmesh_ezdata[grid];
-
+
const typename PLayout::basefab_t& mfab = (*layout_p->getLevelMask(lev))[grid];
-
+
// not callable
// begin amrex_interpolate_cic(pbx.data(), nstride, N, fab.dataPtr(), box.loVect(), box.hiVect(), nComp, plo, dx);
for (int ii = 0; ii < 3; ++ii) {
@@ -490,18 +490,16 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
wxyz_hi[ii] = lxyz[ii] - ijk[ii];
wxyz_lo[ii] = 1.0 - wxyz_hi[ii];
}
-
+
int& i = ijk[0];
int& j = ijk[1];
int& k = ijk[2];
-
+
bool use_coarse = false;
-
+
// AMReX: electrostatic_pic_3d.f90
// use the coarse E near the level boundary
if ( mfab(AmrIntVect_t(i-1, j-1, k-1)) == 1 ) {
-
-
for (int ii = 0; ii < 3; ++ii) {
lxyz[ii] = ( this->R[ip](ii) - cplo[ii] ) * inv_cdx[ii] + 0.5;
ijk[ii] = lxyz[ii];
@@ -510,12 +508,12 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
}
use_coarse = true;
}
-
+
AmrIntVect_t i1(i-1, j-1, k-1);
AmrIntVect_t i3(i-1, j, k-1);
AmrIntVect_t i5(i, j-1, k-1);
AmrIntVect_t i7(i, j, k-1);
-
+
AmrIntVect_t i2(i-1, j-1, k);
AmrIntVect_t i4(i-1, j, k);
AmrIntVect_t i6(i, j-1, k);
@@ -529,7 +527,7 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * cexfab(i4) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * cexfab(i6) +
wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * cexfab(i8);
-
+
pa[ip](1) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * ceyfab(i1) +
wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * ceyfab(i3) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * ceyfab(i5) +
@@ -538,7 +536,7 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * ceyfab(i4) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * ceyfab(i6) +
wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * ceyfab(i8);
-
+
pa[ip](2) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * cezfab(i1) +
wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * cezfab(i3) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * cezfab(i5) +
@@ -556,7 +554,7 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * exfab(i4) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * exfab(i6) +
wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * exfab(i8);
-
+
pa[ip](1) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * eyfab(i1) +
wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * eyfab(i3) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * eyfab(i5) +
@@ -565,7 +563,7 @@ void BoxLibParticle<PLayout>::InterpolateMultiLevelFort(ParticleAttrib<AType> &p
wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * eyfab(i4) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * eyfab(i6) +
wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * eyfab(i8);
-
+
pa[ip](2) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * ezfab(i1) +
wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * ezfab(i3) +
wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * ezfab(i5) +
diff --git a/src/Classic/Algorithms/AmrPartBunch.cpp b/src/Classic/Algorithms/AmrPartBunch.cpp
index 707f90d4eb1c19e52aba6b6ca0e12a0cf122942a..c3dbfa92fb8b4044c1927aa03f4fb037524f021d 100644
--- a/src/Classic/Algorithms/AmrPartBunch.cpp
+++ b/src/Classic/Algorithms/AmrPartBunch.cpp
@@ -20,7 +20,7 @@
//
#include "AmrPartBunch.h"
-#include "Utilities/OpalException.h"
+#include "Utilities/GeneralClassicException.h"
AmrPartBunch::AmrPartBunch(const PartData *ref)
: PartBunchBase<double, 3>(new AmrPartBunch::pbase_t(new AmrLayout_t()), ref)
@@ -40,9 +40,7 @@ AmrPartBunch::AmrPartBunch(const PartData *ref, pbase_t* pbase_p)
amrpbase_mp->initializeAmr();
}
-AmrPartBunch::~AmrPartBunch() {
-
-}
+AmrPartBunch::~AmrPartBunch() { }
AmrPartBunch::pbase_t *AmrPartBunch::getAmrParticleBase() {
@@ -61,7 +59,7 @@ void AmrPartBunch::initialize(FieldLayout_t */*fLayout*/) {
void AmrPartBunch::do_binaryRepart() {
-
+
if ( amrobj_mp && !amrobj_mp->isRefined() ) {
/* In the first call to this function we
* intialize all fine levels
@@ -101,11 +99,11 @@ Vector_t AmrPartBunch::get_hr() const {
void AmrPartBunch::set_meshEnlargement(double dh) {
// set dh_m = dh
PartBunchBase<double, 3>::set_meshEnlargement(dh);
-
+
// update base geometry with new mesh enlargement
AmrLayout_t* layout_p = &amrpbase_mp->getAmrLayout();
layout_p->setBoundingBox(dh);
-
+
// if amrobj_mp != nullptr --> we need to regrid
this->do_binaryRepart();
}
@@ -131,7 +129,7 @@ double AmrPartBunch::getRho(int x, int y, int z) {
FieldLayout_t &AmrPartBunch::getFieldLayout() {
//TODO Implement
- throw OpalException("AmrPartBunch::getFieldLayout() ", "Not yet Implemented.");
+ throw GeneralClassicException("AmrPartBunch::getFieldLayout", "Not yet Implemented.");
return *fieldlayout_m;
}
@@ -142,50 +140,50 @@ void AmrPartBunch::boundp() {
if ( !(R.isDirty() || ID.isDirty() ) && stateOfLastBoundP_ == unit_state_) return; //-DW
stateOfLastBoundP_ = unit_state_;
-
+
if ( amrobj_mp ) {
/* we do an explicit domain mapping of the particles and then
* forbid it during the regrid process, this way it's only
* executed ones --> saves computation
*/
bool isForbidTransform = amrpbase_mp->isForbidTransform();
-
+
if ( !isForbidTransform ) {
this->updateLorentzFactor();
// Lorentz transform + mapping to [-1,1]
amrpbase_mp->domainMapping();
amrpbase_mp->setForbidTransform(true);
}
-
+
this->update();
-
+
if ( !isForbidTransform ) {
amrpbase_mp->setForbidTransform(false);
// map particles back + undo Lorentz transform
amrpbase_mp->domainMapping(true);
}
-
+
} else {
// At this point an amrobj_mp needs already be set
- throw GeneralClassicException("AmrPartBunch::boundp() ",
+ throw GeneralClassicException("AmrPartBunch::boundp",
"AmrObject pointer is not set.");
}
-
+
R.resetDirtyFlag();
-
+
IpplTimings::stopTimer(boundpTimer_m);
}
void AmrPartBunch::computeSelfFields() {
IpplTimings::startTimer(selfFieldTimer_m);
-
+
if ( !fs_m->hasValidSolver() )
- throw OpalException("AmrPartBunch::computeSelfFields() ",
- "No field solver.");
-
+ throw GeneralClassicException("AmrPartBunch::computeSelfFields",
+ "No field solver.");
+
amrobj_mp->computeSelfFields();
-
+
IpplTimings::stopTimer(selfFieldTimer_m);
}
@@ -206,15 +204,15 @@ void AmrPartBunch::computeSelfFields_cycl(double gamma) {
void AmrPartBunch::computeSelfFields_cycl(int bin) {
IpplTimings::startTimer(selfFieldTimer_m);
-
+
/* make sure it is refined in multi-bunch case
*/
if ( !amrobj_mp->isRefined() ) {
amrobj_mp->initFineLevels();
}
-
+
amrobj_mp->computeSelfFields_cycl(bin);
-
+
IpplTimings::stopTimer(selfFieldTimer_m);
}
@@ -226,19 +224,19 @@ void AmrPartBunch::setAmrDomainRatio(const std::vector<double>& ratio) {
void AmrPartBunch::gatherLevelStatistics() {
int nLevel = amrobj_mp->maxLevel() + 1;
-
+
std::unique_ptr<size_t[]> partPerLevel( new size_t[nLevel] );
globalPartPerLevel_m.reset( new size_t[nLevel] );
-
+
for (int i = 0; i < nLevel; ++i)
partPerLevel[i] = globalPartPerLevel_m[i] = 0.0;
-
+
// do not modify LocalNumPerLevel in here!!!
auto& LocalNumPerLevel = amrpbase_mp->getLocalNumPerLevel();
-
+
for (size_t i = 0; i < LocalNumPerLevel.size(); ++i)
partPerLevel[i] = LocalNumPerLevel[i];
-
+
reduce(*partPerLevel.get(),
*globalPartPerLevel_m.get(),
nLevel, std::plus<size_t>());
@@ -256,8 +254,7 @@ void AmrPartBunch::updateLorentzFactor(int bin) {
if ( this->weHaveBins() ) {
gamma = this->getBinGamma(bin);
}
-
-
+
/* At the beginning of simulation the Lorentz factor
* is not yet set since PartBunchBase::calcBeamParameters
* is not yet called. Therefore, we do this ugly workaround.
@@ -266,7 +263,7 @@ void AmrPartBunch::updateLorentzFactor(int bin) {
if ( gamma >= 1.0 ) {
init = false;
}
-
+
if ( init ) {
gamma = 1.0;
}
@@ -289,9 +286,7 @@ void AmrPartBunch::updateLorentzFactor(double gamma) {
}
-void AmrPartBunch::updateFieldContainers_m() {
-
-}
+void AmrPartBunch::updateFieldContainers_m() { }
void AmrPartBunch::updateDomainLength(Vektor<int, 3>& grid) {
grid = amrobj_mp->getBaseLevelGridPoints();
diff --git a/src/Classic/Algorithms/AmrPartBunch.h b/src/Classic/Algorithms/AmrPartBunch.h
index 88845b5076762a8d630a0a950cf3a2d9ea7ba6b1..8a6e232f0ae2810a85af8965db41ffe694582760 100644
--- a/src/Classic/Algorithms/AmrPartBunch.h
+++ b/src/Classic/Algorithms/AmrPartBunch.h
@@ -24,73 +24,71 @@
#include "Algorithms/PartBunchBase.h"
#include "Amr/AmrObject.h"
-class AmrPartBunch : public PartBunchBase<double, 3>
-{
+class AmrPartBunch: public PartBunchBase<double, 3> {
+
public:
typedef AmrParticle_t pbase_t;
-
-public:
- AmrPartBunch(const PartData *ref);
+public:
+ AmrPartBunch(const PartData* ref);
- AmrPartBunch(const PartData *ref, pbase_t* pbase_p);
+ AmrPartBunch(const PartData* ref, pbase_t* pbase_p);
~AmrPartBunch();
- pbase_t *getAmrParticleBase();
+ pbase_t* getAmrParticleBase();
const pbase_t *getAmrParticleBase() const;
void initialize(FieldLayout_t *fLayout);
-
+
// does actually another repartition
void do_binaryRepart();
-
+
Vector_t get_hr() const;
-
+
void set_meshEnlargement(double dh);
-
+
VectorPair_t getEExtrema();
-
+
double getRho(int x, int y, int z);
-
- FieldLayout_t &getFieldLayout();
-
-
+
+ FieldLayout_t& getFieldLayout();
+
void boundp();
-
+
void computeSelfFields();
-
+
void computeSelfFields(int bin);
-
+
void computeSelfFields_cycl(double gamma);
-
+
void computeSelfFields_cycl(int bin);
-
- void setSolver(FieldSolver *fs) {
+
+ void setSolver(FieldSolver* fs) {
PartBunchBase<double, 3>::setSolver(fs);
this->amrobj_mp = fs->getAmrObject();
}
-
+
virtual void setBinCharge(int /*bin*/, double /*q*/) { };
virtual void setBinCharge(int /*bin*/) { };
-
+
/*
* AmrPartBunch only
*/
-
+
const AmrObject* getAmrObject() const {
return this->amrobj_mp;
}
-
+
PoissonSolver *getFieldSolver() {
return fs_m->solver_m;
}
-
+
const PoissonSolver *getFieldSolver() const {
return fs_m->solver_m;
}
-
+
void setBaseLevelMeshSpacing(const Vector_t& hr) {
for (int i = 0; i < 3; ++i)
hr_m[i] = hr[i];
@@ -103,14 +101,13 @@ public:
void setAmrDomainRatio(const std::vector<double>& ratio);
void gatherLevelStatistics();
-
+
/*!
* Only a valid call of root core (core 0)
* @param l is the level
*/
const size_t& getLevelStatistics(int l) const;
-
-
+
/*!
* Update the Lorentz factor before every domainMapping in order
* to have the correct boosted frame for the particle redistribution,
@@ -118,7 +115,7 @@ public:
* @param bin is the energy bin
*/
void updateLorentzFactor(int bin=0);
-
+
/*!
* Update the Lorentz factor before every domainMapping in order
* to have the correct boosted frame for the particle redistribution,
@@ -127,33 +124,32 @@ public:
* @param gamma is the Lorentz factor
*/
void updateLorentzFactor(double gamma);
-
+
//FIXME BCs
void setBCAllPeriodic() {}
void setBCAllOpen() {}
void setBCForDCBeam() {}
-
-
+
+
private:
void updateFieldContainers_m();
-
+
void updateDomainLength(Vektor<int, 3>& grid);
-
+
void updateFields(const Vector_t& hr, const Vector_t& origin);
-
+
private:
-
/* pointer to AMR object that is part
* of solver_m (AmrPoissonSolver) in src/Structure/FieldSolver.h
*/
AmrObject *amrobj_mp;
pbase_t *amrpbase_mp;
-
+
/* We need this due to H5PartWrapper etc, but it's always nullptr.
* Thus, don't use it.
*/
FieldLayout_t* fieldlayout_m;
-
+
std::unique_ptr<size_t[]> globalPartPerLevel_m;
};
diff --git a/src/Classic/Utilities/Util.cpp b/src/Classic/Utilities/Util.cpp
index a3db487292f50ca76da8c5308e5643a39726399b..72e73b14eca7255054d64d4031c6f858739798b4 100644
--- a/src/Classic/Utilities/Util.cpp
+++ b/src/Classic/Utilities/Util.cpp
@@ -215,6 +215,12 @@ namespace Util {
}
}
+ bool isAllDigits(const std::string& str) {
+ return std::all_of(str.begin(),
+ str.end(),
+ [](char c) { return std::isdigit(c); });
+ }
+
KahanAccumulation::KahanAccumulation():
sum(0.0),
correction(0.0)
diff --git a/src/Classic/Utilities/Util.h b/src/Classic/Utilities/Util.h
index 5edf649f77659648f289cedbfe7e58bef014b6d5..a6fdda8035487ea08f51569f6d3392583cf37aa1 100644
--- a/src/Classic/Utilities/Util.h
+++ b/src/Classic/Utilities/Util.h
@@ -220,6 +220,10 @@ namespace Util {
void checkInt(double real, std::string name, double tolerance = 1e-9);
+ /// Check if there are only digits in a string
+ /// from https://stackoverflow.com/questions/19678572/how-to-validate-that-there-are-only-digits-in-a-string
+ bool isAllDigits(const std::string& str);
+
template<class IteratorIn, class IteratorOut>
void toString(IteratorIn first, IteratorIn last, IteratorOut out);
diff --git a/src/Structure/FieldSolver.cpp b/src/Structure/FieldSolver.cpp
index c5c7fed2e3b69c8793eb96e3a5821313c111501b..cbbc8ea07dc98b57dcc03a64306ac35cc856740c 100644
--- a/src/Structure/FieldSolver.cpp
+++ b/src/Structure/FieldSolver.cpp
@@ -45,10 +45,6 @@
#include "Solvers/AMReXSolvers/MLPoissonSolver.h"
#include "Amr/AmrDefs.h"
#include "Algorithms/AmrPartBunch.h"
-
- #include <algorithm>
- #include <cctype>
- #include <functional>
#endif
#ifdef HAVE_AMR_MG_SOLVER
@@ -274,7 +270,7 @@ FieldSolver::FieldSolver():
8);
itsAttr[AMR_TAGGING] = Attributes::makeUpperCaseString("AMR_TAGGING",
- "Refinement criteria [CHARGE_DENSITY | POTENTIAL | EFIELD]",
+ "Refinement criteria [CHARGE_DENSITY | POTENTIAL | EFIELD | MIN_NUM_PARTICLES | MAX_NUM_PARTICLES]",
"CHARGE_DENSITY");
itsAttr[AMR_DENSITY] = Attributes::makeReal("AMR_DENSITY",
@@ -713,15 +709,10 @@ Inform& FieldSolver::printInfo(Inform& os) const {
std::string FieldSolver::getTagging_m() const {
std::string tagging = Attributes::getString(itsAttr[AMR_TAGGING]);
- std::function<bool(const std::string&)> all_digits = [](const std::string& s) {
- // 15. Feb. 2019
- // https://stackoverflow.com/questions/19678572/how-to-validate-that-there-are-only-digits-in-a-string
- return std::all_of(s.begin(), s.end(),
- [](char c) { return std::isdigit(c); });
- };
-
- if ( all_digits(tagging) )
- tagging = AmrObject::enum2string(std::stoi(tagging));
+ // This conversion was introduced to allow numbers for tagging (useful for the sampler).
+ if ( Util::isAllDigits(tagging) ) {
+ tagging = AmrObject::getTaggingString(std::stoi(tagging));
+ }
return tagging;
}
diff --git a/src/Structure/FieldSolver.h b/src/Structure/FieldSolver.h
index 0ecfca7751086c321c769c389ba63e57607edee9..f8939e63ab1eed2be705d49188b5edeef7300714 100644
--- a/src/Structure/FieldSolver.h
+++ b/src/Structure/FieldSolver.h
@@ -98,9 +98,9 @@ public:
void setFieldSolverType();
FieldSolverType getFieldSolverType() const;
- inline Layout_t &getParticleLayout() { return* PL_m; }
+ inline Layout_t& getParticleLayout() { return* PL_m; }
- FieldLayout_t *getFieldLayout() { return FL_m; }
+ FieldLayout_t* getFieldLayout() { return FL_m; }
Inform& printInfo(Inform& os) const;
@@ -111,11 +111,11 @@ public:
bool isAmrSolverType() const;
#ifdef ENABLE_AMR
- AmrObject *getAmrObject() {
+ AmrObject* getAmrObject() {
return itsAmrObject_mp.get();
}
-
- const AmrObject *getAmrObject() const {
+
+ const AmrObject* getAmrObject() const {
return itsAmrObject_mp.get();
}
#endif
@@ -129,9 +129,9 @@ private:
std::string getTagging_m() const;
void initAmrObject_m();
-
+
void initAmrSolver_m();
-
+
std::unique_ptr<AmrObject> itsAmrObject_mp;
#endif