tutorial.math

##########################
# C P I  T U T O R I A L #
##########################
#
# TL;DR. try out the following quick example:
# $ python Cpi.py quick.example
#
# 
###############################################################################
# Some Rules and Suggestions:
# 01. Comment is shown by '#'.
# 02. All commands expect Declare and `` must be ended by ';'.
# 03. Broken line is by '\' followed by enter.
# 04. Commands:
#
#     *  Ccode["COMMAND"]: prints exactly whatever is written in between 
#     quotation marks ,namely COMMAND, in output C file, e.g. 
#     Ccode["printf("Hello World\n");"], prints exactly: 
#     printf("Hello World\n"); in C file.
#
#     *  `COMMAND`: prints exactly whatever is written in between 
#     backquotation (backtick) marks,namely COMMAND, in output C file, 
#     e.g. `printf("Hello World\n");`, prints exactly:
#     printf("Hello World\n"); in C file.
#
#     *  Pcode["COMMAND"]: prints exactly whatever is written in between 
#     quotation marks for python calculation, e.g.
#     Pcode["x = symbols('x')"].
#
#     *  Symm[symmetry relation]: imposes (anti)symmetry on a tensor, e.g
#     Symm[g(i,j) = g(j,i)].
#
#     *	 tensor1(i,-j,...) == tensor2(i,-j,...): 
#     this command assigns tensor2 to tensor1 for all of the components. 
#     It is used when tensor2 is caclulated and one wants to save the result 
#     in tensor1 in C language; thus, one needs to declare tensor1 in 
#     C language first and then issue this command.
#
#     *  Command["COMMAND"]: it executes a Unix command on the file, e.g:
#     Command["sed -i '$d'"] it issues 'sed' command on the output C file.
#     This can become handy for some string manipulations and more tasks.
#
# 05. The order of executing command is the order that they have been 
#     written. 
# 06. Equations are written in python syntax.
# 07. Naming convention of indexed object is as follows:
#     for each negative sign index letter 'D' and for 
#     each positive sign index letter 'U'. For example,
#     g(i,j) expands:  g_UiUj for i and j in 0,...,N; and 
#     x(i,-j) expands: x_UiDj for i and j in 0,...,N.
# 08. Einstein summation convention is denoted by a pair of indices, one
#     with negative sign the other with no sign; e.g Tu(i)*Td(-i) will be 
#     expanded like: 
#     Tu(0)*Td(0) + Tu(1)*Td(1) + Tu(2)*Td(2) + ... . NOTE, '-' sign DOES
#     NOT MEAN covariant index, NEITHER no sign means contravariant index. 
#     One needs to calculate covariant and contravariant tensor separately 
#     with the desired metric. As a result, one can have different metrics 
#     in the calculations, each for specific tensor(s). However, for naming
#     convention of the left hand side (LHS), and only LHS, quantities 
#     in the EQUATIONS, one can abuse the notation and use indices with negative 
#     (no) sign to give LHS object letter 'D' ('U') and improve readability.
#     For example, in the equation below, if one defines x as follows:
#     x(i,-j) = T(i,k)*g(-k,j) then, x expands like: x_UiDj and it is 
#     kept in mind that x(i,j) = x(i,-j) because '-' sign only has 
#     meaning in the context of Einstein summation convention. 
#     In other words, in the example since there is no other tensor for
#     summation, thus '-' sign is meaningless.
# 09. 'KD' and 'EIJK' are reserved to denote Kronecker Delta and 
#     Levi Civita, respectively.
# 10. Try to write tensorial expressions as simple as possible, i.e. instead
#     of writing a gigantic line full of complicated tensorial contractions 
#     and summations, calculate each term separately by defining new variable 
#     and then add them together.
# 11. Do not use '-' sign indices for tensor of rank > 2. I
#     observed (and reported) a small bug in sympy for them.
# 12. It is recommended that you do not use I, E, S, N, C, O, or Q 
#     for variable or symbol names, as those are used in sympy built-ins.
#     Variables prefixed with "CPI_" are preserved, so do not use this
#     prefix as well. Some important built-in variables are:
#
#     CPI_name:  expands to the name of field
#     CPI_index: expands to the index of the field, starting from 0
#     CPI_head:  expands to the root of the field name, e.g., "g_D0D0" => "g"
#
#     for example see the C_macro5 below.
# 13. If you encounter some undeclared indices at the C output file 
#     or some errors during code generating, it means
#     that the sympy could not resolve the expression so you need to 
#     write it more simpler or break it down to more terms each evaluted
#     separately. After all sympy is quite young! 
# 14. Please feel free to contact me at "rashti.alireza" "at" "gmail.com" for
#     bug report or suggestion.
###############################################################################

# A contrive input file to try the taste of Cpi:

# Manifold or grid Dimension (must be grater than 2):
Dimension = 3;

# point on manifold shown by:
Point = ijk;

# If you need different macros add number to their ends
# CPI_name:  expands to the name of field
# CPI_index: expands to the index of the field, starting from 0
# CPI_head:  expands to the root of the field name, e.g., "g_D0D0" => "g"
C_macro1 = GET_FIELD(CPI_name); 
C_macro2 = GET_PAR(CPI_name);
C_macro3 = GET_FUNCTION1(CPI_name);
C_macro4 = GET_FUNCTION2(CPI_name);
# since ';' is a part of C_macro5 itself thus I need trailing ';;'
C_macro5 = double *CPI_name = double *db->v[i_CPI_head+CPI_index];;

