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