Metadata-Version: 2.4
Name: pytdgl3d
Version: 0.3.0
Summary: 3D Time-Dependent Ginzburg-Landau solver
Author-email: Tanvir Hassan <tanvir6307@gmail.com>
License-Expression: MIT
Requires-Python: >=3.14
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=2.4.2
Requires-Dist: scipy>=1.17.0
Requires-Dist: h5py>=3.15.1
Requires-Dist: numba>=0.64.0
Requires-Dist: pint>=0.25.2
Requires-Dist: cloudpickle>=3.1.2
Requires-Dist: matplotlib>=3.10.8
Requires-Dist: Pillow>=12.1.1
Provides-Extra: visualization
Requires-Dist: pyvista>=0.35; extra == "visualization"
Requires-Dist: meshio>=5.0; extra == "visualization"
Provides-Extra: gpu
Requires-Dist: cupy-cuda12x>=14.0.1; extra == "gpu"
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Dynamic: license-file

# 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](#scientific-scope)
- [Main Features](#main-features)
- [Numerical Method](#numerical-method)
- [Package Architecture](#package-architecture)
- [Installation](#installation)
- [Quick Start](#quick-start)
- [Geometry and Device Construction](#geometry-and-device-construction)
- [Solver Usage](#solver-usage)
- [Magnetic Sources](#magnetic-sources)
- [Screening and GPU Acceleration](#screening-and-gpu-acceleration)
- [Terminal Currents and Probe Points](#terminal-currents-and-probe-points)
- [Solution Analysis](#solution-analysis)
- [Visualization](#visualization)
- [HDF5 Output Format](#hdf5-output-format)
- [Testing and Validation](#testing-and-validation)
- [Reproducibility Checklist](#reproducibility-checklist)
- [Known Limitations](#known-limitations)
- [License](#license)
- [Citation](#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

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

The covariant gradient uses the link variable:

```text
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:

```text
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

```text
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:

```text
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

```bash
pip install .
```

### Editable Development Install

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

### Install From requirements.txt

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

```bash
pip install -r requirements.txt
```

### Optional GPU Package

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

```bash
pip install cupy-cuda12x
```

For CUDA 11.x:

```bash
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.

```python
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:

```python
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:

```python
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:

```python
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:

```python
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

```python
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:

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

Plot the surface projection or filled surface projection:

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

### Physical Scales

`Device` exposes useful dimensional scales through Pint:

```python
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` | `|psi|` imposed at terminal sites |
| `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:

```python
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:

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

Callable field:

```python
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`.

```python
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.

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

### Manual Solver Access

For advanced workflows:

```python
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

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

The source uses the symmetric gauge:

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

### Circular Current Loop

```python
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

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

Custom `Parameter` objects can be combined arithmetically:

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

## Screening and GPU Acceleration

Enable magnetic self-field screening:

```python
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:

```python
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.

```python
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:

```python
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:

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

Current-bias feedback:

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

Time-dependent driving:

```python
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:

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

Time-dependent currents:

```python
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:

```python
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:

```python
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:

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

Common properties:

```python
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:

```python
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:

```python
frame = solution.closest_solve_step(target_time=5.0)
```

### Derived Quantities

```python
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:

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

Current density:

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

Regular 3D grid:

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

### Fluxoid Quantization

Circular contour:

```python
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:

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

Rectangular contour:

```python
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:

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

For time series through named paths:

```python
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)`:

```python
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

```python
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

```python
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:

```bash
tdgl3d-visualize
```

Examples:

```bash
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.

```text
/
  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:

```bash
pytest tdgl3d/test
```

Or use the package helper:

```python
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:

```bash
pytest tdgl3d/test/test_fast.py
```

For simulation-level tests, run:

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

## Reproducibility Checklist

For reproducible simulation studies, record:

- Package version and dependency versions:

```python
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:

```bash
pip install pyvista meshio
```

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

```python
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).