# Special argument: it gets used when you want to give a field special
# argument. E.g, there are cases that a tensor components are functions
# and they need few arguments, C_arg does this for you.
C_arg1 = (ijk,lmn);
C_arg2 = (matrix_name,ijk); # note: name will be replaced by the name of
                            # component of the tensor using this C_arg2.
C_arg3 = (U?,D?);#note: it replaces the indices of the tensor in '?' place

# C file headers
`#include <stdio.h>`
`#include <stdlib.h>`
`#include <math.h>`
`#include "mylib.h"`

# Some macros in C to retrieve pre defined quantities in C
`#define GET_FIELD(name) double *const name = Fields[Ind(#name)];`
`#define GET_PAR(name) const int name = GetParameter(#name);`
`#define GET_FUNCTION1(name) Function_T1 *name = GetFunc1(#name);`
`#define GET_FUNCTION2(name) Function_T2 *name = GetFunc2(#name);`
`` # empty line
# main function
`double *example(double **Fields)\n{`

# Declare all of the fields, variables and functions that are known
# in advance and you just want to use their values in your calculation.
# To use them one can either use macro or direct C code or None as below.
# Note how rank,symmetry and type of the quantities are determined.
# Fields are those objects that are called by memory, like pointers 
# to double in C.

`  const double R = GetParameterD("rotation");`

Declare = 
{
 # a parameter for integration
 (obj = Variable,name = R, none);

 # metric
 (obj = Field,name = g(-i,-j), C_macro1);
 
 # derivative of the metric
 (obj = Field,name = dg(-i,-j,-k), C_macro1);

 # Christopher symbols
 (obj = Field,name = GAMMA(i,-j,-k), C_macro1);

 # a field psi
 (obj = Field,name = psi,C_macro1); 

 # a field B
 (obj = Field,name = B(i),C_macro1); 

 # covariant derivative of field psi
 (obj = Field,name = D_psi(-i), C_macro1);

 # partial derivative of field D_psi
 (obj = Field,name = dD_psi(-i,-j), C_macro1);

 # covariant derivative of field B
 (obj = Field,name = DB(i,-j), C_macro5);

 # a parameter 
 (obj = variable,name = A,C_macro2); 

 # a function in C
 (obj = function,name = Cfunc,None);

 # functional fields ex1:
 (obj = field,name = Ftensor1(i,j),C_macro3,C_arg1);

 # functional fields ex2:
 (obj = field,name = Ftensor2(-i,-j),C_macro4,C_arg2);

 # just for reference:
 ###############################################################
 # an old method to declare tensor.
 # metric
 #(obj = Field,name = g, rank = DD, C_macro1); # rank = DD means 
				  	       # second rank tensor with 
					       # covariant indices.
}
# symmetries:
Symm[g(i,j) = g(j,i)];
Symm[dg(i,j,k) = dg(j,i,k)];
Symm[GAMMA(i,j,k) = GAMMA(i,k,j)];

`  double *F;`
`  int nn = totol_nodes();`
`  int ijk;`
``
`  for (ijk = 0; ijk < nn; ++ijk)`
`  {`

# Note: each command ends with ';'
# Note: I haven't used '-' sign for dg  
# due to the sympy bug for tensors of rank 3 and more.

# christoffel connection:
Con(i,-k,-l)    = 1/2*g(i,-m)*(dg(m,k,l)+dg(m,l,k)-dg(k,l,m));
# that's how a symmetry in equation determined.
Symm[Con(i,k,l) = Con(i,l,k)]; 

# covariant derivative of D_psi
DD_psi(i,j)     = dD_psi(i,j) + Con(k,i,j)*D_psi(-k); 
DDpsi           = g(i,j)*DD_psi(-i,-j); # Laplace operator

# conformal longitudinal operator
LB(i,j)	        = DB(i,j) + DB(j,i) - 2/3*g(i,j)*DB(k,-k); 
LB2(i,j)        = EIJK(l,i,j)*LB(-l,m)*B(-m);
Symm[LB2(i,j)   = -LB2(j,i)]; # anti symmetry

# a python code
Pcode["x = symbols('x')"]; # we need this for following integration
Pcode["s = integrate(pi*log(x),(x,1,R))+cos(R)"];

# one might need different declaration here
Declare = 
{
 # a parameter for integration
 (obj = Variable,name = s, none);
}

s0              = s;
F0              = g(-i,-j)*DD_psi(i,j) - s0*g(i,j)*LB2(-i,-j);
F	        = F0*sin(Cfunc(A));
Jac             = Ftensor1(i,j)*Ftensor2(-i,-j) + s0;
`    Jacobian[ijk] = Jac;`

# assigning the value of calculated Con to GAMMA for each point and
# components
GAMMA(i,-j,-k) == Con(i,-j,-k);

# some C code
`    F[ijk] = F;`
`  }`; # end of for
`  return F;`
`}` # end of function

# NOTE: in older version one can use Ccode["..."] instead of `...`