PoissonConvolve3D Class Reference

#include <PoissonConvolve3D.hh>

List of all members.

Detailed Description

class PoissonConvolve3D

A solver for Poisson's equation in 3-D Cartesian coordinates {x,y,z}, using convolution of the Green's function (whch has infinite boundary conditions). The convolution is sped up by using an FFT on a grid twice as large (in each dimension) as the solution (physical) grid.

This class uses the Poisson solver on a rectangular 3-d grid; inside the grid, it calculates E by applying the analytic derivative of the interpolation function to the values of phi at the surrounding 8 grid points. Outside the grid, by default it uses the analytic value of E for a single point charge, with its charge equal to the total charge in the grid, and its position at the mean position of the charge in the grid. More complicated approximations, based on multiple point charges, can be used (not discussed here).

Units are determined by the code using this class. Normally the charge values will not be integers, but rather integers times the positron charge divided by epsilon_0 (the permittivity of free space), the grid sizes will have units of length; the result is that phi() and getE() will then have the appropriate units.

This class solves: del^2 phi = -4*pi*rho On a specified grid in {x,y,z}. rho(x,y,z) is the charge density. The -4*pi and the conversion from charge to charge density is handled internally in putCharge(); they are NOT handled in setRhs() (which is intended primarily for testing).

This class includes an approximation for use outside the grid. The default is simply to use the total charge put into the grid (via putCharge()), and the weighted mean position of that charge. Better approximations can be computed externally as a vector of point charges and used via setApproximation().

Note: array[] contains CHARGE between calls to zeroRhs() and solve(); it contains PHI between calls to solve() and zeroRhs(); status reflects its contents, ans assert()-s are used to ensure accesses are valid..

Accuracy for a 65x65x65 grid of unit size: ... MORE TO COME ...

The CPU time increases with increasing N, and can vary by as much as a factor of 3 for adjacent values. The following are good values: 128, 120, 112, 100, 96, 90, 84, 80, 75, 70, 64, 60, 56, 48, 42, 36, 32. Below 28 the variations are small. The following are particularly SLOW values: 129, 127, 123, 107, 103, 101.

Public Types

enum  PutMode { LOWEST, NEAREST, PRORATED }
 enum PutMode determines how putCharge() places charge into the grid. PRORATED is the default. LOWEST places charge into the grid point lower than (x,y,z). NEAREST places charge into the nearest grid point to (x,y,z). PRORATED places charge into the 8 nearest grid points, pro-rated by proximity to (x,y,z). For 1M particles into a 129^3 grid, NEAREST and LOWER were less than 10% faster than PRORATED, giving slightly more ragged results. More...

Public Member Functions

 PoissonConvolve3D (int _nx, int _ny, int _nz, float _xMin, float _xMax, float _yMin, float _yMax, float _zMin, float _zMax)
 Constructor. nx,ny,nz are the # grid points in the 3 dimensions -- NOTE: these are the size of the PHYSICAL grid, the actual grid is doubled in each dimension. They can be any integers >= 9, but powers of 2 have a performance advantage. Note that 512x512x512 allocates too much memory for a 32-bit program on an Intel processor. xMin,xMax,yMin,yMax,zMin,zMax are the grid limits in 3 dimensions. That implies array[index(0,iy,iz)] is the boundary at x=xMin; array[index(nx-1,iy,iz)] is the boundary at x=xMax; etc. 1.0e-4 is appropriate.
virtual ~PoissonConvolve3D ()
 Destructor.
void updatePosition (float _xMin, float _xMax, float _yMin, float _yMax, float _zMin, float _zMax)
 updatePosition() updates the position of the boundaries. Usually followed by a call to setApproximation(), with the approximation charge positions updated similarly, or by a series of calls to putCharge() and then defaultApproximation().
float phi (float x, float y, float z)
 phi() works for all {x,y,z}, using the solution inside the grid and the approximation outside the grid. Note that phi() is not necessarily continuous or differentiable, but it should be approximately continuous. Do not differentiate this function, use getE() instead (it is MUCH more intelligent about the grid and its bounadries). Intended primarily for testing.
float phiGrid (int ix, int iy, int iz) const
 phiGrid() returns the value of the solution to Poisson's equation (i.e. the electrostatic potential) on grid points. ix,iy,iz MUST be in range: 0..nx-1, 0..ny-1, 0..nz-1 Intended primarily for testing.
void setPhi (int ix, int iy, int iz, float v)
 setPhi() will set a value in array[]. Intended primarily for testing.
void addPhi (int ix, int iy, int iz, float v)
 addPhi() will add v to a value in array[]. Intended primarily for testing.
float rhs (int ix, int iy, int iz) const
 rhs() returns the charge at a grid point of the right-hand-side of Poisson's equation. NOTE: ix,iy,iz MUST be in range: 0..nx-1, 0..ny-1, 0..nz-1
void setRhs (int ix, int iy, int iz, float v)
 setRhs() will set an entry in rhs[]. Normally putCharge() is better. ix,iy,iz MUST be in range: 0..nx-1, 0..ny-1, 0..nz-1 This directly sets the rhs array entry, so you probably need to multiply by -4*pi/(dx*dy*dz) to get the proper R.H.S. for Poisson's equation; putCharge() does that internally. Intended primarily for testing.
bool putCharge (float x, float y, float z, float c)
 putCharge() puts a charge into the rhs grid according to the current putMode. Does nothing if outside the grid. Multiplies by -1.0/(dx*dy*dz) to convert the charge to the R.H.S. of Poisson's equation, which is minus the charge density. NOTE: dividing by epsilon0 is up to the user's choice of units. Accumulates the total charge and mean position of charge actually placed into the grid (for use by defaultApproximation()). Returns true if the charge was placed within the grid, false if not.
