GRIN PART I 1 GRIN GReen’s function INtegrals Irek Szczesniak and Alexander Pletzer Princeton Plasma Physics Laboratory September 2001
GRIN PART I 2 PART I Green’s function method and the GRIN program • Green’s function method • Examples • The GRIN program
GRIN PART I 3 Green’s Function Method The Green’s function is the response to a Dirac excitation. Suppose we want to solve: ∇ 2 ψ = − 4 πρ ( � x ) and we know the Green’s function G satisfying: ∇ 2 G ( � x − � x ′ , � x ) = − 4 πδ ( � x ′ ) which corresponds to the potential due to a point charge. Then the solution is: � ψ ( � G ( � x ′ ) = x ) ρ ( � x ) dV x ′ , � V
GRIN PART I 4 Green’s Second Identity Central to the Green’s function method is Green’s second identity � � dS · { G ( � � x ) − ∇ G ( � x ) } − 4 παψ ( � dV G ( � x ′ , � x ) ∇ ψ ( � x ′ , � x ) ψ ( � x ′ ) = − 4 π x ′ , � x ) ρ ( � x ) S V which is obtained after multiplying ∇ 2 ψ = − 4 πρ ( � x ) by G , integrating over the volume V, and using the definition of G. Here, if � x ′ ∈ V 1 if � x ′ ∈ S 1 where α = 2 if � x ′ / 0 ∈ V ∪ S
GRIN PART I 5 Boundary Conditions The intuitive solution (sum of all source contributions) � ψ ( � x, � x ′ ) = G ( � x ′ ) ρ ( � x ) dV V is a special case of the general solution � � dS · { G ∇ ( � � x ) − ∇ G ( � x ) } − 4 παψ ( � dV G ( � x ′ , � x ) ψ ( � x ′ , � x ) ψ ( � x ′ ) = − 4 π x ′ , � x ) ρ ( � x ) S V which also takes into account the effect of inhomogeneous boundary conditions.
GRIN PART I 6 First Example Given: G ( x, y, x ′ , y ′ ) = − log (( x − x ′ ) 2 + ( y − y ′ ) 2 ) Find the potential ψ due to a source s ( t ) = cos ( t ) distributed on a ring: x = 3 cos ( t ) , t ∈ (0 , 2 π ) y = 3 sin ( t ) Solution: � � 2 π 2 2 dx ( t ) + dy ( t ) y ( x ′ , y ′ ) = G [ x ′ , y ′ , x ( t ) , y ( t )] s ( t ) dt dt dt 0
GRIN PART I 7 Solution to the First Example 10 10 0 5 -10 10 10 0 5 5 -5 0 0 -5 -5 -10 -10 -10 The figure depicts potential ψ . ψ ∼ r inside ψ ∼ 1 r outside
GRIN PART I 8 Second Example Given: U ( θ ) for r = b ψ ( θ ) = 0 for r = a n -n Solve: a b ∇ 2 ψ = 0 in the shaded region ( a < r < b ) to find the normal electric field E on r = a .
GRIN PART I 9 Second Example (continued) Green’s second identity for the example reads: E b ( θ ) − ∂G E a ( θ )+ ∂G � � bdθ {− G ( � x ′ , θ ) � adθ {− G ( � x ′ , θ ) � ∂n 0 }− 4 παψ ( � ∂n U ( θ ) }− x ′ ) = 0 b a x ′ on a : For � bdθ ∂G ( � x a ( θ ′ ) , θ ) � � � x a ( θ ′ ) , θ ) E b ( θ ) + x a ( θ ′ ) , θ ) E a ( θ ) bdθG ( � U ( θ ) = adθG ( � ∂n b b a x ′ on b : For � x b ( θ ′ ) , θ ) bdθ { ∂G ( � − 2 π � � b δ ( θ ′ − θ ) } U ( θ ) = x b ( θ ′ ) , θ ) E b ( θ ) + bdθG ( � ∂n b b � x b ( θ ′ ) , θ ) E a ( θ ) adθG ( � a
GRIN PART I 10 Second Example (continued) Expand � E ( i ) E a = a e i ( θ ) i � E ( i ) E b = b e i ( θ ) i � U ( i ) e i ( θ ) U = i in basis functions e i ( θ ). We then get a coupled linear system of equations: A · E b + B · U = C · E a D · E b + E · U = A · E a After elimination of E b , a linear relation between E a and U can be obtained, and the problem is solved.
GRIN PART I 11 Punch Line To solve the above type of problems, we need the capability to: • compute accurately line integrals involving a kernel with a log singularity, • solve dense linear system of equations.
GRIN PART I 12 GRIN The GRIN code computes integrals of two kinds: � K ( l ′ , l ) ρ ( l ) dl C � � C ′ ν ( l ′ ) K ( l ′ , l ) ρ ( l ) dldl ′ C with K ( l ′ , l ) ∼ log ( l ′ − l ) as l ′ → l if C = C ′ (true for elliptic operators in 2-D).
GRIN PART I 13 GRIN vs. VACUUM VACUUM is a highly Successful code written by M. Chance, which is used to compute the natural boundary conditions in a number of stability codes (PEST, DCON, GATO, ...). Feature VACUUM GRIN portability CRAY, alpha a UNIX max. # of contours hardcoded hardware limited geometry of contours hardcoded almost arbitrary choice of kernel 1 (hardcoded) arbitrary with log singularity basis function Fourier, finite elements user supplied a ongoing work to port to other Unixes
GRIN PART I 14 Green’s functions provided by GRIN There are presently 10 kernel functions ( G or ∂G ∂n ) that come with GRIN. An example is the toroidally averaged Green’s function for the Laplace equation (following Chance [ Phys. Plasmas 4 , 2161 (1997)]): � 1 n + 1 2 , 1 � � � 4 F 2 , n + 1 | p 1 2 2 + 1 n Γ( 1 G n ( R, Z ; R ′ Z ′ ) = 2 + n ) p Γ( n + 1) πR ′ R λ ≡ 1 + ( R ′ − R ) 2 + ( Z ′ − Z ) 2 p ≡ s − 1 λ √ s ≡ s + 1 , , λ 2 − 1 2 R ′ R where F ( a, b, c | p ) is the Gauss hypergeometric function .
GRIN PART I 15 More that comes with GRIN GRIN has vector and matrix classes implemented. The matrix class is interfaced to Lapack. C++ example of matrix operations: Mat Temp = MatMult(kmat_ab, Inverse(kmat_bb)); Vec chi_a = MatMult( MatMult( Inverse(kmat_aa - MatMult(Temp, kmat_ba)), gmat_aa - MatMult(Temp, gmat_ba) ), dchin_a);
GRIN PART I 16 Example of GRIN usage Define the geometry (segments) /* Define segments */ Vec tb(11), ta(5); segment sa, sb; ta.space(0., 1.); tb.space(0., TwoPi); Vec xa(5); xa = 2.; Vec ya(5); ya.space(-0.1, 0.1); sa.load(ta, xa, xb, not_periodic); sb.load(tb, xb, yb, periodic);
GRIN PART I 17 Example of GRIN usage (continued) � Compute dℓ ′ K ( ℓ, ℓ ′ ) α ′ ( ℓ ′ ) /* set observer point */ double tobs = 0.5; /* call procedure */ /* greenLaplaceCartesian = Green’s fct alpha is the basis fct (ext proc) */ double res = greenIntegral( greenLaplaceCartesian, sa, tobs, sb, alpha);
GRIN PART I 18 Example of GRIN usage (continued) � � Compute s dℓα ( ℓ ) s ′ dℓ ′ K ( ℓ, ℓ ′ ) α ′ ( ℓ ′ ) /* beta, alpha are the basis fcts */ double res = twoGreenIntegral( greenLaplaceCartesian, sb, beta, sa, alpha);
GRIN PART I 19 GRIN’s NICHE Reasons to develop GRIN There is a gap between general software as Mathematica and specialized software as VACUUM. This gap is filled in by GRIN. Feature Mathematica 3 GRIN VACUUM Accuracy depends very good very good Flexibility very good good poor Performance poor good very good Friendliness very good good poor Programming from scratch tools provided hard to change
GRIN PART I 20 Inaccurate Results of Mathematica 3 Reasons to develop GRIN GRIN vs Mathematica 0 θ′ −T + ∫ θ′ −T θ′ θ′ +T + ∫ θ′ +T 2 π ∫ 0 + ∫ θ′ −0.1 θ ’ + ∫ θ ’ 2 π ∫ 0 −0.2 ∫ d θ ∂ G(n=2)/ ∂ n 2 π ∫ 0 −0.3 R 0 =1 a=0.95 κ =1 δ =0.95 −0.4 GRIN −0.5 0 1 2 3 4 5 6 θ′ (observer) Discrepancy between GRIN’s and Mathematica 3 results.
GRIN PART I 21 Inaccurate Results of Mathematica 3 (continued) � 2 π The circles represent the value of dθK ( n = 3) returned by the 0 NIntegrate Mathematica function for various observer coordinates. The triangles are those returned by NIntegrate when splitting the interval in two. The crosses (x) are those values returned by NIntegrate after further splitting; these agree much better with the GRIN results obtained using greenIntegral. GRIN’s results are more accurate because of the proper handling of the log singularity.
GRIN PART II 22 PART II Improvements made to GRIN • Accuracy and performance: Green’s functions and 2 F 1 • Portability: building the code on different platforms • Robustness: rely on extensively tested software when possible
GRIN PART II 23 Accuracy and Performance The Green’s functions for the toroidal Laplace equation are expressed in terms of the Gauss hypergeometric series ( 2 F 1 ). As the series converges slowly for some arguments, it is prone to inaccuracy and bad performance.
GRIN PART II 24 The 2 F 1 Function The 2 F 1 function is represented by the Gauss hypergeometric series: ∞ z k ( a ) k ( b ) k � F ( a, b, c | x ) = ( c ) k k ! k =0 where: ( a ) k = a ( a + 1) ... ( a + k − 1) , ( a ) 0 = 1 � �� � k terms The series has the convergence radius r = 1 in the complex plane. For x ∈ (0 , 1) the series converges. For x > = 1 the series diverges. For x → 1 − the series converges slowly.
GRIN PART II 25 The 2 F 1 Function The most important Green’s function (the toroidally averaged Green’s function for the Laplace equation) uses 2 F 1 in the following form: F ( n + 1 2 , 1 2 , n + 1 | x ) for x ∈ (0 , 1) This can be computed using the series: ∞ F ( a, b, a + b | x ) = Γ( a + b ) ( a ) k ( b ) k � [2 ψ ( k + 1) − ( k !) 2 Γ( a )Γ( b ) k =0 ψ ( a + k ) − ψ ( b + k ) − ψ (1 − x )](1 − x ) k for x ∈ (0 , 1) which converges very well for x → 1 − , but slowly for x → 0 + .
GRIN PART II 26 Use 2 F 1 from GSL? GNU Scientific Library (GSL) provides only one general function ( gsl_sf_hyperg_2F1 ) to compute: F ( a, b, c | x ) for x ∈ (0 , 1)
Recommend
More recommend
