<jupyter_start><jupyter_text>IntroductionClassical electromagnetism is most often described using maxwell's equations. Instead, we can also describe it using a Lagrange density and an action which is the spacetime integral over the Lagrange density.The field is represented by a 4-vector in the [spacetime-algebra](https://en.wikipedia.org/wiki/Spacetime_algebra) where the first component is the electric potential and the last three components are the magnetic vector potential. Such as 4-vector is given at every point in spacetime. The Lagrangian density at a spacetime point $X = (t, x, y, z)$ for such a 4-vector field $A(X)$ when speed of light $c = 1$ and without external sources is given by the following equation:$\mathcal{L}(A, X) = \langle \nabla_X A(X) * \widetilde{\nabla_X A}(X) \rangle_0$The [principle of stationary action](https://en.wikipedia.org/wiki/Principle_of_least_action) then states that the classical solution of the field is achieved when the action$S(A) = \int_{X}{\mathcal{L}(A, X) dX}$does not change anymore, that is $\delta S(A) = 0$. GoalBelow we will obtain an entire space-time field configuration $A(X)$ given only some boundary conditions and a function for the action given $A$. We will then use Tensorflow's optimizer tofind a field configuration that makes the action stationary.  Create the spacetime algebraHere we initialize a `tfga.GeometricAlgebra` instance with bases $e_0=e_t, e_1=e_x, e_2=e_y, e_3=e_z$ and corresponding metric $[-1, 1, 1, 1]$. We will use thiswhen calculating the action later as we need the geometric product and reversion operations.<jupyter_code>ga = GeometricAlgebra([-1, 1, 1, 1])<jupyter_output><empty_output><jupyter_text>Calculate the actionNow we create a function which returns the action $S$ given a field configuration $A(X)$ on a discretized spacetime lattice of size $[N, N, N, N]$. We use the following boundary conditions for $A(X)$:$A_{t=-1} = 0, A_{t=N} = 0$$A_{x=-1} = 10 sin(4 * \pi / N * t) e_0, A_{x=N} = -5 e_0$$A_{y=-1} = 0, A_{y=N} = 0$$A_{z=-1} = 0, A_{z=N} = 0$As a reminder, $e_0$ is the electric potential part of the 4-vector, so we have a periodic sine electric potential that changes over time (two periods in total) and amplitude `10` at the lower x boundary and a constant negative electric potential of `-5` at the upper x boundary.<jupyter_code>def get_action(config_a_variable):
    # config_a_variable will be of shape [N, N, N, N, 4].
    # The last axis' values are the e0, e1, e2, e3 parts of the multivector.

    # Finite differences in each direction using padding.
    # Example with zero padding (ie. zeros on the boundary):
    # 1 2 3
    #   1 2 3 0 padded right
    # - 0 1 2 3 padded left
    # = 1 1 1-3 padded right - padded left
    # As spacing we use 1 so we don't need to divide by anything here.
    
    # Also use the boundary conditions in the padded values here.
    # This gets a bit verbose because of the pad syntax esepcially since we only want to pad the
    # first index of the last axis with non-zeros.

    # Create time-varying boundary condition. Start with sine of shape [N].
    # Then reshape to [N, 1, N, N, 1] which we can concatenate with the
    # original values.
    pad_values = 10.0 * tf.sin(2.0 * tf.range(grid_size[0], dtype=tf.float32) * 2.0 * tf.constant(np.pi, dtype=tf.float32) / grid_size[0])
    pad_values = tf.expand_dims(pad_values, axis=-1)
    pad_values = tf.expand_dims(pad_values, axis=-1)
    pad_values = tf.expand_dims(pad_values, axis=-1)
    pad_values = tf.expand_dims(pad_values, axis=-1)
    pad_values = tf.tile(pad_values, [1, 1, grid_size[0], grid_size[0], 1])

    config_left_pad_x = tf.concat([
        tf.concat([pad_values, config_a_variable[..., :1]], axis=1),
        tf.pad(config_a_variable[..., 1:], [[0, 0], [1, 0], [0, 0], [0, 0], [0, 0]]),
    ], axis=-1)

    config_right_pad_x = tf.concat([
        tf.pad(config_a_variable[..., :1], [[0, 0], [0, 1], [0, 0], [0, 0], [0, 0]], constant_values=-5),
        tf.pad(config_a_variable[..., 1:], [[0, 0], [0, 1], [0, 0], [0, 0], [0, 0]]),
    ], axis=-1)

    config_left_pad_y = tf.concat([
        tf.pad(config_a_variable[..., :1], [[0, 0], [0, 0], [1, 0], [0, 0], [0, 0]]),
        tf.pad(config_a_variable[..., 1:], [[0, 0], [0, 0], [1, 0], [0, 0], [0, 0]]),
    ], axis=-1)

    config_dt_a = (
        tf.pad(config_a_variable, [[0, 1], [0, 0], [0, 0], [0, 0], [0, 0]]) -
        tf.pad(config_a_variable, [[1, 0], [0, 0], [0, 0], [0, 0], [0, 0]])
    )
    config_dx_a = config_right_pad_x - config_left_pad_x
    config_dy_a = (
        tf.pad(config_a_variable, [[0, 0], [0, 0], [0, 1], [0, 0], [0, 0]]) -
        config_left_pad_y
    )
    config_dz_a = (
        tf.pad(config_a_variable, [[0, 0], [0, 0], [0, 0], [0, 1], [0, 0]]) -
        tf.pad(config_a_variable, [[0, 0], [0, 0], [0, 0], [1, 0], [0, 0]])
    )

    # Convert to multivectors so we can use GA ops we need in the Lagrangian:
    # the geometric product and reversion.
    config_dt_a = ga.from_tensor_with_kind(config_dt_a, "vector")
    config_dx_a = ga.from_tensor_with_kind(config_dx_a, "vector")
    config_dy_a = ga.from_tensor_with_kind(config_dy_a, "vector")
    config_dz_a = ga.from_tensor_with_kind(config_dz_a, "vector")

    # Sum all the derivatives according to the action / Lagrangian and return a single scalar value
    return (
        tf.reduce_sum(ga.geom_prod(config_dt_a, ga.reversion(config_dt_a))[..., 0]) +
        tf.reduce_sum(ga.geom_prod(config_dx_a, ga.reversion(config_dx_a))[..., 0]) +
        tf.reduce_sum(ga.geom_prod(config_dy_a, ga.reversion(config_dy_a))[..., 0]) +
        tf.reduce_sum(ga.geom_prod(config_dz_a, ga.reversion(config_dz_a))[..., 0])
    )<jupyter_output><empty_output><jupyter_text>Initialize the 4-vector field variable randomly<jupyter_code>grid_size = [16, 16, 16, 16]