bool putCharge (Charge &c)
void setPutMode (PutMode m)
 setPutMode() sets the mode of putCharge(). See PutMode above.
void zeroRhs ()
 zeroRhs() will zero rhs.
G4ThreeVector getE (float x, float y, float z)
 getE() returns the value of the E field at a given point. Uses the solution inside the grid and the approximation outside. It handles grid cells and boundaries intelligently.
G4ThreeVector approxE (double x, double y, double z)
 approxE() evaluates the approximation E field. {x,y,z} can be anywhere, but normally are outside the grid or on its boundary.
float approxPhi (double x, double y, double z)
 approxPhi() evaluates the approximation phi. {x,y,z} can be anywhere, but normally are outside the grid or on its boundary.
G4ThreeVector getPosition (int ix, int iy, int iz)
 getPosition() returns the position corresponding to grid point indexes. Intended primarily for testing.
void setApproximation (std::vector< Charge > v)
 setApproximation() sets the approximation to use outside the grid, and on its boundary. v is a vector of point charges. Updates the boundary values of phi.
void defaultApproximation ()
 defaultApproximation() calls setApproximation() with the total charge and mean position in rhs (from putCharge()). Updates the boundary values of phi, and the interior if interior=true.
void computeGreensFunction ()
 computeGreensFunction() does just that. Must be called whenever the physical size of the grid is changed.
void solve ()
 solve() will solve Poisson's equation
void printRhs ()
 printRhs() prints the rhs. Intended primarily for testing.
void printPhi ()
 printPhi() prints phi. Intended primarily for testing.
int getNx () const
 accessors.
int getNy () const
int getNz () const

Protected Types

enum  Status { UNINITIALIZED, CHARGE, PHI }

Protected Member Functions

int index (int ix, int iy, int iz) const

Protected Attributes

int nx
int ny
int nz
float xMin
float xMax
float yMin
float yMax
float zMin
float zMax
Status status
float dx
float dy
float dz
float * array
fftwf_complex * transf
fftwf_complex * greensFunction
std::vector< Charge > approx
float sumX
float sumY
float sumZ
float totalCharge
float vol
PutMode putMode
int nc
fftwf_plan gfPlan
fftwf_plan plan1
fftwf_plan plan2

Member Enumeration Documentation

enum PoissonConvolve3D::PutMode

enum PutMode determines how putCharge() places charge into the grid. PRORATED is the default. LOWEST places charge into the grid point lower than (x,y,z). NEAREST places charge into the nearest grid point to (x,y,z). PRORATED places charge into the 8 nearest grid points, pro-rated by proximity to (x,y,z). For 1M particles into a 129^3 grid, NEAREST and LOWER were less than 10% faster than PRORATED, giving slightly more ragged results.

Enumerator:
LOWEST 
NEAREST 
PRORATED  
00126 {LOWEST, NEAREST, PRORATED};

enum PoissonConvolve3D::Status [protected]

Enumerator:
UNINITIALIZED 
CHARGE 
PHI  
00128 {UNINITIALIZED, CHARGE, PHI};

Constructor & Destructor Documentation

PoissonConvolve3D::PoissonConvolve3D ( int  _nx,
int  _ny,
int  _nz,
float  _xMin,
float  _xMax,
float  _yMin,
float  _yMax,
float  _zMin,
float  _zMax 
)

Constructor. nx,ny,nz are the # grid points in the 3 dimensions -- NOTE: these are the size of the PHYSICAL grid, the actual grid is doubled in each dimension. They can be any integers >= 9, but powers of 2 have a performance advantage. Note that 512x512x512 allocates too much memory for a 32-bit program on an Intel processor. xMin,xMax,yMin,yMax,zMin,zMax are the grid limits in 3 dimensions. That implies array[index(0,iy,iz)] is the boundary at x=xMin; array[index(nx-1,iy,iz)] is the boundary at x=xMax; etc. 1.0e-4 is appropriate.

References array, gfPlan, greensFunction, nc, nx, ny, nz, plan1, plan2, PRORATED, putMode, status, sumX, sumY, sumZ, totalCharge, transf, UNINITIALIZED, and updatePosition().

00019                                                           : approx() {
00020         nx=_nx; ny=_ny; nz=_nz;
00021         nc = 4*nx*ny*(nz+1); // double size, except last dimension
00022         status = UNINITIALIZED;
00023         array = (float*)fftwf_malloc(sizeof(float)*8*nx*ny*nz);
00024         transf = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nc);
00025         greensFunction = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nc);
00026         sumX = 0.0; sumY = 0.0; sumZ = 0.0; totalCharge = 0.0;
00027         putMode = PRORATED;
00028 
00029         fftwf_set_timelimit(10.0); // otherwise can take up to 226 sec (!)
00030 
00031         gfPlan = fftwf_plan_dft_r2c_3d(2*nx,2*ny,2*nz,array,greensFunction,
00032                                                         FFTW_ESTIMATE);
00033         plan1 = fftwf_plan_dft_r2c_3d(2*nx,2*ny,2*nz,array,transf,
00034                                                         FFTW_ESTIMATE);
00035         plan2 = fftwf_plan_dft_c2r_3d(2*nx,2*ny,2*nz,transf,array,
00036                                                         FFTW_ESTIMATE);
00037 
00038         updatePosition(_xMin,_xMax,_yMin,_yMax,_zMin,_zMax);
00039 }

virtual PoissonConvolve3D::~PoissonConvolve3D (  )  [inline, virtual]

Destructor.

References approx, array, gfPlan, greensFunction, plan1, plan2, and transf.

00164                                      {
00165                 if(array) fftwf_free(array);
00166                 if(transf) fftwf_free(transf);
00167                 if(greensFunction) fftwf_free(greensFunction);
00168                 fftwf_destroy_plan(gfPlan);
00169                 fftwf_destroy_plan(plan1);
00170                 fftwf_destroy_plan(plan2);
00171                 approx.clear();
00172         }

