Generating the pyTDGL logo

In this notebook, we generate the pyTDGL logo (pictured above) by simulating a mesoscopic superconducting film patterned into the letters “py” and subject to a uniform out-of-plane applied field.

# Automatically install tdgl from GitHub only if running in Google Colab
if "google.colab" in str(get_ipython()):
    %pip install --quiet git+https://github.com/loganbvh/py-tdgl.git
    !apt install ffmpeg

%config InlineBackend.figure_formats = {"retina", "png"}
%matplotlib inline

import os
import tempfile
from typing import List, Tuple

os.environ["OPENBLAS_NUM_THREADS"] = "1"

import numpy as np
import matplotlib as mpl
from matplotlib.path import Path
from matplotlib.textpath import TextPath
from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = (5, 4)
plt.rcParams["font.size"] = 11

import tdgl
from tdgl.geometry import box, ensure_unique, close_curve
from tdgl.sources import LinearRamp, ConstantField

USE_TEMPDIR = True
SAVE = True
MAKE_ANIMATIONS = True
image_path = os.path.join(os.pardir, "images")

if USE_TEMPDIR:
    tempdir = tempfile.TemporaryDirectory()
    output_path = tempdir.name
else:
    output_path = "."

Define the letter geometries

# Letters that have loops will require more than one Polygon to define.
loops_uppercase = "ABDOPQR"
loops_lowercase = "abdegopq"

def interp_path(path: Path, points_per_segment: int = 11) -> np.ndarray:
    """Generate polygon coordinates for a given matplotlib Path"""
    # The paths are defined in terms of Bezier curves, so we need to
    # evaluate the curves at many points.
    t = np.linspace(0, 1, points_per_segment)
    segments = [
        bezier(t) for bezier, code in path.iter_bezier() if code != Path.MOVETO
    ]
    points = np.concatenate(segments)
    return ensure_unique(points)

def make_polygons(
    letter: str,
    fontsize: int,
    fontprops: FontProperties,
    resample_points: Tuple[int, int] = (251, 101),
) -> List[tdgl.Polygon]:
    path = TextPath((0, 0), letter, size=fontsize, prop=fontprops)
    if letter in loops_lowercase + loops_uppercase:
        # If the letter has holes in it, split the overall outline and
        # all holes into separate paths.
        jumps = np.where(path.codes == TextPath.MOVETO)[0][1:]
        vertices = np.split(path.vertices, jumps)
        codes = np.split(path.codes, jumps)
        paths = [Path(v[:-1], c[:-1]) for v, c in zip(vertices, codes)]
    else:
        # If the letter doesn't have holes in it, we only need one path.
        paths = [Path(path.vertices[:-1], path.codes[:-1])]
    polygons = [tdgl.Polygon(points=interp_path(p)) for p in paths]
    # Sort the polygons by area in descending order so that the outline
    # will be first.
    polygons = sorted(polygons, key=lambda p: p.area, reverse=True)
    polygons = (
        [polygons[0].resample(resample_points[0]).buffer(2e-3)]
        + [p.resample(resample_points[1]) for p in polygons[1:]]
    )
    # Name the polygons
    for i, p in enumerate(polygons):
        if i == 0:
            p.name = letter
        else:
            p.name = f"{letter}_hole{i}"
    return polygons

Make the Device

fontsize = 10
fontprops = FontProperties(weight="bold")
xi = 0.4
layer = tdgl.Layer(london_lambda=4, coherence_length=xi, thickness=0.1)

p_outer, p_inner  = make_polygons("p", fontsize, fontprops)
y_outer,  = make_polygons("y", fontsize, fontprops)
# Merge the "p" and the "y" into a single "Polygon"
film = p_outer.union(y_outer.translate(dx=5.75)).resample(501).buffer(0)

device = tdgl.Device("py", layer=layer, film=film, holes=[p_inner])

fig, ax = device.draw(legend=False)

Generate the mesh

fig, ax = plt.subplots()
device.make_mesh(max_edge_length=xi / 2, smooth=100)
_ = device.plot(ax=ax, mesh=True, legend=False)

Constructing Voronoi polygons: 100%|█████| 4347/4347 [00:00<00:00, 22715.27it/s]

Solve the model

options = tdgl.SolverOptions(
    solve_time=800,
    field_units="mT",
    output_file=os.path.join(output_path, "py.h5"),
)

# Ramp the applied field from 0 to 1.8 mT between t=0 and t=100, then hold it at 1.8 mT.
applied_vector_potential = (
    LinearRamp(tmin=0, tmax=100)
    * ConstantField(1.8, field_units=options.field_units, length_units=device.length_units)
)

solution = tdgl.solve(
    device,
    options,
    applied_vector_potential=applied_vector_potential,
)

Simulating: 100%|██████████████████████████▉| 800/800 [02:00<00:00,  6.63tau/s ]

Plot the order parameter

fig, ax = solution.plot_order_parameter(figsize=(7, 2))

fig, ax = plt.subplots(figsize=(3, 1), facecolor="none")
ax.set_aspect("equal")
ax.axis("off")

device = solution.device
xs, ys = device.points.T
triangles = device.triangles
psi = solution.tdgl_data.psi
ax.tripcolor(xs, ys, triangles, np.abs(psi), vmin=0, cmap="viridis", shading="gouraud")

if SAVE:
    fig.set_size_inches(3, 1)
    fig.savefig(os.path.join(image_path, "logo-transparent-small.png"), dpi=300)
    fig.set_size_inches(4, 3)
    fig.savefig(os.path.join(image_path, "logo-transparent-large.png"), dpi=300, bbox_inches="tight")

Plot the supercurrent density

fig, ax = solution.plot_currents(
    dataset="supercurrent",
    figsize=(6, 3),
    vmin=0,
    vmax=3,
    streamplot=False,
)

Create an animation

from IPython.display import display, HTML
import h5py
from tdgl.visualization.animate import create_animation

def make_video_from_solution(solution, **kwargs):
    """Generates an HTML5 video from a tdgl.Solution."""
    with tdgl.non_gui_backend():
        with h5py.File(solution.path, "r") as h5file:
            anim = create_animation(h5file, **kwargs)
            video = anim.to_html5_video()
        return HTML(video)

if MAKE_ANIMATIONS:
    video = make_video_from_solution(
        solution,
        quantities=["order_parameter", "phase"],
        dpi=100,
        fps=50,
        title_off=True,
        axes_off=True,
        max_frame=600,
        max_cols=1,
        figure_kwargs=dict(figsize=(4.5, 6)),
    )
    display(video)

if USE_TEMPDIR:
    tempdir.cleanup()

tdgl.version_table()

Software	Version
tdgl	0.5.1; git revision 03859ea [2023-09-21]
Numpy	1.24.3
SciPy	1.10.1
matplotlib	3.7.1
jax	None
numba	0.57.1
IPython	8.14.0
Python	3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]
OS	posix [darwin]
Number of CPUs	Physical: 10, Logical: 10
BLAS Info	OPENBLAS
Thu Sep 21 15:35:01 2023 PDT