pytdgl3d 0.3.0

3D Time-Dependent Ginzburg-Landau solver

tdgl3d

tdgl3d is an open-source Python package for solving three-dimensional time-dependent Ginzburg-Landau (TDGL) problems in superconductors. It represents the superconducting domain by a closed triangular surface mesh, generates a tetrahedral volume mesh, assembles finite-volume operators on mesh vertices and edges, and evolves the complex superconducting order parameter using gauge-invariant link variables.

The package is designed for computational studies of vortex nucleation, vortex motion, magnetic screening, fluxoid quantization, and transport-like boundary conditions in finite three-dimensional geometries. Built-in geometry generators cover common shapes such as cones, cylinders, cuboids, cubes, pyramids, spheres, and ellipsoids, while user-defined closed surface meshes can also be supplied.

Contents

• Scientific Scope
• Main Features
• Numerical Method
• Package Architecture
• Installation
• Quick Start
• Geometry and Device Construction
• Solver Usage
• Magnetic Sources
• Screening and GPU Acceleration
• Terminal Currents and Probe Points
• Solution Analysis
• Visualization
• HDF5 Output Format
• Testing and Validation
• Reproducibility Checklist
• Known Limitations
• License
• Citation

Scientific Scope

tdgl3d solves a generalized TDGL model for a superconducting volume Omega. The primary fields are:

• psi(r, t): complex superconducting order parameter.
• mu(r, t): scalar electric potential.
• A(r, t): vector potential.
• J_s: superconducting current on mesh edges.
• J_n: normal current on mesh edges.

The code uses dimensionless TDGL units internally. Lengths are scaled by the coherence length xi, magnetic fields by Bc2, vector potentials by xi * Bc2, and currents by the corresponding GL current scale exposed through the Device physical-scale helpers.

The intended use cases include:

• 3D Abrikosov vortex nucleation and relaxation.
• Vortex line deformation in non-planar geometries.
• Surface and geometry effects in cones, tips, spheres, and finite-thickness mesoscopic superconductors.
• Fluxoid and winding-number analysis along arbitrary 3D contours.
• Transport-style simulations using terminal volumes and scalar-potential boundary feedback.
• Post-processing of HDF5 TDGL trajectories for visualization and physical observables.

Main Features

• 3D tetrahedral finite-volume mesh support.
• Gauge-invariant edge link variables for vector-potential coupling.
• Adaptive semi-implicit Euler time stepping for TDGL evolution.
• Built-in geometry generators: cone, cylinder, cuboid, box, cube, pyramid, sphere, ellipsoid.
• User-supplied closed triangular surface meshes through Volume.
• Material model through Layer, including xi, London penetration depth, conductivity, u, and gamma.
• Device abstraction for geometry, materials, mesh generation, terminal volumes, probe points, transformations, HDF5 I/O, and physical scales.
• Sparse finite-volume gradient, divergence, and Laplacian operators.
• Optional magnetic self-field screening using a 3D Biot-Savart summation.
• Numba-accelerated CPU screening kernel.
• Optional CuPy CUDA screening kernel.
• Uniform magnetic field and circular current-loop source helpers.
• Spatially and time-dependent parameters using Python callables.
• Seed-solution support for continuation studies.
• HDF5 output containing device, mesh, options, snapshots, and dynamics.
• Solution analysis for order parameter, phase, current density, magnetic moment, magnetic field, vorticity, superfluid velocity, London penetration depth, and fluxoid quantization.
• Matplotlib slice plots, 3D plotting helpers, snapshots, GIF/MP4 animation, live monitoring, vortex detection, and XDMF export.
• Test suite covering geometry, meshing, finite-volume operators, short solver runs, serialization, visualization, screening, and post-processing.

Numerical Method

Mesh and Control Volumes

The computational domain is a closed triangular surface mesh. tdgl3d fills the interior with points and uses SciPy Delaunay tetrahedralization to create a tetrahedral mesh. Tetrahedra with centroids outside the surface are discarded, and boundary vertices are projected back to the original surface.

For each mesh vertex, a barycentric control volume is approximated by distributing each tetrahedron volume equally among its four vertices. The mesh stores:

