GSoC 2017 Application Szymon Mieszczak: Implementation of multiple types of coordinate systems for vectors

Szymon Mieszczak edited this page Apr 3, 2017 · 5 revisions
Clone this wiki locally

Personal information

Name: Szymon Mieszczak

University: Adam Mickiewicz University in Poznań, Poland

Email: szymon.mieszczak@gmail.com

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:

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 CoordSystem3D class.
  • 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 CoordSystem3D class 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) and Point classes. This classes have dependency to Cartesian system. For example class Point has method express_coordinate. We need to take into account different coordinate system.
  • [PR] Rebuild Orienter

Week 7, 8 [July 12th - July 25th]

  • [PR] Look into BassisDependent classes. 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, Vector and BaseScalar will be needed.

Evaluation [July 24th- July 28th]

Week 9, 10 [July 26th - August 8th]

  • [PR] Rebuild every function in function.py file. 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))