From 117eda519580f6400ada1488da2f115792f41f55 Mon Sep 17 00:00:00 2001 From: LeetHappyfeet <138872496+LeetHappyfeet@users.noreply.github.com> Date: Mon, 16 Jun 2025 01:54:32 -0400 Subject: [PATCH] Normalize line endings and fix imports --- Metrics/Get_Alcubierre.py | 131 +++++------ Metrics/Get_AlcubierreComoving.py | 193 ++++++++-------- Metrics/Get_Lentz.py | 179 +++++++-------- Metrics/Get_LentzComoving.py | 171 +++++++------- Metrics/Get_Minkowski.py | 77 +++---- Metrics/build_minkowski_metric.py | 53 ++--- Metrics/getScalars.py | 95 ++++---- Metrics/setMinkowski.py | 81 +++---- Metrics/setMinkowskiThreePlusOne.py | 75 ++++--- Metrics/setMinkowskiThreePlusOne_test.py | 83 +++---- Metrics/threePlusOneBuilder.py | 173 +++++++------- Metrics/verifyTensor.py | 1 + Test1.py | 136 +++++------ Test2.py | 144 ++++++------ analysis/generateUniformField.py | 90 ++++---- analysis/getEulerianTransformationMatrix.py | 102 ++++----- analysis/getEvenPointsOnSphere.py | 65 +++--- analysis/getInnerProduct.py | 82 +++---- analysis/getTrace.py | 62 +++--- analyticalEnergyTensor.py | 98 ++++---- compute_energy_tensor_example.py | 81 +++---- solver/getEnergyTensor.py | 3 +- solver/tools/__init__.py | 1 + solver/tools/met2den2.py | 7 + solver/tools/ricciT.py | 5 + solver/tools/ricciTMem2.py | 5 + solver/tools/takeFiniteDifference1_2.py | 4 + solver/tools/takeFiniteDifference2_2.py | 4 + visualizer.py | 235 ++++++++++---------- 29 files changed, 1243 insertions(+), 1193 deletions(-) diff --git a/Metrics/Get_Alcubierre.py b/Metrics/Get_Alcubierre.py index 3cb3c78..81210d7 100644 --- a/Metrics/Get_Alcubierre.py +++ b/Metrics/Get_Alcubierre.py @@ -1,65 +1,66 @@ -# Filename: Get_Alcubierre.py -import numpy as np -import tensorflow as tf - -from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne -from threePlusOneBuilder import threePlusOneBuilder - -# Define shape function for Alcubierre metric -def shapeFunction_Alcubierre(r, R, sigma): - return tf.exp(-((r - R) / sigma) ** 2) - -# Define function to build the Alcubierre metric -def metricGet_Alcubierre(gridSize, worldCenter, v, R, sigma, gridScale): - # Handle default input arguments - if len(gridScale) < 4: - gridScale += [1] * (4 - len(gridScale)) - - # Assign parameters to metric struct - metric = {} - metric['params'] = { - 'gridSize': gridSize, - 'worldCenter': worldCenter, - 'velocity': v, - 'R': R, - 'sigma': sigma - } - - # Assign quantities to metric struct - metric['type'] = "metric" - metric['name'] = 'Alcubierre' - metric['scaling'] = gridScale - metric['coords'] = "cartesian" - metric['index'] = "covariant" - metric['date'] = tf.strings.format("{}", tf.date) - - # Declare a Minkowski space - alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) - - # Add the Alcubierre modification - for i in range(gridSize[1]): - for j in range(gridSize[2]): - for k in range(gridSize[3]): - # Find grid center x, y, z - x = (i + 1) * gridScale[1] - worldCenter[1] - y = (j + 1) * gridScale[2] - worldCenter[2] - z = (k + 1) * gridScale[3] - worldCenter[3] - - for t in range(gridSize[0]): - # Determine the x offset of the center of the bubble, centered in time - xs = (t + 1) * gridScale[0] - worldCenter[0] * v * tf.constant(299792458.0) - - # Find the radius from the center of the bubble - r = tf.sqrt(tf.pow(x - xs, 2) + tf.pow(y, 2) + tf.pow(z, 2)) - - # Find shape function at this point in r - fs = shapeFunction_Alcubierre(r, R, sigma) - - # Add alcubierre modification to shift vector along x - beta[0][t, i, j, k] = -v * fs - - # Make tensor from the 3+1 functions - metric['tensor'] = threePlusOneBuilder(alpha, beta, gamma) - - return metric - +# Filename: Get_Alcubierre.py +import numpy as np +import tensorflow as tf + +from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne +from threePlusOneBuilder import threePlusOneBuilder + +# Define shape function for Alcubierre metric +def shapeFunction_Alcubierre(r, R, sigma): + return tf.exp(-((r - R) / sigma) ** 2) + +# Define function to build the Alcubierre metric +def metricGet_Alcubierre(gridSize, worldCenter, v, R, sigma, gridScale): + # Handle default input arguments + if len(gridScale) < 4: + gridScale += [1] * (4 - len(gridScale)) + + # Assign parameters to metric struct + metric = {} + metric['params'] = { + 'gridSize': gridSize, + 'worldCenter': worldCenter, + 'velocity': v, + 'R': R, + 'sigma': sigma + } + + # Assign quantities to metric struct + metric['type'] = "metric" + metric['name'] = 'Alcubierre' + metric['scaling'] = gridScale + metric['coords'] = "cartesian" + metric['index'] = "covariant" + metric['date'] = tf.strings.format("{}", tf.date) + + # Declare a Minkowski space + alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) + + # Add the Alcubierre modification + for i in range(gridSize[1]): + for j in range(gridSize[2]): + for k in range(gridSize[3]): + # Find grid center x, y, z + x = (i + 1) * gridScale[1] - worldCenter[1] + y = (j + 1) * gridScale[2] - worldCenter[2] + z = (k + 1) * gridScale[3] - worldCenter[3] + + for t in range(gridSize[0]): + # Determine the x offset of the center of the bubble, centered in time + xs = (t + 1) * gridScale[0] - worldCenter[0] * v * tf.constant(299792458.0) + + # Find the radius from the center of the bubble + r = tf.sqrt(tf.pow(x - xs, 2) + tf.pow(y, 2) + tf.pow(z, 2)) + + # Find shape function at this point in r + fs = shapeFunction_Alcubierre(r, R, sigma) + + # Add alcubierre modification to shift vector along x + beta[0][t, i, j, k] = -v * fs + + # Make tensor from the 3+1 functions + metric['tensor'] = threePlusOneBuilder(alpha, beta, gamma) + + return metric + + diff --git a/Metrics/Get_AlcubierreComoving.py b/Metrics/Get_AlcubierreComoving.py index a23691b..fdab7a5 100644 --- a/Metrics/Get_AlcubierreComoving.py +++ b/Metrics/Get_AlcubierreComoving.py @@ -1,96 +1,97 @@ -import numpy as np -import tensorflow as tf -import logging -import traceback -from Metrics.setMinkowskiThreePlusOne import setMinkowskiThreePlusOne -from Metrics.threePlusOneBuilder import threePlusOneBuilder - -logging.basicConfig(level=logging.INFO, filename='singularities.log', filemode='w') - -tf.compat.v1.disable_eager_execution() - -# Define shape function for Alcubierre metric -def shapeFunction_Alcubierre(r, R, sigma): - return tf.exp(-((r - R) / sigma) ** 2) - -# Define function to build the Alcubierre metric in a Galilean comoving frame -def metricGet_AlcubierreComoving(gridSize, worldCenter, v, R, sigma, gridScale): - try: - # Ensure worldCenter has enough elements - if len(worldCenter) < 4: - raise ValueError("worldCenter must have at least 4 elements.") - - # Initialize TensorFlow session - with tf.compat.v1.Session() as sess: - # Define placeholders for r, R, and sigma - r_placeholder = tf.compat.v1.placeholder(tf.float32, shape=(None,)) - R_placeholder = tf.compat.v1.placeholder(tf.float32) - sigma_placeholder = tf.compat.v1.placeholder(tf.float32) - - # Compute shape function - shape_function = shapeFunction_Alcubierre(r_placeholder, R_placeholder, sigma_placeholder) - - # Initialize variables - sess.run(tf.compat.v1.global_variables_initializer()) - - # Collect all r_value values - r_values = [] - for i in range(gridSize[1]): - for j in range(gridSize[2]): - for k in range(gridSize[3]): - x = (i + 1) * gridScale[1] - worldCenter[1] - y = (j + 1) * gridScale[2] - worldCenter[2] - z = (k + 1) * gridScale[3] - worldCenter[3] - r_value = np.sqrt(x**2 + y**2 + z**2) - r_values.append(r_value) - - # Evaluate shape function for all r_values - fs_values = sess.run(shape_function, feed_dict={r_placeholder: r_values, R_placeholder: R, sigma_placeholder: sigma}) - - # Check if gridSize in time is 1 and raise error if not - if gridSize[0] > 1: - raise ValueError('The time grid is greater than 1, only a size of 1 can be used in comoving') - - # Assign parameters to metric struct - metric = {} - metric['params'] = { - 'gridSize': gridSize, - 'worldCenter': worldCenter, - 'velocity': v, - 'R': R, - 'sigma': sigma - } - - # Assign quantities to metric struct - metric['type'] = "metric" - metric['name'] = 'Alcubierre Comoving' - metric['scaling'] = gridScale - metric['coords'] = "cartesian" - metric['index'] = "covariant" - - # Declare a Minkowski space - alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) - - # Add the Alcubierre modification - t = 0 # only one timeslice is used - idx = 0 # index for fs_values - for i in range(gridSize[1]): - for j in range(gridSize[2]): - for k in range(gridSize[3]): - fs_value = fs_values[idx] - beta[0][t, i, j, k] = v * (1 - fs_value) - idx += 1 - - # Make tensor from the 3+1 functions - metric_tensor = threePlusOneBuilder(alpha, beta, gamma) - - if metric_tensor is None: - raise ValueError("Unable to construct metric tensor.") - - metric['tensor'] = metric_tensor - - return metric - except Exception as e: - logging.error(f"Error occurred while computing the metric: {e}") - logging.error(traceback.format_exc()) - return None +import numpy as np +import tensorflow as tf +import logging +import traceback +from Metrics.setMinkowskiThreePlusOne import setMinkowskiThreePlusOne +from Metrics.threePlusOneBuilder import threePlusOneBuilder + +logging.basicConfig(level=logging.INFO, filename='singularities.log', filemode='w') + +tf.compat.v1.disable_eager_execution() + +# Define shape function for Alcubierre metric +def shapeFunction_Alcubierre(r, R, sigma): + return tf.exp(-((r - R) / sigma) ** 2) + +# Define function to build the Alcubierre metric in a Galilean comoving frame +def metricGet_AlcubierreComoving(gridSize, worldCenter, v, R, sigma, gridScale): + try: + # Ensure worldCenter has enough elements + if len(worldCenter) < 4: + raise ValueError("worldCenter must have at least 4 elements.") + + # Initialize TensorFlow session + with tf.compat.v1.Session() as sess: + # Define placeholders for r, R, and sigma + r_placeholder = tf.compat.v1.placeholder(tf.float32, shape=(None,)) + R_placeholder = tf.compat.v1.placeholder(tf.float32) + sigma_placeholder = tf.compat.v1.placeholder(tf.float32) + + # Compute shape function + shape_function = shapeFunction_Alcubierre(r_placeholder, R_placeholder, sigma_placeholder) + + # Initialize variables + sess.run(tf.compat.v1.global_variables_initializer()) + + # Collect all r_value values + r_values = [] + for i in range(gridSize[1]): + for j in range(gridSize[2]): + for k in range(gridSize[3]): + x = (i + 1) * gridScale[1] - worldCenter[1] + y = (j + 1) * gridScale[2] - worldCenter[2] + z = (k + 1) * gridScale[3] - worldCenter[3] + r_value = np.sqrt(x**2 + y**2 + z**2) + r_values.append(r_value) + + # Evaluate shape function for all r_values + fs_values = sess.run(shape_function, feed_dict={r_placeholder: r_values, R_placeholder: R, sigma_placeholder: sigma}) + + # Check if gridSize in time is 1 and raise error if not + if gridSize[0] > 1: + raise ValueError('The time grid is greater than 1, only a size of 1 can be used in comoving') + + # Assign parameters to metric struct + metric = {} + metric['params'] = { + 'gridSize': gridSize, + 'worldCenter': worldCenter, + 'velocity': v, + 'R': R, + 'sigma': sigma + } + + # Assign quantities to metric struct + metric['type'] = "metric" + metric['name'] = 'Alcubierre Comoving' + metric['scaling'] = gridScale + metric['coords'] = "cartesian" + metric['index'] = "covariant" + + # Declare a Minkowski space + alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) + + # Add the Alcubierre modification + t = 0 # only one timeslice is used + idx = 0 # index for fs_values + for i in range(gridSize[1]): + for j in range(gridSize[2]): + for k in range(gridSize[3]): + fs_value = fs_values[idx] + beta[0][t, i, j, k] = v * (1 - fs_value) + idx += 1 + + # Make tensor from the 3+1 functions + metric_tensor = threePlusOneBuilder(alpha, beta, gamma) + + if metric_tensor is None: + raise ValueError("Unable to construct metric tensor.") + + metric['tensor'] = metric_tensor + + return metric + except Exception as e: + logging.error(f"Error occurred while computing the metric: {e}") + logging.error(traceback.format_exc()) + return None + diff --git a/Metrics/Get_Lentz.py b/Metrics/Get_Lentz.py index f05e163..57ee14c 100644 --- a/Metrics/Get_Lentz.py +++ b/Metrics/Get_Lentz.py @@ -1,89 +1,90 @@ -# Filename: metricGet_Lentz.py - - -import numpy as np -import tensorflow as tf -from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne -from threePlusOneBuilder import threePlusOneBuilder - -# Define function to build the Lentz metric -def metricGet_Lentz(gridSize, worldCenter, v, scale, gridScale): - # Handle default input arguments - if len(scale) < 1: - scale = max(gridSize[1:4]) / 7 - if len(gridScale) < 4: - gridScale += [1] * (4 - len(gridScale)) - - # Assign parameters to metric struct - metric = {} - metric['params'] = { - 'gridSize': gridSize, - 'worldCenter': worldCenter, - 'velocity': v - } - - # Assign quantities to metric struct - metric['type'] = "metric" - metric['name'] = "Lentz" - metric['scaling'] = gridScale - metric['coords'] = "cartesian" - metric['index'] = "covariant" - metric['date'] = tf.strings.format("{}", tf.date) - - # Declare a Minkowski space - alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) - - # Lentz Soliton Terms - for i in range(gridSize[1]): - for j in range(gridSize[2]): - for k in range(gridSize[3]): - x = (i + 1) * gridScale[1] - worldCenter[1] - y = (j + 1) * gridScale[2] - worldCenter[2] - - for t in range(gridSize[0]): - # Determine the x offset of the center of the bubble, centered in time - xs = (t + 1) * gridScale[0] - worldCenter[0] * v * tf.constant(299792458.0) - - xp = x - xs - - # Get Lentz template values - WFX, WFY = getWarpFactorByRegion(xp, y, scale) - - # Assign dxdt term - beta[0][t, i, j, k] = -WFX * v - - # Assign dydt term - beta[1][t, i, j, k] = WFY * v - - # Make tensor from the 3+1 functions - metric['tensor'] = threePlusOneBuilder(alpha, beta, gamma) - - return metric - -# Define function to get Lentz warp factor by region -def getWarpFactorByRegion(xIn, yIn, sizeScale): - x = xIn - y = tf.abs(yIn) - WFX = tf.constant(0, dtype=tf.float32) - WFY = tf.constant(0, dtype=tf.float32) - - # Lentz shift vector template - conditions = [ - (x >= sizeScale) & (x <= 2*sizeScale) & (x-sizeScale >= y), - (x > sizeScale) & (x <= 2*sizeScale) & (x-sizeScale <= y) & (-y+3*sizeScale >= x), - (x > 0) & (x <= sizeScale) & (x+sizeScale > y) & (-y+sizeScale < x), - (x > 0) & (x <= sizeScale) & (x+sizeScale <= y) & (-y+3*sizeScale >= x), - (x > -sizeScale) & (x <= 0) & (-x+sizeScale < y) & (-y+3*sizeScale >= -x), - (x > -sizeScale) & (x <= 0) & (x+sizeScale <= y) & (-y+sizeScale >= x), - (x >= -sizeScale) & (x <= sizeScale) & (x+sizeScale > y) - ] - - choices_WFX = [-2, -1, 0, -0.5, 0.5, 1, 1] - choices_WFY = [0, 1, 1, 0.5, 0.5, 0, 0] - - WFX = tf.case([(conditions[i], lambda: choices_WFX[i]) for i in range(len(conditions))], default=lambda: 1.0) - WFY = tf.case([(conditions[i], lambda: choices_WFY[i]) for i in range(len(conditions))], default=lambda: 0.0) - - WFY = tf.sign(yIn) * WFY - - return WFX, WFY \ No newline at end of file +# Filename: metricGet_Lentz.py + + +import numpy as np +import tensorflow as tf +from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne +from threePlusOneBuilder import threePlusOneBuilder + +# Define function to build the Lentz metric +def metricGet_Lentz(gridSize, worldCenter, v, scale, gridScale): + # Handle default input arguments + if len(scale) < 1: + scale = max(gridSize[1:4]) / 7 + if len(gridScale) < 4: + gridScale += [1] * (4 - len(gridScale)) + + # Assign parameters to metric struct + metric = {} + metric['params'] = { + 'gridSize': gridSize, + 'worldCenter': worldCenter, + 'velocity': v + } + + # Assign quantities to metric struct + metric['type'] = "metric" + metric['name'] = "Lentz" + metric['scaling'] = gridScale + metric['coords'] = "cartesian" + metric['index'] = "covariant" + metric['date'] = tf.strings.format("{}", tf.date) + + # Declare a Minkowski space + alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) + + # Lentz Soliton Terms + for i in range(gridSize[1]): + for j in range(gridSize[2]): + for k in range(gridSize[3]): + x = (i + 1) * gridScale[1] - worldCenter[1] + y = (j + 1) * gridScale[2] - worldCenter[2] + + for t in range(gridSize[0]): + # Determine the x offset of the center of the bubble, centered in time + xs = (t + 1) * gridScale[0] - worldCenter[0] * v * tf.constant(299792458.0) + + xp = x - xs + + # Get Lentz template values + WFX, WFY = getWarpFactorByRegion(xp, y, scale) + + # Assign dxdt term + beta[0][t, i, j, k] = -WFX * v + + # Assign dydt term + beta[1][t, i, j, k] = WFY * v + + # Make tensor from the 3+1 functions + metric['tensor'] = threePlusOneBuilder(alpha, beta, gamma) + + return metric + +# Define function to get Lentz warp factor by region +def getWarpFactorByRegion(xIn, yIn, sizeScale): + x = xIn + y = tf.abs(yIn) + WFX = tf.constant(0, dtype=tf.float32) + WFY = tf.constant(0, dtype=tf.float32) + + # Lentz shift vector template + conditions = [ + (x >= sizeScale) & (x <= 2*sizeScale) & (x-sizeScale >= y), + (x > sizeScale) & (x <= 2*sizeScale) & (x-sizeScale <= y) & (-y+3*sizeScale >= x), + (x > 0) & (x <= sizeScale) & (x+sizeScale > y) & (-y+sizeScale < x), + (x > 0) & (x <= sizeScale) & (x+sizeScale <= y) & (-y+3*sizeScale >= x), + (x > -sizeScale) & (x <= 0) & (-x+sizeScale < y) & (-y+3*sizeScale >= -x), + (x > -sizeScale) & (x <= 0) & (x+sizeScale <= y) & (-y+sizeScale >= x), + (x >= -sizeScale) & (x <= sizeScale) & (x+sizeScale > y) + ] + + choices_WFX = [-2, -1, 0, -0.5, 0.5, 1, 1] + choices_WFY = [0, 1, 1, 0.5, 0.5, 0, 0] + + WFX = tf.case([(conditions[i], lambda: choices_WFX[i]) for i in range(len(conditions))], default=lambda: 1.0) + WFY = tf.case([(conditions[i], lambda: choices_WFY[i]) for i in range(len(conditions))], default=lambda: 0.0) + + WFY = tf.sign(yIn) * WFY + + return WFX, WFY + diff --git a/Metrics/Get_LentzComoving.py b/Metrics/Get_LentzComoving.py index 37357b2..ee10644 100644 --- a/Metrics/Get_LentzComoving.py +++ b/Metrics/Get_LentzComoving.py @@ -1,85 +1,86 @@ -import numpy as np -import tensorflow as tf -from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne -from threePlusOneBuilder import threePlusOneBuilder - -# Define function to build the Lentz metric in Galilean comoving frame -def metricGet_LentzComoving(gridSize, worldCenter, v, scale, gridScale): - # Handle default input arguments - if len(scale) < 1: - scale = max(gridSize[1:4]) / 7 - if len(gridScale) < 4: - gridScale += [1] * (4 - len(gridScale)) - - # Check if gridSize in time is 1 and raise error if not - if gridSize[0] > 1: - raise ValueError('The time grid is greater than 1, only a size of 1 can be used for Lentz') - - # Assign parameters to metric struct - metric = {} - metric['params'] = { - 'gridSize': gridSize, - 'worldCenter': worldCenter, - 'velocity': v - } - - # Assign quantities to metric struct - metric['type'] = "metric" - metric['name'] = "Lentz Comoving" - metric['scaling'] = gridScale - metric['coords'] = "cartesian" - metric['index'] = "covariant" - metric['date'] = tf.strings.format("{}", tf.date) - - # Declare a Minkowski space - alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) - t = 0 # only one timeslice is used - - # Lentz Soliton Terms - for i in range(gridSize[1]): - for j in range(gridSize[2]): - for k in range(gridSize[3]): - x = (i + 1) * gridScale[1] - worldCenter[1] - y = (j + 1) * gridScale[2] - worldCenter[2] - - # Get Lentz template values - WFX, WFY = getWarpFactorByRegion(x, y, scale) - - # Assign dxdt term - beta[0][t, i, j, k] = v * (1 - WFX) - - # Assign dydt term - beta[1][t, i, j, k] = v * WFY - - # Make tensor from the 3+1 functions - metric['tensor'] = threePlusOneBuilder(alpha, beta, gamma) - - return metric - -# Define function to get Lentz warp factor by region -def getWarpFactorByRegion(xIn, yIn, sizeScale): - x = xIn - y = tf.abs(yIn) - WFX = tf.constant(0, dtype=tf.float32) - WFY = tf.constant(0, dtype=tf.float32) - - # Lentz shift vector template - conditions = [ - (x >= sizeScale) & (x <= 2*sizeScale) & (x-sizeScale >= y), - (x > sizeScale) & (x <= 2*sizeScale) & (x-sizeScale <= y) & (-y+3*sizeScale >= x), - (x > 0) & (x <= sizeScale) & (x+sizeScale > y) & (-y+sizeScale < x), - (x > 0) & (x <= sizeScale) & (x+sizeScale <= y) & (-y+3*sizeScale >= x), - (x > -sizeScale) & (x <= 0) & (-x+sizeScale < y) & (-y+3*sizeScale >= -x), - (x > -sizeScale) & (x <= 0) & (x+sizeScale <= y) & (-y+sizeScale >= x), - (x >= -sizeScale) & (x <= sizeScale) & (x+sizeScale > y) - ] - - choices_WFX = [-2, -1, 0, -0.5, 0.5, 1, 1] - choices_WFY = [0, 1, 1, 0.5, 0.5, 0, 0] - - WFX = tf.case([(conditions[i], lambda: choices_WFX[i]) for i in range(len(conditions))], default=lambda: 1.0) - WFY = tf.case([(conditions[i], lambda: choices_WFY[i]) for i in range(len(conditions))], default=lambda: 0.0) - - WFY = tf.sign(yIn) * WFY - - return WFX, WFY \ No newline at end of file +import numpy as np +import tensorflow as tf +from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne +from threePlusOneBuilder import threePlusOneBuilder + +# Define function to build the Lentz metric in Galilean comoving frame +def metricGet_LentzComoving(gridSize, worldCenter, v, scale, gridScale): + # Handle default input arguments + if len(scale) < 1: + scale = max(gridSize[1:4]) / 7 + if len(gridScale) < 4: + gridScale += [1] * (4 - len(gridScale)) + + # Check if gridSize in time is 1 and raise error if not + if gridSize[0] > 1: + raise ValueError('The time grid is greater than 1, only a size of 1 can be used for Lentz') + + # Assign parameters to metric struct + metric = {} + metric['params'] = { + 'gridSize': gridSize, + 'worldCenter': worldCenter, + 'velocity': v + } + + # Assign quantities to metric struct + metric['type'] = "metric" + metric['name'] = "Lentz Comoving" + metric['scaling'] = gridScale + metric['coords'] = "cartesian" + metric['index'] = "covariant" + metric['date'] = tf.strings.format("{}", tf.date) + + # Declare a Minkowski space + alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) + t = 0 # only one timeslice is used + + # Lentz Soliton Terms + for i in range(gridSize[1]): + for j in range(gridSize[2]): + for k in range(gridSize[3]): + x = (i + 1) * gridScale[1] - worldCenter[1] + y = (j + 1) * gridScale[2] - worldCenter[2] + + # Get Lentz template values + WFX, WFY = getWarpFactorByRegion(x, y, scale) + + # Assign dxdt term + beta[0][t, i, j, k] = v * (1 - WFX) + + # Assign dydt term + beta[1][t, i, j, k] = v * WFY + + # Make tensor from the 3+1 functions + metric['tensor'] = threePlusOneBuilder(alpha, beta, gamma) + + return metric + +# Define function to get Lentz warp factor by region +def getWarpFactorByRegion(xIn, yIn, sizeScale): + x = xIn + y = tf.abs(yIn) + WFX = tf.constant(0, dtype=tf.float32) + WFY = tf.constant(0, dtype=tf.float32) + + # Lentz shift vector template + conditions = [ + (x >= sizeScale) & (x <= 2*sizeScale) & (x-sizeScale >= y), + (x > sizeScale) & (x <= 2*sizeScale) & (x-sizeScale <= y) & (-y+3*sizeScale >= x), + (x > 0) & (x <= sizeScale) & (x+sizeScale > y) & (-y+sizeScale < x), + (x > 0) & (x <= sizeScale) & (x+sizeScale <= y) & (-y+3*sizeScale >= x), + (x > -sizeScale) & (x <= 0) & (-x+sizeScale < y) & (-y+3*sizeScale >= -x), + (x > -sizeScale) & (x <= 0) & (x+sizeScale <= y) & (-y+sizeScale >= x), + (x >= -sizeScale) & (x <= sizeScale) & (x+sizeScale > y) + ] + + choices_WFX = [-2, -1, 0, -0.5, 0.5, 1, 1] + choices_WFY = [0, 1, 1, 0.5, 0.5, 0, 0] + + WFX = tf.case([(conditions[i], lambda: choices_WFX[i]) for i in range(len(conditions))], default=lambda: 1.0) + WFY = tf.case([(conditions[i], lambda: choices_WFY[i]) for i in range(len(conditions))], default=lambda: 0.0) + + WFY = tf.sign(yIn) * WFY + + return WFX, WFY + diff --git a/Metrics/Get_Minkowski.py b/Metrics/Get_Minkowski.py index c43c266..7c9ac5e 100644 --- a/Metrics/Get_Minkowski.py +++ b/Metrics/Get_Minkowski.py @@ -1,38 +1,39 @@ -# Filename: metricGet_Minkowski.py - -import numpy as np -import tensorflow as tf - -# Define function to build the Minkowski metric -def metricGet_Minkowski(gridSize, gridScaling): - # Handle default input arguments - if len(gridScaling) < 4: - gridScaling += [1] * (4 - len(gridScaling)) - - # Assign quantities to metric struct - metric = {} - metric['type'] = "metric" - metric['name'] = "Minkowski" - metric['scaling'] = gridScaling - metric['coords'] = "cartesian" - metric['index'] = "covariant" - metric['date'] = tf.strings.format("{}", tf.date) - - # Initialize tensor dictionary - metric['tensor'] = {} - - # dt^2 term - metric['tensor'][(1, 1)] = -tf.ones(gridSize) - - # Non-time diagonal terms - metric['tensor'][(2, 2)] = tf.ones(gridSize) - metric['tensor'][(3, 3)] = tf.ones(gridSize) - metric['tensor'][(4, 4)] = tf.ones(gridSize) - - # Cross terms - for i in range(1, 5): - for j in range(1, 5): - if i != j: - metric['tensor'][(i, j)] = tf.zeros(gridSize) - - return metric \ No newline at end of file +# Filename: metricGet_Minkowski.py + +import numpy as np +import tensorflow as tf + +# Define function to build the Minkowski metric +def metricGet_Minkowski(gridSize, gridScaling): + # Handle default input arguments + if len(gridScaling) < 4: + gridScaling += [1] * (4 - len(gridScaling)) + + # Assign quantities to metric struct + metric = {} + metric['type'] = "metric" + metric['name'] = "Minkowski" + metric['scaling'] = gridScaling + metric['coords'] = "cartesian" + metric['index'] = "covariant" + metric['date'] = tf.strings.format("{}", tf.date) + + # Initialize tensor dictionary + metric['tensor'] = {} + + # dt^2 term + metric['tensor'][(1, 1)] = -tf.ones(gridSize) + + # Non-time diagonal terms + metric['tensor'][(2, 2)] = tf.ones(gridSize) + metric['tensor'][(3, 3)] = tf.ones(gridSize) + metric['tensor'][(4, 4)] = tf.ones(gridSize) + + # Cross terms + for i in range(1, 5): + for j in range(1, 5): + if i != j: + metric['tensor'][(i, j)] = tf.zeros(gridSize) + + return metric + diff --git a/Metrics/build_minkowski_metric.py b/Metrics/build_minkowski_metric.py index b90fb1b..1961324 100644 --- a/Metrics/build_minkowski_metric.py +++ b/Metrics/build_minkowski_metric.py @@ -1,26 +1,27 @@ -# Filename: build_minkowski_metric.py - -import tensorflow as tf - -def build_minkowski_metric(grid_size, grid_scaling=None): - if grid_scaling is None: - grid_scaling = [1, 1, 1, 1] - - metric = {} - metric['type'] = "metric" - metric['name'] = "Minkowski" - metric['scaling'] = grid_scaling - metric['coords'] = "cartesian" - metric['index'] = "covariant" - metric['date'] = tf.strings.format("{}", tf.date) - - # Initialize tensor components - for i in range(1, 5): - for j in range(1, 5): - if i == j: - val = -tf.ones(grid_size) if i == 1 else tf.ones(grid_size) - else: - val = tf.zeros(grid_size) - metric[f"tensor_{i},{j}"] = val - - return metric \ No newline at end of file +# Filename: build_minkowski_metric.py + +import tensorflow as tf + +def build_minkowski_metric(grid_size, grid_scaling=None): + if grid_scaling is None: + grid_scaling = [1, 1, 1, 1] + + metric = {} + metric['type'] = "metric" + metric['name'] = "Minkowski" + metric['scaling'] = grid_scaling + metric['coords'] = "cartesian" + metric['index'] = "covariant" + metric['date'] = tf.strings.format("{}", tf.date) + + # Initialize tensor components + for i in range(1, 5): + for j in range(1, 5): + if i == j: + val = -tf.ones(grid_size) if i == 1 else tf.ones(grid_size) + else: + val = tf.zeros(grid_size) + metric[f"tensor_{i},{j}"] = val + + return metric + diff --git a/Metrics/getScalars.py b/Metrics/getScalars.py index afc2a77..05a0fc3 100644 --- a/Metrics/getScalars.py +++ b/Metrics/getScalars.py @@ -1,47 +1,48 @@ -# Filename: getScalars.py -import numpy as np -import tensorflow as tf -from tensor_utils import c3_inv, change_tensor_index - -def three_plus_one_decomposer(metric): - """ - Finds 3+1 terms from the metric tensor. - - Args: - metric (dict): Dictionary containing the metric tensor. - - Returns: - alpha (numpy.ndarray): 4D array representing the lapse rate. - beta_down (list): List of 3 4D arrays representing the covariant shift vectors. - gamma_down (list): List of 3x3 4D arrays representing the covariant spatial terms. - beta_up (list): List of 3 4D arrays representing the contravariant shift vectors. - gamma_up (list): List of 3x3 4D arrays representing the contravariant spatial terms. - """ - - # Check that the metric is covariant and change index if not - metric = change_tensor_index(metric, "covariant") - - # Covariant shift vector maps to the covariant tensor terms g_0i - beta_down = [metric['tensor'][0][i+1] for i in range(3)] - - # Covariant gamma maps to the covariant tensor terms g_ij - gamma_down = [[metric['tensor'][i+1][j+1] for j in range(3)] for i in range(3)] - - # Transform gamma to contravariant - gamma_up = c3_inv(gamma_down) - - # Find the world grid size - s = np.shape(metric['tensor'][0][0]) - - # Transform beta to contravariant - beta_up = [] - for i in range(3): - temp = np.zeros(s) - for j in range(3): - temp += gamma_up[i][j] * beta_down[j] - beta_up.append(temp) - - # Find lapse using beta and g_00 - alpha = tf.sqrt(beta_up[0] * beta_down[0] + beta_up[1] * beta_down[1] + beta_up[2] * beta_down[2] - metric['tensor'][0][0]) - - return alpha, beta_down, gamma_down, beta_up, gamma_up \ No newline at end of file +# Filename: getScalars.py +import numpy as np +import tensorflow as tf +from tensor_utils import c3_inv, change_tensor_index + +def three_plus_one_decomposer(metric): + """ + Finds 3+1 terms from the metric tensor. + + Args: + metric (dict): Dictionary containing the metric tensor. + + Returns: + alpha (numpy.ndarray): 4D array representing the lapse rate. + beta_down (list): List of 3 4D arrays representing the covariant shift vectors. + gamma_down (list): List of 3x3 4D arrays representing the covariant spatial terms. + beta_up (list): List of 3 4D arrays representing the contravariant shift vectors. + gamma_up (list): List of 3x3 4D arrays representing the contravariant spatial terms. + """ + + # Check that the metric is covariant and change index if not + metric = change_tensor_index(metric, "covariant") + + # Covariant shift vector maps to the covariant tensor terms g_0i + beta_down = [metric['tensor'][0][i+1] for i in range(3)] + + # Covariant gamma maps to the covariant tensor terms g_ij + gamma_down = [[metric['tensor'][i+1][j+1] for j in range(3)] for i in range(3)] + + # Transform gamma to contravariant + gamma_up = c3_inv(gamma_down) + + # Find the world grid size + s = np.shape(metric['tensor'][0][0]) + + # Transform beta to contravariant + beta_up = [] + for i in range(3): + temp = np.zeros(s) + for j in range(3): + temp += gamma_up[i][j] * beta_down[j] + beta_up.append(temp) + + # Find lapse using beta and g_00 + alpha = tf.sqrt(beta_up[0] * beta_down[0] + beta_up[1] * beta_down[1] + beta_up[2] * beta_down[2] - metric['tensor'][0][0]) + + return alpha, beta_down, gamma_down, beta_up, gamma_up + diff --git a/Metrics/setMinkowski.py b/Metrics/setMinkowski.py index 97a09bc..2a36e69 100644 --- a/Metrics/setMinkowski.py +++ b/Metrics/setMinkowski.py @@ -1,40 +1,41 @@ -# Filename: setMinkowski.py -import numpy as np - -def setMinkowski(gridSize): - """ - Builds metric terms for a flat Minkowski space. - - Args: - gridSize (tuple): World size in [t, x, y, z]. - - Returns: - metric (dict): The metric tensor as a dictionary. - """ - - # Initialize the metric tensor as a dictionary - metric = {} - - # dt^2 term - metric[(1, 1)] = -np.ones(gridSize) - - # Non-time diagonal terms - metric[(2, 2)] = np.ones(gridSize) - metric[(3, 3)] = np.ones(gridSize) - metric[(4, 4)] = np.ones(gridSize) - - # Cross terms - metric[(1, 2)] = np.zeros(gridSize) - metric[(2, 1)] = np.zeros(gridSize) - metric[(1, 3)] = np.zeros(gridSize) - metric[(3, 1)] = np.zeros(gridSize) - metric[(2, 3)] = np.zeros(gridSize) - metric[(3, 2)] = np.zeros(gridSize) - metric[(2, 4)] = np.zeros(gridSize) - metric[(4, 2)] = np.zeros(gridSize) - metric[(3, 4)] = np.zeros(gridSize) - metric[(4, 3)] = np.zeros(gridSize) - metric[(1, 4)] = np.zeros(gridSize) - metric[(4, 1)] = np.zeros(gridSize) - - return metric \ No newline at end of file +# Filename: setMinkowski.py +import numpy as np + +def setMinkowski(gridSize): + """ + Builds metric terms for a flat Minkowski space. + + Args: + gridSize (tuple): World size in [t, x, y, z]. + + Returns: + metric (dict): The metric tensor as a dictionary. + """ + + # Initialize the metric tensor as a dictionary + metric = {} + + # dt^2 term + metric[(1, 1)] = -np.ones(gridSize) + + # Non-time diagonal terms + metric[(2, 2)] = np.ones(gridSize) + metric[(3, 3)] = np.ones(gridSize) + metric[(4, 4)] = np.ones(gridSize) + + # Cross terms + metric[(1, 2)] = np.zeros(gridSize) + metric[(2, 1)] = np.zeros(gridSize) + metric[(1, 3)] = np.zeros(gridSize) + metric[(3, 1)] = np.zeros(gridSize) + metric[(2, 3)] = np.zeros(gridSize) + metric[(3, 2)] = np.zeros(gridSize) + metric[(2, 4)] = np.zeros(gridSize) + metric[(4, 2)] = np.zeros(gridSize) + metric[(3, 4)] = np.zeros(gridSize) + metric[(4, 3)] = np.zeros(gridSize) + metric[(1, 4)] = np.zeros(gridSize) + metric[(4, 1)] = np.zeros(gridSize) + + return metric + diff --git a/Metrics/setMinkowskiThreePlusOne.py b/Metrics/setMinkowskiThreePlusOne.py index c8a1d94..de60601 100644 --- a/Metrics/setMinkowskiThreePlusOne.py +++ b/Metrics/setMinkowskiThreePlusOne.py @@ -1,37 +1,38 @@ -# Filename: setMinkowskiThreePlusOne.py -import numpy as np -import logging - -logging.basicConfig(level=logging.INFO) - -def setMinkowskiThreePlusOne(gridSize): - """ - Returns the 3+1 format for flat space. - - Args: - gridSize (tuple): World size in [t, x, y, z]. - - Returns: - alpha (numpy.ndarray): Lapse rate 4D array. - beta (list): Shift vector, 1x3 list of 4D arrays. - gamma (list): Spatial terms, 3x3 list of 4D arrays. - """ - - # Lapse rate alpha - alpha = np.ones(gridSize) - - # Shift vector beta - beta = [np.zeros(gridSize) for _ in range(3)] - - # Spatial terms gamma - gamma = [[np.ones(gridSize) if i == j else np.zeros(gridSize) for j in range(3)] for i in range(3)] - - logging.info(f"Alpha shape: {alpha.shape}") - for i in range(3): - logging.info(f"Beta[{i}] shape: {beta[i].shape}") - for i in range(3): - for j in range(3): - logging.info(f"Gamma[{i}][{j}] shape: {gamma[i][j].shape}") - logging.info(f"Gamma[{i}][{j}] initial values: {gamma[i][j][0,0,0]}") - - return alpha, beta, gamma \ No newline at end of file +# Filename: setMinkowskiThreePlusOne.py +import numpy as np +import logging + +logging.basicConfig(level=logging.INFO) + +def setMinkowskiThreePlusOne(gridSize): + """ + Returns the 3+1 format for flat space. + + Args: + gridSize (tuple): World size in [t, x, y, z]. + + Returns: + alpha (numpy.ndarray): Lapse rate 4D array. + beta (list): Shift vector, 1x3 list of 4D arrays. + gamma (list): Spatial terms, 3x3 list of 4D arrays. + """ + + # Lapse rate alpha + alpha = np.ones(gridSize) + + # Shift vector beta + beta = [np.zeros(gridSize) for _ in range(3)] + + # Spatial terms gamma + gamma = [[np.ones(gridSize) if i == j else np.zeros(gridSize) for j in range(3)] for i in range(3)] + + logging.info(f"Alpha shape: {alpha.shape}") + for i in range(3): + logging.info(f"Beta[{i}] shape: {beta[i].shape}") + for i in range(3): + for j in range(3): + logging.info(f"Gamma[{i}][{j}] shape: {gamma[i][j].shape}") + logging.info(f"Gamma[{i}][{j}] initial values: {gamma[i][j][0,0,0]}") + + return alpha, beta, gamma + diff --git a/Metrics/setMinkowskiThreePlusOne_test.py b/Metrics/setMinkowskiThreePlusOne_test.py index 0b1feae..99297c9 100644 --- a/Metrics/setMinkowskiThreePlusOne_test.py +++ b/Metrics/setMinkowskiThreePlusOne_test.py @@ -1,41 +1,42 @@ -import unittest -from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne -import numpy as np - -class TestSetMinkowskiThreePlusOne(unittest.TestCase): - def test_output_shapes(self): - # Define grid size for testing - gridSize = (2, 3, 3, 3) # Arbitrary grid size - - # Call the function to get alpha, beta, and gamma - alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) - - # Check shapes - self.assertEqual(alpha.shape, gridSize, "Alpha shape mismatch") - self.assertEqual(len(beta), 3, "Number of beta components mismatch") - for i in range(3): - self.assertEqual(beta[i].shape, gridSize, f"Beta[{i}] shape mismatch") - for i in range(3): - for j in range(3): - self.assertEqual(gamma[i][j].shape, gridSize, f"Gamma[{i}][{j}] shape mismatch") - - def test_initial_values(self): - # Define grid size for testing - gridSize = (2, 3, 3, 3) # Arbitrary grid size - - # Call the function to get alpha, beta, and gamma - alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) - - # Check initial values - self.assertTrue(np.all(alpha == 1), "Incorrect initial values for alpha") - for i in range(3): - self.assertTrue(np.all(beta[i] == 0), f"Incorrect initial values for beta[{i}]") - for i in range(3): - for j in range(3): - if i == j: - self.assertTrue(np.all(gamma[i][j] == 1), f"Incorrect initial values for gamma[{i}][{j}] diagonal") - else: - self.assertTrue(np.all(gamma[i][j] == 0), f"Incorrect initial values for gamma[{i}][{j}] off-diagonal") - -if __name__ == '__main__': - unittest.main() \ No newline at end of file +import unittest +from setMinkowskiThreePlusOne import setMinkowskiThreePlusOne +import numpy as np + +class TestSetMinkowskiThreePlusOne(unittest.TestCase): + def test_output_shapes(self): + # Define grid size for testing + gridSize = (2, 3, 3, 3) # Arbitrary grid size + + # Call the function to get alpha, beta, and gamma + alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) + + # Check shapes + self.assertEqual(alpha.shape, gridSize, "Alpha shape mismatch") + self.assertEqual(len(beta), 3, "Number of beta components mismatch") + for i in range(3): + self.assertEqual(beta[i].shape, gridSize, f"Beta[{i}] shape mismatch") + for i in range(3): + for j in range(3): + self.assertEqual(gamma[i][j].shape, gridSize, f"Gamma[{i}][{j}] shape mismatch") + + def test_initial_values(self): + # Define grid size for testing + gridSize = (2, 3, 3, 3) # Arbitrary grid size + + # Call the function to get alpha, beta, and gamma + alpha, beta, gamma = setMinkowskiThreePlusOne(gridSize) + + # Check initial values + self.assertTrue(np.all(alpha == 1), "Incorrect initial values for alpha") + for i in range(3): + self.assertTrue(np.all(beta[i] == 0), f"Incorrect initial values for beta[{i}]") + for i in range(3): + for j in range(3): + if i == j: + self.assertTrue(np.all(gamma[i][j] == 1), f"Incorrect initial values for gamma[{i}][{j}] diagonal") + else: + self.assertTrue(np.all(gamma[i][j] == 0), f"Incorrect initial values for gamma[{i}][{j}] off-diagonal") + +if __name__ == '__main__': + unittest.main() + diff --git a/Metrics/threePlusOneBuilder.py b/Metrics/threePlusOneBuilder.py index f977ef1..e966974 100644 --- a/Metrics/threePlusOneBuilder.py +++ b/Metrics/threePlusOneBuilder.py @@ -1,86 +1,87 @@ -import numpy as np -from Metrics.verifyTensor import verifyTensor -import logging - -logging.basicConfig(level=logging.INFO) - -def is_symmetric(tensor, tol=1e-10): - for (i, j), value in tensor.items(): - if not np.allclose(value, tensor.get((j, i), np.zeros_like(value)), atol=tol): - logging.error(f"Tensor component {(i, j)} is not symmetric with {(j, i)}.") - return False - return True - -def threePlusOneBuilder(alpha, beta, gamma, threshold=1e-10): - """ - Builds the metric tensor given input 3+1 components of alpha, beta, and gamma. - - Args: - alpha (numpy.ndarray): Lapse rate map across spacetime. - beta (list of numpy.ndarray): Covariant shift vector map across spacetime. - gamma (list of list of numpy.ndarray): Covariant spatial term map across spacetime. - threshold (float): Threshold value for the signature criterion. - - Returns: - metricTensor (dict): Metric tensor represented as a dictionary. - """ - - logging.info(f"Shape of alpha: {alpha.shape}") - for i in range(3): - logging.info(f"Shape of beta[{i}]: {beta[i].shape}") - - for i in range(3): - for j in range(3): - logging.info(f"Shape of gamma[{i}][{j}]: {gamma[i][j].shape}") - logging.info(f"gamma[{i}][{j}] sample values: {gamma[i][j][0,0,0]}") - - gamma_up = [] - - for idx, g in enumerate(gamma): - dets = np.linalg.det(g) - logging.info(f"Determinant of gamma[{idx}]: {dets}") - - if np.any(np.abs(dets) < threshold): - logging.warning(f"Singular matrix detected at gamma[{idx}] with determinant {dets}. Using pseudoinverse instead.") - gamma_up.append(np.linalg.pinv(g)) # Use pseudo-inverse for singular matrices - else: - gamma_up.append(np.linalg.inv(g)) - - # Find gridSize - s = alpha.shape - - # Calculate beta_i - beta_up = [np.zeros(s) for _ in range(3)] - for i in range(3): - for j in range(3): - beta_up[i] += gamma_up[i][j] * beta[j] - - # Create time-time component - metricTensor = {(1, 1): -alpha ** 2} - for i in range(3): - metricTensor[(1, 1)] += beta_up[i] * beta[i] - - # Create time-space components - for i in range(2, 5): - metricTensor[(1, i)] = beta[i - 2] - metricTensor[(i, 1)] = beta[i - 2] - - # Create space-space components - for i in range(2, 5): - for j in range(2, 5): - metricTensor[(i, j)] = gamma[i - 2][j - 2] - - # Log the metric tensor before verification - logging.info(f"Constructed metric tensor: {metricTensor}") - - # Verify the symmetry of the metric tensor - if not is_symmetric(metricTensor, threshold): - logging.error("Constructed metric tensor is not symmetric.") - raise ValueError("Constructed metric tensor does not satisfy symmetry criteria.") - - # Verify the metric tensor using external function - if not verifyTensor(metricTensor, threshold): - logging.error(f"Metric tensor failed verification: {metricTensor}") - raise ValueError("Constructed metric tensor does not satisfy signature criteria.") - - return metricTensor \ No newline at end of file +import numpy as np +from Metrics.verifyTensor import verifyTensor +import logging + +logging.basicConfig(level=logging.INFO) + +def is_symmetric(tensor, tol=1e-10): + for (i, j), value in tensor.items(): + if not np.allclose(value, tensor.get((j, i), np.zeros_like(value)), atol=tol): + logging.error(f"Tensor component {(i, j)} is not symmetric with {(j, i)}.") + return False + return True + +def threePlusOneBuilder(alpha, beta, gamma, threshold=1e-10): + """ + Builds the metric tensor given input 3+1 components of alpha, beta, and gamma. + + Args: + alpha (numpy.ndarray): Lapse rate map across spacetime. + beta (list of numpy.ndarray): Covariant shift vector map across spacetime. + gamma (list of list of numpy.ndarray): Covariant spatial term map across spacetime. + threshold (float): Threshold value for the signature criterion. + + Returns: + metricTensor (dict): Metric tensor represented as a dictionary. + """ + + logging.info(f"Shape of alpha: {alpha.shape}") + for i in range(3): + logging.info(f"Shape of beta[{i}]: {beta[i].shape}") + + for i in range(3): + for j in range(3): + logging.info(f"Shape of gamma[{i}][{j}]: {gamma[i][j].shape}") + logging.info(f"gamma[{i}][{j}] sample values: {gamma[i][j][0,0,0]}") + + gamma_up = [] + + for idx, g in enumerate(gamma): + dets = np.linalg.det(g) + logging.info(f"Determinant of gamma[{idx}]: {dets}") + + if np.any(np.abs(dets) < threshold): + logging.warning(f"Singular matrix detected at gamma[{idx}] with determinant {dets}. Using pseudoinverse instead.") + gamma_up.append(np.linalg.pinv(g)) # Use pseudo-inverse for singular matrices + else: + gamma_up.append(np.linalg.inv(g)) + + # Find gridSize + s = alpha.shape + + # Calculate beta_i + beta_up = [np.zeros(s) for _ in range(3)] + for i in range(3): + for j in range(3): + beta_up[i] += gamma_up[i][j] * beta[j] + + # Create time-time component + metricTensor = {(1, 1): -alpha ** 2} + for i in range(3): + metricTensor[(1, 1)] += beta_up[i] * beta[i] + + # Create time-space components + for i in range(2, 5): + metricTensor[(1, i)] = beta[i - 2] + metricTensor[(i, 1)] = beta[i - 2] + + # Create space-space components + for i in range(2, 5): + for j in range(2, 5): + metricTensor[(i, j)] = gamma[i - 2][j - 2] + + # Log the metric tensor before verification + logging.info(f"Constructed metric tensor: {metricTensor}") + + # Verify the symmetry of the metric tensor + if not is_symmetric(metricTensor, threshold): + logging.error("Constructed metric tensor is not symmetric.") + raise ValueError("Constructed metric tensor does not satisfy symmetry criteria.") + + # Verify the metric tensor using external function + if not verifyTensor(metricTensor, threshold): + logging.error(f"Metric tensor failed verification: {metricTensor}") + raise ValueError("Constructed metric tensor does not satisfy signature criteria.") + + return metricTensor + diff --git a/Metrics/verifyTensor.py b/Metrics/verifyTensor.py index faf6839..5244391 100644 --- a/Metrics/verifyTensor.py +++ b/Metrics/verifyTensor.py @@ -25,3 +25,4 @@ def verifyTensor(metricTensor, threshold): return False return True + diff --git a/Test1.py b/Test1.py index 923cbed..60eba44 100644 --- a/Test1.py +++ b/Test1.py @@ -1,68 +1,68 @@ -import sympy as sp -import numpy as np -import matplotlib.pyplot as plt -from analyticalEnergyTensor import analyticalEnergyTensor -import time - -def monitor_progress(step, total_steps, start_time): - elapsed_time = time.time() - start_time - progress = (step / total_steps) * 100 - estimated_total_time = (elapsed_time / step) * total_steps if step != 0 else 0 - estimated_remaining_time = estimated_total_time - elapsed_time - print(f"Step {step}/{total_steps} - {progress:.2f}% complete - Elapsed Time: {elapsed_time:.2f}s - Estimated Remaining Time: {estimated_remaining_time:.2f}s") - -# Define symbolic variables -t, x, y, z = sp.symbols('t x y z') - -# Define an arbitrary metric tensor (for example, a weak gravitational wave) -g = sp.zeros(4, 4) -g[0, 0] = -1 -g[1, 1] = 1 + 0.1 * sp.sin(x) * sp.sin(t) -g[2, 2] = 1 + 0.1 * sp.sin(y) * sp.sin(t) -g[3, 3] = 1 + 0.1 * sp.sin(z) * sp.sin(t) -g[1, 2] = 0.05 * sp.cos(x) * sp.cos(y) * sp.cos(t) -g[2, 1] = g[1, 2] # Symmetric component - -# Define the coordinates -coords = [t, x, y, z] - -# Start time for progress monitoring -start_time = time.time() - -# Compute the stress-energy tensor -T_out = analyticalEnergyTensor(g, coords) - -# Print the resulting stress-energy tensor component-wise -print("Stress-Energy Tensor:") -for i in range(4): - for j in range(4): - print(f"T[{i},{j}] = {T_out[i, j]}") - -# Progress monitoring -total_steps = 1 # Since this is a single calculation, only one step -monitor_progress(1, total_steps, start_time) - -# Define a numerical grid for evaluation -x_vals = np.linspace(-5, 5, 100) -y_vals = np.linspace(-5, 5, 100) -x_grid, y_grid = np.meshgrid(x_vals, y_vals) -t_val = 0 # Fixed time value for visualization -z_val = 0 # Fixed z value for visualization - -# Function to evaluate tensor components numerically -def evaluate_component(component, t_val, x_grid, y_grid, z_val): - component_func = sp.lambdify((t, x, y, z), component, 'numpy') - return component_func(t_val, x_grid, y_grid, z_val) - -# Evaluate a specific component of the tensor for visualization -component = T_out[1, 2] -component_values = evaluate_component(component, t_val, x_grid, y_grid, z_val) - -# Plot the component values -plt.figure(figsize=(10, 8)) -plt.contourf(x_grid, y_grid, component_values, cmap='viridis') -plt.colorbar(label='T[1,1] values') -plt.xlabel('x') -plt.ylabel('y') -plt.title('Stress-Energy Tensor Component T[1,2]') -plt.show() \ No newline at end of file +import sympy as sp +import numpy as np +import matplotlib.pyplot as plt +from analyticalEnergyTensor import analyticalEnergyTensor +import time + +def monitor_progress(step, total_steps, start_time): + elapsed_time = time.time() - start_time + progress = (step / total_steps) * 100 + estimated_total_time = (elapsed_time / step) * total_steps if step != 0 else 0 + estimated_remaining_time = estimated_total_time - elapsed_time + print(f"Step {step}/{total_steps} - {progress:.2f}% complete - Elapsed Time: {elapsed_time:.2f}s - Estimated Remaining Time: {estimated_remaining_time:.2f}s") + +# Define symbolic variables +t, x, y, z = sp.symbols('t x y z') + +# Define an arbitrary metric tensor (for example, a weak gravitational wave) +g = sp.zeros(4, 4) +g[0, 0] = -1 +g[1, 1] = 1 + 0.1 * sp.sin(x) * sp.sin(t) +g[2, 2] = 1 + 0.1 * sp.sin(y) * sp.sin(t) +g[3, 3] = 1 + 0.1 * sp.sin(z) * sp.sin(t) +g[1, 2] = 0.05 * sp.cos(x) * sp.cos(y) * sp.cos(t) +g[2, 1] = g[1, 2] # Symmetric component + +# Define the coordinates +coords = [t, x, y, z] + +# Start time for progress monitoring +start_time = time.time() + +# Compute the stress-energy tensor +T_out = analyticalEnergyTensor(g, coords) + +# Print the resulting stress-energy tensor component-wise +print("Stress-Energy Tensor:") +for i in range(4): + for j in range(4): + print(f"T[{i},{j}] = {T_out[i, j]}") + +# Progress monitoring +total_steps = 1 # Since this is a single calculation, only one step +monitor_progress(1, total_steps, start_time) + +# Define a numerical grid for evaluation +x_vals = np.linspace(-5, 5, 100) +y_vals = np.linspace(-5, 5, 100) +x_grid, y_grid = np.meshgrid(x_vals, y_vals) +t_val = 0 # Fixed time value for visualization +z_val = 0 # Fixed z value for visualization + +# Function to evaluate tensor components numerically +def evaluate_component(component, t_val, x_grid, y_grid, z_val): + component_func = sp.lambdify((t, x, y, z), component, 'numpy') + return component_func(t_val, x_grid, y_grid, z_val) + +# Evaluate a specific component of the tensor for visualization +component = T_out[1, 2] +component_values = evaluate_component(component, t_val, x_grid, y_grid, z_val) + +# Plot the component values +plt.figure(figsize=(10, 8)) +plt.contourf(x_grid, y_grid, component_values, cmap='viridis') +plt.colorbar(label='T[1,1] values') +plt.xlabel('x') +plt.ylabel('y') +plt.title('Stress-Energy Tensor Component T[1,2]') +plt.show() diff --git a/Test2.py b/Test2.py index 8ace59a..7fac7e3 100644 --- a/Test2.py +++ b/Test2.py @@ -1,72 +1,72 @@ -import sympy as sp -import numpy as np -import matplotlib.pyplot as plt -from analyticalEnergyTensor import analyticalEnergyTensor -import time - -def monitor_progress(step, total_steps, start_time): - elapsed_time = time.time() - start_time - progress = (step / total_steps) * 100 - estimated_total_time = (elapsed_time / step) * total_steps if step != 0 else 0 - estimated_remaining_time = estimated_total_time - elapsed_time - print(f"Step {step}/{total_steps} - {progress:.2f}% complete - Elapsed Time: {elapsed_time:.2f}s - Estimated Remaining Time: {estimated_remaining_time:.2f}s") - -# Define symbolic variables -t, x, y, z = sp.symbols('t x y z') - -# Define an arbitrary metric tensor (for example, a weak gravitational wave) -g = sp.zeros(4, 4) -g[0, 0] = -1 -g[1, 1] = 1 + 0.1 * sp.sin(x) * sp.sin(t) -g[2, 2] = 1 + 0.1 * sp.sin(y) * sp.sin(t) -g[3, 3] = 1 + 0.1 * sp.sin(z) * sp.sin(t) -g[1, 2] = 0.05 * sp.cos(x) * sp.cos(y) * sp.cos(t) -g[2, 1] = g[1, 2] # Symmetric component - -# Define the coordinates -coords = [t, x, y, z] - -# Start time for progress monitoring -start_time = time.time() - -# Compute the stress-energy tensor -T_out = analyticalEnergyTensor(g, coords) - -# Print the resulting stress-energy tensor component-wise -print("Stress-Energy Tensor:") -for i in range(4): - for j in range(4): - print(f"T[{i},{j}] = {T_out[i, j]}") - -# Progress monitoring -total_steps = 1 # Since this is a single calculation, only one step -monitor_progress(1, total_steps, start_time) - -# Define a numerical grid for evaluation -x_vals = np.linspace(-5, 5, 100) -y_vals = np.linspace(-5, 5, 100) -x_grid, y_grid = np.meshgrid(x_vals, y_vals) -t_val = 0 # Fixed time value for visualization -z_val = 0 # Fixed z value for visualization - -# Function to evaluate tensor components numerically -def evaluate_component(component, t_val, x_grid, y_grid, z_val): - component_func = sp.lambdify((t, x, y, z), component, 'numpy') - return component_func(t_val, x_grid, y_grid, z_val) - -# Visualize multiple components of the stress-energy tensor -components_to_plot = [(1, 1), (1, 2), (2, 2), (3, 3)] # List of tensor components to visualize -fig, axs = plt.subplots(2, 2, figsize=(12, 10)) - -for ax, (i, j) in zip(axs.flatten(), components_to_plot): - component = T_out[i, j] - component_values = evaluate_component(component, t_val, x_grid, y_grid, z_val) - - c = ax.contourf(x_grid, y_grid, component_values, cmap='viridis') - fig.colorbar(c, ax=ax) - ax.set_title(f'T[{i},{j}]') - ax.set_xlabel('x') - ax.set_ylabel('y') - -plt.tight_layout() -plt.show() \ No newline at end of file +import sympy as sp +import numpy as np +import matplotlib.pyplot as plt +from analyticalEnergyTensor import analyticalEnergyTensor +import time + +def monitor_progress(step, total_steps, start_time): + elapsed_time = time.time() - start_time + progress = (step / total_steps) * 100 + estimated_total_time = (elapsed_time / step) * total_steps if step != 0 else 0 + estimated_remaining_time = estimated_total_time - elapsed_time + print(f"Step {step}/{total_steps} - {progress:.2f}% complete - Elapsed Time: {elapsed_time:.2f}s - Estimated Remaining Time: {estimated_remaining_time:.2f}s") + +# Define symbolic variables +t, x, y, z = sp.symbols('t x y z') + +# Define an arbitrary metric tensor (for example, a weak gravitational wave) +g = sp.zeros(4, 4) +g[0, 0] = -1 +g[1, 1] = 1 + 0.1 * sp.sin(x) * sp.sin(t) +g[2, 2] = 1 + 0.1 * sp.sin(y) * sp.sin(t) +g[3, 3] = 1 + 0.1 * sp.sin(z) * sp.sin(t) +g[1, 2] = 0.05 * sp.cos(x) * sp.cos(y) * sp.cos(t) +g[2, 1] = g[1, 2] # Symmetric component + +# Define the coordinates +coords = [t, x, y, z] + +# Start time for progress monitoring +start_time = time.time() + +# Compute the stress-energy tensor +T_out = analyticalEnergyTensor(g, coords) + +# Print the resulting stress-energy tensor component-wise +print("Stress-Energy Tensor:") +for i in range(4): + for j in range(4): + print(f"T[{i},{j}] = {T_out[i, j]}") + +# Progress monitoring +total_steps = 1 # Since this is a single calculation, only one step +monitor_progress(1, total_steps, start_time) + +# Define a numerical grid for evaluation +x_vals = np.linspace(-5, 5, 100) +y_vals = np.linspace(-5, 5, 100) +x_grid, y_grid = np.meshgrid(x_vals, y_vals) +t_val = 0 # Fixed time value for visualization +z_val = 0 # Fixed z value for visualization + +# Function to evaluate tensor components numerically +def evaluate_component(component, t_val, x_grid, y_grid, z_val): + component_func = sp.lambdify((t, x, y, z), component, 'numpy') + return component_func(t_val, x_grid, y_grid, z_val) + +# Visualize multiple components of the stress-energy tensor +components_to_plot = [(1, 1), (1, 2), (2, 2), (3, 3)] # List of tensor components to visualize +fig, axs = plt.subplots(2, 2, figsize=(12, 10)) + +for ax, (i, j) in zip(axs.flatten(), components_to_plot): + component = T_out[i, j] + component_values = evaluate_component(component, t_val, x_grid, y_grid, z_val) + + c = ax.contourf(x_grid, y_grid, component_values, cmap='viridis') + fig.colorbar(c, ax=ax) + ax.set_title(f'T[{i},{j}]') + ax.set_xlabel('x') + ax.set_ylabel('y') + +plt.tight_layout() +plt.show() diff --git a/analysis/generateUniformField.py b/analysis/generateUniformField.py index 9bb1eaa..523042c 100644 --- a/analysis/generateUniformField.py +++ b/analysis/generateUniformField.py @@ -1,44 +1,46 @@ -# Filename: generateUniformField.py -import tensorflow as tf -import numpy as np - -def generateUniformField(type, numAngularVec, numTimeVec, tryGPU): - - if not (type.lower() == "nulllike" or type.lower() == "timelike"): - raise ValueError('Vector field type not generated, use either: "nulllike", "timelike"') - - if type.lower() == "timelike": - # Generate timelike vectors c^2t^2 > r^2 - bb = tf.linspace(0.0, 1.0, numTimeVec) - VecField = tf.ones((4, numAngularVec, numTimeVec)) - - for jj in range(numTimeVec): - # Build vector field in Cartesian coordinates - spherical_points = tf.stack([tf.ones(numAngularVec), getEvenPointsOnSphere(1.0 - bb[jj], numAngularVec, 1)], axis=0) - VecField[:,:,jj].assign(spherical_points / tf.norm(spherical_points, axis=0)) - - elif type.lower() == "nulllike": - # Build vector field in Cartesian coordinates - spherical_points = getEvenPointsOnSphere(1.0, numAngularVec, 1) - VecField = tf.ones((4, numAngularVec)) - VecField[1:,:] = spherical_points - VecField = VecField / tf.norm(VecField, axis=0) - - # Convert to GPU array if requested - if tryGPU: - VecField = tf.Variable(VecField) - - return VecField - -def getEvenPointsOnSphere(radius, numPoints, dim): - # Function to generate evenly spaced points on a sphere - theta = np.linspace(0, np.pi, numPoints) - phi = np.linspace(0, 2 * np.pi, numPoints) - theta, phi = np.meshgrid(theta, phi) - - x = radius * np.sin(theta) * np.cos(phi) - y = radius * np.sin(theta) * np.sin(phi) - z = radius * np.cos(theta) - - spherical_points = np.stack([x, y, z], axis=-1) - return tf.convert_to_tensor(spherical_points, dtype=tf.float32) +# Filename: generateUniformField.py +import tensorflow as tf +import numpy as np + +def generateUniformField(type, numAngularVec, numTimeVec, tryGPU): + + if not (type.lower() == "nulllike" or type.lower() == "timelike"): + raise ValueError('Vector field type not generated, use either: "nulllike", "timelike"') + + if type.lower() == "timelike": + # Generate timelike vectors c^2t^2 > r^2 + bb = tf.linspace(0.0, 1.0, numTimeVec) + VecField = tf.ones((4, numAngularVec, numTimeVec)) + + for jj in range(numTimeVec): + # Build vector field in Cartesian coordinates + spherical_points = tf.stack([tf.ones(numAngularVec), getEvenPointsOnSphere(1.0 - bb[jj], numAngularVec, 1)], axis=0) + VecField[:,:,jj].assign(spherical_points / tf.norm(spherical_points, axis=0)) + + elif type.lower() == "nulllike": + # Build vector field in Cartesian coordinates + spherical_points = getEvenPointsOnSphere(1.0, numAngularVec, 1) + VecField = tf.ones((4, numAngularVec)) + VecField[1:,:] = spherical_points + VecField = VecField / tf.norm(VecField, axis=0) + + # Convert to GPU array if requested + if tryGPU: + VecField = tf.Variable(VecField) + + return VecField + +def getEvenPointsOnSphere(radius, numPoints, dim): + # Function to generate evenly spaced points on a sphere + theta = np.linspace(0, np.pi, numPoints) + phi = np.linspace(0, 2 * np.pi, numPoints) + theta, phi = np.meshgrid(theta, phi) + + x = radius * np.sin(theta) * np.cos(phi) + y = radius * np.sin(theta) * np.sin(phi) + z = radius * np.cos(theta) + + spherical_points = np.stack([x, y, z], axis=-1) + return tf.convert_to_tensor(spherical_points, dtype=tf.float32) + + diff --git a/analysis/getEulerianTransformationMatrix.py b/analysis/getEulerianTransformationMatrix.py index 37b8c2c..c588565 100644 --- a/analysis/getEulerianTransformationMatrix.py +++ b/analysis/getEulerianTransformationMatrix.py @@ -1,50 +1,52 @@ -# Filename: getEulerianTransformationMatrix.py -import numpy as np - -def get_eulerian_transformation_matrix(g): - if g.ndim == 2: # For simple 4x4 metrics - assert np.array_equal(g, g.T), "g is not symmetric" - - factor0 = g[3, 3] - factor1 = (-g[2, 3] ** 2 + g[2, 2] * factor0) - factor2 = (2 * g[1, 2] * g[1, 3] * g[2, 3] - g[3, 3] * g[1, 2] ** 2 - g[2, 2] * g[1, 3] ** 2 + g[1, 1] * factor1) - factor3 = (-2 * g[3, 3] * g[0, 1] * g[0, 2] * g[1, 2] + 2 * g[0, 2] * g[0, 3] * g[1, 2] * g[2, 3] + - 2 * g[0, 1] * g[0, 2] * g[1, 3] * g[2, 3] + 2 * g[0, 1] * g[0, 3] * g[1, 2] * g[2, 3] - - g[0, 1] ** 2 * g[2, 3] ** 2 - g[0, 2] ** 2 * g[1, 3] ** 2 - g[0, 3] ** 2 * g[1, 2] ** 2 + - g[2, 2] * (-2 * g[0, 1] * g[0, 3] * g[2, 3] + g[3, 3] * g[0, 1] ** 2) + - g[1, 1] * (-2 * g[0, 2] * g[0, 3] * g[2, 3] + g[3, 3] * g[0, 2] ** 2 + g[2, 2] * g[0, 3] ** 2) - - g[0, 0] * factor2) - - M = np.zeros((4, 4)) - M[0, 0] = np.sqrt(factor2 / factor3) - M[1, 0] = (g[0, 1] * g[2, 3] ** 2 + g[0, 2] * g[1, 2] * g[3, 3] - g[0, 2] * g[1, 3] * g[2, 3] - - g[0, 3] * g[1, 2] * g[2, 3] + g[0, 3] * g[1, 3] * g[2, 2] - g[0, 1] * g[2, 2] * g[3, 3]) / np.sqrt( - factor2 * factor3) - M[2, 0] = (g[0, 2] * g[1, 3] ** 2 - g[0, 3] * g[1, 2] * g[1, 3] + g[0, 1] * g[1, 2] * g[3, 3] - - g[0, 1] * g[1, 3] * g[2, 3] - g[0, 2] * g[1, 1] * g[3, 3] + g[0, 3] * g[1, 1] * g[2, 3]) / np.sqrt( - factor2 * factor3) - M[3, 0] = (g[0, 3] * g[1, 2] ** 2 - g[0, 2] * g[1, 2] * g[1, 3] - g[0, 1] * g[1, 2] * g[2, 3] + - g[0, 1] * g[1, 3] * g[2, 2] + g[0, 2] * g[1, 1] * g[2, 3] - g[0, 3] * g[1, 1] * g[2, 2]) / np.sqrt( - factor2 * factor3) - - M[1, 1] = (g[2, 3] ** 2 - g[1, 2] ** 2 - g[1, 3] ** 2 + g[1, 1] * factor0) / np.sqrt(factor0 * factor2) - M[2, 1] = (-g[1, 2] * g[1, 3] + g[1, 1] * g[2, 3]) / np.sqrt(factor0 * factor2) - M[3, 1] = (-g[1, 3] * g[2, 3] + g[1, 2] * factor0) / np.sqrt(factor0 * factor2) - - M[2, 2] = np.sqrt(factor0 / factor1) - M[3, 2] = (-g[2, 3]) / np.sqrt(factor0 * factor1) - - M[3, 3] = np.sqrt(factor0) - - - - return M - -# Example usage: -g = np.array([[1, 0, 0, 0], - [0, -1, 0, 0], - [0, 0, -1, 0], - [0, 0, 0, -1]]) - -M = get_eulerian_transformation_matrix(g) -print(M) +# Filename: getEulerianTransformationMatrix.py +import numpy as np + +def get_eulerian_transformation_matrix(g): + if g.ndim == 2: # For simple 4x4 metrics + assert np.array_equal(g, g.T), "g is not symmetric" + + factor0 = g[3, 3] + factor1 = (-g[2, 3] ** 2 + g[2, 2] * factor0) + factor2 = (2 * g[1, 2] * g[1, 3] * g[2, 3] - g[3, 3] * g[1, 2] ** 2 - g[2, 2] * g[1, 3] ** 2 + g[1, 1] * factor1) + factor3 = (-2 * g[3, 3] * g[0, 1] * g[0, 2] * g[1, 2] + 2 * g[0, 2] * g[0, 3] * g[1, 2] * g[2, 3] + + 2 * g[0, 1] * g[0, 2] * g[1, 3] * g[2, 3] + 2 * g[0, 1] * g[0, 3] * g[1, 2] * g[2, 3] - + g[0, 1] ** 2 * g[2, 3] ** 2 - g[0, 2] ** 2 * g[1, 3] ** 2 - g[0, 3] ** 2 * g[1, 2] ** 2 + + g[2, 2] * (-2 * g[0, 1] * g[0, 3] * g[2, 3] + g[3, 3] * g[0, 1] ** 2) + + g[1, 1] * (-2 * g[0, 2] * g[0, 3] * g[2, 3] + g[3, 3] * g[0, 2] ** 2 + g[2, 2] * g[0, 3] ** 2) - + g[0, 0] * factor2) + + M = np.zeros((4, 4)) + M[0, 0] = np.sqrt(factor2 / factor3) + M[1, 0] = (g[0, 1] * g[2, 3] ** 2 + g[0, 2] * g[1, 2] * g[3, 3] - g[0, 2] * g[1, 3] * g[2, 3] - + g[0, 3] * g[1, 2] * g[2, 3] + g[0, 3] * g[1, 3] * g[2, 2] - g[0, 1] * g[2, 2] * g[3, 3]) / np.sqrt( + factor2 * factor3) + M[2, 0] = (g[0, 2] * g[1, 3] ** 2 - g[0, 3] * g[1, 2] * g[1, 3] + g[0, 1] * g[1, 2] * g[3, 3] - + g[0, 1] * g[1, 3] * g[2, 3] - g[0, 2] * g[1, 1] * g[3, 3] + g[0, 3] * g[1, 1] * g[2, 3]) / np.sqrt( + factor2 * factor3) + M[3, 0] = (g[0, 3] * g[1, 2] ** 2 - g[0, 2] * g[1, 2] * g[1, 3] - g[0, 1] * g[1, 2] * g[2, 3] + + g[0, 1] * g[1, 3] * g[2, 2] + g[0, 2] * g[1, 1] * g[2, 3] - g[0, 3] * g[1, 1] * g[2, 2]) / np.sqrt( + factor2 * factor3) + + M[1, 1] = (g[2, 3] ** 2 - g[1, 2] ** 2 - g[1, 3] ** 2 + g[1, 1] * factor0) / np.sqrt(factor0 * factor2) + M[2, 1] = (-g[1, 2] * g[1, 3] + g[1, 1] * g[2, 3]) / np.sqrt(factor0 * factor2) + M[3, 1] = (-g[1, 3] * g[2, 3] + g[1, 2] * factor0) / np.sqrt(factor0 * factor2) + + M[2, 2] = np.sqrt(factor0 / factor1) + M[3, 2] = (-g[2, 3]) / np.sqrt(factor0 * factor1) + + M[3, 3] = np.sqrt(factor0) + + + + return M + +# Example usage: +g = np.array([[1, 0, 0, 0], + [0, -1, 0, 0], + [0, 0, -1, 0], + [0, 0, 0, -1]]) + +if __name__ == "__main__": + M = get_eulerian_transformation_matrix(g) + print(M) + diff --git a/analysis/getEvenPointsOnSphere.py b/analysis/getEvenPointsOnSphere.py index 485fe69..552e081 100644 --- a/analysis/getEvenPointsOnSphere.py +++ b/analysis/getEvenPointsOnSphere.py @@ -1,32 +1,33 @@ -# Filename: getEvenPointsOnSphere.py -import numpy as np - -def getEvenPointsOnSphere(R, numberOfPoints, useGPU=False): - """ - GETEVENPOINTSONSPHERE: Generate evenly distributed points on a sphere. - - Args: - R (float): Radius of the sphere. - numberOfPoints (int): Number of points to generate. - useGPU (bool): Flag indicating whether to use GPU computation (default is False). - - Returns: - Vector (numpy.ndarray): Array of shape (3, numberOfPoints) containing the generated points. - """ - - if useGPU: - goldenRatio = (1 + 5 ** 0.5) / 2 - Vector = np.zeros((3, numberOfPoints), dtype=np.float32) - else: - goldenRatio = (1 + 5 ** 0.5) / 2 - Vector = np.zeros((3, numberOfPoints), dtype=np.float64) - - for i in range(numberOfPoints): - theta = 2 * np.pi * i / goldenRatio - phi = np.arccos(1 - 2 * (i + 0.5) / numberOfPoints) - - Vector[0, i] = R * np.cos(theta) * np.sin(phi) - Vector[1, i] = R * np.sin(theta) * np.sin(phi) - Vector[2, i] = R * np.cos(phi) - - return Vector +# Filename: getEvenPointsOnSphere.py +import numpy as np + +def getEvenPointsOnSphere(R, numberOfPoints, useGPU=False): + """ + GETEVENPOINTSONSPHERE: Generate evenly distributed points on a sphere. + + Args: + R (float): Radius of the sphere. + numberOfPoints (int): Number of points to generate. + useGPU (bool): Flag indicating whether to use GPU computation (default is False). + + Returns: + Vector (numpy.ndarray): Array of shape (3, numberOfPoints) containing the generated points. + """ + + if useGPU: + goldenRatio = (1 + 5 ** 0.5) / 2 + Vector = np.zeros((3, numberOfPoints), dtype=np.float32) + else: + goldenRatio = (1 + 5 ** 0.5) / 2 + Vector = np.zeros((3, numberOfPoints), dtype=np.float64) + + for i in range(numberOfPoints): + theta = 2 * np.pi * i / goldenRatio + phi = np.arccos(1 - 2 * (i + 0.5) / numberOfPoints) + + Vector[0, i] = R * np.cos(theta) * np.sin(phi) + Vector[1, i] = R * np.sin(theta) * np.sin(phi) + Vector[2, i] = R * np.cos(phi) + + return Vector + diff --git a/analysis/getInnerProduct.py b/analysis/getInnerProduct.py index def57cc..d0096d4 100644 --- a/analysis/getInnerProduct.py +++ b/analysis/getInnerProduct.py @@ -1,40 +1,42 @@ -# Filename: getInnerProduct.py -import numpy as np -from verifyTensor import verifyTensor - -##ensure that you have defined the verifyTensor and c4Inv functions or integrate their functionalities within this code. -##Also, make sure that the input vecA, vecB, and Metric are dictionaries containing the appropriate keys. - -def getInnerProduct(vecA, vecB, Metric): - """ - GETINNERPRODUCT: takes the inner product of two vector fields with their metric - - Args: - vecA (dict): Dictionary containing a vector field represented by a list of arrays. - vecB (dict): Dictionary containing a vector field represented by a list of arrays. - Metric (dict): Dictionary containing a metric tensor represented by a 4x4 array. - - Returns: - innerprod (numpy.ndarray): Inner product of the two vector fields. - """ - - # Verify Metric - if not verifyTensor(Metric, 1): - raise ValueError("Metric is not verified. Please verify metric using verifyTensor(metric).") - - s = np.shape(Metric['tensor'][0][0]) - innerprod = np.zeros(s) - - if vecA['index'] != vecB['index']: - for mu in range(4): - for nu in range(4): - innerprod += vecA['field'][mu] * vecB['field'][nu] - - elif vecA['index'] == vecB['index']: - if vecA['index'] == Metric['index']: - Metric['tensor'] = c4Inv(Metric['tensor']) # flip index - for mu in range(4): - for nu in range(4): - innerprod += vecA['field'][mu] * vecB['field'][nu] * Metric['tensor'][mu][nu] - - return innerprod +# Filename: getInnerProduct.py +import numpy as np +from solver.verifyTensor import verifyTensor +from solver.tools.c4Inv import c4Inv2 as c4Inv + +##ensure that you have defined the verifyTensor and c4Inv functions or integrate their functionalities within this code. +##Also, make sure that the input vecA, vecB, and Metric are dictionaries containing the appropriate keys. + +def getInnerProduct(vecA, vecB, Metric): + """ + GETINNERPRODUCT: takes the inner product of two vector fields with their metric + + Args: + vecA (dict): Dictionary containing a vector field represented by a list of arrays. + vecB (dict): Dictionary containing a vector field represented by a list of arrays. + Metric (dict): Dictionary containing a metric tensor represented by a 4x4 array. + + Returns: + innerprod (numpy.ndarray): Inner product of the two vector fields. + """ + + # Verify Metric + if not verifyTensor(Metric, 1): + raise ValueError("Metric is not verified. Please verify metric using verifyTensor(metric).") + + s = np.shape(Metric['tensor'][0][0]) + innerprod = np.zeros(s) + + if vecA['index'] != vecB['index']: + for mu in range(4): + for nu in range(4): + innerprod += vecA['field'][mu] * vecB['field'][nu] + + elif vecA['index'] == vecB['index']: + if vecA['index'] == Metric['index']: + Metric['tensor'] = c4Inv(Metric['tensor']) # flip index + for mu in range(4): + for nu in range(4): + innerprod += vecA['field'][mu] * vecB['field'][nu] * Metric['tensor'][mu][nu] + + return innerprod + diff --git a/analysis/getTrace.py b/analysis/getTrace.py index 03972ee..36b7e44 100644 --- a/analysis/getTrace.py +++ b/analysis/getTrace.py @@ -1,30 +1,32 @@ -# Filename: getTrace -import numpy as np -from verifyTensor import verifyTensor - -def getTrace(tensor, metric): - """ - GETTRACE: take the trace of a tensor - - Args: - tensor (dict): Dictionary containing a tensor represented by a 4x4 array. - metric (dict): Dictionary containing a metric tensor represented by a 4x4 array. - - Returns: - Trace (numpy.ndarray): Trace of the tensor. - """ - - # Verify metric - if not verifyTensor(metric, 1): - raise ValueError("Metric is not verified. Please verify metric using verifyTensor(metric).") - - Trace = np.zeros_like(metric['tensor'][0][0]) - - if tensor['index'] == metric['index']: - metric['tensor'] = c4Inv(metric['tensor']) - - for a in range(4): - for b in range(4): - Trace += metric['tensor'][a][b] * tensor['tensor'][a][b] - - return Trace +# Filename: getTrace +import numpy as np +from solver.verifyTensor import verifyTensor +from solver.tools.c4Inv import c4Inv2 as c4Inv + +def getTrace(tensor, metric): + """ + GETTRACE: take the trace of a tensor + + Args: + tensor (dict): Dictionary containing a tensor represented by a 4x4 array. + metric (dict): Dictionary containing a metric tensor represented by a 4x4 array. + + Returns: + Trace (numpy.ndarray): Trace of the tensor. + """ + + # Verify metric + if not verifyTensor(metric, 1): + raise ValueError("Metric is not verified. Please verify metric using verifyTensor(metric).") + + Trace = np.zeros_like(metric['tensor'][0][0]) + + if tensor['index'] == metric['index']: + metric['tensor'] = c4Inv(metric['tensor']) + + for a in range(4): + for b in range(4): + Trace += metric['tensor'][a][b] * tensor['tensor'][a][b] + + return Trace + diff --git a/analyticalEnergyTensor.py b/analyticalEnergyTensor.py index c3f66b3..3d50c20 100644 --- a/analyticalEnergyTensor.py +++ b/analyticalEnergyTensor.py @@ -1,49 +1,49 @@ -import sympy as sp - -def analyticalEnergyTensor(g, coords): - # Calculate inverse metric tensor - gi = g.inv() - - # Compute the Christoffel symbols - Christoffel = [[[sp.S(0) for _ in range(4)] for _ in range(4)] for _ in range(4)] - for i in range(4): - for j in range(4): - for k in range(4): - Christoffel[i][j][k] = sum( - 0.5 * gi[i, d] * (g[d, j].diff(coords[k]) + g[d, k].diff(coords[j]) - g[j, k].diff(coords[d])) - for d in range(4) - ) - - # Compute the Ricci tensor - Ricci = sp.zeros(4, 4) - for i in range(4): - for j in range(4): - Ricci[i, j] = sum( - Christoffel[a][i][j].diff(coords[a]) - Christoffel[a][i][a].diff(coords[j]) + - sum( - Christoffel[a][i][j] * Christoffel[b][a][b] - Christoffel[a][i][b] * Christoffel[b][j][a] - for b in range(4) - ) - for a in range(4) - ) - - # Compute the Ricci Scalar - Ricci_Scalar = sum(Ricci[i, j] * gi[i, j] for i in range(4) for j in range(4)) - - # Compute the Stress-Energy Tensor - T = sp.zeros(4, 4) - for i in range(4): - for j in range(4): - T[i, j] = (Ricci[i, j] - 0.5 * Ricci_Scalar * g[i, j]) - - # Convert to Contravariant - T_out = sp.zeros(4, 4) - for i in range(4): - for j in range(4): - T_out[i, j] = sum( - T[a, b] * gi[a, i] * gi[b, j] - for a in range(4) - for b in range(4) - ) - - return T_out \ No newline at end of file +import sympy as sp + +def analyticalEnergyTensor(g, coords): + # Calculate inverse metric tensor + gi = g.inv() + + # Compute the Christoffel symbols + Christoffel = [[[sp.S(0) for _ in range(4)] for _ in range(4)] for _ in range(4)] + for i in range(4): + for j in range(4): + for k in range(4): + Christoffel[i][j][k] = sum( + 0.5 * gi[i, d] * (g[d, j].diff(coords[k]) + g[d, k].diff(coords[j]) - g[j, k].diff(coords[d])) + for d in range(4) + ) + + # Compute the Ricci tensor + Ricci = sp.zeros(4, 4) + for i in range(4): + for j in range(4): + Ricci[i, j] = sum( + Christoffel[a][i][j].diff(coords[a]) - Christoffel[a][i][a].diff(coords[j]) + + sum( + Christoffel[a][i][j] * Christoffel[b][a][b] - Christoffel[a][i][b] * Christoffel[b][j][a] + for b in range(4) + ) + for a in range(4) + ) + + # Compute the Ricci Scalar + Ricci_Scalar = sum(Ricci[i, j] * gi[i, j] for i in range(4) for j in range(4)) + + # Compute the Stress-Energy Tensor + T = sp.zeros(4, 4) + for i in range(4): + for j in range(4): + T[i, j] = (Ricci[i, j] - 0.5 * Ricci_Scalar * g[i, j]) + + # Convert to Contravariant + T_out = sp.zeros(4, 4) + for i in range(4): + for j in range(4): + T_out[i, j] = sum( + T[a, b] * gi[a, i] * gi[b, j] + for a in range(4) + for b in range(4) + ) + + return T_out diff --git a/compute_energy_tensor_example.py b/compute_energy_tensor_example.py index fb2d534..33d917b 100644 --- a/compute_energy_tensor_example.py +++ b/compute_energy_tensor_example.py @@ -1,40 +1,41 @@ -from Metrics.Get_AlcubierreComoving import metricGet_AlcubierreComoving -from solver.getEnergyTensor import getEnergyTensor -from solver.metric_definition import metric -import time - -# Define monitoring function -def monitor_progress(start_time, step, total_steps): - elapsed_time = time.time() - start_time - estimated_total_time = elapsed_time / step * total_steps - estimated_remaining_time = estimated_total_time - elapsed_time - print(f"Step {step}/{total_steps} - Elapsed Time: {elapsed_time:.2f}s - Estimated Remaining Time: {estimated_remaining_time:.2f}s") - -# Define parameters for the Alcubierre metric -gridSize = [1, 100, 100, 100] -worldCenter = [0, 0, 0, 0] # You need to add the fourth element here -v = 0.5 -r_value = 5.0 -R_value = 10.0 -sigma_value = 2.0 -gridScale = [1, 0.1, 0.1, 0.1] - -# Define R and sigma -R = R_value -sigma = sigma_value - -# Check lengths of gridScale and worldCenter -print(len(gridScale)) -print(len(worldCenter)) - -# Get the Alcubierre metric -total_steps = gridSize[1] * gridSize[2] * gridSize[3] -start_time = time.time() -metric = metricGet_AlcubierreComoving(gridSize, worldCenter, v, R, sigma, gridScale) -monitor_progress(start_time, total_steps, total_steps) - -# Compute the stress-energy tensor from the Alcubierre metric -energy_tensor = getEnergyTensor(metric, tryGPU=False, diffOrder='fourth') - -# Print the resulting stress-energy tensor -print(energy_tensor) \ No newline at end of file +from Metrics.Get_AlcubierreComoving import metricGet_AlcubierreComoving +from solver.getEnergyTensor import getEnergyTensor +from solver.metric_definition import metric +import time + +# Define monitoring function +def monitor_progress(start_time, step, total_steps): + elapsed_time = time.time() - start_time + estimated_total_time = elapsed_time / step * total_steps + estimated_remaining_time = estimated_total_time - elapsed_time + print(f"Step {step}/{total_steps} - Elapsed Time: {elapsed_time:.2f}s - Estimated Remaining Time: {estimated_remaining_time:.2f}s") + +# Define parameters for the Alcubierre metric +gridSize = [1, 100, 100, 100] +worldCenter = [0, 0, 0, 0] # You need to add the fourth element here +v = 0.5 +r_value = 5.0 +R_value = 10.0 +sigma_value = 2.0 +gridScale = [1, 0.1, 0.1, 0.1] + +# Define R and sigma +R = R_value +sigma = sigma_value + +# Check lengths of gridScale and worldCenter +print(len(gridScale)) +print(len(worldCenter)) + +# Get the Alcubierre metric +total_steps = gridSize[1] * gridSize[2] * gridSize[3] +start_time = time.time() +metric = metricGet_AlcubierreComoving(gridSize, worldCenter, v, R, sigma, gridScale) +monitor_progress(start_time, total_steps, total_steps) + +# Compute the stress-energy tensor from the Alcubierre metric +energy_tensor = getEnergyTensor(metric, tryGPU=False, diffOrder='fourth') + +# Print the resulting stress-energy tensor +print(energy_tensor) + diff --git a/solver/getEnergyTensor.py b/solver/getEnergyTensor.py index 2101b36..97d7cce 100644 --- a/solver/getEnergyTensor.py +++ b/solver/getEnergyTensor.py @@ -65,4 +65,5 @@ def getEnergyTensor(metric, tryGPU=False, diffOrder='fourth'): 'date': np.datetime64('today') } - return energy \ No newline at end of file + return energy + diff --git a/solver/tools/__init__.py b/solver/tools/__init__.py index 8b13789..139597f 100644 --- a/solver/tools/__init__.py +++ b/solver/tools/__init__.py @@ -1 +1,2 @@ + diff --git a/solver/tools/met2den2.py b/solver/tools/met2den2.py index c300235..ef58b78 100644 --- a/solver/tools/met2den2.py +++ b/solver/tools/met2den2.py @@ -1,3 +1,9 @@ +from solver.tools.c4Inv import c4Inv2 +from solver.tools.ricciT import ricciT2 +from solver.tools.ricciS import ricciS2 +from solver.tools.einT import einT2 +from solver.tools.einE import einE + def met2den2(metricTensor, delta=None, units=None): # Set default values if not provided if delta is None: @@ -25,3 +31,4 @@ def met2den2(metricTensor, delta=None, units=None): # Example usage: # energyDensity = met2den2(metricTensor) + diff --git a/solver/tools/ricciT.py b/solver/tools/ricciT.py index 4e5ed27..3ee5dc0 100644 --- a/solver/tools/ricciT.py +++ b/solver/tools/ricciT.py @@ -1,3 +1,7 @@ +import numpy as np +from solver.tools.takeFiniteDifference1_2 import takeFiniteDifference1_2 +from solver.tools.takeFiniteDifference2_2 import takeFiniteDifference2_2 + def ricciT2(gu, gl, delta): # Initialize Ricci tensor R_munu = [[None for _ in range(4)] for _ in range(4)] @@ -63,3 +67,4 @@ def ricciT2(gu, gl, delta): # Example usage: # R_munu = ricciT2(gu, gl, delta) + diff --git a/solver/tools/ricciTMem2.py b/solver/tools/ricciTMem2.py index 6c27c5c..714d59c 100644 --- a/solver/tools/ricciTMem2.py +++ b/solver/tools/ricciTMem2.py @@ -1,3 +1,7 @@ +import numpy as np +from solver.tools.takeFiniteDifference1_2 import takeFiniteDifference1_2 +from solver.tools.takeFiniteDifference2_2 import takeFiniteDifference2_2 + def ricciTMem2(gu, gl, delta): # Calculate the size of gl{1,1} s = np.shape(gl[0][0]) @@ -51,3 +55,4 @@ def ricciTMem2(gu, gl, delta): R_munu[3][2] = R_munu[2][3] return R_munu + diff --git a/solver/tools/takeFiniteDifference1_2.py b/solver/tools/takeFiniteDifference1_2.py index ca2e7cd..ab4c591 100644 --- a/solver/tools/takeFiniteDifference1_2.py +++ b/solver/tools/takeFiniteDifference1_2.py @@ -1,3 +1,6 @@ +import numpy as np +import torch + def takeFiniteDifference1_2(A, k, delta): s = np.shape(A) @@ -27,3 +30,4 @@ def takeFiniteDifference1_2(A, k, delta): B[:, :, :, -1] = B[:, :, :, -2] return B + diff --git a/solver/tools/takeFiniteDifference2_2.py b/solver/tools/takeFiniteDifference2_2.py index cb76dbb..62828c4 100644 --- a/solver/tools/takeFiniteDifference2_2.py +++ b/solver/tools/takeFiniteDifference2_2.py @@ -1,3 +1,6 @@ +import numpy as np +import torch + def takeFiniteDifference2_2(A, k1, k2, delta): s = np.shape(A) @@ -66,3 +69,4 @@ def takeFiniteDifference2_2(A, k1, k2, delta): A[:, :, x_1, y_1] - A[:, :, x_1, y1] - A[:, :, x1, y_1] + A[:, :, x1, y1]) return B + diff --git a/visualizer.py b/visualizer.py index 06b86a5..a3c179b 100644 --- a/visualizer.py +++ b/visualizer.py @@ -1,117 +1,118 @@ -import tkinter as tk -from tkinter import ttk -import sympy as sp -import numpy as np -import matplotlib.pyplot as plt -from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg -from analyticalEnergyTensor import analyticalEnergyTensor - -# Define a function to create and return the metric tensor -def create_metric_tensor(t_amp, x_amp, y_amp, z_amp, xy_coupling, t, x, y, z): - g = sp.zeros(4, 4) - g[0, 0] = -1 - g[1, 1] = 1 + x_amp * sp.sin(x) * sp.sin(t) - g[2, 2] = 1 + y_amp * sp.sin(y) * sp.sin(t) - g[3, 3] = 1 + z_amp * sp.sin(z) * sp.sin(t) - g[1, 2] = xy_coupling * sp.cos(x) * sp.cos(y) * sp.cos(t) - g[2, 1] = g[1, 2] # Symmetric component - return g - -# Function to evaluate tensor components numerically -def evaluate_component(component, t_val, x_grid, y_grid, z_val, t, x, y, z): - component_func = sp.lambdify((t, x, y, z), component, 'numpy') - return component_func(t_val, x_grid, y_grid, z_val) - -# Function to update the plot -def update_plot(): - # Get parameter values from UI - t_amp = float(t_amp_var.get()) - x_amp = float(x_amp_var.get()) - y_amp = float(y_amp_var.get()) - z_amp = float(z_amp_var.get()) - xy_coupling = float(xy_coupling_var.get()) - - # Define symbolic variables - t, x, y, z = sp.symbols('t x y z') - - # Create metric tensor - g = create_metric_tensor(t_amp, x_amp, y_amp, z_amp, xy_coupling, t, x, y, z) - - # Define the coordinates - coords = [t, x, y, z] - - # Compute the stress-energy tensor - T_out = analyticalEnergyTensor(g, coords) - - # Define a numerical grid for evaluation - x_vals = np.linspace(-5, 5, 100) - y_vals = np.linspace(-5, 5, 100) - x_grid, y_grid = np.meshgrid(x_vals, y_vals) - t_val = 0 # Fixed time value for visualization - z_val = 0 # Fixed z value for visualization - - # Visualize multiple components of the stress-energy tensor - components_to_plot = [(1, 1), (1, 2), (2, 2), (3, 3)] # List of tensor components to visualize - fig, axs = plt.subplots(2, 2, figsize=(12, 10)) - - for ax, (i, j) in zip(axs.flatten(), components_to_plot): - component = T_out[i, j] - component_values = evaluate_component(component, t_val, x_grid, y_grid, z_val, t, x, y, z) - - c = ax.contourf(x_grid, y_grid, component_values, cmap='viridis') - fig.colorbar(c, ax=ax) - ax.set_title(f'T[{i},{j}]') - ax.set_xlabel('x') - ax.set_ylabel('y') - - plt.tight_layout() - - # Clear the old plot - for widget in plot_frame.winfo_children(): - widget.destroy() - - # Embed the new plot in the tkinter window - canvas = FigureCanvasTkAgg(fig, master=plot_frame) - canvas.draw() - canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True) - -# Create the main window -root = tk.Tk() -root.title("Stress-Energy Tensor Visualization") - -# Create input controls -control_frame = ttk.Frame(root) -control_frame.pack(side=tk.LEFT, fill=tk.Y) - -tk.Label(control_frame, text="t amplitude").pack() -t_amp_var = tk.StringVar(value="0.1") -tk.Entry(control_frame, textvariable=t_amp_var).pack() - -tk.Label(control_frame, text="x amplitude").pack() -x_amp_var = tk.StringVar(value="0.1") -tk.Entry(control_frame, textvariable=x_amp_var).pack() - -tk.Label(control_frame, text="y amplitude").pack() -y_amp_var = tk.StringVar(value="0.1") -tk.Entry(control_frame, textvariable=y_amp_var).pack() - -tk.Label(control_frame, text="z amplitude").pack() -z_amp_var = tk.StringVar(value="0.1") -tk.Entry(control_frame, textvariable=z_amp_var).pack() - -tk.Label(control_frame, text="xy coupling").pack() -xy_coupling_var = tk.StringVar(value="0.05") -tk.Entry(control_frame, textvariable=xy_coupling_var).pack() - -update_button = ttk.Button(control_frame, text="Update Plot", command=update_plot) -update_button.pack(pady=10) - -# Create plot area -plot_frame = ttk.Frame(root) -plot_frame.pack(side=tk.RIGHT, fill=tk.BOTH, expand=True) - -# Initial plot -update_plot() - -# Run the tkinter main loop -root.mainloop() \ No newline at end of file +import tkinter as tk +from tkinter import ttk +import sympy as sp +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg +from analyticalEnergyTensor import analyticalEnergyTensor + +# Define a function to create and return the metric tensor +def create_metric_tensor(t_amp, x_amp, y_amp, z_amp, xy_coupling, t, x, y, z): + g = sp.zeros(4, 4) + g[0, 0] = -1 + g[1, 1] = 1 + x_amp * sp.sin(x) * sp.sin(t) + g[2, 2] = 1 + y_amp * sp.sin(y) * sp.sin(t) + g[3, 3] = 1 + z_amp * sp.sin(z) * sp.sin(t) + g[1, 2] = xy_coupling * sp.cos(x) * sp.cos(y) * sp.cos(t) + g[2, 1] = g[1, 2] # Symmetric component + return g + +# Function to evaluate tensor components numerically +def evaluate_component(component, t_val, x_grid, y_grid, z_val, t, x, y, z): + component_func = sp.lambdify((t, x, y, z), component, 'numpy') + return component_func(t_val, x_grid, y_grid, z_val) + +# Function to update the plot +def update_plot(): + # Get parameter values from UI + t_amp = float(t_amp_var.get()) + x_amp = float(x_amp_var.get()) + y_amp = float(y_amp_var.get()) + z_amp = float(z_amp_var.get()) + xy_coupling = float(xy_coupling_var.get()) + + # Define symbolic variables + t, x, y, z = sp.symbols('t x y z') + + # Create metric tensor + g = create_metric_tensor(t_amp, x_amp, y_amp, z_amp, xy_coupling, t, x, y, z) + + # Define the coordinates + coords = [t, x, y, z] + + # Compute the stress-energy tensor + T_out = analyticalEnergyTensor(g, coords) + + # Define a numerical grid for evaluation + x_vals = np.linspace(-5, 5, 100) + y_vals = np.linspace(-5, 5, 100) + x_grid, y_grid = np.meshgrid(x_vals, y_vals) + t_val = 0 # Fixed time value for visualization + z_val = 0 # Fixed z value for visualization + + # Visualize multiple components of the stress-energy tensor + components_to_plot = [(1, 1), (1, 2), (2, 2), (3, 3)] # List of tensor components to visualize + fig, axs = plt.subplots(2, 2, figsize=(12, 10)) + + for ax, (i, j) in zip(axs.flatten(), components_to_plot): + component = T_out[i, j] + component_values = evaluate_component(component, t_val, x_grid, y_grid, z_val, t, x, y, z) + + c = ax.contourf(x_grid, y_grid, component_values, cmap='viridis') + fig.colorbar(c, ax=ax) + ax.set_title(f'T[{i},{j}]') + ax.set_xlabel('x') + ax.set_ylabel('y') + + plt.tight_layout() + + # Clear the old plot + for widget in plot_frame.winfo_children(): + widget.destroy() + + # Embed the new plot in the tkinter window + canvas = FigureCanvasTkAgg(fig, master=plot_frame) + canvas.draw() + canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True) + +# Create the main window +root = tk.Tk() +root.title("Stress-Energy Tensor Visualization") + +# Create input controls +control_frame = ttk.Frame(root) +control_frame.pack(side=tk.LEFT, fill=tk.Y) + +tk.Label(control_frame, text="t amplitude").pack() +t_amp_var = tk.StringVar(value="0.1") +tk.Entry(control_frame, textvariable=t_amp_var).pack() + +tk.Label(control_frame, text="x amplitude").pack() +x_amp_var = tk.StringVar(value="0.1") +tk.Entry(control_frame, textvariable=x_amp_var).pack() + +tk.Label(control_frame, text="y amplitude").pack() +y_amp_var = tk.StringVar(value="0.1") +tk.Entry(control_frame, textvariable=y_amp_var).pack() + +tk.Label(control_frame, text="z amplitude").pack() +z_amp_var = tk.StringVar(value="0.1") +tk.Entry(control_frame, textvariable=z_amp_var).pack() + +tk.Label(control_frame, text="xy coupling").pack() +xy_coupling_var = tk.StringVar(value="0.05") +tk.Entry(control_frame, textvariable=xy_coupling_var).pack() + +update_button = ttk.Button(control_frame, text="Update Plot", command=update_plot) +update_button.pack(pady=10) + +# Create plot area +plot_frame = ttk.Frame(root) +plot_frame.pack(side=tk.RIGHT, fill=tk.BOTH, expand=True) + +# Initial plot +update_plot() + +# Run the tkinter main loop +root.mainloop() +