Member Function Documentation

int PoissonConvolve3D::index ( int  ix,
int  iy,
int  iz 
) const [inline, protected]

References ny, and nz.

Referenced by addPhi(), computeGreensFunction(), getE(), phi(), phiGrid(), printPhi(), printRhs(), putCharge(), rhs(), setPhi(), setRhs(), and zeroRhs().

00144                                                 {
00145                 assert(ix>=0&&ix<2*nx && iy>=0&&iy<2*ny && iz>=0&&iz<2*nz);
00146                 return iz+2*nz*(iy+2*ny*ix);  // C order
00147         }

void PoissonConvolve3D::updatePosition ( float  _xMin,
float  _xMax,
float  _yMin,
float  _yMax,
float  _zMin,
float  _zMax 
) [inline]

updatePosition() updates the position of the boundaries. Usually followed by a call to setApproximation(), with the approximation charge positions updated similarly, or by a series of calls to putCharge() and then defaultApproximation().

References computeGreensFunction(), dx, dy, dz, nx, ny, nz, vol, xMax, xMin, yMax, yMin, zMax, and zMin.

Referenced by Bunch::computeField(), Bunch::init(), and PoissonConvolve3D().

00179                                                                             {
00180                 xMin=_xMin; xMax=_xMax;
00181                 yMin=_yMin; yMax=_yMax;
00182                 zMin=_zMin; zMax=_zMax;
00183                 dx = (xMax-xMin)/(nx-1);
00184                 dy = (yMax-yMin)/(ny-1);
00185                 dz = (zMax-zMin)/(nz-1);
00186                 vol = dx*dy*dz;
00187                 computeGreensFunction();
00188         }

float PoissonConvolve3D::phi ( float  x,
float  y,
float  z 
)

phi() works for all {x,y,z}, using the solution inside the grid and the approximation outside the grid. Note that phi() is not necessarily continuous or differentiable, but it should be approximately continuous. Do not differentiate this function, use getE() instead (it is MUCH more intelligent about the grid and its bounadries). Intended primarily for testing.

References approx, array, dx, dy, dz, index(), M_PI, nx, ny, nz, PHI, status, xMin, yMin, and zMin.

Referenced by approxPhi().

00042 {
00043         assert(status==PHI);
00044         int ix = (int)floor((x-xMin)/dx);
00045         int iy = (int)floor((y-yMin)/dy);
00046         int iz = (int)floor((z-zMin)/dz);
00047         if(ix < 0 || iy < 0 || iz < 0 || ix+1 >= nx || iy+1 >= ny || 
00048                                                                 iz+1 >= nz) {
00049                 // outside grid, use the approximation
00050                 float v=0.0;
00051                 G4ThreeVector pos(x,y,z);
00052                 for(unsigned i=0; i<approx.size(); ++i) {
00053                         G4ThreeVector pt = approx[i].getPosition();
00054                         pt -= pos;
00055                         double r = pt.mag();
00056                         v += approx[i].c / r / (4.0*M_PI);
00057                 }
00058                 return v;
00059         }
00060 
00061         // within grid, interpolate linearly among 8 surrounding grid points
00062         float fx = (x - (xMin+ix*dx))/dx;
00063         assert(-0.01 <= fx && fx <= 1.01);
00064         float fy = (y - (yMin+iy*dy))/dy;
00065         assert(-0.01 <= fy && fy <= 1.01);
00066         float fz = (z - (zMin+iz*dz))/dz;
00067         assert(-0.01 <= fz && fz <= 1.01);
00068         float gx=1.0-fx, gy=1.0-fy, gz=1.0-fz;
00069 
00070         float v =
00071                 gx*gy*gz*array[index(ix,iy,iz)] +
00072                 fx*gy*gz*array[index(ix+1,iy,iz)] +
00073                 gx*fy*gz*array[index(ix,iy+1,iz)] +
00074                 fx*fy*gz*array[index(ix+1,iy+1,iz)] +
00075                 gx*gy*fz*array[index(ix,iy,iz+1)] +
00076                 fx*gy*fz*array[index(ix+1,iy,iz+1)] +
00077                 gx*fy*fz*array[index(ix,iy+1,iz+1)] +
00078                 fx*fy*fz*array[index(ix+1,iy+1,iz+1)];
00079 
00080         return v;
00081 }

float PoissonConvolve3D::phiGrid ( int  ix,
int  iy,
int  iz 
) const [inline]

phiGrid() returns the value of the solution to Poisson's equation (i.e. the electrostatic potential) on grid points. ix,iy,iz MUST be in range: 0..nx-1, 0..ny-1, 0..nz-1 Intended primarily for testing.

References array, index(), PHI, and status.

00204                 { assert(status==PHI);  return array[index(ix,iy,iz)]; }

void PoissonConvolve3D::setPhi ( int  ix,
int  iy,
int  iz,
float  v 
) [inline]

setPhi() will set a value in array[]. Intended primarily for testing.

References array, index(), PHI, and status.

00209                 { assert(status==PHI);  array[index(ix,iy,iz)] = v; }

void PoissonConvolve3D::addPhi ( int  ix,
int  iy,
int  iz,
float  v 
) [inline]

addPhi() will add v to a value in array[]. Intended primarily for testing.

References array, index(), PHI, and status.

00214                 { assert(status==PHI);  array[index(ix,iy,iz)] += v; }

float PoissonConvolve3D::rhs ( int  ix,
int  iy,
int  iz 
) const [inline]

rhs() returns the charge at a grid point of the right-hand-side of Poisson's equation. NOTE: ix,iy,iz MUST be in range: 0..nx-1, 0..ny-1, 0..nz-1

