Solver

Simulating the dynamics of a tdgl.Device for a given applied magnetic vector potential and set of bias currents is as simple as calling the tdgl.solve() function. The solver implements the finite-volume and implicit Euler methods described in detail in Theoretical Background. The behavior of the solver is determined an instance of tdgl.SolverOptions.

The applied vector potential can be specified as a scalar (indicating the vector potential associated with a uniform magnetic field), a function with signature func(x, y, z) -> [Ax, Ay, Az], or a tdgl.Parameter. The physical units for the applied vector potential are field_units * device.length_units.

The bias or terminal currents (if any) can be specified as a dictionary like terminal_currents = {terminal_name: current}, where current is a float in units of the specified current_units. For time-dependent applied currents, one can provide a function with signature terminal_currents(time: float) -> {terminal_name: current}, where time is the dimensionless time. In either case, the sum of all terminal currents must be zero at every time step and every terminal in the device must be included in the dictionary to ensure current conservation.

tdgl.solve(device, options, applied_vector_potential=0, terminal_currents=None, disorder_epsilon=1, seed_solution=None)

Solve a TDGL model.

Parameters:

• device (Device) – The tdgl.Device to solve.

• options (SolverOptions) – An instance tdgl.SolverOptions specifying the solver parameters.

• applied_vector_potential (Union[Callable, float]) – A function or tdgl.Parameter that computes the applied vector potential as a function of position (x, y, z), or of position and time (x, y, z, *, t). If a float B is given, the applied vector potential will be that of a uniform magnetic field with strength B field_units.

• terminal_currents (Union[Callable, Dict[str, float], None]) – A dict of {terminal_name: current} or a callable with signature func(time: float) -> {terminal_name: current}, where current is a float in units of current_units and time is the dimensionless time.

• disorder_epsilon (Union[float, Callable]) – A float in range [-1, 1], or a callable with signature disorder_epsilon(r: Tuple[float, float]) -> epsilon, where epsilon is a float in range [-1, 1]. Setting \(\epsilon(\mathbf{r})=T_c(\mathbf{r})/T_c - 1 < 1\) suppresses the critical temperature at position \(\mathbf{r}\), which can be used to model inhomogeneity.

• seed_solution (Optional[Solution]) – A tdgl.Solution instance to use as the initial state for the simulation.

Return type:

Optional[Solution]

Returns:

A tdgl.Solution instance. Returns None if the simulation was cancelled during the thermalization stage.

class tdgl.SolverOptions(solve_time, skip_time=0.0, dt_init=1e-06, dt_max=0.1, adaptive=True, adaptive_window=10, max_solve_retries=10, adaptive_time_step_multiplier=0.25, output_file=None, terminal_psi=0.0, gpu=False, sparse_solver=SparseSolver.SUPERLU, pause_on_interrupt=True, save_every=100, progress_interval=0, monitor=False, monitor_update_interval=1.0, field_units='mT', current_units='uA', include_screening=False, max_iterations_per_step=1000, screening_tolerance=0.001, screening_step_size=0.1, screening_step_drag=0.5)

Options for the TDGL solver.

Parameters:

• solve_time (float) – Total simulation time, after any thermalization.

• skip_time (float) – Amount of ‘thermalization’ time to simulate before recording data.

• dt_init (float) – Initial time step.

• dt_max (float) – Maximum adaptive time step.

• adaptive (bool) – Whether to use an adpative time step. Setting dt_init = dt_max is equivalent to setting adaptive = False.

• adaptive_window (int) – Number of most recent solve steps to consider when computing the time step adaptively.

• max_solve_retries (int) – The maximum number of times to reduce the time step in a given solve iteration before giving up.

• adaptive_time_step_multiplier (float) – The factor by which to multiple the time step dt for each adaptive solve retry.

• terminal_psi (Union[float, complex, None]) – Fixed value for the order parameter in current terminals.

• output_file (Optional[str]) – Path to an HDF5 file in which to save the data. If the file name already exists, a unique name will be generated. If output_file is None, the solver results will not be saved to disk.

• gpu (bool) – Use the GPU via CuPy. This option requires a GPU and the CuPy Python package, which can be installed via pip.

