SymFields: An Open Source Symbolic Fields Analysis Tool for General Curvilinear Coordinates in Python
SSymFields: An Open Source Symbolic Fields Analysis Tool forGeneral Curvilinear Coordinates in Python
CHU NanInstitute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, ChinaDecember 22, [email protected]
Abstract
An open source symbolic tool for vector fields analysis ’SymFields’ is developed in Python. The Sym-Fields module is constructed upon Python symbolic module sympy, which could only conduct scalerfield analysis. With SymFields module, you can conduct vector analysis for general curvilinear coordi-nates regardless whether it is orthogonal or not. In SymFields, the differential operators based on metrictensor are normalized to real physical values, which means your can use real physical value of the vectorfields as inputs. This could greatly free the physicists from the tedious calculation under complicatedcoordinates.
The plasma physics MHD theory is a combination of fluid equations and Maxwell’s equations. TheMHD equations involve the solution to multiple vector fields: electric field (cid:126)E ( (cid:126)R ) , magnetic field (cid:126)B ( (cid:126)R ) ,displacement vector (cid:126)ξ ( (cid:126)R ) , etc. Thus, the derivation in MHD theory usually becomes very complicatednot because of physics issue but due to tremendous long terms in equations. To avoid this kind of tediouswork, symbolic calculation is developed in many programming languages. For example, the GeneralVector Analysis (GVA) toolbox for Mathematica software developed by Prof. H. Qin, which couldconduct symbolic fields analysis for general coordinates with elegant expression after simplification [1].Matlab and Python also have some elementary symbolic calculation functions. However, Mathematicaand Matlab are commercial software with high price, which is hardly affordable to me. In open sourcePython language, there is already a fundamental symbolic calculation module sympy. However, it doesnot provide the general analysis functions for vector fields. The SageMath project is another usefulopen source symbolic calculation software ( ). It provides manyuseful functions for symbolic calculation in vector fields. Although it depends on many modules inpython, its source codes is not compatible with Python. Since we do not wish to get involved with a newlanguage, we would like to develop an open source symbolic fields analysis tool solely within Python. a r X i v : . [ c s . S C ] D ec herefore, we developed the SymFields module in Python for fields analysis in general coordinates.It could not only deal with the vector analysis in commonly seeing orthogonal coordinates, but alsocapable to analyse fields in non-orthogonal coordinates. In this paper, the second section talks aboutfields analysis related mathematics for general coordinates (orthogonal and non-orthogonal). The thirdsection discusses the realization of field analysis symbolic calculation in SymFields module in python.The fourth section give benchmark examples to use SymFields in both orthogonal and non-orthogonalcoordinates. The last section is summary. The SymFields module is available on github at ( https://github.com/DocNan/SymFields ) under GNU General Public License v3.0. The multiple choose of coordinates in Mathematics can sometime make the formulas in fields analysiscomplicated. In orthogonal coordinates, the vector analysis can be simplified due to the orthogonality.The differential field operators can be easily expressed with Lame coefficients in orthogonal coordinates[2]. However, for the more general non-orthogonal coordinates, one must use metric tensor to expressthese operators. Therefore, we shall start our discussion from the general coordinates.
In general R curvilinear coordinates, we can define two sets of basis vectors. First we pick a point P inCartesian coordinate with location vector: (cid:126)R = ( x, y, z ) = x (cid:126)e x + y (cid:126)e y + z (cid:126)e z . Suppose for the same pointP in general curvilinear coordinate, its coordinate is: ( ξ , ξ , ξ ) , where ξ i ( (cid:126)R ) = ξ i ( x, y, z ) is a functionof the Cartesian coordinates. The contra-variant vector basis is defined as: (cid:126)g = ∇ ξ (cid:126)g = ∇ ξ (cid:126)g = ∇ ξ (1)The reciprocal covariant basis vector is defined as: (cid:126)g = ∂ (cid:126)R∂ξ (cid:126)g = ∂ (cid:126)R∂ξ (cid:126)g = ∂ (cid:126)R∂ξ (2)They can be converted through the simple definition that: (cid:126)g i = (cid:126)g j × (cid:126)g k (cid:126)g i · ( (cid:126)g j × (cid:126)g k ) , (cid:126)g i = (cid:126)g j × (cid:126)g k (cid:126)g i · ( (cid:126)g j × (cid:126)g k ) (3)With the above definition we can easily find the orthogonality between the vector basis as: (cid:126)g i · (cid:126)g j = δ ij = (cid:40) i (cid:54) = j )1 ( i = j ) (4)2ith the two set of basis vectors, any vector filed (cid:126)A = (cid:126)A ( (cid:126)R ) can also be expressed in two equal ways: (cid:126)A = (cid:88) i A i (cid:126)g j = (cid:88) i A i (cid:126)g i (5)But we must pay attention that unlike the Cartesian coordinates, the basis vector of general curvilinearcoordinates: (cid:126)g i and (cid:126)g j are not unit vector: | (cid:126)g i | (cid:54) = 1 (cid:54) = | (cid:126)g j | . Thus the vector components A i and A i undercontra and covariant vector basis are not real physical value, instead they are called the contra- andcovariant components of this vector field [3]. The contra variant metric tensor is defined as: g ij = (cid:126)g i · (cid:126)g j = (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g (cid:126)g · (cid:126)g = g g g g g g g g g (6)Similarly, the covariant metric tensor is defined as: g ij = (cid:126)g i · (cid:126)g j . Due to the position exchange propertyof dot product between vectors ( (cid:126)g i · (cid:126)g j = (cid:126)g j · (cid:126)g i ), we can easily find that the transposition of the metrictensor equals to itself as: g ij = g ji = g Tij . The contra and covariant basis vector share one import relationas the product of the two metric tensors is unit matrix: g ij g ij = I . The Jacobian of a coordinate represent the element volume under this coordinate, it is defined as: J = (cid:126)g · ( (cid:126)g × (cid:126)g ) = (cid:113) | g ij | = ∂ ( x, y, z ) ∂ ( ξ , ξ , ξ ) = ∂ (cid:126)R∂ξ · ( ∂ (cid:126)R∂ξ × ∂ (cid:126)R∂ξ ) (7)Due to the reciprocal relation between contra and covariant basis vectors, their Jacobian also hasrelationship as: J (cid:48) = (cid:126)g · ( (cid:126)g × (cid:126)g ) = (cid:112) | g ij | = 1 (cid:112) | g ij | = 1 J = ∂ ( ξ , ξ , ξ ) ∂ ( x, y, z ) = ∇ ξ · ( ∇ ξ × ∇ ξ ) (8) Suppose we have a scaler field: U = U ( ξ , ξ , ξ ) and a vector field: (cid:126)A = A (cid:126)g + A (cid:126)g + A (cid:126)g , where thecontra-variant component A j = A j ( ξ , ξ , ξ ) is a scaler field, then the differential operators in generalcurvilinear coordinates are [3, 4]:* Nabla operator: ∇ ∼ ∂∂ξ i ∇ ξ i = (cid:126)g i ∂∂ξ i (9)3 Gradient: ∇ U = (cid:88) i ∂U∂ξ i ∇ ξ i = (cid:88) i ∂U∂ξ i (cid:126)g i = (cid:88) i,j g ij ∂U∂ξ i (cid:126)g j (10)* Divergence: ∇ · (cid:126)A = (cid:88) i J ∂∂ξ i ( J A i ) (11)* Curl: ∇ × (cid:126)A = 1 J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126)g (cid:126)g (cid:126)g ∂∂ξ ∂∂ξ ∂∂ξ A A A (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 J (cid:15) ijk ∂A k ∂ξ j (cid:126)g i (12)* Laplacian: Since Laplacian is the combination of gradient and divergence, and can be calculateddirectly with ∆ (cid:126)A = ∇ (cid:126)A = ∇ · ∇ U , we will not give the specific expression to it here. A type of common error frequently occurs in using the differential operators with metric tensor. For anexample, in cylinder coordinate, the Jacobian is: J = H r H φ H z = r . According to formula (11), thedivergence is calculated via: ∇ · (cid:126)A = (cid:80) i J ∂∂ξ i ( J A i )= r ( ∂∂r ( r A r ) + ∂∂θ ( r A φ ) + ∂∂z ( r A z ))= r ( r A r ) + ∂A φ ∂φ + ∂A z ∂z ( W rong ) (13)However, this expression for divergence operator in cylinder coordinate is wrong, the correct one shallbe: ∇ · (cid:126)A = r ∂∂r ( rA r ) + r ∂A φ ∂φ + ∂A z ∂z . The reason that cause this error is that the differential operatorsgiven in the above section use non-unity basis vectors. If we want to really use them in physics, weneed to make normalization and convert the basis vectors to unity vectors to get the real meaningfulphysical components values. Which will be critical for the realization of the SymFields we will betalking about in the next section. Suppose the unity contra and covariant basis vector as: (cid:126)e i = (cid:126)g i | (cid:126)g i | = (cid:126)g i H i , (cid:126)e i = (cid:126)g i | (cid:126)g i | = H i cosθ i (cid:126)g i . Where we uses the relation that: (cid:126)g i · (cid:126)g i = δ ii = 1 ⇔ | (cid:126)g i || (cid:126)g i | cosθ i = | (cid:126)g i | H i cosθ i = 1 (14)Among them H i = | (cid:126)g i | = | ∂ (cid:126)R∂ξ i | is the Lame coefficient in ξ i coordinate direction [2], θ i is the angle be-tween the vectors (cid:126)g i and (cid:126)g i . Here we make the physical value normalization to the differential operatorswith unity covariant basis vector (cid:126)e i .Suppose we have a vector field: (cid:126)A = (cid:80) i A i (cid:126)g i = (cid:80) i A i (cid:126)g i , its normalization shall be: (cid:126)A = (cid:80) i ¯ A i (cid:126)e i = (cid:80) i ¯ A i (cid:126)e i , where: ¯ A i = A i H i , ¯ A i = A i | (cid:126)g i | . Thus the contra and covariant components of the fieldfunction for differential operator can be easily replaced with physical value as: (cid:40) A i = ¯ A i H i A i = ¯ A i | (cid:126)g i | = ¯ A i |∇ ξ i | = ¯ A i H i cosθ i (15)4 Dot product normalization: (cid:126)A · (cid:126)B = (cid:88) ij ( ¯ A i (cid:126)e i ) · ( ¯ B j (cid:126)e j ) = (cid:88) ij ¯ A i ¯ B j (cid:126)e i · (cid:126)e j = (cid:88) ij ¯ A i ¯ B j g ij H i H j (16)Where (cid:126)e i · (cid:126)e j = (cid:126)g i H i · (cid:126)g j H j = g ij H i H j = e ij = cosθ ij , θ ij is the angle between covariant basis vectors (cid:126)g i and (cid:126)g j .* Cross product normalization: (cid:126)A × (cid:126)B = (cid:80) ij A i (cid:126)g i × B j (cid:126)g j = (cid:80) ijk A i B j (cid:15) ijk J (cid:126)g k = (cid:80) ijk A i B j (cid:15) ijk J H k (cid:126)e k = J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126)g (cid:126)g (cid:126)g A A A B B B (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) H (cid:126)e H (cid:126)e H (cid:126)e g i ¯ A i /H i g i ¯ A i /H i g i ¯ A i /H i g i ¯ B j /H j g j ¯ B j /H j g j ¯ B j /H j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (17)Where (cid:15) ijk is Levi–Civita symbol and A i = g ij A j [5].* Gradient normalization: ∇ U ( ξ , ξ , ξ ) = (cid:88) i,j g ij ∂U∂ξ i (cid:126)g j = (cid:88) i,j g ij ∂U∂ξ i H j (cid:126)e j (18)Notice that for scaler field function, U = U ( (cid:126)ξ ) is already a normalized physical value. Thus we onlyneed to convert the basis vectors of gradient to unity ones.* Divergence normalization: ∇ · (cid:126)A = (cid:88) i J ∂∂ξ i ( J A i ) = (cid:88) i J ∂∂ξ i ( J ¯ A i H i ) (19)* Curl normalization: ∇ × (cid:126)A = 1 J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126)g (cid:126)g (cid:126)g ∂∂ξ ∂∂ξ ∂∂ξ A A A (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) H (cid:126)e H (cid:126)e H (cid:126)e ∂∂ξ ∂∂ξ ∂∂ξ g i ¯ A i /H i g i ¯ A i /H i g i ¯ A i /H i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (20)Where A i = g ij A j = g ij ¯ A j /H j , also notice that Einstein’s summation assumption is used here. In orthogonal coordinates, the Jacobian is the product of all Lame coefficients: J = H H H . Thecontra-variant metric tensor shall be: ( g ij ) = /H /H
00 0 1 /H = ( δ ij H i H j ) (21)Where δ ij is delta function, when i = j , δ ii = 1 , otherwise it takes 0 value. And the angle between (cid:126)g i and (cid:126)g i is: θ i = 0 ⇔ (cid:126)g i (cid:107) (cid:126)g i ⇔ cosθ i = 1 ⇔ H i | (cid:126)g i | = 1 (22)5hus the general normalized differential operators can be rewritten as:* Gradient: ∇ U = (cid:88) i,j g ij ∂U∂ξ i H i (cid:126)e i = (cid:88) i ∂U/∂ξ i H i (cid:126)e i (23)* Divergence: ∇ · (cid:126)A = (cid:88) i J ∂∂ξ i ( J ¯ A i H i ) = 1 H H H (cid:88) i ∂∂ξ i ( ¯ A i H j H k ) (24)* Curl: ∇ × (cid:126)A = 1 J (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) H (cid:126)e H (cid:126)e H (cid:126) ∂∂ξ ∂∂ξ ∂∂ξ g i ¯ A i /H i g i ¯ A i /H i g i ¯ A i /H i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 H H H (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) H (cid:126)e H (cid:126)e H (cid:126)e ∂∂ξ ∂∂ξ ∂∂ξ H ¯ A H ¯ A H ¯ A (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (25)The above expression of differential operators in orthogonal coordinates are the same as the ones intextbooks [2]. This also proves the correctness of our normalization. To make the calculation inputs in a minimum state, we only need to set the coordinates ( ξ , ξ , ξ ) andcontra-variant metric tensor g ij as inputs for Python symbolic calculation. The related Lame coefficientsand covariant metric tensor could all be calculate from the two inputs. Since the covariant metric tensor is the inverse matrix of contra-variant metric tensor. It is very easyto use sympy module to calculate the inverse matrix with Gaussian elimination method: g ij = ( g ij ) − . Whether the coordinates is orthogonal or not, we can always get the Lame coefficients from thecovariant metric tensor as: H i = | (cid:126)g i | = √ g ii and similarly we have: | (cid:126)g i | = |∇ ξ i | = (cid:112) g ii . For the cross product and curl operator, if we set the vector field components as: (cid:126)A = A i (cid:126)g i . Then aftercalculation we can get the results with covariant basis as: (cid:126)g j . However, to make an agreement duringnormalization, we choose the contra-variant components of field ¯ A j by default. Thus we have to convertthe covariant components A i to normalized contra-variant components. This could also be achieved withthe help of co-variant metric tensor as: A i = g ij A j = g ij ¯ A j /H j . Following the formulas in section 2, we developed several functions in SymFields module to realizethe symbolic field analysis. They are: 6 e t r i c ( )
After import * from SymFields module, you can use the above functions to realize fields analysis forgeneral curvilinear coordinates, regardless whether it is orthogonal or not.
In any curvilinear coordinates, the curl of gradient for a scaler field will always be zero: ∇ × ( ∇ U ) = 0 .Thus we can use this rule to test our code. The benchmark test code is presented here: i m p o r t sympyfrom S y m F i e l d s i m p o r t * The output results are:
I n : g r a dOut :I n : c u r l g r a dOut :
Similarly, the divergence of curl operator for a vector field will also be zero: ∇ · ( ∇ × (cid:126)A ) = 0 . Therelated code is presented here: 7 r = sympy . F u n c t i o n ( ’ A r ’ )A p h i = sympy . F u n c t i o n ( ’ A phi ’ )A z = sympy . F u n c t i o n ( ’ A z ’ )A r = A r ( r , p h i , z )A p h i = A p h i ( r , p h i , z )A z = A z ( r , p h i , z )A = [ A r , A phi , A z ]c u r l = C u r l (A, X, c o o r d i n a t e = ’ C y l i n d e r ’ )d i v c u r l = Div ( c u r l , X, c o o r d i n a t e = ’ C y l i n d e r ’ )d i v c u r l . d o i t ( )I n : c u r lOut :I n : d i v c u r l . d o i t ( )Out : We also calculated the Laplacian operator with: ∇ U = ∇ · ∇ U . The detailed codes goes below: L a p l a c i a n = Div ( g r a d , X, c o o r d i n a t e = ’ C y l i n d e r ’ , e v a l u a t i o n =1 )I n : L a p l a c i a nOut :
After manually simplification, you will find it is the correct expression of Laplacian operator in cylin-der coordinate as the formula in textbooks. 8 .1.4 Build-in coordinates in SymFields module
Figure 1: Build-in coordinates in SymFields module. From left to right, Cartesian, Cylinder coordinate,Sphere and Toroidal coordinates.In SymFields module, we have already constructed the metric tensors for several frequently used or-thogonal coordinates. They are: Cartesian, Cylinder, Sphere and Toroidal coordinates as shown in Fig.. To use these build in coordinates, you just need to set string value in the optional input like this:functionA(coordinate=’Cylinder’). For other coordinates, you can first get the mapping function rela-tion between this curvilinear coordinates ( ξ , ξ , ξ ) and Cartesian coordinates (x, y, z). Then you cancalculate the related metric with function Metric() in SymFields module. Finally, you can conduct allthe rest differential operations with the calculated metric tensor as inputs for the operator functions inSymFields module. Figure 2: Non-orthogonal z axis shifted cylinder coordinate.The above subsection, we have tested the SymFields module in an orthogonal (cylinder) coordinates,now we will benchmark it under non-orthogonal coordinates by given an α angle shift to the originalZ coordinate of (orthogonal) cylinder coordinate. As shown in Fig. 2, the z (cid:48) axis is shifted with angle α to z axis in the y-o-z plane. The mapping from original cylinder coordinates to the shifted cylinder9oordinates shall be: x = r (cid:48) cosφ (cid:48) y = r (cid:48) sinφ (cid:48) + z (cid:48) sinαz = z (cid:48) cosα (26)With SymFields, we can easily get its covariant metric tensor by: Let’s further check the curl of gradient and divergence of curl in this non-orthogonal coordinate. Tomake the expressions more simple, we value the shift angle as: α = π/ .* curl of gradient a l p h a = p i / 3Xi = [ r 2 , p h i 2 , z 2 ]x = r 2 * sympy . c o s ( p h i 2 )y = r 2 * sympy . s i n ( p h i 2 ) + z 2 * sympy . s i n ( a l p h a )z = z 2 * sympy . c o s ( a l p h a )R = [ x , y , z ]M co = M e t r i c ( Xi=Xi , R=R , c o o r d i n a t e = ’ s h i f t e d c y l i n d e r ’ , c o n t r a =0 , e v a l u a t i o n = 1 )M co = sympy . s i m p l i f y ( M co )M c o n t r a = M co . i n v ( method = ’GE ’ )U2 = sympy . F u n c t i o n ( ’ U2 ’ )U2 = U2 ( r 2 , p h i 2 , z 2 )g r a d 2 = Grad ( U2 , Xi , c o o r d i n a t e = ’ s h i f t e d c y l i n d e r ’ , m e t r i c = M c o n t r a , e v a l u a t i o n = 1 )c u r l g r a d 2 = C u r l ( g r a d 2 , Xi , c o o r d i n a t e = ’ s h i f t e d c y l i n d e r ’ , m e t r i c = M c o n t r a , e v a l u a t i o n = 1 )I n : g r a d 2Out : I n : sympy . s i m p l i f y ( c u r l g r a d 2 [ 0 ] )I n : sympy . s i m p l i f y ( c u r l g r a d 2 [ 1 ] )I n : sympy . s i m p l i f y ( c u r l g r a d 2 [ 2 ] )Out :
We can find that all 3 components of the curl of gradient is 0.Then let’s check the divergence of curl.
A r2 = sympy . F u n c t i o n ( ’ A r2 ’ )A p h i 2 = sympy . F u n c t i o n ( ’ A phi2 ’ )A z2 = sympy . F u n c t i o n ( ’ A z2 ’ )A r2 = A r2 ( r 2 , p h i 2 , z 2 )A p h i 2 = A p h i 2 ( r 2 , p h i 2 , z 2 )A z2 = A z2 ( r 2 , p h i 2 , z 2 )A2 = [ A r2 , A phi2 , A z2 ]c u r l 2 = C u r l ( A2 , Xi , c o o r d i n a t e = ’ s h i f t e d c y l i n d e r ’ , m e t r i c = M c o n t r a , e v a l u a t i o n = 1 )d i v c u r l 2 = Div ( c u r l 2 , Xi , c o o r d i n a t e = ’ s h i f t e d c y l i n d e r ’ , m e t r i c = M c o n t r a , e v a l u a t i o n = 1 )I n : sympy . s i m p l i f y ( c u r l 2 )Out :
We can also find the expression of curl in this α = π/ shifted cylinder coordinate is also very longand complicated. However, we can also verify that the divergence of curl remains zero. I n : sympy . s i m p l i f y ( d i v c u r l 2 )Out :
In this paper, we report the development of an open source symbolic calculation tool for vector fieldanalysis in Python. The SymFields module is constructed upon Python symbolic module sympy, whichcould only conduct scaler field analysis. Within SymFields module, the vector fields operation aredefined upon the metric tensor of a general curvilinear coordinates. Which means you can conductvector analysis for any curvilinear coordinates regardless whether it is orthogonal or not. Four orthogonal11oordinates: Cartesian, Cylinder, Sphere and Toroidal are set at build-in coordinates. You can extend itto any other coordinates by providing a new metric tensor. In SymFields, the differential operators basedon metric tensor are normalized to real physical values, which means your can use real physical valueof the vector fields as inputs. Thus could greatly free the physicists from the tedious calculation undercomplicated coordinates.
References [1] H. Qin, W. M. Tang, and G. Rewoldt. Symbolic vector analysis in plasma physics.
ComputerPhysics Communications , 116(1):107 – 120, 1999.[2] Introduction to Advanced Mathematics (printed in Chinese), University of Science and Technologyof China Press, by group author of Mathematics teaching and research team, 2008 Hefei, China.[3] V. I. Piercey. Lame and metric coefficients for curvilinear coordinates in R . Lecture notes, Univer-sity of Arizona.[4] William Denis D’haeseleer, William Nicholas Guy Hitchon, James D. Callen, and J. Leon Shohet. Flux Coordinates and Magnetic Field Structure . Springer Berlin Heidelberg, 1991.[5] Leonid P Lebedev, Michael J Cloud, and Victor A Eremeyev.