References array, CHARGE, index(), M_PI, and status.

00220                 { assert(status==CHARGE);
00221                   return 4.0*M_PI*array[index(ix,iy,iz)]; }

void PoissonConvolve3D::setRhs ( int  ix,
int  iy,
int  iz,
float  v 
) [inline]

setRhs() will set an entry in rhs[]. Normally putCharge() is better. ix,iy,iz MUST be in range: 0..nx-1, 0..ny-1, 0..nz-1 This directly sets the rhs array entry, so you probably need to multiply by -4*pi/(dx*dy*dz) to get the proper R.H.S. for Poisson's equation; putCharge() does that internally. Intended primarily for testing.

References array, CHARGE, index(), and status.

00230                 { assert(status==CHARGE); array[index(ix,iy,iz)] = v; }

bool PoissonConvolve3D::putCharge ( float  x,
float  y,
float  z,
float  c 
)

putCharge() puts a charge into the rhs grid according to the current putMode. Does nothing if outside the grid. Multiplies by -1.0/(dx*dy*dz) to convert the charge to the R.H.S. of Poisson's equation, which is minus the charge density. NOTE: dividing by epsilon0 is up to the user's choice of units. Accumulates the total charge and mean position of charge actually placed into the grid (for use by defaultApproximation()). Returns true if the charge was placed within the grid, false if not.

References array, CHARGE, dx, dy, dz, index(), LOWEST, M_PI, NEAREST, nx, ny, nz, PRORATED, putMode, status, sumX, sumY, sumZ, totalCharge, xMin, yMin, and zMin.

Referenced by Bunch::computeField(), putCharge(), Bunch::testBfield(), and Bunch::testEfield().

00084 {
00085         assert(status==CHARGE);
00086         int ix = (int)floor((x-xMin)/dx);
00087         int iy = (int)floor((y-yMin)/dy);
00088         int iz = (int)floor((z-zMin)/dz);
00089         // (omit boundaries)
00090         if(ix<=0 || iy<=0 || iz<=0 || ix+1>=nx || iy+1>=ny || iz+1>=nz) {
00091                 return false;
00092         }
00093 
00094         sumX += c*x;
00095         sumY += c*y;
00096         sumZ += c*z;
00097         totalCharge += c;
00098 
00099         c /= 4.0*M_PI;
00100 
00101         switch(putMode) {
00102         case LOWEST:
00103                 array[index(ix,iy,iz)] += c;
00104                 break;
00105         case NEAREST:
00106                 if(x > xMin+ix*dx+dx/2.0) ++ix; // could be upper boundary
00107                 if(y > yMin+iy*dy+dy/2.0) ++iy; // but is guaranteed to 
00108                 if(z > zMin+iz*dz+dz/2.0) ++iz; // be in range.
00109                 array[index(ix,iy,iz)] += c;
00110                 break;
00111         case PRORATED:
00112                 { float fx = (x - (xMin+ix*dx))/dx;
00113                   assert(-0.01 <= fx && fx <= 1.01);
00114                   float fy = (y - (yMin+iy*dy))/dy;
00115                   assert(-0.01 <= fy && fy <= 1.01);
00116                   float fz = (z - (zMin+iz*dz))/dz;
00117                   assert(-0.01 <= fz && fz <= 1.01);
00118                   float gx=1.0-fx, gy=1.0-fy, gz=1.0-fz;
00119                   array[index(ix,iy,iz)] += gx*gy*gz*c;
00120                   array[index(ix+1,iy,iz)] += fx*gy*gz*c;
00121                   array[index(ix,iy+1,iz)] += gx*fy*gz*c;
00122                   array[index(ix+1,iy+1,iz)] += fx*fy*gz*c;
00123                   array[index(ix,iy,iz+1)] += gx*gy*fz*c;
00124                   array[index(ix+1,iy,iz+1)] += fx*gy*fz*c;
00125                   array[index(ix,iy+1,iz+1)] += gx*fy*fz*c;
00126                   array[index(ix+1,iy+1,iz+1)] += fx*fy*fz*c;
00127                 }
00128                 break;
00129         }
00130 
00131         return true;
00132 }

bool PoissonConvolve3D::putCharge ( Charge &  c  )  [inline]

References Charge::c, putCharge(), Charge::x, Charge::y, and Charge::z.

00241 { return putCharge(c.x,c.y,c.z,c.c); }

void PoissonConvolve3D::setPutMode ( PutMode  m  )  [inline]

setPutMode() sets the mode of putCharge(). See PutMode above.

References putMode.

00244 { putMode = m; }

void PoissonConvolve3D::zeroRhs (  ) 

zeroRhs() will zero rhs.

References array, CHARGE, index(), nx, ny, nz, status, sumX, sumY, sumZ, and totalCharge.

Referenced by Bunch::computeField(), Bunch::testBfield(), and Bunch::testEfield().

00135 {
00136         sumX = sumY = sumZ = totalCharge = 0;
00137         for(int i=0; i<2*nx; ++i) {
00138                 for(int j=0; j<2*ny; ++j) {
00139                         for(int k=0; k<2*nz; ++k) {
00140                                 array[index(i,j,k)] = 0.0;
00141                         }
00142                 }
00143         }
00144 
00145         // indicate CHARGE is now in array[]
00146         status = CHARGE;
00147 }

G4ThreeVector PoissonConvolve3D::getE ( float  x,
float  y,
float  z 
)

getE() returns the value of the E field at a given point. Uses the solution inside the grid and the approximation outside. It handles grid cells and boundaries intelligently.

References approxE(), array, dx, dy, dz, index(), nx, ny, nz, PHI, status, xMax, xMin, yMax, yMin, zMax, and zMin.

Referenced by Bunch::addFieldValue(), Bunch::testBfield(), and Bunch::testEfield().

