```python
from sympy import init_printing; init_printing();
from sympy import symbols, sin, cos, tan, acos, sqrt, Eq, Function, Array
from einsteinpy.symbolic import MetricTensor, ChristoffelSymbols, RiemannCurvatureTensor

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all" # display all expression in one cell instead of the last one
```

## Cylindrical Coordinates

Transform Relation from $(r, \theta, \phi)$ to $(x, y, z)$:

\begin{equation}
\left\{\begin{matrix}
x =& R &\sin\phi \\ 
y =& R &\sin\phi \\ 
z =& Z &
\end{matrix}\right.
\end{equation}

Transform Relation from $(x, y, z)$ to $(r, \theta, \phi)$:
\begin{equation}
\left\{\begin{matrix}
 & R^2 =& x^2+y^2 & \\ 
 & Z =& z & \\ 
\tan & \phi  =& y/x & 
\end{matrix}\right.
\end{equation}

Many programming languages have `atan2(y, x)` to avoid the function `atan(y/x)` which is not smart enough.


```python
R = symbols('R', positive=True); Z, phi = symbols('Z phi', real=True)
x, y, z = [f(R, Z, phi) for f in symbols('x, y, z', real=True, cls=Function)]
```


```python
eq_list = [
    Eq(x, R * cos(phi)),
    Eq(y, R * sin(phi)),
    Eq(z, Z)]
for eq in eq_list: eq
diffed_eq_list = []
for eq in eq_list:
    for var in [R, Z, phi]:
        diffed_eq_list.append(
            Eq(eq.lhs.diff(var),
               eq.rhs.diff(var)))# .subs(x**2+y**2+z**2, r**2).simplify() )
```


```python
[eq for eq in diffed_eq_list[0:3]]
[eq for eq in diffed_eq_list[3:6]]
[eq for eq in diffed_eq_list[6:9]]
```

\begin{equation}
\vec{r}_{1} = \frac{\partial \vec{r}}{\partial R} = \frac{\partial (x\hat{x}+y\hat{y}+z\hat{z})}{\partial R},\quad
\vec{r}_{2} = \frac{\partial \vec{r}}{\partial Z} = \frac{\partial (x\hat{x}+y\hat{y}+z\hat{z})}{\partial Z},\quad
\vec{r}_{3} = \frac{\partial \vec{r}}{\partial \phi} = \frac{\partial (x\hat{x}+y\hat{y}+z\hat{z})}{\partial \phi}
\end{equation}


```python
r_i = [None, None, None] # The three components would be 
for i, var in enumerate([R, Z, phi]):
    r_i_temp = Array([x.diff(var), y.diff(var), z.diff(var)])
    for eq in diffed_eq_list:
        r_i_temp = r_i_temp.subs(eq.lhs, eq.rhs)
    r_i[i] = r_i_temp; r_i[i]
    
```


```python
from sympy import trigsimp
# define the metric for 3d cylindrical coordinates
metric = [[0 for i in range(3)] for i in range(3)]
for i in range(3):
    for j in range(3):
        metric[i][j] = (lambda vec1, vec2: 
                        vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2])(r_i[i], r_i[j])
        metric[i][j] = trigsimp(metric[i][j])
        # metric[i][j]
# creating metric object
m_obj = MetricTensor(metric, (R, Z, phi), config='ll'); m_obj.tensor()
g = m_obj.tensor().tomatrix().det(); g
```


```python
ch = ChristoffelSymbols.from_metric(m_obj)
ch.config
```


```python
# Calculating Riemann Tensor from Christoffel Symbols
rm1 = RiemannCurvatureTensor.from_christoffels(ch)
rm1.tensor()
```


```python
# Calculating Riemann Tensor from Metric Tensor
rm2 = RiemannCurvatureTensor.from_metric(m_obj)
rm2.tensor()
```

## Construct $\rho$ and $\vec{J}$ Field


```python
from sympy import exp, Matrix, sqrt, trigsimp
from einsteinpy.symbolic.tensor import BaseRelativityTensor
from einsteinpy.symbolic.vector import GenericVector

```


```python
# define the Levi-Civita tensor for 3d cylindrical coordinates
levi_civita_arr_lll = [[[0 for k in range(3)] for j in range(3)] for i in range(3)]
for i in range(3):
    for j in range(3):
        for k in range(3):
            levi_civita_arr_lll[i][j][k] = Matrix([r_i[i], r_i[j], r_i[k]]).det() / sqrt(g)
            levi_civita_arr_lll[i][j][k] = trigsimp(levi_civita_arr_lll[i][j][k])
```


```python
levi_civita = BaseRelativityTensor(
    levi_civita_arr_lll,
    (R,Z,phi), config='lll', parent_metric=m_obj)
levi_civita.tensor()
```


```python
density_df = BaseRelativityTensor(
    arr=exp(-R**2),
    syms=(R,Z,phi), config='', parent_metric=m_obj); density_df.arr

sqrt_g_df = BaseRelativityTensor(
    arr=sqrt(g),
    syms=(R,Z,phi), config='', parent_metric=m_obj); sqrt_g_df.arr

nabla_density_df = GenericVector(
    arr=[exp(-R**2).diff(var) for var in [R,Z,phi]],
    syms=(R,Z,phi), config='l', parent_metric=m_obj)
nabla_density_df

B_comp_list = [f(R,Z,phi) for f in symbols("B^1:4", real=True, cls=Function)]
B_df = GenericVector(
    arr=B_comp_list,
    syms=(R,Z,phi), config='u', parent_metric=m_obj)
B_df.tensor()
```


```python
J_df = GenericVector([0, exp(-R**2), exp(-R**2)], (R,Z,phi), config='u', parent_metric=m_obj)
J_df
# j_df = j_df.change_config('l')
J2_df = GenericVector([0,exp(-R**2),0], (R,Z,phi), config='u', parent_metric=m_obj)
```


```python
J_df == J2_df
J_df.tensor() == j2_df.tensor()
```


```python
from einsteinpy.symbolic.tensor import tensor_product
tensor_product(density_df, J_df)
```


```python
def tensor_cross(tensor1, tensor2, i=None, j=None):
    from einsteinpy.symbolic.tensor import tensor_product as prod
    middle_tensor = \
        prod( 
            prod(
                prod(levi_civita, sqrt_g_df), 
                tensor1, i=0,j=0
                ), 
            tensor2, i=0,j=0 
            )
#     middle_tensor.change_config('llluu')
    return middle_tensor
```


```python

```


```python
# nabla_density == j x B
nabla_density_df.tensor(), nabla_density_df.config
tensor_cross(J_df, B_df).tensor() # j x B
tensor_cross(J_df, B_df).config
```


```python
ideal_MHD_eq = Eq(
    nabla_density_df.tensor(), 
    tensor_cross(J_df, B_df).tensor())
```


```python
def divide_Array_Eq(eq):
    from sympy import Eq
    assert(eq.lhs.shape == eq.rhs.shape)
    eq_shape = eq.lhs.shape
    arr_order = len(eq_shape) 
    if arr_order > 1: # Not yet tested for high order tensor
        eq_list = []
        for i in range(eq_shape[0]):
            eq_list.append(
                divide_Array_Eq(Eq(eq.lhs[i], eq.rhs[i]))
            )
        return eq_list
    elif arr_order == 1:
        return [Eq(eq.lhs[i], eq.rhs[i]) for i in range(eq_shape[0])]
    else:# arr_order == 0
        return eq
```


```python
from sympy import solve, nonlinsolve
B_sol = solve(
    divide_Array_Eq(ideal_MHD_eq), 
    B_df.tensor().tolist())
B_sol
```


```python
B_comp_list
B_comp_N_list = B_comp_list
for i in range(len(B_comp_N_list)):
    B_comp_N_list[i] = B_comp_N_list[i].subs(B_sol)
    B_comp_N_list[i] = B_comp_N_list[i].subs({B_comp_list[1]: 4, B_comp_list[2]: 6})
B_comp_N_list



B_N_df = GenericVector(
    arr=B_comp_N_list,
    syms=(R,Z,phi), config='u', parent_metric=m_obj)
B_N_df.tensor()
```


```python
from scipy.integrate import ode as scipy_ode
import ray; ray.init()
```


```python




```


```python
from scipy.integrate import ode as _scipy_ode

# @ray.remote
class FieldTracer(_scipy_ode): # FieldTracer is a kind of ode integrator. 
    def __init__(self, vec_field):
        from sympy.utilities.lambdify import lambdify
        assert(vec_field.config == 'u')
        A__1, A__2, A__3 = vec_field.tensor().tolist()
        A__1 = lambdify(vec_field.syms, A__1)
        A__2 = lambdify(vec_field.syms, A__2)
        A__3 = lambdify(vec_field.syms, A__3)
        def dfdt(t, x):
            du1 = A__1(*x)
            du2 = A__2(*x)
            du3 = A__3(*x)
            return [du1, du2, du3]
        _scipy_ode.__init__(self, dfdt)
        
    def integrate_until(self, t1, dt):
        import numpy as np
        num_of_t = int(t1 / dt); num_of_t
        u_arr = np.empty((num_of_t, 3)); u_arr[0, :] = self._y
        t_arr = np.arange(num_of_t) * dt
        i = 0
        while self.successful() and i < num_of_t-1:
            i += 1
            u_arr[i, :] = self.integrate(self.t+dt)
        return u_arr, t_arr
    
    def set_f_params(self, *args, **kwargs):
        raise DeprecationWarning("The set_f_params method has been removed from the FieldTracer class.")
    
    def set_jac_params(self, *args, **kwargs):
        raise DeprecationWarning("The set_jac_params method has been removed from the FieldTracer class.")
    
# tracers = [FieldTracer.remote(B_N_df)for i in range(16)]
t = FieldTracer(B_N_df).set_integrator('vode', method='bdf')
# [t.set_integrator.remote('vode', method='bdf') for t in tracers]

u0, t0 = (10.0, 0.0, 0.0), 0.0
t.set_initial_value(u0, t0)
# [t.set_initial_value.remote(u0, t0) for t in tracers]
RZphi_arr, t_arr = t.integrate_until(10, 0.05)
# ray.get([t.integrate_until.remote(10, 0.05) for t in tracers])




```


```python
from sympy.utilities.lambdify import lambdify
import numpy as np
B__1, B__2, B__3 = B_N_df.tensor().tolist()
# B__1_arr = lambdify(B_N_df.syms, B__1, modules='numpy')(RZphi_arr[:, 0], RZphi_arr[:, 1], RZphi_arr[:, 2])
# B__2_arr = lambdify(B_N_df.syms, B__2, modules='numpy')(RZphi_arr[:, 0], RZphi_arr[:, 1], RZphi_arr[:, 2])
# B__3_arr = lambdify(B_N_df.syms, B__3, modules='numpy')(RZphi_arr[:, 0], RZphi_arr[:, 1], RZphi_arr[:, 2])
B__1_arr = np.zeros(len(RZphi_arr[:,0]))
B__2_arr = np.ones(len(RZphi_arr[:,0])) * 4
B__3_arr = np.ones(len(RZphi_arr[:,0])) * 6

```


```python
r_i_arr = [0 for i in range(3)]
for i in range(3):
    r_i_arr[i] = lambdify(B_N_df.syms, r_i[i], modules='numpy')(
        RZphi_arr[:, 0], RZphi_arr[:, 1], RZphi_arr[:, 2])
    r_i_arr[i][2] = np.zeros_like(r_i_arr[i][1])


B_arr = B__1_arr[:, None] * np.stack(r_i_arr[0], axis=-1) + \
        B__2_arr[:, None] * np.stack(r_i_arr[1], axis=-1) + \
        B__3_arr[:, None] * np.stack(r_i_arr[2], axis=-1)

```


```python
import plotly.graph_objects as go
import numpy as np

fig = go.Figure(data=go.Cone(
    x=RZphi_arr[:, 0] * np.cos(RZphi_arr[:, 2]), 
    y=RZphi_arr[:, 0] * np.sin(RZphi_arr[:, 2]), 
    z=RZphi_arr[:, 1],
    u=B_arr[:, 0],
    v=B_arr[:, 1],
    w=B_arr[:, 2],
    colorscale='Blues',
    sizemode="absolute"
#     marker=dict(
#         size=4,
#         color=p_arr,
#         colorscale='Viridis',
#     ),
#     line=dict(
#         color=p_arr,
#         width=2
#     )
))

fig.update_layout(
    scene=dict(
        aspectratio = dict( x=1, y=1, z=1.0 ),
        aspectmode = 'data'))
# fig.show()

```


```python

```


```python

```

b78d21fd91957d9c1b8c2034253a1e97470e8f00
ipynb
Jupyter Notebook
nb/from_cylindrical_to_metric_Christoffel_Riemann_tensor.ipynb
jiaxin1996/silkpy
7720d47b33b731d9e11e67d99c8574514b8f177b
MIT
nb/from_cylindrical_to_metric_Christoffel_Riemann_tensor.ipynb
jiaxin1996/silkpy
7720d47b33b731d9e11e67d99c8574514b8f177b
MIT
nb/from_cylindrical_to_metric_Christoffel_Riemann_tensor.ipynb
jiaxin1996/silkpy
7720d47b33b731d9e11e67d99c8574514b8f177b
MIT
Qwen/Qwen-72B
1. YES
2. YES