Field cooling in TDGL? about py-tdgl HOT 1 CLOSED

import numpy as np
import tdgl
from tdgl.geometry import circle

# Define the device
film_radius = 1
hole_radius = film_radius / 3
layer = tdgl.Layer(london_lambda=0.1, coherence_length=0.05, thickness=0.01)
film = tdgl.Polygon("film", points=circle(film_radius)).resample(201)
hole = tdgl.Polygon("hole", points=circle(hole_radius)).resample(201)
device = tdgl.Device(
    "ring",
    layer=layer,
    film=film,
    holes=[hole],
    length_units="um",
)

# Make the mesh
device.make_mesh(max_edge_length=0.025)

# Setup the field cooling protocol:
# 1. Use disorder_epsilon to create a temporary normal metal channel connecting the
#    film exterior and the hole.
# 2. Apply a uniform magnetic field sufficient to trap a single flux quantum in the hole.
# 3. Remove the normal metal channel, then turn off the magnetic field.

options = tdgl.SolverOptions(
    output_file="ring.h5",
    solve_time=250,
    field_units="mT",
    current_units="uA",
    dt_max=0.01,
)

def epsilon(r, *, t):
    """Sets T_c = 0 for a portion of the film until time t = 100."""
    x, y = r  # r = (x, y) is the position within the film
    # t is time in units of \tau_0
    if t < 100 and y > 0 and np.abs(x) < 0.1:
        return -1
    return 1


def scale_func(x, y, z, *, t):
    """A scale factor that turns off the applied field after time t = 200."""
    if t < 200:
        return 1
    return 0

# Apply a uniform magnetic field equal to 1.5 * Phi_0 / film_area
# to trap one Phi_0 in the hole
film_area = np.pi * (film_radius * tdgl.ureg("um"))**2
applied_field = 1.5 * (tdgl.ureg("Phi_0") / film_area).to("mT")
print(f"Applied field = {applied_field}")

# Combine the scale factor and uniform magnetic field
# to define the applied magnetic vector potential.
applied_vector_potential = (
    # Time-dependent scale factor
    tdgl.sources.Scale(scale_func)
    # Uniform applied magnetic field
    * tdgl.sources.ConstantField(
        applied_field.to(options.field_units).magnitude,
        length_units=device.length_units,
        field_units=options.field_units,
    )
)

# Solve the model
solution = tdgl.solve(
    device,
    options,
    applied_vector_potential=applied_vector_potential,
    disorder_epsilon=epsilon,
)