• sites: vertex coordinates in dimensionless units.
• elements: tetrahedral connectivity.
• boundary_indices: mesh sites on the boundary.
• volumes: per-site control volumes.
• edge_mesh: unique edges, edge centers, edge directions, edge lengths, approximate dual-face areas, and boundary-edge indices.

Finite-Volume Operators

The finite-volume operators are assembled on the edge mesh:

• build_gradient: edge gradient from site values to edge values.
• build_divergence: site divergence from edge values.
• build_laplacian: finite-volume Laplacian assembled from edge weights.
• MeshOperators: cached divergence, scalar Laplacian, covariant gradient, and covariant Laplacian.

The edge weights use edge lengths, site control volumes, and approximate dual-face areas. This gives a practical Delaunay-Voronoi-style discretization for 3D exploratory calculations.

Gauge-Invariant Link Variables

The vector potential is evaluated at edge midpoints. For each oriented edge i -> j, the code computes a link variable

U_ij = exp(-i A_ij . (r_j - r_i))

The covariant gradient uses the link variable:

grad_A(psi)_ij = (U_ij psi_j - psi_i) / L_ij

This preserves gauge-invariant phase coupling on the mesh edges.

Time Stepping

The TDGL update uses a semi-implicit Euler step. The nonlinear term in |psi|^2 is handled through a quadratic solve for the next order-parameter magnitude. Adaptive stepping reduces the time step on failed updates and uses a sliding-window estimate of recent |psi|^2 changes to adjust the next step.

At each successful step, the solver:

1. Updates time-dependent disorder or pair-breaking fields.
2. Advances psi.
3. Solves for scalar potential mu.
4. Computes supercurrent and normal current.
5. Applies terminal-current feedback when terminals are configured.
6. Updates the induced vector potential if screening is enabled.
7. Periodically applies a Coulomb-gauge projection to the induced field.
8. Saves snapshots according to save_every.

Screening

When include_screening=True, the induced vector potential is updated from the site current density using a Biot-Savart-type volume summation:

A_induced(r_e) ~ sum_j J(r_j) V_j / |r_e - r_j|

The screening update uses a Polyak heavy-ball iteration controlled in the current implementation by screening_step_size and screening_drag. screening_tolerance is stored in the solver options as a screening convergence parameter. The CPU implementation uses Numba. The optional GPU path uses a CuPy raw CUDA kernel.

Package Architecture

tdgl3d/
  __init__.py              Public API exports
  about.py                 Version and dependency table helpers
  em.py                    Unit registry, vector potentials, Biot-Savart tools
  fluxoid.py               3D fluxoid contours and fluxoid integration
  geometry.py              Surface mesh generators
  parameter.py             Callable spatial/time-dependent parameters
  testing.py               Test runner helper

  device/
    layer.py               Material parameters
    volume.py              Closed triangular surface mesh container
    meshing.py             Interior fill and tetrahedral mesh generation
    device.py              Device abstraction and HDF5 I/O

  finite_volume/
    mesh.py                Tetrahedral mesh and mesh quality metrics
    edge_mesh.py           Edge topology and edge geometry
    operators.py           Gradient, divergence, Laplacian, link variables
    util.py                Mesh topology and geometry utilities

  solver/
    options.py             SolverOptions dataclass
    solve.py               Public solve() function and terminal validation
    solver.py              TDGLSolver implementation
    runner.py              Simulation loop and HDF5 writer
    screening.py           Numba/CuPy induced-vector-potential kernels

  solution/
    data.py                TDGLData, DynamicsData, HDF5 data helpers
    solution.py            Solution object and analysis methods
    plot_solution.py       Solution-bound Matplotlib plotting helpers

  sources/
    constant.py            Uniform-field source
    loop.py                Circular current-loop vector potential
    scaling.py             LinearRamp and Scale parameters

  visualization/
    common.py              Quantity enum, plot defaults, utility helpers
    io.py                  HDF5 plotting data extraction
    snapshot.py            Snapshot generation
    animate.py             Animation generation
    interactive.py         Interactive plotting helpers
    convert.py             XDMF export
    monitor.py             Live monitoring
    plots3d.py             3D plotting helpers
    vortex.py              Vortex detection and vortex-tube rendering
    analysis.py            Supercurrent and multi-slice analysis panels
    physics.py             Voltage, V-I, free-energy, and Lorentz plots

Installation

Requirements

