Skip to content

GitLab

S
src
Commit e16c2223 authored by Zheng Dawei's avatar Zheng Dawei
Browse files
Options

update efficient Poisson solver for the first phrase, to be debugged and tested with examples next

parent de0c7799
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 2586 additions and 521 deletions
ippl/src/FFT/FFT.h
View file @ e16c2223
This diff is collapsed.
ippl/src/FFT/FFT.hpp
View file @ e16c2223
This diff is collapsed.
ippl/src/FFT/FFTBase.h
View file @ e16c2223
......@@ -38,12 +38,13 @@ std::ostream& operator<<(std::ostream&, const FFTBase<Dim,T>&);
inline
std::string getTransformType(unsigned int i)
{
static const char* transformTypeString_g[4] = { "complex-to-complex FFT",
static const char* transformTypeString_g[5] = { "complex-to-complex FFT",
"real-to-complex FFT",
"sine transform",
"cosine transform" };
"cosine transform",
"implicitly zeo paded complex-to-complex FFT" };
return std::string(transformTypeString_g[i % 4]);
return std::string(transformTypeString_g[i % 5]);
}
/*!
......@@ -64,8 +65,7 @@ public:
typedef NDIndex<Dim> Domain_t; // domain type
// Enumeration of transform types, used by derived FFT classes
enum FFT_e { ccFFT, rcFFT, sineFFT, cosineFFT };
enum FFT_e { ccFFT, rcFFT, sineFFT, cosineFFT, impFFT};
// Type used for performing 1D FFTs
#if defined(IPPL_USE_SCSL_FFT)
typedef SCSL<T> InternalFFT_t;
......
ippl/src/FFT/fftpack.F
View file @ e16c2223
......@@ -21,7 +21,64 @@ c Cray machines don't support double precision type. Use real.
#define FLOAT DOUBLE PRECISION
#endif
c---------------------------DOUBLE PRECISION----------------------------
subroutine impcffti(n, wsave, zeta)
c----------------------------TJW
implicit FLOAT (a-h,o-z)
implicit integer (i-n)
c----------------------------TJW
dimension wsave(*), zeta(*)
data pi /3.14159265358979/
if (n .eq. 1) return
iw1 = n+n+1
iw2 = iw1+n+n
call cffti1 (n,wsave(iw1),wsave(iw2))
theta=pi/dble(n)
do i=0,n-1
omega = dble(i)*theta
zeta(2*i+1)= cos(omega)
zeta(2*i+2)= -sin(omega)
end do
return
end
subroutine impcfftb(n, c, d, wsave, zeta)
c----------------------------TJW
implicit FLOAT (a-h,o-z)
implicit integer (i-n)
c----------------------------TJW
dimension c(*), d (*), wsave(*), zeta(*)
call cfftb (n,c,wsave)
call cfftb (n,d,wsave)
do i = 1, n
c(2*i-1)=c(2*i-1)+d(2*i-1)*zeta(2*i-1)+d(2*i)*zeta(2*i)
c(2*i) =c(2*i)+d(2*i)*zeta(2*i-1)-d(2*i-1)*zeta(2*i)
end do
return
end
subroutine impcfftf(n, c, d, wsave, zeta)
c----------------------------TJW
implicit FLOAT (a-h,o-z)
implicit integer (i-n)
c----------------------------TJW
dimension c(*), d (*), wsave(*), zeta(*)
do i = 1,n
d(2*i-1)=c(2*i-1)*zeta(2*i-1)-c(2*i)*zeta(2*i)
d(2*i) =c(2*i)*zeta(2*i-1)+c(2*i-1)*zeta(2*i)
end do
call cfftf (n,c,wsave)
call cfftf (n,d,wsave)
return
end
subroutine rffti (n,wsave)
c----------------------------TJW
implicit FLOAT (a-h,o-z)
......
ippl/src/FFT/fftpack_FFT.h
View file @ e16c2223
......@@ -36,6 +36,9 @@
#define sint_ SINT
#define costi_ COSTI
#define cost_ COST
#define impcffti_ IMPCFFTI
#define impcfftf_ IMPCFFTF
#define impcfftb_ IMPCFFTB
#define fcffti_ FCFFTI
#define fcfftf_ FCFFTF
#define fcfftb_ FCFFTB
......@@ -60,6 +63,9 @@
#define sint_ sint
#define costi_ costi
#define cost_ cost
#define impcffti_ impcffti
#define impcfftf_ impcfftf
#define impcfftb_ impcfftb
#define fcffti_ fcffti
#define fcfftf_ fcfftf
#define fcfftb_ fcfftb
......@@ -87,6 +93,10 @@ extern "C" {
// double-precision cosine transform
void costi_(int& n, double& wsave);
void cost_(int& n, double& r, double& wsave);
//double-precision implicitly zero padded CC FFT
void impcffti_(int& n, double& wsave, double& zeta);
void impcfftf_(int& n, double& r, double& s, double& wsave, double& zeta);
void impcfftb_(int& n, double& r, double& s, double& wsave, double& zeta);
// single-precision CC FFT
void fcffti_(int& n, float& wsave);
void fcfftf_(int& n, float& r, float& wsave);
......@@ -147,6 +157,9 @@ public:
static void rrffti(int n, double* wsave) { sinti_(n, *wsave); }
static void rrsinti(int n, double* wsave) { sinti_(n, *wsave); }
static void rrcosti(int n, double* wsave) { costi_(n, *wsave); }
static void impcffti(int n, double* wsave, double* zeta) {
impcffti_(n,*wsave,*zeta);
}
// forward and backward CC FFT
static void ccfftf(int n, double* r, double* wsave) {cfftf_(n, *r, *wsave);}
static void ccfftb(int n, double* r, double* wsave) {cfftb_(n, *r, *wsave);}
......@@ -159,7 +172,10 @@ public:
static void rrsint(int n, double* r, double* wsave) { sint_(n, *r, *wsave); }
// cosine transform
static void rrcost(int n, double* r, double* wsave) { cost_(n, *r, *wsave); }
// imp CC FFT
static void impcfftf(int n, double *r, double *s, double* wsave, double *zeta) { impcfftf_(n,*r, *s, *wsave, *zeta); }
static void impcfftb(int n, double *r, double *s, double* wsave, double *zeta) { impcfftb_(n,*r, *s, *wsave, *zeta); }
};
......@@ -191,6 +207,11 @@ public:
// invoke FFT on complex data for given dimension and direction
void callFFT(unsigned transformDim, int direction, Complex_t* data);
// invoke FFT on complex data for given dimension and direction
void callFFT(unsigned transformDim, int direction, Complex_t* data,
Complex_t* data2);
// invoke FFT on real data for given dimension and direction
void callFFT(unsigned transformDim, int direction, T* data);
......@@ -200,6 +221,7 @@ private:
int* transformType_m; // transform type for each dimension
int* axisLength_m; // length of each transform dimension
T** trig_m; // trigonometric tables
T** trig_zeta; // trigonometric tables for IMPCFFT
};
......@@ -211,9 +233,12 @@ template <class T>
inline
FFTPACK<T>::~FFTPACK(void) {
// delete storage
for (unsigned d=0; d<numTransformDims_m; ++d)
delete [] trig_m[d];
for (unsigned d=0; d<numTransformDims_m; ++d) {
delete [] trig_m[d];
delete [] trig_zeta[d];
}
delete [] trig_m;
delete [] trig_zeta;
delete [] transformType_m;
delete [] axisLength_m;
}
......@@ -236,6 +261,7 @@ FFTPACK<T>::setup(unsigned numTransformDims, const int* transformTypes,
// allocate and initialize trig table
trig_m = new T*[numTransformDims_m];
trig_zeta = new T*[numTransformDims_m];
for (d=0; d<numTransformDims_m; ++d) {
switch (transformType_m[d]) {
case 0: // CC FFT
......@@ -254,6 +280,11 @@ FFTPACK<T>::setup(unsigned numTransformDims, const int* transformTypes,
trig_m[d] = new T[static_cast<int>(3 * axisLength_m[d]) + 15];
FFTPACK_wrap<T>::rrcosti(axisLength_m[d], trig_m[d]);
break;
case 4: // Imp FFT
trig_m[d] = new T[static_cast<int>(4 * axisLength_m[d]) + 15];
trig_zeta[d] = new T[static_cast<int>(2 * axisLength_m[d])];
FFTPACK_wrap<T>::impcffti(axisLength_m[d], trig_m[d], trig_zeta[d]);
break;
default:
ERRORMSG("Unknown transform type requested!!" << endl);
break;
......@@ -262,6 +293,43 @@ FFTPACK<T>::setup(unsigned numTransformDims, const int* transformTypes,
return;
}
// invoke FFT on complex data for given dimension and direction
template <class T>
inline void
FFTPACK<T>::callFFT(unsigned transformDim, int direction,
FFTPACK<T>::Complex_t* data,
FFTPACK<T>::Complex_t* data2) {
// check transform dimension and direction arguments
PAssert(transformDim<numTransformDims_m);
PAssert(direction==+1 || direction==-1);
// cast complex number pointer to T* for calling Fortran routines
T* rdata = reinterpret_cast<T*>(data);
T* rdata2 = reinterpret_cast<T*>(data2);
// branch on transform type for this dimension
switch (transformType_m[transformDim]) {
case 4: // CC FFT
if (direction == +1) {
// call forward complex-to-complex FFT
FFTPACK_wrap<T>::impcfftf(axisLength_m[transformDim], rdata,
rdata2, trig_m[transformDim],
trig_zeta[transformDim]);
}
else {
// call backward complex-to-complex FFT
FFTPACK_wrap<T>::impcfftb(axisLength_m[transformDim], rdata,
rdata2, trig_m[transformDim],
trig_zeta[transformDim]);
}
break;
default:
ERRORMSG("This settings are only used for IMP CC FFT" << endl);
break;
}
return;
}
// invoke FFT on complex data for given dimension and direction
template <class T>
......@@ -319,6 +387,9 @@ FFTPACK<T>::callFFT(unsigned transformDim, int direction,
case 3: // Cosine transform
ERRORMSG("Input for real-to-real FFT Cossine transform should be real!!" << endl);
break;
case 4: // Implicitly zero padded FFT transform
ERRORMSG("Input for implicitly zero padeed FFT should have two inputs!!" << endl);
break;
default:
ERRORMSG("Unknown transform type requested!!" << endl);
break;
......@@ -354,6 +425,9 @@ FFTPACK<T>::callFFT(unsigned transformDim, int direction, T* data) {
FFTPACK_wrap<T>::rrcost(axisLength_m[transformDim], data,
trig_m[transformDim]);
break;
case 4: // IMP CC FFT
ERRORMSG("Input for IMP complex-to-complex FFT should be complex!!" << endl);
break;
default:
ERRORMSG("Unknown transform type requested!!" << endl);
break;
......
ippl/src/Field/CompressedBrickIterator.hpp
View file @ e16c2223
......@@ -301,14 +301,17 @@ permute(const CompressedBrickIterator<T,D1>& iter,
{
// Compare with every Index in perm.
for (d2=0; d2<D2; ++d2)
{
// If they match, check the next one.
if ( perm[d2].sameBase( current[d1] ) )
{
goto FoundIt;
}
}
// Didn't find it.
// Make sure the length is 1.
PAssert( current[d1].length() == 1 );
FoundIt:
;
}
......@@ -332,8 +335,7 @@ permute(const CompressedBrickIterator<T,D1>& iter,
// The size of the loop comes from the permuted NDIndex.
permute.SetCount(d2,perm[d2].length());
// Set the counters to zero.
permute.ResetCounter(d2);
// Set the stride to zero in case we don't find a match below.
permute.ResetCounter(d2); // Set the stride to zero in case we don't find a match below.
permute.SetStride(d2,0);
// Check each Index in current to find a match.
for (d1=0; d1<D1; ++d1)
......
ippl/test/FFT/CMakeLists.txt
View file @ e16c2223
......@@ -34,6 +34,7 @@ add_executable (TestRC TestRC.cpp)
#add_executable (TestRCGPU TestRCGPU.cpp)
#add_executable (TestRCMIC TestRCMIC.cpp)
add_executable (TestFFTCos TestFFTCos.cpp)
add_executable (TestImpRC TestImpRC.cpp)
target_link_libraries (TestFFT ${IPPL_LIBS})
......@@ -41,5 +42,6 @@ target_link_libraries (TestRC ${IPPL_LIBS})
#target_link_libraries (TestRCGPU ${IPPL_LIBS})
#target_link_libraries (TestRCMIC ${IPPL_LIBS})
target_link_libraries (TestFFTCos ${IPPL_LIBS})
target_link_libraries (TestImpRC ${IPPL_LIBS})
ippl/test/FFT/TestFFTCos.cpp
View file @ e16c2223
......@@ -10,8 +10,7 @@
*
***************************************************************************/
#include "Ippl.h"
#include <complex>
#include <complex>
using namespace std;
int main(int argc, char *argv[])
......@@ -49,19 +48,30 @@ int main(int argc, char *argv[])
ngrid[d] = atoi(argv[1]);
NDIndex<D> ndiStandard;
NDIndex<D> ndiStandardG01;
for (unsigned int d=0; d<D; d++)
ndiStandard[d] = Index(ngrid[d]);
ndiStandardG01[0] = Index(ngrid[2]);
ndiStandardG01[1] = Index(ngrid[0]-1);
ndiStandardG01[2] = Index(ngrid[1]);
msg << ndiStandard << ndiStandard[0].getBase()<<endl;
msg << ndiStandardG01 <<ndiStandardG01[0].getBase()<<endl;
msg << ndiStandard << endl;
// all parallel layout, standard domain, normal axis order
FieldLayout<D> layoutPPStan(ndiStandard,allParallel);
FieldLayout<D> layoutG0(ndiStandardG01,allParallel);
FieldLayout<D> layoutG1(ndiStandardG01,allParallel);
msg << layoutPPStan << endl;
BareField<double,D> RField(layoutPPStan);
BareField<double,D> RField_save(layoutPPStan);
BareField<double,D> diffField(layoutPPStan);
BareField<double,D> G0Field(layoutG0);
BareField<double,D> G1Field(layoutG1);
double kx;
kx = 1.0; // wavenumbers
......@@ -95,19 +105,25 @@ int main(int argc, char *argv[])
cosTransformDims[d] = true;
// ToDo
FFT<CosTransform,D,double> cosfft2(ndiStandard, cosTransformDims, compressTemps);
FFT<CosTransform,D,double> cosfft3(ndiStandard, ndiStandardG01, compressTemps);
// FFT<SineTransform,D,double> sinefft2(ndiStandard, cosTransformDims, compressTemps);
msg << &cosfft2 << endl;
//cosfft2.write(&cout);
// msg << &sinefft2 << endl;
msg << "In-place cosine transform using all-parallel layout ..." << endl;
IpplTimings::startTimer(cosFFTTimer);
// ToDo
cosfft2.transform(-1, RField);
cosfft2.transform(1, RField);
// cosfft2.transform(-1, RField);
// cosfft2.transform(1, RField);
cosfft3.transform(-1, RField,G0Field,G1Field);
msg << G0Field<<endl;
//sinefft2.transform(-1, RField);
//sinefft2.transform( 1, RField);
IpplTimings::stopTimer(cosFFTTimer);
......
ippl/test/FFT/TestFFTCos.cpp.TrivalTest 0 → 100644
View file @ e16c2223
/***************************************************************************
*
* The IPPL Framework
*
* Filename: TestFFTCos.cpp
*
* Usage: TestFFTCos #GridPoints
*
* Visit http://people.web.psi.ch/adelmann/ for more details
*
***************************************************************************/
#include "Ippl.h"
#include <complex>
using namespace std;
int main(int argc, char *argv[])
{
Ippl ippl(argc,argv);
static IpplTimings::TimerRef mainprgTimer = IpplTimings::getTimer("mainprg");
static IpplTimings::TimerRef cosFFTTimer = IpplTimings::getTimer("cosFFT");
IpplTimings::startTimer(mainprgTimer);
Inform msg("TestFFTCos ",0);
Inform msg2all("TestFFTCos ",INFORM_ALL_NODES);
//msg << "HEllo " << endl;
//msg2all << "HEllo " << endl;
const unsigned D=3U;
unsigned ngrid[D]; // grid sizes
double realDiff;
// Various counters, constants, etc:
double pi = acos(-1.0);
double twopi = 2.0*pi;
// Layout information:
e_dim_tag allParallel[D]; // Specifies SERIAL, PARALLEL dims
for (unsigned int d=0; d<D; d++)
allParallel[d] = PARALLEL;
// Compression of temporaries:
bool compressTemps = false;
for (unsigned int d=0; d<D; d++)
ngrid[d] = atoi(argv[1]);
NDIndex<D> ndiStandard;
for (unsigned int d=0; d<D; d++)
ndiStandard[d] = Index(ngrid[d]);
msg << ndiStandard << endl;
// all parallel layout, standard domain, normal axis order
FieldLayout<D> layoutPPStan(ndiStandard,allParallel);
msg << layoutPPStan << endl;
BareField<double,D> RField(layoutPPStan);
BareField<double,D> RField_save(layoutPPStan);
BareField<double,D> diffField(layoutPPStan);
double kx;
kx = 1.0; // wavenumbers
double xfact, yfact, zfact,ky,kz;
ky = 2.0; kz = 3.0; // wavenumbers
xfact = pi/(ngrid[0] + 1.0);
yfact = 2.0*twopi/(ngrid[1]);
zfact = 2.0*twopi/(ngrid[2]);
RField[ndiStandard[0]][ndiStandard[1]][ndiStandard[2]] =
1.0 * ( 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 ) ) +
0.0 * (-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 ) );
RField_save = RField;
bool cosTransformDims[D];
for (unsigned int d=0; d<D; ++d)
cosTransformDims[d] = true;
// ToDo
FFT<CosTransform,D,double> cosfft2(ndiStandard, cosTransformDims, compressTemps);
// FFT<SineTransform,D,double> sinefft2(ndiStandard, cosTransformDims, compressTemps);
msg << &cosfft2 << endl;
//cosfft2.write(&cout);
// msg << &sinefft2 << endl;
msg << "In-place cosine transform using all-parallel layout ..." << endl;
IpplTimings::startTimer(cosFFTTimer);
// ToDo
cosfft2.transform(-1, RField);
cosfft2.transform(1, RField);
//sinefft2.transform(-1, RField);
//sinefft2.transform( 1, RField);
IpplTimings::stopTimer(cosFFTTimer);
diffField = Abs(RField - RField_save);
realDiff = max(diffField);
msg << "fabs(realDiff) = " << fabs(realDiff) << endl;
IpplTimings::stopTimer(mainprgTimer);
IpplTimings::print();
IpplTimings::print(string(argv[0])+string(".timing"));
return 0;
}
/***************************************************************************
* $RCSfile: TestFFT.cpp.org,v $ $Author: adelmann $
* $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:36 $
***************************************************************************/
ippl/test/FFT/TestImpRC.cpp 0 → 100644
View file @ e16c2223
This diff is collapsed.
src/Classic/Algorithms/PBunchDefs.h
View file @ e16c2223
......@@ -47,6 +47,7 @@ typedef Field<int, 3, Mesh_t, Center_t> IField_t;
typedef Field<dcomplex, 3, Mesh_t, Center_t> CxField_t;
typedef FFT<RCTransform, 3, double> FFT_t;
typedef FFT<SineTransform, 3, double> SINE_t;
typedef FFT<CosTransform, 3, double> COS_t;
typedef FFT<CCTransform, 3, double> FFTC_t;
namespace ParticleType {
......
src/Classic/Algorithms/PartBunch.cpp
View file @ e16c2223
......@@ -653,6 +653,7 @@ size_t PartBunch::calcNumPartsOutside(Vector_t x) {
void PartBunch::computeSelfFields(int binNumber) {
IpplTimings::startTimer(selfFieldTimer_m);
INFOMSG("*** DAWEI IS TESTING line 656 ***" << endl);
/// Set initial charge density to zero. Create image charge
/// potential field.
......@@ -988,6 +989,7 @@ void PartBunch::resizeMesh() {
void PartBunch::computeSelfFields() {
IpplTimings::startTimer(selfFieldTimer_m);
INFOMSG("*** DAWEI IS TESTING line 992 ***" << endl);
rho_m = 0.0;
eg_m = Vector_t(0.0);
......@@ -1255,6 +1257,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
void PartBunch::computeSelfFields_cycl(double gamma) {
IpplTimings::startTimer(selfFieldTimer_m);
INFOMSG("*** DAWEI IS TESTING line 1260 ***" << endl);
/// set initial charge density to zero.
rho_m = 0.0;
......@@ -1308,7 +1311,9 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
/// calculate Possion equation (without coefficient: -1/(eps))
IpplTimings::startTimer(compPotenTimer_m);
INFOMSG("*** DAWEI IS TESTING line 1315 ***"<<fs_m->getFieldSolverType()<< endl);
fs_m->solver_m->computePotential(rho_m, hr_scaled);
INFOMSG("*** DAWEI IS TESTING line 1316 ***" << endl);
IpplTimings::stopTimer(compPotenTimer_m);
......@@ -1450,6 +1455,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
*/
void PartBunch::computeSelfFields_cycl(int bin) {
IpplTimings::startTimer(selfFieldTimer_m);
INFOMSG("*** DAWEI IS TESTING line 1455 ***" << endl);
/// set initial charge dentsity to zero.
rho_m = 0.0;
......
src/Solvers/EFFTPoissonSolver.cpp
View file @ e16c2223
This diff is collapsed.
src/Solvers/EFFTPoissonSolver.h
View file @ e16c2223
......@@ -46,7 +46,7 @@ public:
// constructor and destructor
EFFTPoissonSolver(PartBunch &bunch, std::string greensFuntion);
EFFTPoissonSolver(Mesh_t *mesh, FieldLayout_t *fl, std::string greensFunction, std::string bcz);
EFFTPoissonSolver(Mesh_t *mesh, FieldLayout_t *fl, std::string greensFunction, std::string bcz);
~EFFTPoissonSolver();
......@@ -56,22 +56,34 @@ public:
// given a charge-density field rho and a set of mesh spacings hr,
// compute the scalar potential with image charges at -z
void computePotential(Field_t &rho, Vector_t hr, double zshift);
void computePotential(Field_t &rho, Vector_t hr, double zshift); //invalid
// given a charge-density field rho and a set of mesh spacings hr,
// compute the scalar potential in open space
void computePotential(Field_t &rho, Vector_t hr);
// the multiply function for the efft Poisson solver
void multiplyrhogrn();
void submultiplyrhogrn(
CxField_t &rho0, CxField_t &rho1,
CxField_t &imgrho0, CxField_t &imgrho1,
Field_t &grn, Field_t &grnexk0, Field_t &grnexk1);
// compute the green's function for a Poisson problem and put it in in grntm_m
// uses grnIField_m to eliminate excess calculation in greenFunction()
// given mesh information in nr and hr
void greensFunction();
// void greensFunction(); invalid
/// compute the integrated Green function as described in <A HREF="http://prst-ab.aps.org/abstract/PRSTAB/v9/i4/e044204">Three-dimensional quasistatic model for high brightness beam dynamics simulation</A> by Qiang et al.
void integratedGreensFunction();
// void integratedGreensFunction(); invalid
/// efficient integrated Green function
void effntintegratedGreensFunction();
/// compute the shifted integrated Green function as described in <A HREF="http://prst-ab.aps.org/abstract/PRSTAB/v9/i4/e044204">Three-dimensional quasistatic model for high brightness beam dynamics simulation</A> by Qiang et al.
void shiftedIntGreensFunction(double zshift);
// void shiftedIntGreensFunction(double zshift); not valid
double getXRangeMin() {return 1.0;}
double getXRangeMax() {return 1.0;}
......@@ -84,10 +96,10 @@ public:
Inform &print(Inform &os) const;
private:
void mirrorRhoField() FIELDASSIGNOPTIMIZATION;
void mirrorRhoField(Field_t & ggrn2);// FIELDASSIGNOPTIMIZATION;
// void mirrorRhoField() FIELDASSIGNOPTIMIZATION;
// void mirrorRhoField(Field_t & ggrn2);// FIELDASSIGNOPTIMIZATION;
// rho2_m is the charge-density field with mesh doubled in each dimension
// rho2_m is the charge-density field with mesh doubled in the first dimension
Field_t rho2_m;
// real field with layout of complex field: domain3_m
......@@ -96,11 +108,21 @@ private:
// rho2tr_m is the Fourier transformed charge-density field
// domain3_m and mesh3_ are used
CxField_t rho2tr_m;
CxField_t rho2tr_m01;
CxField_t rho2tr_m10;
CxField_t rho2tr_m11;
CxField_t imgrho2tr_m;
CxField_t imgrho2tr_m01;
CxField_t imgrho2tr_m10;
CxField_t imgrho2tr_m11;
Field_t igrn_kji;
Field_t igrn_kije;
Field_t igrn_kijo;
// grntr_m is the Fourier transformed Green's function
// domain3_m and mesh3_ are used
CxField_t grntr_m;
Field_t igrnexk0;//temp field, same size as original rho_m, with kij order
Field_t igrnexk1;//temp field, same size as original rho_m
#ifdef OPAL_DKS
//pointer for Fourier transformed Green's function on GPU
......@@ -122,35 +144,47 @@ private:
// the FFT object
std::unique_ptr<FFT_t> fft_m;
std::unique_ptr<COS_t> cost_m;
// mesh and layout objects for rho_m
Mesh_t *mesh_m;
FieldLayout_t *layout_m;
// mesh and layout objects for rho2_m
// mesh and layout objects for rho2_m, double sized in the first dimension
std::unique_ptr<Mesh_t> mesh2_m;
std::unique_ptr<FieldLayout_t> layout2_m;
//
//after the fft transform for rho2_m
std::unique_ptr<Mesh_t> mesh3_m;
std::unique_ptr<FieldLayout_t> layout3_m;
// mesh and layout for integrated greens function
std::unique_ptr<Mesh_t> mesh4_m;
std::unique_ptr<FieldLayout_t> layout4_m;
std::unique_ptr<Mesh_t> mesh4_kji;
std::unique_ptr<FieldLayout_t> layout4_kji;
// std::unique_ptr<Mesh_t> mesh4_jik;
// std::unique_ptr<FieldLayout_t> layout4_jik;
std::unique_ptr<Mesh_t> mesh4_kij;
std::unique_ptr<FieldLayout_t> layout4_kij;
// tmp
Field_t tmpgreen;
Field_t tmpgreen_kji;
// domains for the various fields
NDIndex<3> domain_m; // original domain, gridsize
// mesh and gridsize defined outside of FFT class, given as
// parameter to the constructor (mesh and layout object).
NDIndex<3> domain2_m; // doubled gridsize (2*Nx,2*Ny,2*Nz)
NDIndex<3> domain2_m; // doubled gridsize (2*Nx,Ny,Nz)
NDIndex<3> domain3_m; // field for the complex values of the RC transformation
NDIndex<3> domain4_m;
NDIndex<3> domain4_m; // field for green function
NDIndex<3> domain4_kji; // field for green function
//NDIndex<3> domain4_jik; // field for green function temp
NDIndex<3> domain4_kij; // field for green function
// (2*Nx,Ny,2*Nz)
NDIndex<3> domainFFTConstruct_m;
NDIndex<3> domainCOSTConstruct_jik;