implementations.problem_classes.Quench module

class Quench(Cv=1000.0, K=1000.0, u_thresh=0.03, u_max=0.06, Q_max=1.0, leak_range=(0.45, 0.55), leak_type='linear', leak_transition='step', order=2, stencil_type='center', bc='neumann-zero', nvars=128, newton_tol=1e-08, newton_maxiter=99, lintol=1e-08, liniter=99, direct_solver=True, inexact_linear_ratio=None, min_lintol=1e-12, reference_sol_type='scipy')[source]

Bases: ptype

This is a toy problem [1] to emulate a magnet that has been cooled to temperatures where superconductivity is possible. However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment sufficiently, there will be a runaway effect heating up the entire magnet. This effect has actually lead to huge magnets being destroyed at CERN in the past and hence warrants investigation.

The model we use is a 1d heat equation with Neumann-zero boundary conditions, meaning this magnet is totally insulated from its environment except for the leak. We add a non-linear term that heats parts of the domain that exceed a certain temperature threshold as well as the leak itself.

The problem is discretised with finite difference in space and treated fully-implicitly.

Parameters:

  • Cv (float, optional) – Volumetric heat capacity.

  • K (float, optional) – Thermal conductivity.

  • u_thresh (float, optional) – Threshold for temperature.

  • u_max (float, optional) – Maximum temperature.

  • Q_max (float, optional) – Maximum heat source power density.

  • leak_range (tuple of float) – Range of the leak.

  • leak_type (str, optional) – Type of leak, choose between 'linear' or 'exponential'.

  • leak_transition (str, optional) – Indicates how the heat in the leak propagates, choose between 'step' and 'Gaussian'.

  • order (int, optional) – Order of the finite difference discretization.

  • stencil_type (str, optional) – Type of stencil for finite differences.

  • bc (str, optional) – Type of boundary conditions. Default is 'neumann-zero'.

  • nvars (int, optional) – Spatial resolution.

  • newton_tol (float, optional) – Tolerance for Newton to terminate.

  • newton_maxiter (int, optional) – Maximum number of Newton iterations to be done.

  • lintol (float, optional) – Tolerance for linear solver to be done.

  • liniter (int, optional) – Maximum number of linear iterations inside the Newton solver.

  • direct_solver (bool, optional) – Indicates if a direct solver should be used.

  • inexact_linear_ratio (float, optional) – Ratio of tolerance of linear solver to the Newton residual, overrides lintol

  • min_lintol (float, optional) – Minimal tolerance for the linear solver

  • reference_sol_type (str, optional) – Indicates which method should be used to compute a reference solution. Choose between 'scipy', 'SDC', or 'DIRK'.

A

FD discretization matrix of the ND grad operator.

Type:

sparse matrix (CSC)

Id

Identity matrix of the same dimension as A.

Type:

sparse matrix (CSC)

dx

Distance between two spatial nodes.

Type:

float

xv

Spatial grid values.

Type:

np.1darray

leak

Indicates the leak.

Type:

np.1darray of bool

References

dtype_f

alias of mesh

dtype_u

alias of mesh

eval_f(u, t)[source]

Evaluate the full right-hand side.

Parameters:

  • u (dtype_u) – Current values of the numerical solution.

  • t (float) – Current time at which the numerical solution is computed.

Returns:

f – The right-hand side of the problem.

Return type:

dtype_f

eval_f_non_linear(u, t)[source]

Get the non-linear part of f.

Parameters:

  • u (dtype_u) – Current values of the numerical solution:

  • t (float) – Current time at which the numerical solution is computed.

Returns:

me – The non-linear part of the right-hand side.

Return type:

dtype_u

get_non_linear_Jacobian(u)[source]

Evaluate the non-linear part of the Jacobian only.

Parameters:

u (dtype_u) – Current values of the numerical solution.

Returns:

The derivative of the non-linear part of the solution w.r.t. to the solution.

Return type:

scipy.sparse.csc

solve_system(rhs, factor, u0, t)[source]

Simple Newton solver for \((I - factor \cdot f)(\vec{u}) = \vec{rhs}\).

Parameters:

  • rhs (dtype_f) – Right-hand side.

  • factor (float) – Abbrev. for the local stepsize (or any other factor required).

  • u0 (dtype_u) – Initial guess for the iterative solver.

  • t (float) – Current time (e.g. for time-dependent BCs).

Returns:

u – The solution as mesh.

Return type:

dtype_u

u_exact(t, u_init=None, t_init=None)[source]

Routine to compute the exact solution at time \(t\).

Parameters:

t (float) – Time of the exact solution.

Returns:

me – The exact solution.

Return type:

dtype_u

class QuenchIMEX(Cv=1000.0, K=1000.0, u_thresh=0.03, u_max=0.06, Q_max=1.0, leak_range=(0.45, 0.55), leak_type='linear', leak_transition='step', order=2, stencil_type='center', bc='neumann-zero', nvars=128, newton_tol=1e-08, newton_maxiter=99, lintol=1e-08, liniter=99, direct_solver=True, inexact_linear_ratio=None, min_lintol=1e-12, reference_sol_type='scipy')[source]

Bases: Quench

This is a toy problem [1] to emulate a magnet that has been cooled to temperatures where superconductivity is possible. However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment sufficiently, there will be a runaway effect heating up the entire magnet. This effect has actually lead to huge magnets being destroyed at CERN in the past and hence warrants investigation.

The model we use is a 1d heat equation with Neumann-zero boundary conditions, meaning this magnet is totally insulated from its environment except for the leak. We add a non-linear term that heats parts of the domain that exceed a certain temperature threshold as well as the leak itself.

The problem is discretised with finite difference in space and treated semi-implicitly.

dtype_f

alias of imex_mesh

eval_f(u, t)[source]

Routine to evaluate the right-hand side of the problem.

Parameters:

  • u (dtype_u) – Current values of the numerical solution.

  • t (float) – Current time of the numerical solution is computed.

Returns:

f – The right-hand side of the problem.

Return type:

dtype_f

solve_system(rhs, factor, u0, t)[source]

Simple linear solver for \((I - factor \cdot f_{expl})(\vec{u}) = \vec{rhs}\).

Parameters:

  • rhs (dtype_f) – Right-hand side for the linear system.

  • factor (float) – Abbrev. for the local stepsize (or any other factor required).

  • u0 (dtype_u) – Initial guess for the iterative solver.

  • t (float) – Current time (e.g. for time-dependent BCs).

Returns:

me – The solution as mesh.

Return type:

dtype_u

u_exact(t, u_init=None, t_init=None)[source]

Routine to compute the exact solution at time \(t\).

Parameters:

t (float) – Time of the exact solution.

Returns:

me – The exact solution.

Return type:

dtype_u