The package metadata targets:

Python >= 3.14

Core Python dependencies:

• NumPy
• SciPy
• h5py
• Numba
• Pint
• cloudpickle
• Matplotlib
• Pillow

Optional dependencies:

• CuPy for CUDA GPU screening.
• PyVista for 3D visualization.
• meshio for XDMF export.
• pytest for development and testing.

Install From a Local Checkout

pip install .

Editable Development Install

pip install -e ".[dev,visualization]"

Install From requirements.txt

requirements.txt includes the optional packages listed in this repository:

pip install -r requirements.txt

Optional GPU Package

Install a CuPy package matching your CUDA version. For CUDA 12.x:

pip install cupy-cuda12x

For CUDA 11.x:

pip install cupy-cuda11x

Quick Start

This example builds a cone-shaped superconductor, generates a tetrahedral mesh, solves a short TDGL run with a uniform dimensionless Bz field, and plots the final order-parameter magnitude on a cross-section.

import matplotlib.pyplot as plt
import tdgl3d

layer = tdgl3d.Layer(
    coherence_length=1.0,
    london_lambda=2.0,
)

vertices, facets = tdgl3d.cone(
    radius_bottom=2.0,
    radius_top=0.3,
    height=3.0,
    n_radial=24,
    n_height=8,
)

volume = tdgl3d.Volume("cone", vertices, facets)

device = tdgl3d.Device(
    name="cone_device",
    layer=layer,
    volume=volume,
    length_units="um",
)

device.make_mesh(max_edge_length=0.9, min_points=150)
print(device.mesh_stats())

options = tdgl3d.SolverOptions(
    solve_time=2.0,
    dt_init=1e-3,
    dt_max=5e-2,
    adaptive=True,
    save_every=50,
    output_file="cone_solution.h5",
)

solution = tdgl3d.solve(
    device=device,
    options=options,
    applied_vector_potential=0.3,
    seed_noise=0.01,
)

fig, ax = solution.plot_order_parameter(slice_axis="z")
plt.show()

Geometry and Device Construction

Built-In Geometry Generators

Each geometry function returns:

vertices, facets = tdgl3d.geometry_function(...)

where vertices has shape (N, 3) and facets has shape (M, 3).

Function	Purpose
cone(radius_bottom, radius_top, height)	Cone or truncated cone aligned with z
cylinder(radius, height)	Cylinder, implemented as a cone special case
cuboid(lx, ly, lz)	Rectangular box
box(lx, ly, lz)	Alias-style box generator
cube(side)	Cube
pyramid(...)	Pyramid with polygonal or rectangular base
sphere(radius)	Spherical surface mesh
ellipsoid(a, b, c)	Axis-aligned ellipsoid

Example:

vertices, facets = tdgl3d.sphere(radius=2.0, n_phi=30, n_theta=60)
tdgl3d.close_surface(vertices, facets)
volume = tdgl3d.Volume("sphere", vertices, facets)

Custom Closed Surface Meshes

You may construct a Volume from your own closed triangular surface mesh:

volume = tdgl3d.Volume(
    name="custom",
    vertices=vertices_array,
    facets=triangles_array,
)

The surface must be watertight. A closed triangular mesh should have each surface edge shared by exactly two triangles. Use close_surface() to check simple watertightness before meshing.

Device Transformations

Volume and Device support basic transformations:

device_shifted = device.translate(dx=1.0, dy=0.0, dz=0.5)
device_scaled = device.scale(sx=1.0, sy=1.0, sz=2.0)
device_rotated = device.rotate(angle=45.0, axis="z")

Transformed devices do not keep the old mesh; regenerate the mesh after changing geometry.

Mesh Generation

device.make_mesh(max_edge_length=1.0, min_points=300)

Parameters:

• max_edge_length: target mesh spacing in the device length units.
• min_points: minimum interior point target used during point filling.

Inspect the generated mesh:

stats = device.mesh_stats()
quality = device.mesh.quality_metrics()
print(stats)
print(quality["summary"])

Plot the surface projection or filled surface projection:

device.plot(mesh=True)
device.draw(alpha=0.5)

Physical Scales

Device exposes useful dimensional scales through Pint:

print(device.coherence_length)
print(device.Bc2)
print(device.A0)
print(device.K0)
print(device.tau0(conductivity=1e7))
print(device.V0(conductivity=1e7))

