Commit 527134a5 authored by snuverink_j's avatar snuverink_j
Browse files

fix unused argument warnings in IPPL tests; fix spacing

parent 76c9e561
// -*- C++ -*-
/***************************************************************************
*
* The IPPL Framework
*
* This program was prepared by PSI.
* All rights in the program are reserved by PSI.
* Neither PSI nor the author(s)
* makes any warranty, express or implied, or assumes any liability or
* responsibility for the use of this software
*
* Visit http://www.acl.lanl.gov/POOMS for more details
*
***************************************************************************/
// -*- C++ -*-
/***************************************************************************
*
* The IPPL Framework
*
*
* Visit http://people.web.psi.ch/adelmann/ for more details
*
***************************************************************************/
// TestFFT.cpp , Tim Williams 1/27/1997
// Updated by Julian Cummings, 3/31/98
......@@ -40,10 +15,10 @@ enum InterPolT {NGP,CIC};
enum BCT {OOO,OOP,PPP,DDD,DDO,DDP}; // OOO == all dim, open BC
enum TestCases {test1};
bool Configure(int argc, char *argv[], InterPolT *interPol,
unsigned int *nx, unsigned int *ny, unsigned int *nz,
TestCases *test2Perform,
int *serialDim, unsigned int *processes, unsigned int *nLoop)
bool Configure(int argc, char *argv[], InterPolT *interPol,
unsigned int *nx, unsigned int *ny, unsigned int *nz,
TestCases *test2Perform,
int *serialDim, unsigned int *processes, unsigned int *nLoop)
{
Inform msg("Configure ");
......@@ -71,7 +46,7 @@ bool Configure(int argc, char *argv[], InterPolT *interPol,
*nLoop = atoi(argv[++i]);
} else if (s == "-Decomp") {
*serialDim = atoi(argv[++i]);
}
}
else {
errmsg << "Illegal format for or unknown option '" << s.c_str() << "'.";
errmsg << endl;
......@@ -113,25 +88,25 @@ int main(int argc, char *argv[])
const unsigned D=3U;
testmsg << "Dimensionality: D= " << D << " P= " << processes;
testmsg << " nx= " << nx << " ny= " << ny << " nz= " << nz << endl;
unsigned ngrid[D]; // grid sizes
// Used in evaluating correctness of results:
double realDiff;
double realDiff = -1.0;
// Various counters, constants, etc:
unsigned int d;
/*int tag = */ Ippl::Comm->next_tag(IPPL_APP_TAG0);
double pi = acos(-1.0);
double twopi = 2.0*pi;
// Timer:
Timer timer;
// Layout information:
unsigned vnodes; // number of vnodes; input at cmd line
e_dim_tag allParallel[D]; // Specifies SERIAL, PARALLEL dims
for (d=0; d<D; d++)
for (d=0; d<D; d++)
allParallel[d] = PARALLEL;
allParallel[0] = SERIAL;
......@@ -147,18 +122,18 @@ int main(int argc, char *argv[])
testmsg << "In-place CC transform using all-parallel layout ..." << endl;
// ================BEGIN INTERACTION LOOP====================================
//------------------------complex<-->complex-------------------------------
// Complex test Fields
// create standard domain
NDIndex<D> ndiStandard;
for (d=0; d<D; d++)
for (d=0; d<D; d++)
ndiStandard[d] = Index(ngrid[d]);
// create new domain with axes permuted to match FFT output
NDIndex<D> ndiPermuted;
ndiPermuted[0] = ndiStandard[D-1];
for (d=1; d<D; d++)
for (d=1; d<D; d++)
ndiPermuted[d] = ndiStandard[d-1];
// create half-size domain for RC transform along zeroth axis
......@@ -167,7 +142,7 @@ int main(int argc, char *argv[])
// create new domain with axes permuted to match FFT output
NDIndex<D> ndiPermuted0h;
ndiPermuted0h[0] = ndiStandard0h[D-1];
for (d=1; d<D; d++)
for (d=1; d<D; d++)
ndiPermuted0h[d] = ndiStandard0h[d-1];
// create half-size domain for sine transform along zeroth axis
......@@ -177,21 +152,21 @@ int main(int argc, char *argv[])
// create new domain with axes permuted to match FFT output
NDIndex<D> ndiPermuted1h;
ndiPermuted1h[0] = ndiStandard1h[D-1];
for (d=1; d<D; d++)
for (d=1; d<D; d++)
ndiPermuted1h[d] = ndiStandard1h[d-1];
// all parallel layout, standard domain, normal axis order
FieldLayout<D> layoutPPStan(ndiStandard,allParallel,vnodes);
FieldLayout<D> layoutPPStan0h(ndiStandard0h,allParallel,vnodes);
FieldLayout<D> layoutPPStan1h(ndiStandard1h,allParallel,vnodes);
// create test Fields for complex-to-complex FFT
BareField<std::complex<double>,D> CFieldPPStan(layoutPPStan);
BareField<std::complex<double>,D> CFieldPPStan_save(layoutPPStan);
BareField<double,D> diffFieldPPStan(layoutPPStan);
// For calling FieldDebug functions from debugger, set up output format:
setFormat(4,3);
......@@ -205,22 +180,22 @@ int main(int argc, char *argv[])
yfact = 2.0*twopi/(ngrid[1]);
zfact = 2.0*twopi/(ngrid[2]);
kx = 1.0; ky = 2.0; kz = 3.0; // wavenumbers
CFieldPPStan[ndiStandard[0]][ndiStandard[1]][ndiStandard[2]] =
CFieldPPStan[ndiStandard[0]][ndiStandard[1]][ndiStandard[2]] =
sfact * ( sin( (ndiStandard[0]+1) * kx * xfact +
ndiStandard[1] * ky * yfact +
ndiStandard[2] * kz * zfact ) +
sin( (ndiStandard[0]+1) * kx * xfact -
sin( (ndiStandard[0]+1) * kx * xfact -
ndiStandard[1] * ky * yfact -
ndiStandard[2] * kz * zfact ) ) +
ndiStandard[2] * kz * zfact ) ) +
cfact * (-cos( (ndiStandard[0]+1) * kx * xfact +
ndiStandard[1] * ky * yfact +
ndiStandard[2] * kz * zfact ) +
cos( (ndiStandard[0]+1) * kx * xfact -
ndiStandard[2] * kz * zfact ) +
cos( (ndiStandard[0]+1) * kx * xfact -
ndiStandard[1] * ky * yfact -
ndiStandard[2] * kz * zfact ) );
CFieldPPStan_save = CFieldPPStan; // Save values for checking later
// CC FFT tests
// Instantiate complex<-->complex FFT object
// Transform in all directions
......@@ -246,16 +221,4 @@ int main(int argc, char *argv[])
testmsg << " nx= " << nx << " ny= " << ny << " nz= " << nz;
testmsg << " ||d||= " << fabs(realDiff) << endl;
return 0;
}
/***************************************************************************
* $RCSfile: TestFFT-SPP.cpp,v $ $Author: adelmann $
* $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:36 $
***************************************************************************/
/***************************************************************************
* $RCSfile: addheaderfooter,v $ $Author: adelmann $
* $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:17 $
* IPPL_VERSION_ID: $Id: addheaderfooter,v 1.1.1.1 2003/01/23 07:40:17 adelmann Exp $
***************************************************************************/
}
\ No newline at end of file
// -*- C++ -*-
/***************************************************************************
*
* The IPPL Framework
*
* This program was prepared by PSI.
* All rights in the program are reserved by PSI.
* Neither PSI nor the author(s)
* makes any warranty, express or implied, or assumes any liability or
* responsibility for the use of this software
*
* Visit http://www.acl.lanl.gov/POOMS for more details
*
***************************************************************************/
// -*- C++ -*-
/***************************************************************************
*
* The IPPL Framework
*
*
* Visit http://people.web.psi.ch/adelmann/ for more details
*
***************************************************************************/
// TestFFT.cpp , Tim Williams 1/27/1997
// Updated by Julian Cummings, 3/31/98
......@@ -35,13 +10,13 @@ using namespace std;
int main(int argc, char *argv[])
{
Ippl ippl(argc,argv);
Inform testmsg(NULL,0);
const unsigned D=3U;
testmsg << "%%%%%%% Dimensionality: D = " << D << " %%%%%%%" << endl;
unsigned ngrid[D]; // grid sizes
// Used in evaluating correctness of results:
......@@ -76,21 +51,21 @@ int main(int argc, char *argv[])
if( Ippl::Comm->myNode() == Parent ) {
bool vnodesOK = false;
while (!vnodesOK) {
compressTemps = false;
constInput = false;
for (d=0; d<D; d++)
ngrid[d] = 32;
// now broadcast data to other nodes
Message *mess = new Message();
putMessage( *mess, vnodes );
if (vnodes > 0) {
putMessage(*mess, compressTemps);
putMessage(*mess, constInput);
for (d=0; d<D; d++)
putMessage( *mess, ngrid[d] );
}
Ippl::Comm->broadcast_all(mess, tag);
compressTemps = false;
constInput = false;
for (d=0; d<D; d++)
ngrid[d] = 32;
// now broadcast data to other nodes
Message *mess = new Message();
putMessage( *mess, vnodes );
if (vnodes > 0) {
putMessage(*mess, compressTemps);
putMessage(*mess, constInput);
for (d=0; d<D; d++)
putMessage( *mess, ngrid[d] );
}
Ippl::Comm->broadcast_all(mess, tag);
}
}
// now each node recieves the data
......@@ -101,13 +76,13 @@ int main(int argc, char *argv[])
if (vnodes <= 0) break;
getMessage(*mess, compressTemps);
getMessage(*mess, constInput);
for (d=0; d<D; d++)
for (d=0; d<D; d++)
getMessage( *mess, ngrid[d] );
delete mess;
//------------------------complex<-->complex-------------------------------
// Complex test Fields
// create standard domain
NDIndex<D> ndiStandard;
for (d=0; d<D; d++)
......@@ -198,27 +173,27 @@ int main(int argc, char *argv[])
// Rather more complete test functions (sine or cosine mode):
std::complex<double> sfact(1.0,0.0); // (1,0) for sine mode; (0,0) for cosine mode
std::complex<double> cfact(0.0,0.0); // (0,0) for sine mode; (1,0) for cosine mode
double xfact, kx, yfact, ky, zfact, kz;
xfact = pi/(ngrid[0] + 1.0);
yfact = 2.0*twopi/(ngrid[1]);
zfact = 2.0*twopi/(ngrid[2]);
kx = 1.0; ky = 2.0; kz = 3.0; // wavenumbers
CFieldPPStan[ndiStandard[0]][ndiStandard[1]][ndiStandard[2]] =
CFieldPPStan[ndiStandard[0]][ndiStandard[1]][ndiStandard[2]] =
sfact * ( sin( (ndiStandard[0]+1) * kx * xfact +
ndiStandard[1] * ky * yfact +
ndiStandard[2] * kz * zfact ) +
sin( (ndiStandard[0]+1) * kx * xfact -
ndiStandard[1] * ky * yfact -
ndiStandard[2] * kz * zfact ) ) +
ndiStandard[1] * ky * yfact +
ndiStandard[2] * kz * zfact ) +
sin( (ndiStandard[0]+1) * kx * xfact -
ndiStandard[1] * ky * yfact -
ndiStandard[2] * kz * zfact ) ) +
cfact * (-cos( (ndiStandard[0]+1) * kx * xfact +
ndiStandard[1] * ky * yfact +
ndiStandard[2] * kz * zfact ) +
cos( (ndiStandard[0]+1) * kx * xfact -
ndiStandard[1] * ky * yfact -
ndiStandard[2] * kz * zfact ) );
ndiStandard[2] * kz * zfact ) +
cos( (ndiStandard[0]+1) * kx * xfact -
ndiStandard[1] * ky * yfact -
ndiStandard[2] * kz * zfact ) );
CFieldPPStan_save = CFieldPPStan; // Save values for checking later
// RC FFT tests
RFieldPPStan = real(CFieldPPStan_save);
CFieldPPStan0h = std::complex<double>(0.0,0.0);
......@@ -236,7 +211,7 @@ int main(int argc, char *argv[])
timer.start();
// Test real<-->complex transform (simple test: forward then inverse
// transform, see if we get back original values.
cout << "PE " << pe << " about to invoke forward FFT::transform()" << endl;
testmsg << "Forward transform ..." << endl;
rcfft.transform("forward", RFieldPPStan, CFieldPPStan0h, constInput);
......@@ -260,14 +235,14 @@ int main(int argc, char *argv[])
timer.start();
// Test real<-->complex transform (simple test: forward then inverse
// transform, see if we get back original values.
testmsg << "Forward transform ..." << endl;
rcfft.transform("forward", RFieldSPStan, CFieldSPStan0h, constInput);
testmsg << "Inverse transform ..." << endl;
rcfft.transform("inverse", CFieldSPStan0h, RFieldSPStan, constInput);
timer.stop();
diffFieldSPStan = Abs(RFieldSPStan - RFieldSPStan_save);
......@@ -286,13 +261,13 @@ int main(int argc, char *argv[])
timer.start();
// Test real<-->complex transform (simple test: forward then inverse
// transform, see if we get back original values.
testmsg << "Forward transform ..." << endl;
rcfft.transform("forward", RFieldSPStan, CFieldSPPerm0h, constInput);
testmsg << "Inverse transform ..." << endl;
rcfft.transform("inverse", CFieldSPPerm0h, RFieldSPStan, constInput);
timer.stop();
diffFieldSPStan = Abs(RFieldSPStan - RFieldSPStan_save);
......@@ -303,17 +278,7 @@ int main(int argc, char *argv[])
}
testmsg << "CPU time used = " << timer.cpu_time() << " secs." << endl;
timer.clear();
}
}
testmsg << "test is correct: " << correct << endl;
return 0;
}
/***************************************************************************
* $RCSfile: TestRC.cpp,v $ $Author: adelmann $
* $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:36 $
***************************************************************************/
/***************************************************************************
* $RCSfile: addheaderfooter,v $ $Author: adelmann $
* $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:17 $
* IPPL_VERSION_ID: $Id: addheaderfooter,v 1.1.1.1 2003/01/23 07:40:17 adelmann Exp $
***************************************************************************/
}
\ No newline at end of file
......@@ -5,221 +5,221 @@
template<typename Particles>
void createParticleDistribution(Particles & P, std::string distribution, unsigned count, double qi, Vektor<double,3> extend_l,Vektor<double,3> extend_r, Vektor<double,3> source, double sphere_radius=1.,unsigned n_dummies=0 ) {
enum DistType {UNIFORM, RANDOM, EVEN, MANUAL, LINE,TWOSTREAM, LANDAU};
DistType type = UNIFORM;
if(distribution == std::string("uniform"))
{
type = UNIFORM;
}
else if(distribution == std::string("random"))
{
type = RANDOM;
}
else if(distribution == std::string("even"))
{
type = EVEN;
}
else if(distribution == std::string("manual"))
{
type = MANUAL;
}
else if(distribution == std::string("line"))
{
type = LINE;
}
else if(distribution == std::string("twostream"))
{
type = TWOSTREAM;
}
else if(distribution == std::string("landau"))
{
type = LANDAU;
}
else {
std::cout << "Unrecognized distribution " << distribution << std::endl;
}
if (P->singleInitNode())
{
size_t index = 0;
if (type == UNIFORM) {
std::cout << "SAMPLING UNIFORM DISTRIBUTION" << std::endl;
double density = count*3./(sphere_radius*sphere_radius*sphere_radius*4*3.14159);
double pdist = std::pow(1./density,1/3.);
//SET TO FIXED DIST:
pdist = 4.*sphere_radius/15/4;
for (double x=-2.*sphere_radius; x<2*sphere_radius; x+=pdist){
for (double y=-2.*sphere_radius; y<2*sphere_radius; y+=pdist){
for (double z=-2.*sphere_radius; z<2*sphere_radius; z+=pdist){
Vektor<double, 3> pos(x,y,z);
//msg << "checking particle pos = " << pos << endl;
if (dot(source-pos,source-pos)<sphere_radius*sphere_radius){
//msg << " Particle accepted " << endl;
P->create(1);
P->Q[index]=qi;
P->R[index++]=pos;
}
}
}
}
if (n_dummies>0) {
n_dummies=0;
for (double x=-2.*sphere_radius; x<2*sphere_radius; x+=pdist){
for (double y=-2.*sphere_radius; y<2*sphere_radius; y+=pdist){
for (double z=-2.*sphere_radius; z<2*sphere_radius; z+=pdist){
Vektor<double, 3> pos(x,y,z);
//msg << "checking particle pos = " << pos << endl;
if (dot(source-pos,source-pos)>sphere_radius*sphere_radius){
//msg << " Particle accepted " << endl;
P->create(1);
P->Q[index]=0;
P->R[index++]=pos;
n_dummies++;
}
}
}
}
}
}
if (type == RANDOM) {
std::cout << "SAMPLING RANDOM DISTRIBUTION" << std::endl;
//create n particles within a sphere of radius R around source and charge 1/n
P->create(count);
std::default_random_engine generator;
//uniform distributed between 0 and 1
std::uniform_real_distribution<double> unidistribution(0,1);
//std::uniform_real_distribution<double> unidistribution(source[0]-sphere_radius, source[0]+sphere_radius);
auto uni = std::bind(unidistribution, generator);
//normal distribution
std::normal_distribution<double> normdistribution(0,1.0);
auto normal = std::bind(normdistribution, generator);
//use rejection method to place particles uniformly in sphere
double density = count*3./(sphere_radius*sphere_radius*sphere_radius*4*3.14159);
double pdist = std::pow(1./density,1/3.);
for (unsigned i = 0; i<count; ++i) {
Vektor<double, 3> X(normal(),normal(),normal());
double U = uni();
Vektor<double, 3> pos = source + sphere_radius*std::pow(U,1./3.)/std::sqrt(dot(X,X))*X;
//BEGIN reject close particles
bool reject = 0;
for (unsigned k=0; k<index; ++k) {
if (sqrt(dot(P->R[k]-pos,P->R[k]-pos))<0.5*pdist)
reject = 1;
}
if (!reject){
Vektor<double, 3> vel(1.,1.,1.); //normal velocity distribution in z dir
//Vektor<double, 3> vel(0,0,0); //normal velocity distribution in z dir
P->Q[index]=-qi;
P->m[index]=1;
P->total_charge+=qi;
P->v[index]=vel;
P->R[index++]=pos;
}
else {
i--;
}
}
//place particles with 0 charge outside the sphere
P->create(n_dummies);
std::uniform_real_distribution<double> o_unidistribution(-2.*sphere_radius,2.*sphere_radius);
auto o_uni = std::bind(o_unidistribution, generator);
for (unsigned i = 0; i<n_dummies; ++i) {
//Vektor<double, Dim> X(normal(),normal(),normal());
bool check_out_of_sphere = 1;
while(check_out_of_sphere){
Vektor<double, 3> pos(o_uni(),o_uni(),o_uni());
if (dot(source-pos,source-pos)>sphere_radius*sphere_radius){
std::cout << "dummy placed " << std::endl;
//P->Q[index]=1./count;
P->Q[index]=0;
P->R[index++]=pos;
check_out_of_sphere=0;
}
}
}
}
if (type == EVEN) {
std::cout << "SAMPLING EVEN DISTRIBUTION" << std::endl;
P->create(count);
std::default_random_engine generator;
//uniform distributed between 0 and 1
std::uniform_real_distribution<double> unidistribution(extend_l[0],extend_r[0]);
std::normal_distribution<double> normdistribution(0,1.0);
auto normal = std::bind(normdistribution, generator);
//std::uniform_real_distribution<double> unidistribution(source[0]-sphere_radius, source[0]+sphere_radius);
auto uni = std::bind(unidistribution, generator);
for (unsigned i = 0; i<count; ++i) {
Vektor<double, 3> X(uni(),uni(),uni());
Vektor<double, 3> pos = X;
Vektor<double, 3> vel(normal(),normal(),normal());
//Vektor<double, 3> vel(0,0,0);
P->Q[index]=qi;
P->R[index++]=pos;
}
}
if (type == MANUAL) {
std::cout << "SAMPLING MANUAL DISTRIBUTION" << std::endl;
for (unsigned i = 0; i<count; ++i) {
std::cout << "please give coordinates of particle " << i << std::endl;
double x,y,z;
std::cin >> x; std::cin >> y; std::cin >>z;
Vektor<double, 3> pos(0,0,0);
pos[0]=x; pos[1]=y; pos[2]=z;
std::cout << "pos = " << pos << std::endl;
P->create(1);
P->Q[index]=qi;
Vektor<double, 3> vel(0,0,0.1); //normal velocity distribution in z dir
P->v[index]=vel;
P->R[index++]=pos;
}
}
if (type == LINE) {
std::cout << "SAMPLING LINE DISTRIBUTION" << std::endl;
P->create(count);
std::default_random_engine generator;
//uniform distributed between 0 and 1
std::uniform_real_distribution<double> unidistribution(extend_l[0],extend_r[0]);
std::normal_distribution<double> normaldist(1,.5);
auto uni = std::bind(unidistribution, generator);
auto normal = std::bind(normaldist, generator);
for (unsigned i = 0; i<count; ++i) {
Vektor<double, 3> X(0,0,uni()); //place all particles on z axis