Numerical Tensor Computation

from itensorpy.numerical import NumericalMetric import numpy as np # Define a grid of r values from 2.1M to 20M # (starting outside event horizon for Schwarzschild) r_values = np.linspace(2.1, 20.0, 100) theta_values = np.linspace(0, np.pi, 50) # Create numerical Schwarzschild metric evaluated at each grid point schwarzschild_field = NumericalMetric.schwarzschild( r_values=r_values, theta_values=theta_values, M=1.0 ) # Examine metric component g_tt at a specific point print(f"g_tt at r=3M, theta=π/2: {schwarzschild_field.component(0, 0, r=3.0, theta=np.pi/2)}") # Should output approximately -0.6667 (equals -(1-2M/r) with M=1, r=3)

from itensorpy import Metric from itensorpy.numerical import NumericalTensor import numpy as np # First create symbolic wormhole metric with parameter alpha from sympy import symbols, exp, sin t, r, theta, phi = symbols('t r theta phi') alpha = symbols('alpha', positive=True) # Define metric components g_tt = -exp(2*alpha) g_rr = 1 g_thth = r**2 + alpha**2 g_phph = (r**2 + alpha**2) * sin(theta)**2 # Create symbolic metric symbolic_wormhole = Metric( coords=[t, r, theta, phi], components=[ [g_tt, 0, 0, 0], [0, g_rr, 0, 0], [0, 0, g_thth, 0], [0, 0, 0, g_phph] ] ) # Define coordinate grid r_vals = np.linspace(-10, 10, 100) # Wormhole connects two regions theta_vals = np.linspace(0, np.pi, 50) # Convert to numerical with alpha=1 numerical_wormhole = NumericalTensor.from_symbolic( symbolic_wormhole, grid={ r: r_vals, theta: theta_vals }, params={alpha: 1.0} )

import numpy as np from itensorpy.numerical import NumericalMetric # Create a grid for a 2D metric x = np.linspace(-5, 5, 100) y = np.linspace(-5, 5, 100) X, Y = np.meshgrid(x, y) # Define a simple curved metric (e.g., with a Gaussian "bump") sigma = 2.0 gaussian_bump = np.exp(-(X**2 + Y**2)/(2*sigma**2)) # Define components (3×3 metric: t, x, y) g_tt = -np.ones_like(X) # t component is flat g_xx = 1 + gaussian_bump # x component has a bump g_yy = 1 + gaussian_bump # y component has a bump # Create metric from components curved_2d = NumericalMetric( components=[ [g_tt, np.zeros_like(X), np.zeros_like(X)], [np.zeros_like(X), g_xx, np.zeros_like(X)], [np.zeros_like(X), np.zeros_like(X), g_yy] ], coordinates=['t', 'x', 'y'], coordinate_values={'x': x, 'y': y} )

# Continuing with the Schwarzschild metric field # Calculate the Christoffel symbols numerically christoffel = schwarzschild_field.compute_christoffel_symbols() # Calculate the Riemann curvature tensor riemann = schwarzschild_field.compute_riemann_tensor() # Calculate Kretschmann scalar (RiemannⁿᵐᵏˡRᵢⱼₖₗ) at each point kretschmann = schwarzschild_field.compute_kretschmann_scalar() # The Kretschmann scalar for Schwarzschild is 48M²/r⁶ # Plot it to verify import matplotlib.pyplot as plt plt.figure(figsize=(10, 6)) # Extract values at equator (theta = π/2) equatorial_indices = len(theta_values) // 2 # Middle of theta array k_values = kretschmann[equatorial_indices, :] plt.plot(r_values, k_values, 'b-', label='Numerical') plt.plot(r_values, 48/(r_values**6), 'r--', label='Analytical (48M²/r⁶)') plt.xlabel('Radius (r/M)') plt.ylabel('Kretschmann Scalar (M⁻⁴)') plt.legend() plt.yscale('log') plt.title('Kretschmann Scalar for Schwarzschild Metric') plt.tight_layout() plt.savefig('kretschmann.png', dpi=300)

from itensorpy.numerical import GeodesicSolver import numpy as np import matplotlib.pyplot as plt # Create a solver for Schwarzschild spacetime with M=1 solver = GeodesicSolver.schwarzschild(M=1.0) # Initial conditions for a photon (null geodesic) # Starting from r=20M, θ=π/2 (equatorial plane), φ=0, with specific initial direction r0 = 20.0 theta0 = np.pi/2 phi0 = 0.0 t0 = 0.0 # Initial time # Impact parameter b = 5.2M (critical impact parameter ≈ 5.196M for photon sphere) b = 5.2 # For photons, we normally parameterize with an affine parameter # where the energy at infinity is 1 initial_conditions = { 't': t0, 'r': r0, 'theta': theta0, 'phi': phi0, 'dt_dlambda': 1.0, # Energy at infinity normalized to 1 'dr_dlambda': -np.sqrt(1.0 - (b/r0)**2 * (1.0 - 2.0/r0)), # Incoming radial velocity 'dtheta_dlambda': 0.0, # Motion in equatorial plane 'dphi_dlambda': b/r0**2 # Angular momentum / r² } # Solve the geodesic equations solution = solver.solve( initial_conditions=initial_conditions, lambda_span=[0, 100], # Affine parameter range method='RK45', # Adaptive Runge-Kutta events=[solver.event_horizon_crossing], # Stop if horizon is crossed dense_output=True # Allow interpolation between steps ) # Extract the trajectory t_vals = solution.y[0] r_vals = solution.y[1] phi_vals = solution.y[3] # y[2] is theta, fixed at π/2 # Convert to Cartesian for plotting x_vals = r_vals * np.cos(phi_vals) y_vals = r_vals * np.sin(phi_vals) # Plot trajectory plt.figure(figsize=(10, 10)) # Plot black hole theta = np.linspace(0, 2*np.pi, 100) plt.fill(2*np.cos(theta), 2*np.sin(theta), 'k') # Event horizon at r=2M plt.plot(3*np.cos(theta), 3*np.sin(theta), 'k--') # Photon sphere at r=3M # Plot geodesic plt.plot(x_vals, y_vals, 'r-', label=f'Photon geodesic (b={b}M)') plt.axis('equal') plt.grid(True, alpha=0.3) plt.xlabel('x (M)') plt.ylabel('y (M)') plt.title('Light Deflection in Schwarzschild Spacetime') plt.legend() plt.savefig('light_deflection.png', dpi=300)

