[capmatxyf] [cmsetxy] [fxyfin] [fxygen] [fxyinit] [fxyvers]

#include top.h @(#) File FXY.M, version $Revision: 3.8 $, $Date: 2011/08/27 00:35:50 $ # Copyright (c) 1990-1998, The Regents of the University of California. # All rights reserved. See LEGAL.LLNL for full text and disclaimer. written by David P. Grote It is the source file for the package FXY of the PIC code WARPxy, but it may be useful by itself. It contains: Capacity matrix solver which works in xy space

[fxyinit]subroutinefxyinituse FXYversion Called at first reference to package (not nec. a "run" etc.).callfxyvers (STDOUT)returnend

subroutinefxyvers(iout) use FXYversion integer(ISZ):: iout Echoes code version,etc. to output files when they're createdcallprintpkgversion(iout,"Fieldsolver FXY",versfxy)returnend

subroutinefxygen()returnend

[capmatxyf]subroutinefxyfin() use FXYversion Deallocates dynamic storage used by test drivercallgfree("FXYvars")returnend= Capacity matrix solver = =

[vpxy]subroutinecmsetxy(phi,kxsq,kysq,attx,atty,filt,xlen,ylen,nx,ny, & nxguardphi,nyguardphi,dx,dy, & xmmin,ymmin,scrtch,phisave,xywork,l2symtry,l4symtry) use CapMatxy integer(ISZ):: nx,ny,nxguardphi,nyguardphi real(kind=8):: phi(-nxguardphi:nx+nxguardphi, & -nyguardphi:ny+nyguardphi) real(kind=8):: kxsq(0:nx-1),kysq(0:ny-1) real(kind=8):: attx(0:nx-1),atty(0:ny-1) real(kind=8):: filt(5,2) real(kind=8):: scrtch(*) real(kind=8):: phisave(-nxguardphi:nx+nxguardphi, & -nyguardphi:ny+nyguardphi) real(kind=8):: xywork(0:nx,0:ny) real(kind=8):: xlen,ylen,dx,dy,xmmin,ymmin logical(ISZ):: l2symtry,l4symtry Sets up capacity matrix for xy fieldsolver. Note that the capacity matrix is calculated using field solves which do not do any filtering. (Experience has shown that when filtering is used, the inverse capacity matrix is very ill-conditioned and so the capacity matrix becomes garbled.) real(kind=8):: wx,wy,tt1,tt2,cx,cy integer(ISZ):: ix,iy,i,j,info --- If the number of conductor points is zero, print a warning --- and return without doing anythingif (ncxy == 0) thencallremark("NOTICE: the capacity matrix field solver is turned on")callremark(" but no conductors have been defined.")returnendif--- Make a copy of phi since it will be scrambled below. phisave = phi find capacity matrix --- find inverse capacity matrixdoi=1,ncxy --- start with zero phi call zeroarry(phi,(nx+1)*(ny+1)) phi = 0. --- set phi with 1 Coulomb at conductor points (charge is scaled by --- 2 or 4 at a symmetric boundary) ix = int((xcond(i)-xmmin)/dx) wx = (xcond(i)-xmmin)/dx - ix iy = int((ycond(i)-ymmin)/dy) wy = (ycond(i)-ymmin)/dy - iy cx = 1. cy = 1.if (l2symtry.and.iy == 0) cy = 2.if (l4symtry.and.ix == 0) cx = 2.if (l4symtry.and.iy == 0) cy = 2. phi(ix ,iy ) = (1.-wx)*(1.-wy)*cx*cy phi(ix+1,iy ) = wx *(1.-wy)*cy phi(ix ,iy+1) = (1.-wx)* wy *cx phi(ix+1,iy+1) = wx * wy --- Solve for fields, piece-by-piece so that the filtering can be --- skipped.callvpois2d(2,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)callvpois2d(4,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)callvpois2d(5,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry) --- fill matrix with phidoj=1,ncxy ix = int((xcond(j)-xmmin)/dx) wx = (xcond(j)-xmmin)/dx - ix iy = int((ycond(j)-ymmin)/dy) wy = (ycond(j)-ymmin)/dy - iy cmatxy(j,i) = phi(ix ,iy )*(1.-wx)*(1.-wy) + & phi(ix+1,iy )* wx *(1.-wy) + & phi(ix ,iy+1)*(1.-wx)* wy + & phi(ix+1,iy+1)* wx * wyenddoenddo--- invert to get capacity matrix (checking for singular matrix) --- These routines don't seem to give as good a result as those below. call ssifa(cmatxy,ncxy,ncxy,kpvtxy,info) if (info .ne. 0) then call kaboom("cmsetxy: ERROR: Capacity matrix is singular - it is likely that & two or more conductor points are within the same grid cell.") endif call ssidi(cmatxy,ncxy,ncxy,kpvtxy,tt1,tt2,xywork,001) --- These routines are from LAPACK.#ifWORDSIZE == 64callssytrf("U",ncxy,cmatxy,ncxy,kpvtxy,xywork,(1+nx)*(1+ny),info)#else#ifdefCYGWINcalldsytrf_("U",ncxy,cmatxy,ncxy,kpvtxy,xywork,(1+nx)*(1+ny),info)#elsecalldsytrf("U",ncxy,cmatxy,ncxy,kpvtxy,xywork,(1+nx)*(1+ny),info)#endif#endifif (info.ne.0) thenprint*,infocallkaboom("cmsetxy: ERROR: Capacity matrix is singular - it is likely that & two or more conductor points are within the same grid cell.")endif--- pcond is passed as workspace.#ifWORDSIZE == 64callssytri("U",ncxy,cmatxy,ncxy,kpvtxy,pcond,info)#else#ifdefCYGWINcalldsytri_("U",ncxy,cmatxy,ncxy,kpvtxy,pcond,info)#elsecalldsytri("U",ncxy,cmatxy,ncxy,kpvtxy,pcond,info)#endif#endif--- Fill in lower half (which is the same as the upper half)doj=2,ncxydoi=1,j-1 cmatxy(j,i) = cmatxy(i,j)enddoenddo--- Return the original saved phi. phi = phisavereturnend

