1 -- 
    2 --      Patricia Fasel
    3 --      Los Alamos National Laboratory
    4 --      1990 August
    5 --
    6 module Potential (potential) where
    7 
    8 import  PicType
    9 import  Consts
   10 import  Utils
   11 import Array--1.3
   12 
   13 -- Given charge density matrix, rho
   14 -- Compute new electrostatic potential phi' where del2(phi') = rho
   15 -- phi from the previous timestep is used as the initial value
   16 -- assume:      phi = phi' + error
   17 -- then:        d_phi = Laplacian(phi) = Laplacian(phi') + Laplacian(error)
   18 --           d_error = d_phi - rho = Laplacian(error)
   19 --          error' = InvLaplacian(d_error) = InvLaplacian(Laplacian(error))
   20 --          phi' = phi - error'
   21 
   22 potential :: Phi -> Rho -> Indx -> Indx -> Phi
   23 potential phi rho depth nIter
   24     | nIter == 0        = phi
   25     | otherwise         = potential phi' rho depth (nIter-1)
   26                           where
   27                               phi' = vCycle rho phi nCell depth
   28 
   29 
   30 -- vCycle is a multigrid laplacian inverse
   31 -- Given d_phi, find phi where Laplacian(phi) = d_phi
   32 -- Algorithm is to invert d_phi on a course mesh and interpolate to get phi
   33 
   34 vCycle :: Phi -> Rho -> Indx -> Indx -> Phi
   35 vCycle phi rho n depth =
   36         if (depth == 0)
   37             then relax phi' rho n
   38             else correct phi' eCoarse n nHalf
   39         where
   40             nHalf = n `div` 2
   41             nHalf' = nHalf-1
   42             phi' = relax phi rho n
   43             rho' = residual phi' rho n
   44             rCoarse = coarseMesh rho' n
   45             eZero = array ((0,0), (nHalf',nHalf')) 
   46                         [((i,j), 0.0) | i<-[0..nHalf'], j<-[0..nHalf']]
   47             eCoarse = vCycle eZero rCoarse nHalf (depth-1)
   48 
   49 
   50 -- laplacian operator
   51 -- mesh configuration where e=(i,j) position, b + d + f + h - 4e
   52 -- a   b   c
   53 -- d   e   f
   54 -- g   h   i
   55  
   56 laplacianOp :: Mesh -> Range -> Value
   57 laplacianOp mesh [iLo, i, iHi, jLo, j, jHi] =
   58         -(mesh!(iLo,j)+mesh!(i,jHi)+mesh!(i,jLo)+mesh!(iHi,j)-4*mesh!(i,j))
   59 
   60 
   61 -- subtract laplacian of mesh from mesh'
   62 -- residual = mesh' - Laplacian(mesh)
   63 
   64 residual :: Phi -> Rho -> Indx -> Rho
   65 residual mesh mesh' n =
   66         applyOpToMesh (residualOp mesh') mesh n
   67         where
   68             residualOp mesh' mesh [iLo, i, iHi, jLo, j, jHi] =
   69                 mesh'!(i,j) - laplacianOp mesh [iLo, i, iHi, jLo, j, jHi]
   70 
   71 
   72 relax :: Phi -> Rho -> Indx -> Phi
   73 relax mesh mesh' n =
   74         applyOpToMesh (relaxOp mesh') mesh n
   75         where
   76             relaxOp mesh' mesh [iLo, i, iHi, jLo, j, jHi] =
   77                 0.25 * mesh'!(i,j) + 
   78                 0.25 * (mesh!(iLo,j)+mesh!(i,jLo)+mesh!(i,jHi)+mesh!(iHi,j))
   79 
   80 correct :: Phi -> Mesh -> Indx -> Indx -> Phi
   81 correct phi eCoarse n' nHalf =
   82         array ((0,0), (n,n))
   83         [((i,j) , phi!(i,j) + eFine!(i,j)) | i <- [0..n], j <- [0..n]]
   84         where
   85             eFine = fineMesh eCoarse nHalf
   86             n = n'-1