• sparse_solver (Union[SparseSolver, str]) – One of "superlu", "umfpack", "pardiso", or "cupy". "umfpack" requires suitesparse, which can be installed via conda, and scikit-umfpack, which can be installed via pip. "pardiso" requires an Intel CPU and the pypardiso package, which can be installed via pip or conda. "cupy" requires a GPU and the CuPy Python package, which can be installed via pip.

• field_units (str) – The units for magnetic fields.

• current_units (str) – The units for currents.

• pause_on_interrupt (bool) – Pause the simulation in the event of a KeyboardInterrupt.

• save_every (int) – Save interval in units of solve steps.

• progress_interval (int) – Minimum number of solve steps between progress bar updates.

• monitor (bool) – Plot data in real time as the simulation is running.

• monitor_update_interval (float) – The monitor update interval in seconds.

• include_screening (bool) – Whether to include screening in the simulation.

• max_iterations_per_step (int) – The maximum number of screening iterations per solve step.

• screening_tolerance (float) – Relative tolerance for the induced vector potential, used to evaluate convergence of the screening calculation within a single time step.

• screening_step_size (float) – Step size \(\alpha\) for Polyak’s method.

• screening_step_drag (float) – Drag parameter \(\beta\) for Polyak’s method.

class tdgl.TDGLSolver(device, options, applied_vector_potential=0.0, terminal_currents=None, disorder_epsilon=1.0, seed_solution=None)

Solver for a TDGL model.

An instance of tdgl.TDGLSolver is created and executed by calling tdgl.solve().

Parameters:

• device (Device) – The tdgl.Device to solve.

• options (SolverOptions) – An instance tdgl.SolverOptions specifying the solver parameters.

• applied_vector_potential (Union[Callable, float]) – A function or tdgl.Parameter that computes the applied vector potential as a function of position (x, y, z), or of position and time (x, y, z, *, t). If a float B is given, the applied vector potential will be that of a uniform magnetic field with strength B field_units. If the applied vector potential is time-dependent, this argument must be a tdgl.Parameter.

• terminal_currents (Union[Callable, Dict[str, float], None]) – A dict of {terminal_name: current} or a callable with signature func(time: float) -> {terminal_name: current}, where current is a float in units of current_units and time is the dimensionless time.

• disorder_epsilon (Union[Callable, float]) – A float <= 1, or a function that returns \(\epsilon\leq 1\) as a function of position r=(x, y) or position and time (x, y, *, t). Setting \(\epsilon(\mathbf{r}, t)=T_c/T - 1 < 1\) suppresses the order parameter at position \(\mathbf{r}=(x, y)\), which can be used to model inhomogeneity.

• seed_solution (Optional[Solution]) – A tdgl.Solution instance to use as the initial state for the simulation.

update_mu_boundary(time)

Computes the terminal current density for a given time step and updates the scalar potential boundary conditions accordingly.

Parameters:

time (float) – The current value of the dimensionless time.

Return type:

None

update_applied_vector_potential(time)

Evaluates the time-dependent vector potential.

Parameters:

time (float) – The current value of the dimensionless time.

Return type:

ndarray

Returns:

The new value of the applied vector potential.

update_epsilon(time)

Evaluates the time-dependent disorder parameter \(\epsilon\).

Parameters:

time (float) – The current value of the dimensionless time.

Return type:

ndarray

Returns:

The new value of \(\epsilon\)

static solve_for_psi_squared(*, psi, abs_sq_psi, mu, epsilon, gamma, u, dt, psi_laplacian)

Solves for \(\psi^{n+1}\) and \(|\psi^{n+1}|^2\) given \(\psi^n\) and \(\mu^n\).

Parameters:

• psi (ndarray) – The current value of the order parameter, \(\psi^n\)

• abs_sq_psi (ndarray) – The current value of the superfluid density, \(|\psi^n|^2\)

• mu (ndarray) – The current value of the electric potential, \(\mu^n\)

• epsilon (ndarray) – The disorder parameter, \(\epsilon\)

• gamma (float) – The inelastic scattering parameter, \(\gamma\).

• u (float) – The ratio of relaxation times for the order parameter, \(u\)

• dt (float) – The time step

• psi_laplacian (spmatrix) – The covariant Laplacian for the order parameter

Return type:

Optional[Tuple[ndarray, ndarray]]

Returns:

None if the calculation failed to converge, otherwise the new order parameter \(\psi^{n+1}\) and superfluid density \(|\psi^{n+1}|^2\).