config_a_variable = tf.Variable(tf.random.normal([*grid_size, 4], seed=0))<jupyter_output><empty_output><jupyter_text>Optimize the 4-vector field variable to make the action stationaryIn order to make the action stationary we use a loss function that is minimal when the action is stationary (ie. the gradient of the action with respect to the field configuration is 0).We use the mean-squared error to create such a loss function, although other functions such as the absolute value would work too.We use Tensorflow's Adam optimizer to find a field configuration which minimizes the loss.<jupyter_code>optimizer = tf.optimizers.Adam(0.01)

@tf.function
def train_step(config_a_variable):
    # Principle of stationary action:
    # Minimize the distance of gradient of the action to zero with respect to our field
    with tf.GradientTape() as tape_outer:
        tape_outer.watch(config_a_variable)
        with tf.GradientTape() as tape:
            tape.watch(config_a_variable)
            loss = get_action(config_a_variable)

        grads = tape.gradient(loss, [config_a_variable])
        grads_mse = tf.reduce_mean(tf.square(grads))
    grads2 = tape_outer.gradient(grads_mse, [config_a_variable])
    optimizer.apply_gradients(zip(grads2, [config_a_variable]))

for i in range(3000):
    train_step(config_a_variable)<jupyter_output><empty_output><jupyter_text>Extract and visualize the optimized electric fieldNow we can take the result, that is the $A$ at every spacetime point and visualize it. Obviously we can't visualize a 4 dimensional 4-vector field. However we can look atindividual 2D slices of the electric potential field, which is the first component of the 4-vector, where the other two coordinates take on a specific value.<jupyter_code># Plot electric potential slices. We are not plotting the boundaries here.