from itensorpy.numerical import NumericalTensor, NumericalMetric import numpy as np class BrillWaveMetric(NumericalMetric): """ Implements a Brill wave initial data for numerical relativity with amplitude A and width σ. This represents a gravitational wave packet in axisymmetric form. """ def __init__(self, rho_values, z_values, amplitude=1.0, width=1.0): """ Parameters: rho_values: Array of cylindrical radial coordinates z_values: Array of z coordinates amplitude: Wave amplitude (A) width: Wave width parameter (σ) """ self.rho_values = rho_values self.z_values = z_values self.A = amplitude self.sigma = width # Create mesh grid RHO, Z = np.meshgrid(rho_values, z_values) # Calculate metric components in cylindrical coordinates # using Brill wave ansatz: # ds² = Ψ⁴(dρ² + dz² + ρ²dφ²) # where Ψ is the conformal factor # Define q function (wave shape) r_squared = RHO**2 + Z**2 q = self.A * np.exp(-r_squared/(self.sigma**2)) * RHO**2 # Solve for conformal factor Ψ (simplified here - # real implementation would solve elliptic PDE) psi = 1.0 + q/8.0 psi_to_4 = psi**4 # Define metric components g_rhorho = psi_to_4 g_zz = psi_to_4 g_phiphi = psi_to_4 * RHO**2 # Build the 3×3 spatial metric (rho, z, phi) components = [ [g_rhorho, np.zeros_like(RHO), np.zeros_like(RHO)], [np.zeros_like(RHO), g_zz, np.zeros_like(RHO)], [np.zeros_like(RHO), np.zeros_like(RHO), g_phiphi] ] # Initialize the metric super().__init__( components=components, coordinates=['rho', 'z', 'phi'], coordinate_values={'rho': rho_values, 'z': z_values} ) def compute_ADM_mass(self): """ Compute the ADM mass of the Brill wave by surface integral. """ # Implementation details would go here # Typically involves computing asymptotic behavior of metric pass # Usage rho = np.linspace(0.01, 10, 100) # Avoid rho=0 for numerical stability z = np.linspace(-10, 10, 100) # Create the metric brill = BrillWaveMetric(rho, z, amplitude=2.0, width=1.5) # Compute derived quantities ricci_scalar = brill.compute_scalar_curvature() # Visualize the Ricci scalar import matplotlib.pyplot as plt from matplotlib import cm RHO, Z = np.meshgrid(rho, z) plt.figure(figsize=(10, 8)) plt.pcolormesh(RHO, Z, ricci_scalar, cmap=cm.viridis, shading='auto') plt.colorbar(label='Ricci scalar') plt.xlabel('ρ') plt.ylabel('z') plt.title('Ricci Scalar for Brill Wave (A=2.0, σ=1.5)') plt.tight_layout() plt.savefig('brill_wave_ricci.png', dpi=300)

from itensorpy.numerical import NumericalMetric import numpy as np # Check for CUDA support from itensorpy.utils import cuda_available print(f"CUDA available: {cuda_available()}") # Create a high-resolution grid r = np.linspace(2.1, 100.0, 500) theta = np.linspace(0, np.pi, 500) # Create metric kerr = NumericalMetric.kerr( r_values=r, theta_values=theta, M=1.0, a=0.9 ) # Compute the Riemann tensor with GPU acceleration riemann_gpu = kerr.compute_riemann_tensor(use_gpu=True) print(f"Riemann tensor shape: {riemann_gpu.shape}") # For comparison, time the CPU version on a smaller grid import time r_small = np.linspace(2.1, 100.0, 100) theta_small = np.linspace(0, np.pi, 100) kerr_small = NumericalMetric.kerr(r_values=r_small, theta_values=theta_small, M=1.0, a=0.9) start = time.time() riemann_cpu = kerr_small.compute_riemann_tensor(use_gpu=False) cpu_time = time.time() - start start = time.time() riemann_gpu_small = kerr_small.compute_riemann_tensor(use_gpu=True) gpu_time = time.time() - start print(f"CPU time: {cpu_time:.3f}s, GPU time: {gpu_time:.3f}s") print(f"GPU speedup: {cpu_time/gpu_time:.1f}x")