00150 {
00151         assert(status==PHI);
00152         int ix = (int)floor((x-xMin)/dx);
00153         int iy = (int)floor((y-yMin)/dy);
00154         int iz = (int)floor((z-zMin)/dz);
00155         if(ix < 0 || iy < 0 || iz < 0 || ix+1 >= nx || iy+1 >= ny || 
00156                                                                 iz+1 >= nz) {
00157                 return approxE(x,y,z);
00158         }
00159 
00160 #ifdef SIMPLE_LINEAR
00161         // Has the problem that the derivative is constant in each box.
00162         // variables used in mathematica notebook "interpolationOn3DGrid.nb"
00163         // -- it differentiates the 8-point interpolation formula.
00164         double x0=xMin+ix*dx, y0=yMin+iy*dy, z0=zMin+iz*dz;
00165         double x1=x0+dx, y1=y0+dy, z1=z0+dz;
00166         double v000=array[index(ix,iy,iz)];
00167         double v100=array[index(ix+1,iy,iz)];
00168         double v010=array[index(ix,iy+1,iz)];
00169         double v110=array[index(ix+1,iy+1,iz)];
00170         double v001=array[index(ix,iy,iz+1)];
00171         double v101=array[index(ix+1,iy,iz+1)];
00172         double v011=array[index(ix,iy+1,iz+1)];
00173         double v111=array[index(ix+1,iy+1,iz+1)];
00174 
00175         // E = - grad array[]
00176         double Ex=-(1.0/(dx*dy*dz))*(dy*(dz*(-v000+v100)+(v000-v001-v100+v101)
00177                 *(z-z0))+(y-y0)*(dz*(v000-v010-v100+v110)-
00178                 (v000-v001-v010+v011-v100+v101+v110-v111)*(z-z0)));
00179         double Ey=-(1.0/(dx*dy*dz))*(dx*(dz*(-v000+v010)+(v000-v001-v010+v011)
00180                 *(z-z0))+(x-x0)*(dz*(v000-v010-v100+v110)-
00181                 (v000-v001-v010+v011-v100+v101+v110-v111)*(z-z0)));
00182         double Ez=-(1.0/(dx*dy*dz))*(dx*(dy*(-v000+v001)+(v000-v001-v010+v011)
00183                 *(y-y0))+(x-x0)*(dy*(v000-v001-v100+v101)-
00184                 (v000-v001-v010+v011-v100+v101+v110-v111)*(y-y0)));
00185 #else
00186         // linear interpolation between derivatives in adjacent cells.
00187         // Associate the derivative in each grid box with the center of the box;
00188         // interpolate the derivative between the current box and the closer
00189         // adjacent one, or the boundary if this box is an edge.
00190         float fx = (x - (xMin+ix*dx))/dx;
00191         assert(-0.01 <= fx && fx <= 1.01);
00192         float fy = (y - (yMin+iy*dy))/dy;
00193         assert(-0.01 <= fy && fy <= 1.01);
00194         float fz = (z - (zMin+iz*dz))/dz;
00195         assert(-0.01 <= fz && fz <= 1.01);
00196         float gx=1.0-fx, gy=1.0-fy, gz=1.0-fz;
00197 
00198         double v000=array[index(ix,iy,iz)];
00199         double v100=array[index(ix+1,iy,iz)];
00200         double v010=array[index(ix,iy+1,iz)];
00201         double v110=array[index(ix+1,iy+1,iz)];
00202         double v001=array[index(ix,iy,iz+1)];
00203         double v101=array[index(ix+1,iy,iz+1)];
00204         double v011=array[index(ix,iy+1,iz+1)];
00205         double v111=array[index(ix+1,iy+1,iz+1)];
00206 
00207         // partial d/dx
00208         double v0=gy*gz*v000 + fy*gz*v010 + gy*fz*v001 + fy*fz*v011;
00209         double v1=gy*gz*v100 + fy*gz*v110 + gy*fz*v101 + fy*fz*v111;
00210         double Ex=-(v1-v0)/dx; // for this grid box
00211         if(fx > 0.5) {
00212                 if(ix+1 == nx-1) {
00213                         // interpolate between center of box and xMax
00214                         double Ex2=approxE(xMax,y,z).x();
00215                         double f=2.0*fx-1.0;
00216                         Ex = (1.0-f)*Ex + f*Ex2;
00217                 } else {
00218                         // interpolate between this box and next higher one
00219                         double v200=array[index(ix+2,iy,iz)];
00220                         double v210=array[index(ix+2,iy+1,iz)];
00221                         double v201=array[index(ix+2,iy,iz+1)];
00222                         double v211=array[index(ix+2,iy+1,iz+1)];
00223                         double v2=gy*gz*v200 + fy*gz*v210 + gy*fz*v201 + 
00224                                                                 fy*fz*v211;
00225                         double f=fx-0.5;
00226                         Ex = (1.0-f)*Ex - f*(v2-v1)/dx;
00227                 }
00228         } else {
00229                 if(ix == 0) {
00230                         // interpolate between center of box and xMin
00231                         double Ex2=approxE(xMin,y,z).x();
00232                         double f=fx*2.0;
00233                         Ex = (1.0-f)*Ex2 + f*Ex;
00234                 } else {
00235                         // interpolate between this box and next lower one
00236                         double v200=array[index(ix-1,iy,iz)];
00237                         double v210=array[index(ix-1,iy+1,iz)];
00238                         double v201=array[index(ix-1,iy,iz+1)];
00239                         double v211=array[index(ix-1,iy+1,iz+1)];
00240                         double v2=gy*gz*v200 + fy*gz*v210 + gy*fz*v201 + 
00241                                                                 fy*fz*v211;
00242                         double f=fx+0.5;
00243                         Ex = -(1.0-f)*(v0-v2)/dx + f*Ex;
00244                 }
00245         }
00246 
00247         // partial d/dy
00248         v0=gx*gz*v000 + fx*gz*v100 + gx*fz*v001 + fx*fz*v101;
00249         v1=gx*gz*v010 + fx*gz*v110 + gx*fz*v011 + fx*fz*v111;
00250         double Ey=-(v1-v0)/dy; // for this grid box
00251         if(fy > 0.5) {
00252                 if(iy+1 == ny-1) {
00253                         // interpolate between center of box and yMax
00254                         double Ey2=approxE(x,yMax,z).y();
00255                         double f=2.0*fy-1.0;
00256                         Ey = (1.0-f)*Ey + f*Ey2;
00257                 } else {
00258                         // interpolate between this box and next higher one
00259                         double v020=array[index(ix,iy+2,iz)];
00260                         double v120=array[index(ix+1,iy+2,iz)];
00261                         double v021=array[index(ix,iy+2,iz+1)];
00262                         double v121=array[index(ix+1,iy+2,iz+1)];
00263                         double v2=gx*gz*v020 + fx*gz*v120 + gx*fz*v021 + 
00264                                                                 fx*fz*v121;
00265                         double f=fy-0.5;
00266                         Ey = (1.0-f)*Ey - f*(v2-v1)/dy;
00267                 }
00268         } else {
00269                 if(iy == 0) {
00270                         // interpolate between center of box and yMin
00271                         double Ey2=approxE(x,yMin,z).y();
00272                         double f=fy*2.0;
00273                         Ey = (1.0-f)*Ey2 + f*Ey;
00274                 } else {
00275                         // interpolate between this box and next lower one
00276                         double v020=array[index(ix,iy-1,iz)];
00277                         double v120=array[index(ix+1,iy-1,iz)];
00278                         double v021=array[index(ix,iy-1,iz+1)];
00279                         double v121=array[index(ix+1,iy-1,iz+1)];
00280                         double v2=gx*gz*v020 + fx*gz*v120 + gx*fz*v021 + 
00281                                                                 fx*fz*v121;
00282                         double f=fy+0.5;
00283                         Ey = -(1.0-f)*(v0-v2)/dy + f*Ey;
00284                 }
00285         }
00286 
00287         // partial d/dz
00288         v0=gx*gy*v000 + fx*gy*v100 + gx*fy*v010 + fx*fy*v110;
00289         v1=gx*gy*v001 + fx*gy*v101 + gx*fy*v011 + fx*fy*v111;
00290         double Ez=-(v1-v0)/dz; // for this grid box
00291         if(fz > 0.5) {
00292                 if(iz+1 == nz-1) {
00293                         // interpolate between center of box and zMax
00294                         double Ez2=approxE(x,y,zMax).z();
00295                         double f=2.0*fz-1.0;
00296                         Ez = (1.0-f)*Ez + f*Ez2;
00297                 } else {
00298                         // interpolate between this box and next higher one
00299                         double v002=array[index(ix,iy,iz+2)];
00300                         double v102=array[index(ix+1,iy,iz+2)];
00301                         double v012=array[index(ix,iy+1,iz+2)];
00302                         double v112=array[index(ix+1,iy+1,iz+2)];
00303                         double v2=gx*gy*v002 + fx*gy*v102 + gx*fy*v012 + 
00304                                                                 fx*fy*v112;
00305                         double f=fz-0.5;
00306                         Ez = (1.0-f)*Ez - f*(v2-v1)/dz;
00307                 }
00308         } else {
00309                 if(iz == 0) {
00310                         // interpolate between center of box and zMin
00311                         double Ez2=approxE(x,y,zMin).z();
00312                         double f=fz*2.0;
00313                         Ez = (1.0-f)*Ez2 + f*Ez;
00314                 } else {
00315                         // interpolate between this box and next lower one
00316                         double v002=array[index(ix,iy,iz-1)];
00317                         double v102=array[index(ix+1,iy,iz-1)];
00318                         double v012=array[index(ix,iy+1,iz-1)];
00319                         double v112=array[index(ix+1,iy+1,iz-1)];
00320                         double v2=gx*gy*v002 + fx*gy*v102 + gx*fy*v012 + 
00321                                                                 fx*fy*v112;
00322                         double f=fz+0.5;
00323                         Ez = -(1.0-f)*(v0-v2)/dz + f*Ez;
00324                 }
00325         }
00326 #endif
00327 
00328         return G4ThreeVector(Ex,Ey,Ez);
00329 }

