Example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1cac1684",
   "metadata": {},
   "source": [
    "This notebook demonstrates how ebweyl can be used to compute the electric and magnetic parts of the Weyl tensor of a given spacetime and classify it according to the Petrov types.\n",
    "\n",
    "In this process we will compute the spatial Ricci tensor and scalar as well as the Weyl scalar for an arbitrary tetrad base and the invariant scalars: I, J, K, L , N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5765ddeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import ebweyl as ebw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f80ddf2",
   "metadata": {},
   "source": [
    "The documentation on the ebweyl module is accessible with the help() function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "025cb700",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function ricci_tensor_down3 in module ebweyl:\n",
      "\n",
      "ricci_tensor_down3(self, Gudd3)\n",
      "    Compute spatial Ricci tensor with both indices down.\n",
      "    \n",
      "    Parameters : \n",
      "        Gudd3 : (3, 3, 3, Nx, Ny, Nz) array_like \n",
      "                Spatial Christoffel symbol with one indice up and two down\n",
      "                You can get this from Weyl.christoffel_symbol_udd3()\n",
      "            \n",
      "    Returns : \n",
      "        (3, 3, Nx, Ny, Nz) array_like\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(ebw.Weyl.ricci_tensor_down3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95592162",
   "metadata": {},
   "source": [
    "# Define the Spacetime\n",
    "\n",
    "For this example we use the A.Harvey and T.Tsoubelis vacuum Bianchi IV plane wave homogeneous spacetime (see page 191 of 'Dynamical Systems in Cosmology' by J.Wainwright and G.F.R.Ellis). With Cartesian coordinates the metric takes the form:\n",
    "\n",
    "$$\n",
    "g_{\\alpha\\beta} = \n",
    "\\begin{pmatrix}\n",
    "   -1 &  0  & 0 & 0 \\\\\n",
    "    0 & t^2 & 0 & 0 \\\\\n",
    "    0 &  0  & te^{x} & te^{x}(x+\\log{t}) \\\\\n",
    "    0 &  0  & te^{x}(x+\\log{t}) & te^{x}((x+\\log{t})^2+1)\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "This is an analytical spacetime but the code presented here works for any numerical spacetime, the user needs the metric, extrinsic curvature and stress-energy tensor. \n",
    "\n",
    "For this metric the extrinsic curvature is: $K_{\\alpha\\beta} = \\frac{-1}{2}\\partial_{t}(g_{\\alpha\\beta})$.\n",
    "\n",
    "This metric is a solution to vacuum, so the stress-energy tensor is: $T_{\\alpha\\beta} = 0$.\n",
    "\n",
    "Let's create arrays for these tensors. \n",
    "\n",
    "#### Create a data grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6b39298d",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 10  # Data grid size\n",
    "# If you are running this on your own laptop and don't want the \n",
    "# computation to take too much time, you should make N smaller. \n",
    "# (But keep in mind that the accuracy increases as N increases)\n",
    "\n",
    "# Number of data points per side\n",
    "Nx = 64\n",
    "Ny = 72\n",
    "Nz = 80\n",
    "# Elementary grid size\n",
    "dx = L / Nx\n",
    "dy = L / Ny\n",
    "dz = L / Nz\n",
    "\n",
    "# Cartesian coordinates\n",
    "x, y, z = np.meshgrid(np.arange(- L / 2, L / 2, dx), \n",
    "                      np.arange(- L / 2, L / 2, dy), \n",
    "                      np.arange(- L / 2, L / 2, dz), indexing='ij')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32f211fd",
   "metadata": {},
   "source": [
    "#### Create metric tensor array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6f6064de",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 1.5  # An arbitrary time\n",
    "B = (x+np.log(t))  # Function from referenced book\n",
    "zeros = np.zeros(np.shape(x))\n",
    "ones = np.ones(np.shape(x))\n",
    "\n",
    "gdown4 = np.array([[-ones, zeros, zeros, zeros],\n",
    "                   [zeros, (t*t)*ones, zeros, zeros],\n",
    "                   [zeros, zeros, t*np.exp(x), t*np.exp(x)*B],\n",
    "                   [zeros, zeros, t*np.exp(x)*B, \n",
    "                    t*np.exp(x)*(B*B+1)]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1942aca",
   "metadata": {},
   "source": [
    "#### Create extrinsic curvature tensor array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0d0457f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "dtB = 1/t  # Time derivative of B function\n",
    "# Time derivative of metric\n",
    "dtgdown3 = np.array([[2*t*ones, zeros, zeros],\n",
    "                     [zeros, np.exp(x), \n",
    "                      np.exp(x)*B + t*np.exp(x)*dtB],\n",
    "                     [zeros, np.exp(x)*B + t*np.exp(x)*dtB, \n",
    "                      np.exp(x)*(B*B+1) + t*np.exp(x)*(2*dtB*B)]])\n",
    "\n",
    "Kdown3 = (-1/2)*dtgdown3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d1445f3",
   "metadata": {},
   "source": [
    "#### Create stress-energy tensor array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2668049f",
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa = 8*np.pi  # Einstein's constant with G = c = 1\n",
    "Tdown4 = np.zeros(np.shape(gdown4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa2151dd",
   "metadata": {},
   "source": [
    "# Electric $E_{\\alpha\\beta}$ and magnetic $B_{\\alpha\\beta}$ parts of the Weyl tensor\n",
    "\n",
    "#### Define the FiniteDifference class\n",
    "\n",
    "Default values are : periodic_boundary=True, fd_order6=False\n",
    "\n",
    "But periodic boundaries would be inapropriate for this spacetime, the code will then use a combination of forward, centered and backward finite difference schemes instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "836b7b87",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6th order finite difference schemes are defined\n"
     ]
    }
   ],
   "source": [
    "FD = ebw.FiniteDifference([dx, dy, dz], [Nx, Ny, Nz], \n",
    "                          periodic_boundary=False, \n",
    "                          fd_order=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82615902",
   "metadata": {},
   "source": [
    "#### Define the Weyl class\n",
    "\n",
    "The FD class is given to the Weyl class and we won't use it again here, but you can use it on it's own to compute spatial derives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "25c92250",
   "metadata": {},
   "outputs": [],
   "source": [
    "EBW = ebw.Weyl(FD, gdown4, Kdown3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc9e54f6",
   "metadata": {},
   "source": [
    "This will automatically compute the standard terms of the 3+1 formulation, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e0467866",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lapse =  1.0\n"
     ]
    }
   ],
   "source": [
    "print('lapse = ', np.mean(EBW.alpha))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35c41818",
   "metadata": {},
   "source": [
    "Notice that here I'm showing the mean. This is because np.shape(EBW.alpha) = (N, N, N) as the lapse can depend on space.\n",
    "\n",
    "#### Compute Spatial Ricci tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2a34f638",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gudd3 = EBW.christoffel_symbol_udd3()\n",
    "RicciTdown3 = EBW.ricci_tensor_down3(Gudd3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d431878",
   "metadata": {},
   "source": [
    "This then easily provides the spatial Ricci scalar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1e96e68c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R3 =  -0.8888922876460006\n",
      "Analytical answer:  -0.8888888888888888\n"
     ]
    }
   ],
   "source": [
    "RicciS3 = EBW.trace_rank2tensor3(RicciTdown3)\n",
    "print('R3 = ', np.mean(FD.cutoffmask(RicciS3)))\n",
    "print('Analytical answer: ', -2/t**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ee3048b",
   "metadata": {},
   "source": [
    "Notice that here I'm applying the FD.cutoffmask function before taking the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "45d551a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method cutoffmask in module ebweyl:\n",
      "\n",
      "cutoffmask(f) method of ebweyl.FiniteDifference instance\n",
      "    Remove points affected by the boundary condition.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(FD.cutoffmask)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e969f7a7",
   "metadata": {},
   "source": [
    "The backward and forward finite difference schemes slightly underperform compared to the centered one, so that is why I remove those before taking the mean. I wouldn't bother doing this if I were using periodic boundary conditions.\n",
    "\n",
    "#### Compute $E_{\\alpha\\beta}$ and $B_{\\alpha\\beta}$ projected along the normal to the hypersurface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "96deeab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Endown3 = EBW.eweyl_n_tensor_down3(RicciTdown3, kappa, Tdown4)\n",
    "Bndown3 = EBW.bweyl_n_tensor_down3(Gudd3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bd01a54",
   "metadata": {},
   "source": [
    "These are both covariant tensors with both indices being spatial, hence the 3 at the end. \n",
    "\n",
    "If I want the indices to include time as well, I can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8a501db3",
   "metadata": {},
   "outputs": [],
   "source": [
    "Endown4 = EBW.ebweyl_n_3D_to_4D(Endown3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac124a59",
   "metadata": {},
   "source": [
    "Lets look at their norm : $|E| = \\sqrt{E_{\\alpha\\beta}E^{\\alpha\\beta}} = \\sqrt{E_{ij}E^{ij}}$. Here Greek indices are spacetime and Latin indices are only space.\n",
    "\n",
    "As this spacetime has a plane wave we have $|E|=|B|$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "282f4e6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "En_norm = EBW.norm_rank2tensor3(Endown3)\n",
    "# or = EBW.norm_rank2tensor4(Endown4)\n",
    "Bn_norm = EBW.norm_rank2tensor3(Bndown3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f8facaec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E^2 =  0.3142708322346237\n",
      "B^2 =  0.31426956479598134\n",
      "Analytical answer:  0.31426968052735443\n"
     ]
    }
   ],
   "source": [
    "print('E^2 = ', np.mean(FD.cutoffmask(En_norm)))\n",
    "print('B^2 = ', np.mean(FD.cutoffmask(Bn_norm)))\n",
    "print('Analytical answer: ', 1/(np.sqrt(2)*t**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68b1f947",
   "metadata": {},
   "source": [
    "### $E_{\\alpha\\beta}$ and $B_{\\alpha\\beta}$ along a different vector\n",
    "\n",
    "With the terms projected along the normal to the hypersurface we can construct the Weyl tensor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2eb5d8f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Cdown4 = EBW.weyl_tensor_down4(Endown3, Bndown3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "711d25b7",
   "metadata": {},
   "source": [
    "This can then be projected along any time-like vector. We chose the 4-velocity, for this metric it simply corresponds to the normal to the hypersurface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "aa9b5c24",
   "metadata": {},
   "outputs": [],
   "source": [
    "uup4 = EBW.nup4  # spacetime indices up"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc1dcfc",
   "metadata": {},
   "source": [
    "Then the new $E_{\\alpha\\beta}$ and $B_{\\alpha\\beta}$ are obtained as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e4febe82",
   "metadata": {},
   "outputs": [],
   "source": [
    "Eudown4 = EBW.eweyl_u_tensor_down4(Cdown4, uup4)\n",
    "Budown4 = EBW.bweyl_u_tensor_down4(Cdown4, uup4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e42c3a8d",
   "metadata": {},
   "source": [
    "Because the 4-velocity here corresponds to the normal to the hypersurface, I demonstrate here that $|E|$ and $|B|$ are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1b935beb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|E| =  0.3142708322346237\n",
      "|B| =  0.3142695647959811\n",
      "Analytical answer:  0.31426968052735443\n"
     ]
    }
   ],
   "source": [
    "Eu_norm = EBW.norm_rank2tensor4(Eudown4)\n",
    "Bu_norm = EBW.norm_rank2tensor4(Budown4)\n",
    "print('|E| = ', np.mean(FD.cutoffmask(Eu_norm)))\n",
    "print('|B| = ', np.mean(FD.cutoffmask(Bu_norm)))\n",
    "print('Analytical answer: ', 1/(np.sqrt(2)*t**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84ac837a",
   "metadata": {},
   "source": [
    "# Petrov type\n",
    "\n",
    "To find the Petrov type, I construct the Weyl scalars $\\Psi$ and then compute the scalar invariants, (J, I, L, K, N). Those are then used to determine the Petrov type according to Figure 9.1 of \"Exact Solutions to Einstein's Field Equations\" 2003, 2nd edition, by H.Stephani, D.Kramer, M.MacCallum, C.Hoenselaers and E.Herlt.\n",
    "\n",
    "In order to compute the $\\Psi$s I need the full Weyl tensor provided in cell [17]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d3602112",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Psi_0 =  8.962147991068151e-07  + i( 1.1107501262046878e-08 )\n",
      "Psi_1 =  3.345964485086464e-17  + i( -3.4175973710854855e-15 )\n",
      "Psi_2 =  -1.0845286764132787e-07  + i( 2.445126924179975e-07 )\n",
      "Psi_3 =  3.345964485086464e-17  + i( -3.4175973710854855e-15 )\n",
      "Psi_4 =  8.962147991068151e-07  + i( 1.1107501262046878e-08 )\n"
     ]
    }
   ],
   "source": [
    "Psis = EBW.weyl_psi_scalars(Cdown4, EBW.nup4)\n",
    "print('Psi_0 = ', np.mean(FD.cutoffmask(np.real(Psis[0]))), ' + i(', np.mean(FD.cutoffmask(np.imag(Psis[0]))), ')')\n",
    "print('Psi_1 = ', np.mean(FD.cutoffmask(np.real(Psis[1]))), ' + i(', np.mean(FD.cutoffmask(np.imag(Psis[1]))), ')')\n",
    "print('Psi_2 = ', np.mean(FD.cutoffmask(np.real(Psis[2]))), ' + i(', np.mean(FD.cutoffmask(np.imag(Psis[2]))), ')')\n",
    "print('Psi_3 = ', np.mean(FD.cutoffmask(np.real(Psis[3]))), ' + i(', np.mean(FD.cutoffmask(np.imag(Psis[3]))), ')')\n",
    "print('Psi_4 = ', np.mean(FD.cutoffmask(np.real(Psis[4]))), ' + i(', np.mean(FD.cutoffmask(np.imag(Psis[4]))), ')')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d2164a7",
   "metadata": {},
   "source": [
    "These scalars are computed for an arbitrary null vector base, in order to use these for the Petrov classification we would need a frame rotation for these to be along the principal null directions. Instead we use the scalar invariants:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "586e2e56",
   "metadata": {},
   "outputs": [],
   "source": [
    "invars = EBW.invariant_scalars(Psis)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33a5115f",
   "metadata": {},
   "source": [
    "If $I^3 - 27 J^2 = 0$ it means that this spacetime is special, either of type II, D, III, N, or O,\n",
    "\n",
    "otherwise it is of type I, the most general type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "5f21958b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0515536073234914e-06"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val = invars['I']**3 - 27*invars['J']**2\n",
    "np.mean(FD.cutoffmask(val**(1/6)).real)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f73095d3",
   "metadata": {},
   "source": [
    "Then it is a special type. \n",
    "Note: \n",
    "- here I put the value to the power of 1/6, that is to show a value of the same order of magnitude as the Weyl tensor. \n",
    "- here we are only considering the real part, to do this properly the imaginary part also should be considered.\n",
    "- again I'm showing the mean, that is because this spacetime is homogeneous.\n",
    "- this value will not be exactly zero, most importantly it will tend towards zero as we increase the resolution. One also would need to determine the numerical error with 2 other results obtained with different resolutions.\n",
    "\n",
    "If $I = J = 0$ then the spacetime is either of type III, N, or O, otherwise it is of type II or D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "2c45be35",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I^{1/2} =  8.162565306745274e-07\n",
      "J^{1/3} =  4.7101243316840303e-07\n"
     ]
    }
   ],
   "source": [
    "print('I^{1/2} = ', np.mean(FD.cutoffmask(invars['I']**(1/2)).real))\n",
    "print('J^{1/3} = ', np.mean(FD.cutoffmask(invars['J']**(1/3)).real))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a925d162",
   "metadata": {},
   "source": [
    "Then it is of type III, N, or O.\n",
    "\n",
    "(Note, if it was of type II or D we would check $K$ and $N$ to distinguish them: if $K = N = 0$ then the spacetime is of type D, otherwise it is on type II.)\n",
    "\n",
    "Next if $K = L = 0$ then the spacetime is of type N or O, otherwise it is on type III."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2e8509c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K^{1/3} =  1.5817571581080633e-09\n",
      "L^{1/2} =  2.644143200466357e-07\n",
      "N^{1/4} =  7.885218047976581e-07\n"
     ]
    }
   ],
   "source": [
    "print('K^{1/3} = ', np.mean(FD.cutoffmask(invars['K']**(1/3)).real))\n",
    "print('L^{1/2} = ', np.mean(FD.cutoffmask(invars['L']**(1/2)).real))\n",
    "print('N^{1/4} = ', np.mean(FD.cutoffmask(invars['N']**(1/4)).real))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "834923e0",
   "metadata": {},
   "source": [
    "Then it is of type N, or O.\n",
    "\n",
    "Next if $|B| = |E| = 0$ then the spacetime is of type O, it is conformally flat, otherwise it is of type N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3dd41d9b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|B| =  0.3142695647959811\n",
      "|E| =  0.3142708322346237\n"
     ]
    }
   ],
   "source": [
    "print('|B| = ', np.mean(FD.cutoffmask(Bu_norm)))\n",
    "print('|E| = ', np.mean(FD.cutoffmask(Eu_norm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdeb20a6",
   "metadata": {},
   "source": [
    "As we have shown previously $|B| = |E| \\neq 0$ therefore this spacetime is of Petrov type N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dfdcf60",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9464190",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}