adaptive_euler_step(step, psi, abs_sq_psi, mu, epsilon, dt)

Updates the order parameter and time step in an adaptive Euler step.

Parameters:

• step (int) – The solve step index, \(n\)

• psi (ndarray) – The current value of the order parameter, \(\psi^n\)

• abs_sq_psi (ndarray) – The current value of the superfluid density, \(|\psi^n|^2\)

• mu (ndarray) – The current value of the electric potential, \(\mu^n\)

• epsilon (ndarray) – The disorder parameter \(\epsilon^n\)

• dt (float) – The tentative time step, which will be updated

Return type:

Tuple[ndarray, ndarray, float]

Returns:

\(\psi^{n+1}\), \(|\psi^{n+1}|^2\), and \(\Delta t^{n}\).

solve_for_observables(psi, dA_dt)

Solves for the scalar potential \(\mu\), the supercurrent density \(\mathbf{J}_s\), and the normal current density \(\mathbf{J}_n\).

Parameters:

• psi (ndarray) – The order parameter.

• dA_dt (Union[float, ndarray]) – The time-derivative of the vector potential.

Return type:

Tuple[ndarray, ndarray, ndarray]

Returns:

\(\mu\), \(\mathbf{J}_s\), and \(\mathbf{J}_n\)

get_induced_vector_potential(current_density, A_induced_vals, velocity)

Computes a new value of the induced vector potential based on Polyak’s method.

Parameters:

• current_density (ndarray) – The total current density \(\mathbf{J}_s + \mathbf{J}_n\)

• A_induced_vals (List[ndarray]) – A running list of the induced vector potential for previous iterations of Polyak’s method.

• velocity (List[ndarray]) – A running list of the “velocities” for previous iterations of Polyak’s method.

Return type:

Tuple[ndarray, float]

Returns:

A new value for the induced vector potential, and the relative error in the induced vector potential between this iteration of Polyak’s method and the previous iteration.

update(state, running_state, dt, *, psi, mu, supercurrent, normal_current, induced_vector_potential, applied_vector_potential=None, epsilon=None)

This method is called at each time step to update the state of the system.

Parameters:

• state (Dict[str, Real]) – The solver state, i.e., the solve step, time, and time step

• running_state (RunningState) – A container for scalar data that is saved at each time step

• dt (float) – The time step for the previous solve step

• psi (ndarray) – The order parameter

• mu (ndarray) – The scalar potential

• supercurrent (ndarray) – The supercurrent density

• normal_current (ndarray) – The normal current density

• induced_vector_potential (ndarray) – The induced vector potential

• applied_vector_potential (Optional[ndarray]) – The applied vector potential. This will be None in the case of a time-independent vector potential.

• epsilon (Optional[ndarray]) – The disorder parameter epsilon. This will be None in the case of a time-independent epsilon.

Return type:

SolverResult

Returns:

A tdgl.SolverResult instance for the solve step.

solve()

Runs the solver.

Return type:

Optional[Solution]

Returns:

A tdgl.Solution instance. Returns None if the simulation was cancelled during the thermalization stage.

enum tdgl.solver.options.SparseSolver(value)

Supported sparse linear solvers.

Valid values are as follows:

SUPERLU: str = <SparseSolver.SUPERLU: 'superlu'>

UMFPACK: str = <SparseSolver.UMFPACK: 'umfpack'>

PARDISO: str = <SparseSolver.PARDISO: 'pardiso'>

CUPY: str = <SparseSolver.CUPY: 'cupy'>

class tdgl.Parameter(func, time_dependent=False, **kwargs)

A callable object that computes a scalar or vector quantity as a function of position coordinates x, y (and optionally z and time t).

Addition, subtraction, multiplication, and division between multiple Parameters and/or numbers is supported. The result of any of these operations is a CompositeParameter object.

Parameters:

• func (Callable) – A callable/function that actually calculates the parameter’s value. The function must take x, y (and optionally z) as the first and only positional arguments, and all other arguments must be keyword arguments. Therefore func should have a signature like func(x, y, z, a=1, b=2, c=True), func(x, y, *, a, b, c), func(x, y, z, *, a, b, c), or func(x, y, z, *, a, b=None, c=3). For time-dependent Parameters, func must also take time t as a keyword-only argument.

• time_dependent (bool) – Specifies that func is a function of time t.

• kwargs – Keyword arguments for func.