plt.figure(figsize=(7, 7))
plt.imshow(config_a_variable[..., 0, 0, 0])
plt.colorbar()
plt.title("Electric potential in TX plane Y=0, Z=0")
plt.xlabel("X")
plt.ylabel("T")
plt.show()

plt.figure(figsize=(7, 7))
plt.imshow(config_a_variable[..., 0, 5, :, :, 0])
plt.colorbar()
plt.title("Electric potential in YZ plane T=0, X=5")
plt.xlabel("Z")
plt.ylabel("Y")
plt.show()

plt.figure(figsize=(7, 7))
plt.imshow(config_a_variable[..., 2, :, :, 0, 0])
plt.colorbar()
plt.title("Electric potential in XY plane T=2, Z=0")
plt.xlabel("Y")
plt.ylabel("X")
plt.show()<jupyter_output><empty_output><jupyter_text>In the first figure we can see the potential close to X=0 (where we applied the sine boundary condition) changing over time.The second figure shows the YZ slice at T=0, X=5 where the potential is almost constant but we still have a radial symmetry.The last figure shows the XY slice at T=2, Z=0 where the potential takes its maximum value around X=0 if we look at the first figure. We can also see that on upper boundary of X that we have a negative potential as we applied a constant negative electric potential for boundary condition there.We can also visualize the XY slices over time in a video. For this I saved the XY slices at all times and converted them to a webm using ffmpeg. Here we can see the electric potential close to X=0 changing over time as we expected from the boundary condition. (Direct link: [em_output/electric_potential.webm](https://raw.githubusercontent.com/RobinKa/tfga/master/notebooks/em_output/electric_potential.webm))<jupyter_code>Video("./em_output/electric_potential.webm")<jupyter_output><empty_output><jupyter_text>Next we can look at the electric vector field corresponding to the electric potential: $E = -\nabla_{x,y,z} \langle A(X) \rangle_{e0} - \nabla_t \langle A(X) \rangle_{e1,e2,e3}$<jupyter_code>def draw_electric_field_xy(t, z):
    # Extract XY slice of electric potential [T=t, X, Y, Z=0, 0]
    electric_potential = config_a_variable[t, :, :, z, 0]
    magnetic_potential_t = config_a_variable[t, :, :, z, 1:]
    magnetic_potential_t2 = config_a_variable[t+1, :, :, z, 1:]

    # The electric field can be obtained from the 4-vector potential:
    # E = - (d/dx, d/dy, d/dz) <A>_e0 - d/dt <A>_e1,e2,e3
    # We can use finite differences again to approximate the derivatives.
    # We also need to get rid of the last element of the respective other axis,
    # since we couldn't calculate the last finite difference as that would
    # require using the boundary condition (which is possible, but would require extra code).

    # Start with -(d/dx, d/dy, d/dz) <A>_e0
    ex = -(electric_potential[1:, :-1] - electric_potential[:-1, :-1])
    ey = -(electric_potential[:-1, 1:] - electric_potential[:-1, :-1])
    
    # Calculate d/dt <A>_e1,e2,e3 and add it to the previous calculation
    dt_mag_a = -(magnetic_potential_t2[-1, :-1] - magnetic_potential_t[:-1, :-1])

    ex += dt_mag_a[..., 0]
    ey += dt_mag_a[..., 1]

    ys, xs = np.meshgrid(np.arange(ex.shape[0]), np.arange(ex.shape[1]))

    plt.figure(figsize=(10, 10))
    plt.quiver(ys, xs, ey, ex, scale=10, scale_units="inches")
    plt.xlabel("Y")
    plt.ylabel("X")
    plt.title("Electric field XY at T=%d, Z=%d" % (t, z))

draw_electric_field_xy(t=2, z=0)
plt.show()<jupyter_output><empty_output><jupyter_text>And again I made a video showing all the time slices. (Direct link: [em_output/electric_field.webm](https://raw.githubusercontent.com/RobinKa/tfga/master/notebooks/em_output/electric_field.webm))<jupyter_code>Video("./em_output/electric_field.webm")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.127874071453838, 'fraction_numerical': 0.026353024407499117, 'mean_word_length': 4.08955895589559, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '453714', 'license': 'MIT', 'n_tokens_mistral': 4103, 'n_tokens_neox': 3801, 'n_words': 1281, 'path': 'notebooks/em.ipynb', 'repo_name': 'RAufschlaeger/tfga', 'size': '11308'}