Solver Usage

SolverOptions

SolverOptions controls integration, output, screening, sparse solver choice, GPU use, and monitoring.

Option	Default	Meaning
solve_time	100.0	Simulation time after optional skip period
skip_time	0.0	Initial time not saved to snapshots
dt_init	1e-4	Initial time step
dt_max	1e-2	Maximum adaptive time step
adaptive	True	Enable adaptive time stepping
adaptive_window	10	Window for time-step adaptation
adaptive_time_step_multiplier	0.25	Retry multiplier after failed step
max_solve_retries	10	Maximum retries per step
output_file	None	HDF5 output path; temporary file if unset
field_units	"mT"	Metadata field units
current_units	"uA"	Metadata current units
save_every	100	Save every N solver steps
include_screening	False	Enable magnetic self-field iteration
screening_step_size	0.1	Polyak screening step size
screening_tolerance	1e-3	Stored screening convergence parameter
screening_drag	0.5	Polyak heavy-ball drag
terminal_psi	0.0	`
sparse_solver	"superlu"	Sparse backend name
gpu	False	Use CuPy screening kernel
monitor	False	Enable SWMR-compatible live monitoring
monitor_update_interval	2.0	Monitor refresh interval
pause_on_interrupt	True	Save current state on keyboard interrupt

Example:

options = tdgl3d.SolverOptions(
    solve_time=10.0,
    skip_time=1.0,
    dt_init=1e-3,
    dt_max=5e-2,
    adaptive=True,
    save_every=100,
    output_file="run.h5",
)

Applied Vector Potential

applied_vector_potential may be:

• None or 0: no applied vector potential.
• float: interpreted as a uniform Bz field.
• numpy.ndarray with shape (num_edges, 3): value at edge midpoints.
• Callable or Parameter: evaluated as A(x, y, z).

Uniform field:

solution = tdgl3d.solve(device, options, applied_vector_potential=0.3)

Callable field:

def A_custom(x, y, z):
    import numpy as np
    A = np.zeros((len(x), 3))
    A[:, 0] = -0.5 * y
    A[:, 1] = 0.5 * x
    return A

solution = tdgl3d.solve(device, options, applied_vector_potential=A_custom)

Disorder or Pair-Breaking Parameter

disorder_epsilon can be a scalar, array, spatial callable, or time-dependent callable. Time-dependent functions must accept keyword-only t.

def epsilon(x, y, z, *, t):
    import numpy as np
    value = max(1.0 - 0.1 * t, 0.5)
    return value * np.ones_like(x)

solution = tdgl3d.solve(
    device,
    options,
    disorder_epsilon=epsilon,
)

Initial Conditions

Initial-condition controls:

• seed_noise: small complex noise added to the initial state.
• field_cooling=True: initialize a moderate-amplitude random-phase state to help vortex nucleation.
• seed_solution: initialize from a previous Solution; different meshes are handled by interpolation.

solution2 = tdgl3d.solve(
    device,
    options,
    applied_vector_potential=0.4,
    seed_solution=solution1,
)

Manual Solver Access

For advanced workflows:

from tdgl3d.solver import TDGLSolver

solver = TDGLSolver(
    device=device,
    options=options,
    applied_vector_potential=0.3,
)

for _ in range(100):
    dt, converged = solver.update()

Most users should prefer tdgl3d.solve().

Magnetic Sources

Uniform Field Source

A = tdgl3d.sources.ConstantField(Bz=0.5)
solution = tdgl3d.solve(device, options, applied_vector_potential=A)

The source uses the symmetric gauge:

A = (Bz / 2) (-y, x, 0)

Circular Current Loop

loop = tdgl3d.sources.CurrentLoop(
    current=1.0,
    radius=5.0,
    center=(0.0, 0.0, 0.0),
)

solution = tdgl3d.solve(device, options, applied_vector_potential=loop)

Time-Dependent Scaling

ramp = tdgl3d.sources.LinearRamp(
    tmin=0.0,
    tmax=10.0,
    vmin=0.0,
    vmax=1.0,
)

Custom Parameter objects can be combined arithmetically:

param = 0.5 * tdgl3d.sources.Scale(2.0)

Screening and GPU Acceleration

Enable magnetic self-field screening:

options = tdgl3d.SolverOptions(
    solve_time=5.0,
    include_screening=True,
    screening_step_size=0.1,
    screening_drag=0.5,
    gpu=False,
)

Use the GPU kernel when CuPy is installed:

options.gpu = True

Notes:

• The screening calculation is typically the most expensive part of a run.
• The CPU path is Numba-accelerated and may compile on first use.
• The GPU path requires a compatible CUDA runtime and a matching CuPy package.
• Very large meshes can exceed GPU memory because the kernel evaluates many source-site and target-edge interactions.

Terminal Currents and Probe Points

Terminal Volumes

Terminals are represented as Volume objects associated with a Device. The solver identifies mesh sites belonging to those terminal regions and uses them for scalar-potential boundary feedback and optional psi suppression.

terminal_left = tdgl3d.Volume("left", left_vertices, left_facets)
terminal_right = tdgl3d.Volume("right", right_vertices, right_facets)

device = tdgl3d.Device(
    "wire",
    layer=layer,
    volume=body,
    terminals=[terminal_left, terminal_right],
)

Inspect terminal metadata:

print(device.terminal_info())

Terminal Driving

The solver accepts terminal_currents in tdgl3d.solve() as either a dictionary or a callable returning a dictionary. The implemented terminal update distinguishes two value forms:

• Scalar value: direct scalar-potential bias at that terminal.
• {"current": I}: current-bias feedback that adjusts terminal mu to drive the measured normal current toward I.

Voltage-bias style input:

solution = tdgl3d.solve(
    device,
    options,
    terminal_currents={"left": 0.5, "right": -0.5},
)

Current-bias feedback:

solution = tdgl3d.solve(
    device,
    options,
    terminal_currents={
        "left": {"current": 1.0},
        "right": {"current": -1.0},
    },
)

Time-dependent driving:

def drive(t):
    return {
        "left": {"current": 1.0},
        "right": {"current": -1.0},
    }

solution = tdgl3d.solve(
    device,
    options,
    terminal_currents=drive,
)

Current Conservation Check

Static currents:

tdgl3d.validate_terminal_currents(
    {"left": 1.0, "right": -1.0},
    terminal_names=["left", "right"],
)

Time-dependent currents:

def currents(t):
    return {"left": 1.0, "right": -1.0}

tdgl3d.validate_terminal_currents(
    currents,
    terminal_names=["left", "right"],
    options=options,
)

validate_terminal_currents() currently validates flat scalar dictionaries or callables returning flat scalar dictionaries. It checks that terminal names are known and that the scalar values sum to zero. For nested current-bias dictionaries such as {"left": {"current": 1.0}}, check conservation manually before passing them to solve().

The sum over all terminal currents must be zero.

Probe Points

Probe points are physical coordinates used to record scalar potential during the run:

device = tdgl3d.Device(
    "device_with_probes",
    layer=layer,
    volume=volume,
    probe_points=[
        [-1.0, 0.0, 0.5],
        [1.0, 0.0, 0.5],
    ],
)

After solving:

voltage = solution.dynamics.voltage(i=0, j=1)
mean_v = solution.dynamics.mean_voltage(i=0, j=1, tmin=1.0)

Solution Analysis

Load a solution:

solution = tdgl3d.Solution.from_hdf5("run.h5")

Common properties:

psi = solution.order_parameter
abs_psi = solution.order_parameter_magnitude
phase = solution.phase
mu = solution.scalar_potential
Js = solution.supercurrent
Jn = solution.normal_current
times = solution.times

Select data by saved frame:

data_range = solution.data_range
data0 = solution.get_data_at_step(data_range[0])
last = solution.get_data_at_step(data_range[1])

Find the closest saved frame to a target simulation time:

frame = solution.closest_solve_step(target_time=5.0)

Derived Quantities

magnetic_moment = solution.magnetic_moment()
magnetic_field = solution.magnetic_field()
free_energy = solution.free_energy()
superfluid_velocity = solution.superfluid_velocity()
lambda_eff = solution.london_penetration_depth()
vorticity = solution.vorticity

Interpolation

Order parameter:

positions = solution.mesh.sites[:10]
psi_at_positions = solution.interp_order_parameter(positions)

Current density:

J = solution.interp_current_density(
    positions,
    dataset=None,       # None = supercurrent + normal current
    method="nearest",
)

Regular 3D grid:

xg, yg, zg, Jgrid = solution.grid_current_density(
    grid_shape=(40, 40, 40),
    method="nearest",
)

Fluxoid Quantization

Circular contour:

contour = tdgl3d.make_fluxoid_contour_circle(
    center=(0.0, 0.0, 1.0),
    radius=0.5,
    normal=(0.0, 0.0, 1.0),
    n_points=64,
)

fluxoid = solution.fluxoid(contour, kappa=solution.device.kappa)
print(fluxoid.flux_part, fluxoid.supercurrent_part)

Validation helper:

result = solution.validate_flux_quantization(
    contour,
    kappa=solution.device.kappa,
    tolerance=0.1,
)
print(result)

Rectangular contour:

contour = tdgl3d.make_fluxoid_contour_rectangular(
    center=(0.0, 0.0, 1.0),
    width=1.0,
    height=1.0,
    normal=(0.0, 0.0, 1.0),
    n_points_per_side=16,
)

Current Through a Surface

For a sampled planar surface:

current = solution.current_through_surface(
    surface_points=points_on_surface,
    surface_normal=[1.0, 0.0, 0.0],
)

For time series through named paths:

from tdgl3d.solution.data import get_current_through_paths

paths = {"center": points_on_surface}
currents = get_current_through_paths("run.h5", paths)

Visualization

Solution Plot Shortcuts

Each plotting method returns (fig, ax):

solution.plot_order_parameter(slice_axis="z")
solution.plot_phase(slice_axis="z")
solution.plot_scalar_potential(slice_axis="z")
solution.plot_supercurrent(slice_axis="z")
solution.plot_currents(slice_axis="z")
solution.plot_vorticity(slice_axis="z")
solution.plot_vector_potential(slice_axis="z")

Slice options:

• slice_axis: "x", "y", or "z".
• slice_value: coordinate of the slice. If omitted, the midpoint is used.

Snapshot Generation

from tdgl3d.visualization import generate_snapshots

figures = generate_snapshots(
    input_path="run.h5",
    steps=[0, 5, 10],
    quantities=["order_parameter", "phase", "supercurrent"],
    slice_axis="z",
)

Animation

from tdgl3d.visualization import create_animation

create_animation(
    input_file="run.h5",
    output_file="order_parameter.gif",
    quantity="order_parameter",
    fps=10,
    ngrid=80,
)

For MP4 output, Matplotlib needs a working ffmpeg writer.

Command-Line Interface

The package defines the console script:

tdgl3d-visualize

Examples:

tdgl3d-visualize -i run.h5 snapshot --steps 0 5 10
tdgl3d-visualize -i run.h5 interactive
tdgl3d-visualize -i run.h5 -o solution.xdmf convert
tdgl3d-visualize -i run.h5 monitor

3D and Vortex Utilities

The visualization package also includes:

• plot_3d_order_parameter
• plot_3d_scalar_potential
• plot_3d_vortex_cores
• plot_3d_volume
• plot_3d_isosurface
• detect_vortex_positions
• build_core_tubes
• make_vortex_slice_gif
• make_vortex_3d_gif
• plot_supercurrent_analysis
• plot_multi_z_slices
• plot_voltage_vs_time
• plot_vi_curve
• plot_free_energy_density
• plot_total_free_energy_vs_time

HDF5 Output Format

The solver writes an HDF5 file containing the full state required for post-processing.

/
  attrs:
    solver
    version
    field_units
    current_units
    total_steps
    total_saves

  device/
    attrs:
      name
      length_units
    layer/
      attrs:
        coherence_length
        london_lambda
        thickness
        conductivity
        u
        gamma
    volume/
      attrs:
        name
      vertices
      facets
    holes/
    terminals/
    probe_points

  mesh/
    sites
    elements
    boundary_indices
    volumes
    edge_mesh/
      edges
      edge_centers
      directions
      edge_lengths
      dual_face_areas
      boundary_edge_indices

  options/
    attrs:
      solve_time
      skip_time
      dt_init
      dt_max
      adaptive
      save_every
      include_screening
      gpu

  data/
    0/
      attrs:
        dt
        time
        step
      psi
      mu
      supercurrent
      normal_current
      A_applied
      A_induced
    1/
    ...

  dynamics/
    time
    dt
    step
    free_energy
    mu

psi is stored as two columns: real and imaginary parts. The Solution and TDGLData loaders reconstruct it as a complex array.

Because the device and mesh are embedded in the solution file, most analysis and visualization functions do not require the original Python script that created the simulation.

Testing and Validation

Run the test suite:

pytest tdgl3d/test

Or use the package helper:

import tdgl3d
tdgl3d.testing.run()

The tests exercise:

• Geometry generation and closed-volume containment checks.
• Tetrahedral mesh generation and finite-volume topology.
• Mesh quality metrics.
• Gradient, divergence, Laplacian, and link-variable construction.
• Short TDGL solves on cone and box geometries.
• Seed-solution continuation and interpolation across meshes.
• Time-dependent pair-breaking/disorder functions.
• Fluxoid contour generation and fluxoid computation.
• Numba screening kernels and optional CPU/GPU consistency.
• HDF5 device and solution round trips.
• Solver option storage.
• Terminal current validation.
• Solution analysis methods.
• Snapshot, animation, monitor, and CLI import paths.
• Plotting helpers and visualization utilities.

For a quick smoke test, run:

pytest tdgl3d/test/test_fast.py

For simulation-level tests, run:

pytest tdgl3d/test/test_cone.py
pytest tdgl3d/test/test_features.py

Reproducibility Checklist

For reproducible simulation studies, record:

• Package version and dependency versions:

import tdgl3d
print(tdgl3d.version_table(verbose=True))

• Python version and operating system.
• Geometry generator parameters or the exact custom surface mesh.
• Material parameters in Layer.
• Device length units.
• Mesh generation parameters: max_edge_length, min_points.
• Solver options.
• Applied vector potential definition.
• Disorder or pair-breaking function.
• Random-noise amplitude and random seed behavior where relevant.
• Whether screening and GPU acceleration were enabled.
• Output HDF5 file.

The HDF5 solution stores most simulation metadata, including the embedded device, mesh, options, saved TDGL snapshots, and dynamics.

Troubleshooting

Mesh Generation Is Slow

Reduce min_points, increase max_edge_length, or start with a simpler surface mesh. Very fine surface meshes and small edge-length targets can create large Delaunay problems.

Solver Diverges or Reports NaN Values

Try:

• Decrease dt_init.
• Decrease dt_max.
• Increase mesh resolution.
• Reduce abrupt changes in time-dependent fields.
• Use seed_noise=0.0 for stable Meissner-like initialization.
• Use a seed solution for continuation in field/current sweeps.

No Data Appears in the Output File

Check save_every, solve_time, and skip_time. Snapshots are saved when simulation time has passed skip_time and the solver step is divisible by save_every.

GPU Screening Fails

Check that:

• CuPy is installed.
• The CuPy package matches the CUDA runtime.
• A CUDA-capable GPU is visible.
• The mesh is small enough for available GPU memory.

Set gpu=False to use the Numba CPU path.

Interactive or 3D Visualization Fails

Install optional visualization dependencies:

pip install pyvista meshio

For headless servers, use Matplotlib's non-GUI backend:

with tdgl3d.non_gui_backend():
    figures = tdgl3d.generate_snapshots("run.h5")

Known Limitations

• The mesher is intended for accessible Python workflows and exploratory simulations. It is not a replacement for specialized constrained tetrahedral mesh generators for very complex CAD geometries.
• Holes are represented in the Device data model, but robust hole-aware constrained meshing requires careful geometry preparation.
• Dual-face areas are approximate and based on local tetrahedral geometry.
• Screening uses a direct Biot-Savart summation, which can be expensive for large meshes.
• Current-through-surface helpers use nearest-edge interpolation and should be interpreted as practical post-processing estimates.
• The solver operates in dimensionless TDGL units internally; dimensional interpretation requires consistent material parameters and length units.
• Sparse solver options currently validate several backend names, but the implemented solve path primarily uses SciPy sparse direct solvers.

License

MIT License.

Citation

If you use our work, please cite:

Tanvir Hassan and A. Hasnat, tdgl3d: Three-Dimensional Time-Dependent Ginzburg-Landau Solver in Python (manuscript submitted).