G4ThreeVector PoissonConvolve3D::approxE ( double  x,
double  y,
double  z 
)

approxE() evaluates the approximation E field. {x,y,z} can be anywhere, but normally are outside the grid or on its boundary.

References approx, dx, dy, dz, E, and M_PI.

Referenced by Bunch::addFieldValue(), and getE().

00332 {
00333         double tol = ((dx<dy?dy:dx)<dz ? (dx<dy?dy:dx) : dz) * 0.05;
00334         G4ThreeVector v(0.0,0.0,0.0);
00335         G4ThreeVector pos(x,y,z);
00336         for(unsigned i=0; i<approx.size(); ++i) {
00337                 G4ThreeVector p = approx[i].getPosition();
00338                 p -= pos; // note - sign
00339                 double r=p.mag();
00340                 if(r < tol) continue; // should never happen
00341                 double denom=4.0*M_PI*r*r*r;
00342                 double c=approx[i].getCharge();
00343                 G4ThreeVector E(c*p.x()/denom,c*p.y()/denom,c*p.z()/denom);
00344                 v -= E; // note - sign above
00345         }
00346         return v;
00347 }

float PoissonConvolve3D::approxPhi ( double  x,
double  y,
double  z 
)

approxPhi() evaluates the approximation phi. {x,y,z} can be anywhere, but normally are outside the grid or on its boundary.