subroutinecapmatxyf(iwhich,phi,kxsq,kysq,attx,atty,filt,xlen,ylen,nx,ny, & nxguardphi,nyguardphi, & dx,dy,xmmin,ymmin,scrtch,phisave,xywork, & l2symtry,l4symtry) use CapMatxy integer(ISZ):: iwhich,nx,ny,nxguardphi,nyguardphi real(kind=8):: phi(-nxguardphi:nx+nxguardphi, & -nyguardphi:ny+nyguardphi) real(kind=8):: kxsq(0:nx-1),kysq(0:ny-1) real(kind=8):: attx(0:nx-1),atty(0:ny-1) real(kind=8):: filt(5,2) real(kind=8):: scrtch(*) real(kind=8):: phisave(-nxguardphi:nx+nxguardphi, & -nyguardphi:ny+nyguardphi) real(kind=8):: xywork(0:nx,0:ny) real(kind=8):: xlen,ylen,dx,dy,xmmin,ymmin logical(ISZ):: l2symtry,l4symtry Applies the capacity matrix for xy code electrostatic field solver. Multiplies phi by cap to get induced rho. Assumes that beam's rho has been set, and copied into phi. When filtering is being done, care is taken so that only the beam's charge density is filtered and not the induced charge on the conductor. Note that fieldsolve is set up now so that the rho array is not used, since phi (which holds rho initially) is saved after the first FFT. This allows the capacity matrix to be used in bends. real(kind=8):: wx,wy,cx,cy integer(ISZ):: ix,iy,i,j --- Initialize the arrays for poisson solve and the capacity matrix.if (iwhich == 1.or.iwhich == 0) thencallvpois2d(1,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)callcmsetxy(phi,kxsq,kysq,attx,atty,filt,xlen,ylen,nx,ny, & nxguardphi,nyguardphi,dx,dy, & xmmin,ymmin,scrtch,phisave,xywork,l2symtry,l4symtry)elseif(iwhich > 1) then--- Make call to do the specialized action requested.callvpois2d(iwhich,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)endif--- If a full field solve was not desired, return.if (iwhich > 0)return--- Do the first field solve (note that rho must have been copied into --- phi for this to work) --- This is done piece-by-piece so that when filtering is done, the --- filtered source can be saved. --- First transform phi (actually the source) and filter it.callvpois2d(2,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)callvpois2d(3,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry) --- Now, save the transformed and possibly filtered source. phisave = phi --- Finish the first field solve: divide by k-squared and untransform.callvpois2d(4,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)callvpois2d(5,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry) --- Extract phi error from FFT solution and zero qconddoi=1,ncxy ix = int((xcond(i)-xmmin)/dx) wx = (xcond(i)-xmmin)/dx - ix iy = int((ycond(i)-ymmin)/dy) wy = (ycond(i)-ymmin)/dy - iy pcond(i) = vcond(i) - (phi(ix ,iy )*(1.-wx)*(1.-wy) + & phi(ix+1,iy )* wx *(1.-wy) + & phi(ix ,iy+1)*(1.-wx)* wy + & phi(ix+1,iy+1)* wx * wy) qcond(i) = 0.enddo--- Multiply pcond by capacity matrix to get induced chargedoj=1,ncxydoi=1,ncxy qcond(i) = qcond(i) + pcond(j)*cmatxy(i,j)enddoenddo--- Zero out phi (the source is added back in later). call zeroarry(phi,(nx+1)*(ny+1)) phi = 0. --- Deposit induced charge onto griddoi=1,ncxy ix = int((xcond(i)-xmmin)/dx) wx = (xcond(i)-xmmin)/dx - ix iy = int((ycond(i)-ymmin)/dy) wy = (ycond(i)-ymmin)/dy - iy cx = 1. cy = 1.if (l2symtry.and.iy == 0) cy = 2.if (l4symtry.and.ix == 0) cx = 2.if (l4symtry.and.iy == 0) cy = 2. phi(ix ,iy ) = phi(ix ,iy ) + qcond(i)*(1.-wx)*(1.-wy)*cx*cy phi(ix+1,iy ) = phi(ix+1,iy ) + qcond(i)* wx *(1.-wy)*cy phi(ix ,iy+1) = phi(ix ,iy+1) + qcond(i)*(1.-wx)* wy *cx phi(ix+1,iy+1) = phi(ix+1,iy+1) + qcond(i)* wx * wyenddo--- Redo the field solve, again piece-by-piece so that the already --- transformed and filtered source may be added and so that the --- filtering can be skipped. Note that the initial transform is --- only applied to the induced charge.callvpois2d(2,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry) --- Sum the induced and transformed beam charge.doi=0,nxdoj=0,ny phi(i,j) = phisave(i,j) + phi(i,j)enddoenddo--- Complete the field solve.callvpois2d(4,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)callvpois2d(5,phi,phi,kxsq,kysq,attx,atty,filt, & xlen,ylen,nx,ny,nxguardphi,nyguardphi, & scrtch,xywork,0,l2symtry,l4symtry)returnend