fxy.F



[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
 

      subroutine fxyinit
      use FXYversion

   Called at first reference to package (not nec. a "run" etc.).


      call fxyvers (STDOUT)

      return
      end

[fxyinit]
      subroutine fxyvers (iout)
      use FXYversion
      integer(ISZ):: iout
   Echoes code version,etc. to output files when they're created
      call printpkgversion(iout,"Fieldsolver FXY",versfxy)
      return
      end

      subroutine fxygen()

      return
      end

      subroutine fxyfin()
      use FXYversion

   Deallocates dynamic storage used by test driver


      call gfree ("FXYvars")

      return
      end
                                                                           =
  Capacity matrix solver                                                   =
                                                                           =

[capmatxyf]
      subroutine cmsetxy(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 anything
      if (ncxy == 0) then
        call remark("NOTICE: the capacity matrix field solver is turned on")
        call remark("        but no conductors have been defined.")
        return
      endif

      --- Make a copy of phi since it will be scrambled below.
      phisave = phi

  find capacity matrix

      --- find inverse capacity matrix
      do i=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.
        call vpois2d(2,phi,phi,kxsq,kysq,attx,atty,filt,
     &               xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &               scrtch,xywork,0,l2symtry,l4symtry)
        call vpois2d(4,phi,phi,kxsq,kysq,attx,atty,filt,
     &               xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &               scrtch,xywork,0,l2symtry,l4symtry)
        call vpois2d(5,phi,phi,kxsq,kysq,attx,atty,filt,
     &               xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &               scrtch,xywork,0,l2symtry,l4symtry)

        --- fill matrix with phi
        do j=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 *    wy
        enddo

      enddo

      --- 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.
#if WORDSIZE == 64
      call ssytrf("U",ncxy,cmatxy,ncxy,kpvtxy,xywork,(1+nx)*(1+ny),info)
#else
#ifdef CYGWIN
      call dsytrf_("U",ncxy,cmatxy,ncxy,kpvtxy,xywork,(1+nx)*(1+ny),info)
#else
      call dsytrf("U",ncxy,cmatxy,ncxy,kpvtxy,xywork,(1+nx)*(1+ny),info)
#endif
#endif
      if (info .ne. 0) then
        print*,info
        call kaboom("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.
#if WORDSIZE == 64
      call ssytri("U",ncxy,cmatxy,ncxy,kpvtxy,pcond,info)
#else
#ifdef CYGWIN
      call dsytri_("U",ncxy,cmatxy,ncxy,kpvtxy,pcond,info)
#else
      call dsytri("U",ncxy,cmatxy,ncxy,kpvtxy,pcond,info)
#endif
#endif

      --- Fill in lower half (which is the same as the upper half)
      do j=2,ncxy
        do i=1,j-1
          cmatxy(j,i) = cmatxy(i,j)
        enddo
      enddo

      --- Return the original saved phi.
      phi = phisave

      return
      end

[vpxy]
      subroutine capmatxyf(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) then

        call vpois2d(1,phi,phi,kxsq,kysq,attx,atty,filt,
     &               xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &               scrtch,xywork,0,l2symtry,l4symtry)

        call cmsetxy(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.
        call vpois2d(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.
      call vpois2d(2,phi,phi,kxsq,kysq,attx,atty,filt,
     &             xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &             scrtch,xywork,0,l2symtry,l4symtry)
      call vpois2d(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.
      call vpois2d(4,phi,phi,kxsq,kysq,attx,atty,filt,
     &             xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &             scrtch,xywork,0,l2symtry,l4symtry)
      call vpois2d(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 qcond
      do i=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 charge
      do j=1,ncxy
        do i=1,ncxy
          qcond(i) = qcond(i) + pcond(j)*cmatxy(i,j)
        enddo
      enddo

      --- Zero out phi (the source is added back in later).
      call zeroarry(phi,(nx+1)*(ny+1))
      phi = 0.

      --- Deposit induced charge onto grid
      do i=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 *    wy
      enddo

      --- 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.
      call vpois2d(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.
      do i=0,nx
        do j=0,ny
          phi(i,j) = phisave(i,j) + phi(i,j)
        enddo
      enddo

      --- Complete the field solve.
      call vpois2d(4,phi,phi,kxsq,kysq,attx,atty,filt,
     &             xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &             scrtch,xywork,0,l2symtry,l4symtry)
      call vpois2d(5,phi,phi,kxsq,kysq,attx,atty,filt,
     &             xlen,ylen,nx,ny,nxguardphi,nyguardphi,
     &             scrtch,xywork,0,l2symtry,l4symtry)

      return
      end