References approx, dx, dy, dz, M_PI, and phi().

00350 {
00351         double tol = ((dx<dy?dy:dx)<dz ? (dx<dy?dy:dx) : dz) * 0.05;
00352         double phi=0.0;
00353         G4ThreeVector pos(x,y,z);
00354         for(unsigned i=0; i<approx.size(); ++i) {
00355                 G4ThreeVector p = approx[i].getPosition();
00356                 p -= pos; // note - sign
00357                 double r=p.mag();
00358                 if(r < tol) r = tol; // should never happen
00359                 double c=approx[i].getCharge();
00360                 phi += c/r;
00361         }
00362         return phi/(4.0*M_PI);
00363 }

G4ThreeVector PoissonConvolve3D::getPosition ( int  ix,
int  iy,
int  iz 
) [inline]

getPosition() returns the position corresponding to grid point indexes. Intended primarily for testing.

References dx, dy, dz, xMin, yMin, and zMin.

00267                 { return G4ThreeVector(xMin+ix*dx,yMin+iy*dy,zMin+iz*dz); }

void PoissonConvolve3D::setApproximation ( std::vector< Charge >  v  )  [inline]

setApproximation() sets the approximation to use outside the grid, and on its boundary. v is a vector of point charges. Updates the boundary values of phi.

References approx.

Referenced by Bunch::computeField(), and Bunch::testBfield().

00273                 { approx = v; }

void PoissonConvolve3D::defaultApproximation (  ) 

defaultApproximation() calls setApproximation() with the total charge and mean position in rhs (from putCharge()). Updates the boundary values of phi, and the interior if interior=true.

References approx, sumX, sumY, sumZ, and totalCharge.

Referenced by Bunch::testEfield().

00366 {
00367         Charge c(sumX/totalCharge,sumY/totalCharge,sumZ/totalCharge,
00368                                                                 totalCharge);
00369         approx.clear();
00370         approx.push_back(c);
00371 }

void PoissonConvolve3D::computeGreensFunction (  ) 

computeGreensFunction() does just that. Must be called whenever the physical size of the grid is changed.

References array, dx, dy, dz, gfPlan, index(), nx, ny, nz, status, and UNINITIALIZED.

Referenced by updatePosition().

00374 {
00375         // fill in the physical region
00376         // include the 1/N overall normalization
00377         // (Need one more sample than usual, for the symmetry.)
00378         float norm=0.125f/(float)(nx*ny*nz);
00379         for(int ix=0; ix<=nx; ++ix) {
00380                 float x2=ix*dx*ix*dx;
00381                 for(int iy=0; iy<=ny; ++iy) {
00382                         float y2=iy*dy*iy*dy;
00383                         for(int iz=0; iz<=nz; ++iz) {
00384                                 float z2=iz*dz*iz*dz;
00385                                 float r=sqrtf(x2+y2+z2);
00386                                 if(ix==0 && iy==0 && iz==0)
00387                                         r = (dx+dy+dz)/3.0;
00388                                 array[index(ix,iy,iz)] = norm/r;
00389                         }
00390                 }
00391         }
00392 
00393         // now fill in the un-physical region via symmetry
00394         for(int ix=0; ix<2*nx; ++ix) {
00395                 int jx=ix;
00396                 if(jx > nx) jx=2*nx-jx;
00397                 assert(jx>=0 && jx<=nx);
00398                 for(int iy=0; iy<2*ny; ++iy) {
00399                         int jy=iy;
00400                         if(jy > ny) jy=2*ny-jy;
00401                         assert(jy>=0 && jy<=ny);
00402                         for(int iz=0; iz<2*nz; ++iz) {
00403                                 if(ix<=nx && iy<=ny && iz<=nz) continue;
00404                                 int jz=iz;
00405                                 if(jz > nz) jz=2*nz-jz;
00406                                 assert(jz>=0 && jz<=nz);
00407                                 array[index(ix,iy,iz)] = array[index(jx,jy,jz)];
00408                         }
00409                 }
00410         }
00411 
00412         // compute the FFT of the Green's function
00413         fftwf_execute(gfPlan);
00414 
00415         // indicate we clobbered array[]
00416         status = UNINITIALIZED;
00417 }

void PoissonConvolve3D::solve (  ) 

solve() will solve Poisson's equation

References CHARGE, greensFunction, nc, PHI, plan1, plan2, status, and transf.

Referenced by Bunch::computeField(), Bunch::testBfield(), and Bunch::testEfield().

00420 {
00421         assert(status==CHARGE);
00422 
00423         // compute the FFT of array[], which contains the charge,
00424         // into transf[].
00425         fftwf_execute(plan1);
00426 
00427         // multiply by the FFT of the Green's function
00428         for(int i=0; i<nc; ++i) {
00429                 double tmp = transf[i][0]*greensFunction[i][0] -
00430                                         transf[i][1]*greensFunction[i][1];
00431                 transf[i][1] = transf[i][1]*greensFunction[i][0] +
00432                                         transf[i][0]*greensFunction[i][1];
00433                 transf[i][0] = tmp;
00434         }
00435 
00436         // inverse transform transf[] to get phi into array[]
00437         fftwf_execute(plan2); // transf[] is overwritten, but don't care
00438 
00439         // indicate a valid phi in array[]
00440         status = PHI;
00441 }

void PoissonConvolve3D::printRhs (  ) 

printRhs() prints the rhs. Intended primarily for testing.

References array, CHARGE, index(), nx, ny, nz, and status.

00444 {
00445         assert(status==CHARGE);
00446         for(int ix=0; ix<nx; ++ix) {
00447                 for(int iy=0; iy<ny; ++iy) {
00448                         for(int iz=0; iz<nz; ++iz) {
00449                                 int k=index(ix,iy,iz);
00450                                 if(array[k] == 0.0) continue;
00451                                 printf("rhs[%d] = %.3f  ix=%d iy=%d iz=%d\n",
00452                                         k,array[k],ix,iy,iz);
00453                         }
00454                 }
00455         }
00456 }

