Commit d48062ee authored by snuverink_j's avatar snuverink_j
Browse files

Merge branch 'master' into cof

parents 8820018f 15ff33df
......@@ -164,7 +164,7 @@ if (BUILD_OPAL_UNIT_TESTS)
endif ()
option (ENABLE_AMR "Enable AMReX based AMR solver" OFF)
IF (ENABLE_AMR)
if (ENABLE_AMR)
message ("Enable AMR_SOLVER " ${ENABLE_AMR})
enable_language (Fortran)
......@@ -179,7 +179,6 @@ IF (ENABLE_AMR)
add_definitions(${AMREX_DEFINES})
add_definitions(-DENABLE_AMR)
add_compile_options (-Wno-unused-variable -Wno-unused-but-set-variable -Wno-maybe-uninitialized)
endif ()
option (ENABLE_AMR_MG_SOLVER "Enable AMR MG solver" OFF)
......@@ -215,6 +214,13 @@ if (ENABLE_SAAMG_SOLVER OR ENABLE_AMR_MG_SOLVER)
endif ()
endif ()
option (ENABLE_OPAL_FEL "Enable OPAL FEL" OFF)
if (ENABLE_OPAL_FEL)
message ("Enable OPAL FEL: " ${ENABLE_OPAL_FEL})
find_package (MITHRA MODULE REQUIRED)
add_DEFINITIONS (-DOPAL_FEL)
endif()
option (DBG_SCALARFIELD "Enable dump of scalar field rho_m" OFF)
if (DBG_SCALARFIELD)
message ("\nWrite scalar rho_m field is enabled ")
......
#
# Find MITHRA package
# https://github.com/aryafallahi/mithra
#
# MITHRA_INCLUDE_DIR
# MITHRA_FOUND
find_path (MITHRA_INCLUDE_DIR mithra/classes.hh
HINTS $ENV{MITHRA_INCLUDE_PATH} $ENV{MITHRA_INCLUDE_DIR} $ENV{MITHRA_PREFIX}/include $ENV{MITHRA_DIR}/include $ENV{MITHRA}/include
PATHS ENV C_INCLUDE_PATH CPLUS_INCLUDE_PATH
)
if (MITHRA_INCLUDE_DIR)
set (MITHRA_FOUND "YES")
endif ()
if (MITHRA_FOUND)
if (NOT MITHRA_FIND_QUIETLY)
message (STATUS "Found MITHRA include dir: ${MITHRA_INCLUDE_DIR}")
endif ()
else (MITHRA_FOUND)
if (MITHRA_FIND_REQUIRED)
message (FATAL_ERROR "Could not find MITHRA!")
endif (MITHRA_FIND_REQUIRED)
endif (MITHRA_FOUND)
......@@ -32,11 +32,13 @@ template<class PLayout>
class AmrParticleBase : public IpplParticleBase<PLayout> {
public:
typedef typename PLayout::ParticlePos_t ParticlePos_t;
typedef typename PLayout::ParticleIndex_t ParticleIndex_t;
typedef typename PLayout::SingleParticlePos_t SingleParticlePos_t;
typedef typename PLayout::AmrField_t AmrField_t;
typedef typename PLayout::AmrFieldContainer_t AmrFieldContainer_t;
typedef typename PLayout::ParticlePos_t ParticlePos_t;
typedef typename PLayout::ParticleIndex_t ParticleIndex_t;
typedef typename PLayout::SingleParticlePos_t SingleParticlePos_t;
typedef typename PLayout::AmrField_t AmrField_t;
typedef typename PLayout::AmrVectorField_t AmrVectorField_t;
typedef typename PLayout::AmrScalarFieldContainer_t AmrScalarFieldContainer_t;
typedef typename PLayout::AmrVectorFieldContainer_t AmrVectorFieldContainer_t;
typedef long SortListIndex_t;
typedef std::vector<SortListIndex_t> SortList_t;
......@@ -125,7 +127,8 @@ public:
inline bool isForbidTransform() const;
/*!
* Linear mapping to AMReX computation domain [-1, 1]^3. All dimensions
* Linear mapping to AMReX computation domain [-1, 1]^3 including the Lorentz
* transform. All dimensions
* are mapped by the same scaling factor.
* The potential and electric field need to be scaled afterwards appropriately.
* @param PData is the particle data
......@@ -140,6 +143,15 @@ public:
*/
inline const double& getScalingFactor() const;
void setLorentzFactor(const Vector_t& lorentzFactor);
// void lorentzTransform(bool inverse = false);
private:
void getLocalBounds_m(Vector_t &rmin, Vector_t &rmax);
void getGlobalBounds_m(Vector_t &rmin, Vector_t &rmax);
protected:
IpplTimings::TimerRef updateParticlesTimer_m;
IpplTimings::TimerRef sortParticlesTimer_m;
......@@ -153,6 +165,16 @@ protected:
*/
double scale_m;
/*!
* Lorentz factor used for the domain mapping.
* Is updated in AmrBoxLib
*
*/
Vector_t lorentzFactor_m;
// bool isLorentzTransformed_m;
private:
ParticleLevelCounter_t LocalNumPerLevel_m;
};
......
......@@ -7,6 +7,8 @@
template<class PLayout>
AmrParticleBase<PLayout>::AmrParticleBase() : forbidTransform_m(false),
scale_m(1.0),
lorentzFactor_m(1.0, 1.0, 1.0),
// isLorentzTransformed_m(false),
LocalNumPerLevel_m()
{
updateParticlesTimer_m = IpplTimings::getTimer("AMR update particles");
......@@ -20,6 +22,8 @@ AmrParticleBase<PLayout>::AmrParticleBase(PLayout* layout)
: IpplParticleBase<PLayout>(layout),
forbidTransform_m(false),
scale_m(1.0),
lorentzFactor_m(1.0, 1.0, 1.0),
// isLorentzTransformed_m(false),
LocalNumPerLevel_m()
{
updateParticlesTimer_m = IpplTimings::getTimer("AMR update particles");
......@@ -94,18 +98,32 @@ void AmrParticleBase<PLayout>::performDestroy(bool updateLocalNum) {
template<class PLayout>
void AmrParticleBase<PLayout>::create(size_t M) {
// size_t localnum = LocalNumPerLevel_m[0];
// particles are created at the coarsest level
LocalNumPerLevel_m[0] += M;
IpplParticleBase<PLayout>::create(M);
// for (size_t i = localnum; i < LocalNumPerLevel_m[0]; ++i) {
// this->Grid[i] = 0;
// this->Level[i] = 0;
// }
}
template<class PLayout>
void AmrParticleBase<PLayout>::createWithID(unsigned id) {
// size_t localnum = LocalNumPerLevel_m[0];
++LocalNumPerLevel_m[0];
IpplParticleBase<PLayout>::createWithID(id);
// this->Grid[localnum] = 0;
// this->Level[localnum] = 0;
}
......@@ -214,10 +232,35 @@ const double& AmrParticleBase<PLayout>::domainMapping(bool inverse) {
double scale = scale_m;
Vector_t gamma = lorentzFactor_m;
if ( !inverse ) {
Vector_t rmin, rmax;
bounds(this->R, rmin, rmax);
// if ( !this->DestroyList.empty() ) {
// this->performDestroy(true);
// }
Vector_t rmin = Vector_t(0.0, 0.0, 0.0);
Vector_t rmax = Vector_t(0.0, 0.0, 0.0);
getGlobalBounds_m(rmin, rmax);
/* in case of 1 particle, the bunch is rotated
* transformed to the local frame such that this
* particle lies on the origin (0, 0, 0).
*/
if ( this->getTotalNum() == 1 ||
(rmin == Vector_t(0.0, 0.0, 0.0) && rmax == Vector_t( 0.0, 0.0, 0.0)) ) {
rmin = Vector_t(-1.0, -1.0, -1.0);
rmax = Vector_t( 1.0, 1.0, 1.0);
}
/* Lorentz transfomration factor
* is not equal 1.0 only in longitudinal
* direction
*/
rmin *= gamma;
rmax *= gamma;
Vector_t tmp = Vector_t(std::max( std::abs(rmin[0]), std::abs(rmax[0]) ),
std::max( std::abs(rmin[1]), std::abs(rmax[1]) ),
std::max( std::abs(rmin[2]), std::abs(rmax[2]) )
......@@ -225,13 +268,23 @@ const double& AmrParticleBase<PLayout>::domainMapping(bool inverse) {
scale = std::max( tmp[0], tmp[1] );
scale = std::max( scale, tmp[2] );
} else {
// inverse Lorentz transform
gamma = 1.0 / gamma;
}
if ( std::isnan(scale) || std::isinf(scale) ) {
if ( !Ippl::myNode() )
throw IpplException("AmrParticleBase::domainMapping()",
"Scale factor is Nan or Inf");
}
Vector_t vscale = Vector_t(scale, scale, scale);
for (unsigned int i = 0; i < this->getLocalNum(); ++i)
this->R[i] /= vscale;
// Lorentz transform + mapping to [-1, 1]
for (unsigned int i = 0; i < this->getLocalNum(); ++i) {
this->R[i] = this->R[i] * gamma / vscale;
}
scale_m = 1.0 / scale;
......@@ -246,4 +299,70 @@ const double& AmrParticleBase<PLayout>::getScalingFactor() const {
return scale_m;
}
template<class PLayout>
void AmrParticleBase<PLayout>::setLorentzFactor(const Vector_t& lorentzFactor) {
lorentzFactor_m = lorentzFactor;
}
template<class PLayout>
void AmrParticleBase<PLayout>::getLocalBounds_m(Vector_t &rmin, Vector_t &rmax) {
const size_t localNum = this->getLocalNum();
if (localNum == 0) {
double max = 1e10;
rmin = Vector_t( max, max, max);
rmax = Vector_t(-max, -max, -max);
return;
}
rmin = this->R[0];
rmax = this->R[0];
for (size_t i = 1; i < localNum; ++ i) {
for (unsigned short d = 0; d < 3u; ++ d) {
if (rmin(d) > this->R[i](d)) rmin(d) = this->R[i](d);
else if (rmax(d) < this->R[i](d)) rmax(d) = this->R[i](d);
}
}
}
template<class PLayout>
void AmrParticleBase<PLayout>::getGlobalBounds_m(Vector_t &rmin, Vector_t &rmax) {
this->getLocalBounds_m(rmin, rmax);
double min[6];
for (unsigned int i = 0; i < 3; ++i) {
min[2*i] = rmin[i];
min[2*i + 1] = -rmax[i];
}
allreduce(min, 6, std::less<double>());
for (unsigned int i = 0; i < 3; ++i) {
rmin[i] = min[2*i];
rmax[i] = -min[2*i + 1];
}
}
// template<class PLayout>
// void AmrParticleBase<PLayout>::lorentzTransform(bool inverse) {
//
// if ( isLorentzTransformed_m && !inverse ) {
// return;
// }
//
// isLorentzTransformed_m = true;
//
// Vector_t gamma = lorentzFactor_m;
//
// if ( inverse ) {
// gamma = 1.0 / gamma;
// isLorentzTransformed_m = false;
// }
//
// for (std::size_t i = 0; i < this->getLocalNum(); ++i)
// this->R[i] *= gamma;
// }
#endif
......@@ -16,6 +16,10 @@ public:
return descr_.c_str();
}
virtual const std::string& where() const {
return meth_;
}
private:
std::string descr_;
......
#ifndef ABSTRACT_SOLVER_H
#define ABSTRACT_SOLVER_H
#include "AmrOpal.h"
#include "Amr/AmrDefs.h"
class AbstractSolver
{
public:
typedef std::unique_ptr<AmrOpal> amropal_p;
typedef amr::AmrFieldContainer_t AmrFieldContainer_t;
virtual ~AbstractSolver() {};
virtual
void solve(amropal_p& amropal,
AmrFieldContainer_t &rho,
AmrFieldContainer_t &phi,
AmrFieldContainer_t &efield,
unsigned short baseLevel,
unsigned short finestLevel,
bool prevAsGuess = true) = 0;
};
#endif
#include "AmrOpal.h"
#include "AmrOpal_F.h"
#include <AMReX_PlotFileUtil.H>
// #include <MultiFabUtil.H>
......@@ -9,7 +8,7 @@
// AmrOpal::AmrOpal() { }
AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Array<int>& n_cell_in, int coord,
AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Vector<int>& n_cell_in, int coord,
#ifdef IPPL_AMR
PartBunchAmr<amrplayout_t>* bunch)
#else
......@@ -26,7 +25,7 @@ AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Array<
nCharge_m(1.0e-15),
minNumPart_m(1),
maxNumPart_m(1),
meshScaling_m(Vector_t(1.0, 1.0, 1.0))
meshScaling_m(Vector_t(D_DECL(1.0, 1.0, 1.0)))
{
initBaseLevel();
......@@ -35,14 +34,14 @@ AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Array<
nChargePerCell_m[0]->setVal(0.0, 1);
}
AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Array<int>& n_cell_in, int coord)
AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Vector<int>& n_cell_in, int coord)
: AmrMesh(rb, max_level_in, n_cell_in, coord),
tagging_m(kChargeDensity),
scaling_m(0.75),
nCharge_m(1.0e-15),
minNumPart_m(1),
maxNumPart_m(1),
meshScaling_m(Vector_t(1.0, 1.0, 1.0))
meshScaling_m(Vector_t(D_DECL(1.0, 1.0, 1.0)))
{
finest_level = 0;
......@@ -55,7 +54,7 @@ AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Array<
MakeNewLevel(0, 0.0, ba, dm);
}
AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Array<int>& n_cell_in, int coord,
AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Vector<int>& n_cell_in, int coord,
const std::vector<int>& refratio)
: AmrMesh(rb, max_level_in, n_cell_in, coord, refratio),
tagging_m(kChargeDensity),
......@@ -63,7 +62,7 @@ AmrOpal::AmrOpal(const amrex::RealBox* rb, int max_level_in, const amrex::Array<
nCharge_m(1.0e-15),
minNumPart_m(1),
maxNumPart_m(1),
meshScaling_m(Vector_t(1.0, 1.0, 1.0))
meshScaling_m(Vector_t(D_DECL(1.0, 1.0, 1.0)))
{
finest_level = 0;
......@@ -277,16 +276,16 @@ void AmrOpal::writePlotFile(std::string filename, int step) {
chargeOnGrid[i]->setVal(0.0);
}
amrex::Array<std::string> varnames(1, "rho");
amrex::Vector<std::string> varnames(1, "rho");
amrex::Array<const amrex::MultiFab*> tmp(finest_level + 1);
amrex::Vector<const amrex::MultiFab*> tmp(finest_level + 1);
for (/*unsigned*/ int i = 0; i < finest_level + 1; ++i) {
tmp[i] = chargeOnGrid[i].get();
}
const auto& mf = tmp;
amrex::Array<int> istep(finest_level+1, step);
amrex::Vector<int> istep(finest_level+1, step);
amrex::WriteMultiLevelPlotfile(filename, finest_level + 1, mf, varnames,
Geom(), 0.0, istep, refRatio());
......@@ -345,7 +344,7 @@ void
AmrOpal::regrid (int lbase, amrex::Real time)
{
int new_finest = 0;
amrex::Array<amrex::BoxArray> new_grids(finest_level+2);
amrex::Vector<amrex::BoxArray> new_grids(finest_level+2);
MakeNewGrids(lbase, time, new_finest, new_grids);
......@@ -428,7 +427,7 @@ AmrOpal::MakeNewLevel (int lev, amrex::Real time,
SetBoxArray(lev, new_grids);
SetDistributionMap(lev, new_dmap);
nChargePerCell_m[lev] = std::unique_ptr<amrex::MultiFab>(new amrex::MultiFab(new_grids, dmap[lev], 1, 1));
nChargePerCell_m[lev] = std::unique_ptr<amrex::MultiFab>(new amrex::MultiFab(new_grids, new_dmap, 1, 1));
}
void AmrOpal::ClearLevel(int lev) {
......@@ -480,30 +479,31 @@ void AmrOpal::tagForChargeDensity_m(int lev, amrex::TagBoxArray& tags, amrex::Re
#pragma omp parallel
#endif
{
amrex::Array<int> itags;
for (amrex::MFIter mfi(*nChargePerCell_m[lev],false/*true*/); mfi.isValid(); ++mfi) {
for (amrex::MFIter mfi(*nChargePerCell_m[lev], false/*true*/); mfi.isValid(); ++mfi) {
const amrex::Box& tilebx = mfi.validbox();//mfi.tilebox();
amrex::TagBox& tagfab = tags[mfi];
amrex::FArrayBox& fab = (*nChargePerCell_m[lev])[mfi];
// We cannot pass tagfab to Fortran becuase it is BaseFab<char>.
// So we are going to get a temporary integer array.
tagfab.get_itags(itags, tilebx);
// data pointer and index space
int* tptr = itags.dataPtr();
const int* tlo = tilebx.loVect();
const int* thi = tilebx.hiVect();
state_error(tptr, ARLIM_3D(tlo), ARLIM_3D(thi),
BL_TO_FORTRAN_3D((*nChargePerCell_m[lev])[mfi]),
&tagval, &clearval,
ARLIM_3D(tilebx.loVect()), ARLIM_3D(tilebx.hiVect()),
ZFILL(dx), ZFILL(prob_lo), &scale, &nCharge_m);
//
// Now update the tags in the TagBox.
//
tagfab.tags_and_untags(itags, tilebx);
for (int i = tlo[0]; i <= thi[0]; ++i) {
for (int j = tlo[1]; j <= thi[1]; ++j) {
#if AMREX_SPACEDIM == 3
for (int k = tlo[2]; k <= thi[2]; ++k) {
#endif
amrex::IntVect iv(D_DECL(i,j,k));
if ( std::abs( fab(iv) ) >= nCharge_m )
tagfab(iv) = tagval;
else
tagfab(iv) = clearval;
#if AMREX_SPACEDIM == 3
}
#endif
}
}
}
}
}
......@@ -596,30 +596,31 @@ void AmrOpal::tagForPotentialStrength_m(int lev, amrex::TagBoxArray& tags, amrex
#pragma omp parallel
#endif
{
amrex::Array<int> itags;
for (amrex::MFIter mfi(*phi[lev],false/*true*/); mfi.isValid(); ++mfi) {
for (amrex::MFIter mfi(*phi[lev], false/*true*/); mfi.isValid(); ++mfi) {
const amrex::Box& tilebx = mfi.validbox();//mfi.tilebox();
amrex::TagBox& tagfab = tags[mfi];
amrex::FArrayBox& fab = (*phi[lev])[mfi];
// We cannot pass tagfab to Fortran becuase it is BaseFab<char>.
// So we are going to get a temporary integer array.
tagfab.get_itags(itags, tilebx);
// data pointer and index space
int* tptr = itags.dataPtr();
const int* tlo = tilebx.loVect();
const int* thi = tilebx.hiVect();
tag_potential_strength(tptr, ARLIM_3D(tlo), ARLIM_3D(thi),
BL_TO_FORTRAN_3D((*phi[lev])[mfi]),
&tagval, &clearval,
ARLIM_3D(tilebx.loVect()), ARLIM_3D(tilebx.hiVect()),
ZFILL(dx), ZFILL(prob_lo), &scale, &value);
//
// Now update the tags in the TagBox.
//
tagfab.tags_and_untags(itags, tilebx);
for (int i = tlo[0]; i <= thi[0]; ++i) {
for (int j = tlo[1]; j <= thi[1]; ++j) {
#if AMREX_SPACEDIM == 3
for (int k = tlo[2]; k <= thi[2]; ++k) {
#endif
amrex::IntVect iv(D_DECL(i,j,k));
if ( std::abs( fab(iv) ) >= value )
tagfab(iv) = tagval;
else
tagfab(iv) = clearval;
#if AMREX_SPACEDIM == 3
}
#endif
}
}
}
}
......@@ -701,14 +702,14 @@ void AmrOpal::tagForEfieldStrength_m(int lev, amrex::TagBoxArray& tags, amrex::R
const int clearval = amrex::TagBox::CLEAR;
const int tagval = amrex::TagBox::SET;
amrex::Real min[3] = {0.0, 0.0, 0.0};
amrex::Real max[3] = {0.0, 0.0, 0.0};
for (int i = 0; i < 3; ++i) {
amrex::Real min[AMREX_SPACEDIM] = { D_DECL(0.0, 0.0, 0.0) };
amrex::Real max[AMREX_SPACEDIM] = { D_DECL(0.0, 0.0, 0.0) };
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
max[i] = scaling_m * grad_phi[lev]->max(i);
min[i] = scaling_m * grad_phi[lev]->min(i);
}
amrex::Real efield[3] = {0.0, 0.0, 0.0};
for (int i = 0; i < 3; ++i)
amrex::Real efield[AMREX_SPACEDIM] = { D_DECL(0.0, 0.0, 0.0) };
for (int i = 0; i < AMREX_SPACEDIM; ++i)
efield[i] = std::max( std::abs(min[i]), std::abs(max[i]) );
......@@ -733,17 +734,23 @@ void AmrOpal::tagForEfieldStrength_m(int lev, amrex::TagBoxArray& tags, amrex::R
for (int i = tlo[0]; i <= thi[0]; ++i) {
for (int j = tlo[1]; j <= thi[1]; ++j) {
#if AMREX_SPACEDIM == 3
for (int k = tlo[2]; k <= thi[2]; ++k) {
amrex::IntVect iv(i,j,k);
#endif
amrex::IntVect iv(D_DECL(i,j,k));
if (std::abs(fab(iv, 0)) >= efield[0])
tagfab(iv) = tagval;
else if (std::abs(fab(iv, 1)) >= efield[1])
tagfab(iv) = tagval;
#if AMREX_SPACEDIM == 3
else if (std::abs(fab(iv, 2)) >= efield[2])
tagfab(iv) = tagval;
#endif
else
tagfab(iv) = clearval;
#if AMREX_SPACEDIM == 3
}
#endif
}
}
......@@ -774,7 +781,7 @@ void AmrOpal::tagForMomentum_m(int lev, amrex::TagBoxArray& tags, amrex::Real sc
size_t lBegin = LocalNumPerLevel.begin(lev);
size_t lEnd = LocalNumPerLevel.end(lev);
Vector_t pmax = Vector_t(0.0, 0.0, 0.0);
Vector_t pmax = Vector_t(D_DECL(0.0, 0.0, 0.0));
for (size_t i = lBegin; i < lEnd; ++i) {
const Vector_t& tmp = bunch_m->P[i];
pmax = ( dot(tmp, tmp) > dot(pmax, pmax) ) ? tmp : pmax;
......@@ -810,13 +817,17 @@ void AmrOpal::tagForMomentum_m(int lev, amrex::TagBoxArray& tags, amrex::Real sc
for (int i = tlo[0]; i <= thi[0]; ++i) {
for (int j = tlo[1]; j <= thi[1]; ++j) {
#if AMREX_SPACEDIM == 3
for (int k = tlo[2]; k <= thi[2]; ++k) {
amrex::IntVect iv(i, j, k);
#endif
amrex::IntVect iv(D_DECL(i, j, k));
if ( cells.find(iv) != cells.end() )
tagfab(iv) = tagval;
else
tagfab(iv) = clearval;
#if AMREX_SPACEDIM == 3
}
#endif