GSoC 2017 Application Szymon Mieszczak: Implementation of multiple types of coordinate systems for vectors - sympy/sympy GitHub Wiki
Personal information
Name: Szymon Mieszczak
University: Adam Mickiewicz University in Poznań, Poland
Email: [email protected]
Github: szymag
Short bio
I'm from Poland, Lower Silesia Voivodeship. I have obtained my bachelor’s degree in Physics at Adam Mickiewicz University in Poznań, Poland. Now, I study physics at graduate level at the same University. My bachelor's thesis and my future master's thesis deal with spin waves.
Programming details
I'm working on PC with Linux Mint 18.1. I really like the comfort coming from
using big monitor and separate keyboard. I think work is then much more
efficient. My favorite Python’s IDE is PyCharm. I really like it because it
integrates everything in one place. As an editor it allows me to personalize
for example highlighting of syntax, show me usage of variables or check if my
code is consistent with PEP 8 standard. As a GUI it offers me a terminal, the
Python shell, git log and TODO list into one window. What's more it has
built-in debbuger and profiler.
I have learned C++ in High School and at beginning of my studies. During my studies I have been also using Mathematica, Matlab, LabView, on different courses, for different projects. I have been learning Python since 2014, when my friend has created a course about it. I felt in love with this language. After that, I have been using Python for almost every task, from plotting, to doing research.
Within my studies I am working on (1, 2) spin waves in magnonic (quasi)crystal. Magnonic crystal (4) is a magnetic metamaterial with spatial periodic arrangement of material parameter. I'm writing a program to calculate and visualize physical properties such as dispersion relation, dynamic magnetization profile etc. Short description: I need to transform a differential equation (Landau-Lifshitz) into an eigenproblem using Bloch theory. I achieve that, by expanding material parameter and searching dynamic magnetization excitation (which is one of the sought solution) into Fourier series with period equal to elementary unit cell size. After solving eigen problem I get eigenvalue, which are allowed frequency of spin waves and eigenvector, which are profile of dynamic components of magnetization vector. In my program I mainly use NumPy package. I also combined matrix operation with multiprocessing package from standard library. I was really surprised, how this combination work smoothly and efficiently. I would like to develop further my program, both from programing and physical sides. I hope that my program will be as useful as equivalent program in photonics (3).
I have theoretical background but also some experience in using Git. At the early beginning of developing my program I started work with it. I've read several materials and I have finished one of popular tutorials (5).
SymPy Involvement
Involvement
My involvement to SymPy mainly focus on geometry package. I just look around
in code and found several TODO there. After I went deeper, I've noticed many
things which needed improvements. I've also found and fixed several bugs and
I've added a few minor improvements in several other places. I also improved
coverage in this package, cleaned up tests, using good programming practices.
I'm looking there further, because I'm still finding things that should be
improved.
Pull Requests:
- Repair is_tangent method in Ellipse class [Open]
- Polygon intersection [Open]
- Simplification in Ellipse intersection method [Merged]
- Clean up test_line [Merged]
- [Parabola intersection]https://github.com/sympy/sympy/pull/12347 [Merged]
- Simplification
intersection()inEllipseclass [Merged] - Improve plane class [Merged]
- Use linsolve to finding intersections in Plane class. Realization of TODO task [Merged]
- Improve geometry coverage [Merged]
The Project
Overview
I would like to generalize package Vector to support any type of orthogonal
curvilinear coordinates (link).
Many problems have symmetry different than Cartesian Coordinates and it is more convenient to use those coordinates that convey symmetry of problem. For example, a suitable coordinate system for a system of ionized hydrogen element are prolate spheroidal coordinates, which corresponds to symmetry of electrical potential generated by two nucleus with some distance between them. My improvement aims at filling this gap.
During my studies I had vector analysis several times, both on a mathematics
course and physics, like classical and relativistic mechanics or
electrodynamics. My knowledge suits a material contained in books “Mathematical Method for Physicist” G.B. Arfken and “Mathematical Methods for Scientist and Engineers” D. A McQuarrie. Another qualification I posses is an experience in
building complex projects which I got by building my project, that has
similar size and complexity to the Vector package. I know that some work has
been done by @jbbskinny and @Upabjojr. I agree with them, that it
should be done by creating general class and every operation on vector or
translation between coordinates should be done by Lamé coefficient.
@jbbskinny did a nice job by adding these Lamé coefficients and a method
transforming from one coordinate system to another to the CoordSystem3D
class. I think that I could use that in my project.
Every operation in vector analysis can be extend in terms of Lamé coefficient
(at least in orthogonal curvilinear coordinate system). My idea to solve this
problem is to define a class CoordSystem3D (following to @jbbskinny)
which will be handling everything that now is implemented in the
CoordSysCartesian, but it will have also a possibility to set int the constructor
a type of orthogonal coordinate system. I would like to add (like
@jbbskinny) several built-in system, but allow the User to introduce
additional ones. As a default, I will put Cartesian one.
Initialization should look like that:
from sympy.vector import CoordSystem3D
s = CoordSystem3D(‘s’, coord=spherical) # Initialization spherical coordinate
c = CoordSystem3D(‘c’) # Initialization Cartesian coordinate
t = CoordSystem3D('t', eq1=x*cos(y),
eq2=x*sin(y), eq3=z) # Initialization of cylindrical coordinate
# using transformation equations
I would like to arrange this problem using Factory Pattern,
where the Client will be the CoordSystem3D class , the Factory from the
diagram under this link will be some class creating proper instance of a class
inheriting from LameCoeff, which will act as the Product. It allows us to preserve
Open/Closed principle,
because CoordSystem3D and LameCoeff will work independently. Each of them
shouldn't know about how the other one works and it will be easy to add
a new coordinate system, by implementing new class or, simply, passing custom
Lamé coefficients.
In the Attachment section I'm putting a code, which demonstrates this concept.
Theory
Introduction of a new coordinate systems involve reconstruction of every method
which is subordinate from coordinate system. Let's look at example, how to
calculate gradient, divergence and curl in spherical polar coordinate
system, but applying Lamé coefficient. The code below works properly and can be
copied and run. Enclosed result are calculated, using this function.
We need to import necessary package and create variables:
from __future__ import division
from sympy import *
from sympy import sin, cos, Matrix
from sympy.vector import CoordSysCartesian
from sympy.abc import rho, phi, theta
x, y, z = symbols('x y z', real=True)
r = symbols('r', real=True, positive=True)
Transformation equations from Cartesian to spherical polar coordinates are the following:
x = r * sin(phi) * cos(rho)
y = r * sin(rho) * sin(phi)
z = r * cos(phi)
Now we define scale factor (Lamé coefficient):
h1 = simplify(sqrt(diff(x, r)**2 + diff(y, r)**2 + diff(z, r)**2))
h2 = simplify(sqrt(diff(x, rho)**2 + diff(y, rho)**2 + diff(z, rho)**2))
h3 = simplify(sqrt(diff(x, phi)**2 + diff(y, phi)**2 + diff(z, phi)**2))
We are ready to define function, which calculate gradient, curl and divergence.
Gradient:
def gradient(scalar):
s = CoordSysCartesian('spherical') # forget, that it is Cartesian coordinate,
# we just need unit vector.
grad = s.i / h1 * diff(scalar, r) + s.j / h2 * diff(scalar, phi) + \
s.k / h3 * diff(scalar, rho)
return simplify(grad)
Example taken from book "Mathematical Methods for Scientists and Engineers", Donald McQuarrie:
gradient(r**2*sin(phi))
(2⋅r⋅sin(φ)) spherical_i + ⎛ r⋅cos(φ) ⎞ spherical_j
⎜────────────⎟
⎜ _________⎟
⎜ ╱ 2 ⎟
⎝╲╱ sin (φ) ⎠
Divergence:
def divergence(vector):
div = 1 / (h1*h2*h3) * (diff(h2 * h3 * vector[0], r) \
+ diff(h1 * h3 * vector[1], rho)\
+ diff(h1 * h2 * vector[2], phi))
return div
Example taken from book "Mathematical Methods for Scientists and Engineers", Donald McQuarrie:
divergence([r, 0, 0])
3
Curl:
def curl(vector):
s = CoordSysCartesian('spherical') # forget,that it is Cartesian coordinate
curl = 1 / (h1*h2*h3) * (
s.i * h1 * (diff(vector[2] * h3, phi) - diff(vector[1] * h2, rho)) \
+ s.j * (diff(vector[0] * h1, rho) - diff(vector[2] * h3, r))\
+ s.k * (diff(vector[1] * h2, r) - diff(vector[0] * h1, phi)))
return curl
Example taken from book "Mathematical Methods for Scientists and Engineers", Donald McQuarrie:
curl([r, 0, 0])
0 # vector zero
Fix up
I would like to improve the existing structure. Now, when we want for example calculate a divergence, we need to do the following (taken from documentation):
from sympy.vector import CoordSysCartesian
R = CoordSysCartesian('R')
v1 = R.x*R.y*R.z * (R.i+R.j+R.k)
divergence(v1, R)
We manually need to put information to the divergence function that this
calculation is made in R coordinate system. From the user's perspective it is
not necessary, because, if we create vector in one coordinate system, we
can’t set different coordinate system. My plan here is to remove this
argument from every method they have and forward information about coordinate
system with vector. From @Upabjojr
(PR #12417)
I know that the Vector has information about its coordinate system. We can get
it using args so it is possible to remove explicit setting coordinate in
a function.
Overloading operators
following the @brombo suggestion,
I would like to introduce overload operator * and ^ for differential operation.
Operation will look this:
gradient(f) = del * f
divergence(V) = del * V
curl(V) = del^V
Where f is a scalar field and V is a vector field.
In Python overloading operator is really easy task, but we need to create good structure
for nabla operator, to handle this future. I treat this future as extra functionality,
which will be added when work will be according to plan.
Possible extension:
Using Lamé coefficient, instead of hard typing every coordinates separately, it would be also beneficial, if we would like to extend our vector package to be able to do vector integration. Infinitesimal volume in a curvilinear orthogonal system can be formulated with the Jacobian determinant, which looks as following:
⎡∂ ∂ ∂ ⎤
⎢──q₁ ──q₁ ──q₁⎥
⎢∂x ∂y ∂z ⎥
⎢ ⎥
⎢∂ ∂ ∂ ⎥
⎢──q₂ ──q₂ ──q₂⎥
⎢∂x ∂y ∂z ⎥
⎢ ⎥
⎢∂ ∂ ∂ ⎥
⎢──q₃ ──q₃ ──q₃⎥
⎣∂x ∂y ∂z ⎦
But, for an orthogonal coordinate system, we known that the Jacobian can be reduced to a product of Lamé coefficient, so our differential element are:
dxdydz = h₁h₂h₃dq₁dq₂dq₃
With that fact, we can simply add vector integration in different orthogonal curvilinear system. I would like to handle this feature after successfully finishing my GSoC project.
Additional information
Before the beginning of the coding period I have still studies but I will be spending as much time as possible to get prepared to work. In that time, I would like prepare majority of tests for new classes and methods. Fortunately, I'm having holiday during the whole June, so nothing will inhibit me from working. During holiday (June-September) I'm going to work exclusively on my GSoC's project, 40+ hours per week. I'm thinking about spending several weeks in Sudetes to work in a nice surroundings, but still working.
In Timeline section I indicated place where I'm thinking of creating PR.
I want to emphasize, that every subtask will have they own PR, and I will
send to the SymPy's repo only commits with a complete realization of some task.
I know that the maintainers of the SymPy are volunteers and I want to save
their time.
Timeline
Week 1,2 [May 30th - June 13th]
- Create from the beginning or adapt some elements from @jbbskinny PR
CoordSystem3Dclass. - This stage is the longest and the most difficult.
- At the early beginning or even before GSoC starts, I will create PR, which will be updated.
Week 3,4 [June 14th - June 27th]
- Finishing and polishing previous part. I think ready
CoordSystem3Dclass with documentation and fully tested is nice thing to evaluation. - [PR] Create structure for nabla operator, to introduce overloading
*and^.
Evaluation [June 26th- June 30th]
Week 5, 6 [June 28th - July 11th]
- [PR] Rebuild all
Vector(Vector,BaseVector,VectorAdd,VectorMul,VectorZero) andPointclasses. This classes have dependency to Cartesian system. For example classPointhas methodexpress_coordinate. We need to take into account different coordinate system. - [PR] Rebuild
Orienter
Week 7, 8 [July 12th - July 25th]
- [PR] Look into
BassisDependentclasses. These classes have overridden method, but I'm still not sure, if we shouldn't change something here. It's better to leave some time and look at that problem when previous task will be ready. - [PR] I would like to change unit vector symbol depending on the coordinate system.
For example, for Cartesian unit vector, common symbol are j, j, k, but for spherical
coordinate symbol unit vectors are r, phi, theta. The changes in
CoordSystem3D,VectorandBaseScalarwill be needed.
Evaluation [July 24th- July 28th]
Week 9, 10 [July 26th - August 8th]
- [PR] Rebuild every function in
function.pyfile. This part, in my opinion, will be the easiest, but also the most responsible. We need to create tests for every function here in different coordinate system. Here I will use test-driven development to be sure, that every method is absolutely correct.
Week 11, 12 [August 9th - August 22th]
- Polishing every previous part. I would like spend here the most possible time on adding examples to documentation. My ideal example here is Wolfram Documentation, where everything what is needed, user can pull out from documentation. I think it's really instructive.
Attachment
I'm presenting here basic concept of introducing the Lamé coefficient.
from __future__ import division
from sympy import *
from sympy import sin, cos
x, y, z = symbols('x y z', real=True)
class LameCoef:
def __init__(self, *args, **kwargs):
pass
def trasformation_equations(self):
pass
def get_lame_coefficient(self):
return self.h1, self.h2, self.h3
def h1(self):
pass
def h2(self):
pass
def h3(self):
pass
class CartesianCoef(LameCoef): # Lame coefficient for Cartesian coordinate system
def trasformation_equations(self):
return 1, 1, 1
def h1(self):
return 1
def h2(self):
return 1
def h3(self):
return 1
class SphericalCoef(LameCoef): # Lame coefficient for Spherical coordinate system
def trasformation_equations(self):
return x * sin(y) * cos(z), x * sin(y) * sin(z), x * cos(y)
def h1(self):
return 1
def h2(self):
return x
def h3(self):
return x * cos(y)
class CurvilinearCoef(LameCoef): # General class for any curvilinear system
def __init__(self, *args, **kwargs):
LameCoef.__init__(self, *args, **kwargs)
try:
self.eq1, self.eq2, self.eq3 = kwargs['eq1'], \
kwargs['eq2'], \
kwargs['eq3']
except KeyError:
raise ValueError('Wrong set of parameters')
def trasformation_equations(self):
return self.eq1, self.eq2, self.eq3
def h1(self):
return simplify(sqrt(diff(self.eq1, x) ** 2 +
diff(self.eq2, x) ** 2 +
diff(self.eq3, x) ** 2))
def h2(self):
return simplify(sqrt(diff(self.eq1, y) ** 2 +
diff(self.eq2, y) ** 2 +
diff(self.eq3, y) ** 2))
def h3(self):
return simplify(sqrt(diff(self.eq1, z) ** 2 +
diff(self.eq2, z) ** 2 +
diff(self.eq3, z) ** 2))
class CoeffProvider:
coordinates_mapping = {
'cartesian': CartesianCoef,
'spherical': SphericalCoef
}
def get_coefficients(self, *args, **kwargs):
if len(args) == 1:
return CoeffProvider.coordinates_mapping[args[0]]()
elif len(kwargs) == 3:
self.eq1, self.eq2, self.eq3 = kwargs['eq1'], kwargs['eq2'], kwargs['eq3']
return CurvilinearCoef(**kwargs)
else:
raise ValueError('Wrong set of parameters')
class CoordSys3D:
def __init__(self, name, *args, **kwargs):
self.name = name
if len(args) == 1:
self.coeff = CoeffProvider().get_coefficients(args[0])
else:
try:
self.coeff = CoeffProvider().get_coefficients(eq1=kwargs['eq1'],
eq2=kwargs['eq2'],
eq3=kwargs['eq3'])
except KeyError:
raise ValueError('Wrong set of parameters')
# Imitation of unit vector
@property
def q1(self):
return symbols('q1')
@property
def q2(self):
return symbols('q2')
@property
def q3(self):
return symbols('q3')
def print_coefficients(self):
print(self.coeff.h1(), self.coeff.h2(), self.coeff.h3())
def gradient(scalar, coordinate):
a = 1 / coordinate.coeff.h1() * diff(scalar, x)
b = 1 / coordinate.coeff.h2() * diff(scalar, y)
c = 1 / coordinate.coeff.h3() * diff(scalar, z)
return simplify([coordinate.q1 * a, coordinate.q2 * b, coordinate.q3 * c])
if __name__ == '__main__':
spherical_1 = CoordSys3D('sph1', 'spherical')
spherical_1.print_coefficients()
cartesian_1 = CoordSys3D('cart', 'cartesian')
cartesian_1.print_coefficients()
spherical_2 = CoordSys3D('sph2', eq1=x * sin(y) * cos(z),
eq2=x * sin(y) * sin(z),
eq3=x * cos(y))
spherical_2.print_coefficients()
cylindrical = CoordSys3D('cyl', eq1=x * cos(y), eq2=x * sin(y), eq3=z)
cylindrical.print_coefficients()
print(gradient(x ** 2 * sin(z), spherical_1))
print(gradient(z * cos(y) / x, cylindrical))