void PoissonConvolve3D::printPhi (  ) 

printPhi() prints phi. Intended primarily for testing.

References array, index(), nx, ny, nz, PHI, and status.

00459 {
00460         assert(status==PHI);
00461         for(int ix=0; ix<nx; ++ix) {
00462                 for(int iy=0; iy<ny; ++iy) {
00463                         for(int iz=0; iz<nz; ++iz) {
00464                                 int k=index(ix,iy,iz);
00465                                 printf("phi[%d] = %.3f  ix=%d iy=%d iz=%d\n",
00466                                         k,array[k],ix,iy,iz);
00467                         }
00468                 }
00469         }
00470 }

int PoissonConvolve3D::getNx (  )  const [inline]

accessors.

References nx.

00297 { return nx; }

int PoissonConvolve3D::getNy (  )  const [inline]

References ny.

00298 { return ny; }

int PoissonConvolve3D::getNz (  )  const [inline]

References nz.

00299 { return nz; }

Member Data Documentation

int PoissonConvolve3D::nx [protected]

Referenced by computeGreensFunction(), getE(), getNx(), phi(), PoissonConvolve3D(), printPhi(), printRhs(), putCharge(), updatePosition(), and zeroRhs().

int PoissonConvolve3D::ny [protected]

Referenced by computeGreensFunction(), getE(), getNy(), index(), phi(), PoissonConvolve3D(), printPhi(), printRhs(), putCharge(), updatePosition(), and zeroRhs().

int PoissonConvolve3D::nz [protected]

Referenced by computeGreensFunction(), getE(), getNz(), index(), phi(), PoissonConvolve3D(), printPhi(), printRhs(), putCharge(), updatePosition(), and zeroRhs().

float PoissonConvolve3D::xMin [protected]

Referenced by getE(), getPosition(), phi(), putCharge(), and updatePosition().

float PoissonConvolve3D::xMax [protected]

Referenced by getE(), and updatePosition().

float PoissonConvolve3D::yMin [protected]

Referenced by getE(), getPosition(), phi(), putCharge(), and updatePosition().

float PoissonConvolve3D::yMax [protected]

Referenced by getE(), and updatePosition().

float PoissonConvolve3D::zMin [protected]

Referenced by getE(), getPosition(), phi(), putCharge(), and updatePosition().

float PoissonConvolve3D::zMax [protected]

Referenced by getE(), and updatePosition().

Status PoissonConvolve3D::status [protected]

Referenced by addPhi(), computeGreensFunction(), getE(), phi(), phiGrid(), PoissonConvolve3D(), printPhi(), printRhs(), putCharge(), rhs(), setPhi(), setRhs(), solve(), and zeroRhs().

float PoissonConvolve3D::dx [protected]

Referenced by approxE(), approxPhi(), computeGreensFunction(), getE(), getPosition(), phi(), putCharge(), and updatePosition().

float PoissonConvolve3D::dy [protected]

Referenced by approxE(), approxPhi(), computeGreensFunction(), getE(), getPosition(), phi(), putCharge(), and updatePosition().

float PoissonConvolve3D::dz [protected]

Referenced by approxE(), approxPhi(), computeGreensFunction(), getE(), getPosition(), phi(), putCharge(), and updatePosition().

float* PoissonConvolve3D::array [protected]

Referenced by addPhi(), computeGreensFunction(), getE(), phi(), phiGrid(), PoissonConvolve3D(), printPhi(), printRhs(), putCharge(), rhs(), setPhi(), setRhs(), zeroRhs(), and ~PoissonConvolve3D().

fftwf_complex* PoissonConvolve3D::transf [protected]

Referenced by PoissonConvolve3D(), solve(), and ~PoissonConvolve3D().

fftwf_complex* PoissonConvolve3D::greensFunction [protected]

Referenced by PoissonConvolve3D(), solve(), and ~PoissonConvolve3D().

std::vector<Charge> PoissonConvolve3D::approx [protected]

Referenced by approxE(), approxPhi(), defaultApproximation(), phi(), setApproximation(), and ~PoissonConvolve3D().

float PoissonConvolve3D::sumX [protected]

Referenced by defaultApproximation(), PoissonConvolve3D(), putCharge(), and zeroRhs().

float PoissonConvolve3D::sumY [protected]

Referenced by defaultApproximation(), PoissonConvolve3D(), putCharge(), and zeroRhs().

float PoissonConvolve3D::sumZ [protected]

Referenced by defaultApproximation(), PoissonConvolve3D(), putCharge(), and zeroRhs().

float PoissonConvolve3D::totalCharge [protected]

Referenced by defaultApproximation(), PoissonConvolve3D(), putCharge(), and zeroRhs().

float PoissonConvolve3D::vol [protected]

Referenced by updatePosition().

PutMode PoissonConvolve3D::putMode [protected]

Referenced by PoissonConvolve3D(), putCharge(), and setPutMode().

int PoissonConvolve3D::nc [protected]

Referenced by PoissonConvolve3D(), and solve().

fftwf_plan PoissonConvolve3D::gfPlan [protected]

Referenced by computeGreensFunction(), PoissonConvolve3D(), and ~PoissonConvolve3D().

fftwf_plan PoissonConvolve3D::plan1 [protected]

Referenced by PoissonConvolve3D(), solve(), and ~PoissonConvolve3D().

fftwf_plan PoissonConvolve3D::plan2 [protected]

Referenced by PoissonConvolve3D(), solve(), and ~PoissonConvolve3D().

The documentation for this class was generated from the following files:
g4beamline