`pyloggrid.LogGrid.Framework`

Solver for the log grids.

Module Contents

Classes

`_CustIntegrator`	Abstract integrator class
`ViscDopri`	DOPRI5-based solver made specifically for log grids.
`ETD4RK`	Based on Cox & Matthews (2002). Exponential time differencing for stiff systems. Journal of Computational Physics, 176(2), 430-455.
`ETD35`	Based on https://github.com/whalenpt/rkstiff
`Solver`	Generic top-level object to handle solving log-grid equations

Functions

`fields_to_1D`(→ numpy.ndarray)	Convert N-dimensional fields to 1D arrays
`D1_to_fields`() → dict[str, numpy.ndarray])	Reverse of `fields_to_1D()`

fields_to_1D(fields: dict[str, ndarray], field_names: list[str]) → ndarray

Convert N-dimensional fields to 1D arrays

Parameters:

field_names – [name1, …] name of the fields (order is important, as it it used to convert back from 1D)
fields – N-dimensional named fields, grid-shaped

Returns:

1D ndarray containing all the appended 1D fields, as ordered by field_names

D1_to_fields(linear: ndarray, field_names: list[str], field_shape: int, Ellipsis) → dict[str, ndarray]

Reverse of fields_to_1D()

Converts one 1D array to a dict of the encoded fields in their original shape

Parameters:

linear – 1D ndarray of fields
field_names – [name1, …] name of the fields (order is important)
field_shape – shape of the grid / fields

Returns:

the fields with their full shape

class _CustIntegrator(equation_nl: Callable[[float, ndarray], ndarray], equation_l: Callable[[float, ndarray], ndarray], save_step: Callable[[float, float, ndarray], bool], init_t: float, solver_params: dict, y0: ndarray)

Abstract integrator class

select_initial_step(t0: float, y0: ndarray, order: int, rtol: float = 1e-06, atol: float = 1e-06) → float

Empirically select a good initial step. The algorithm is described in [1].

Parameters:

t0 – Initial value of the independent variable.
y0 – Initial value of the dependent variable.
order – Error estimator order. It means that the error controlled by the algorithm is proportional to step_size ** (order + 1).
rtol – Desired relative tolerance.
atol – Desired absolute tolerance.

Returns:

h_abs = Absolute value of the suggested initial step.

References

abstract solve()

class ViscDopri(equation_nl: Callable[[float, ndarray], ndarray], equation_l: Callable[[float, ndarray], ndarray], save_step: Callable[[float, float, ndarray], bool], init_t: float, solver_params: dict, y0: ndarray, dt_params: dict[str, float | None])

Bases: _CustIntegrator

DOPRI5-based solver made specifically for log grids.

A nonlinear update step is decoupled from the linear update step. Special care is given to handling numerical (float-induced) random numerical errors which yield huge gradients. Based on scipy’s DOPRI5 implementation.

Parameters:

equation_nl – linear-less update step
equation_l – linear update step
init_t – initial time
y0 – initial 1D array
solver_params – = {"atol"[absolute tolerance]: 1e-6, "rtol"[relative tolerance]: 1e-4, "exact_viscous_splitting"[if True, use exact viscous splitting]: False}
dt_params – {dt0 = initial time step, dtmin, dtmax}

rk_step(fun: Callable[[float, ndarray], ndarray], t: float, y: ndarray, h: float, A: ndarray, B: ndarray, C: ndarray, K: ndarray) → tuple[ndarray, ndarray]

Perform a single Runge-Kutta step. This function computes a prediction of an explicit Runge-Kutta method and also estimates the error of a less accurate method. Notation for Butcher tableau is as in [2].

Parameters:

fun – Right-hand side of the system.
t – Current time.
y – Current state.
h – Step to use.
A – Coefficients for combining previous RK stages to compute the next stage. For explicit methods the coefficients at and above the main diagonal are zeros.
B – Coefficients for combining RK stages for computing the final prediction.
C – Coefficients for incrementing time for consecutive RK stages. The value for the first stage is always zero.
K – Storage array for putting RK stages here. Stages are stored in rows. The last row is a linear combination of the previous rows with coefficients

Returns: tuple (y_new, f_new) where y_new is the solution at t + h computed with a higher accuracy, and f_new is the derivative fun(t + h, y_new).

References

rkdp(force_step: bool = False) → tuple[ndarray, float, float]: Performs one RKDP step

solve() → None

Main solving loop.

The solver runs until an interruption is fired by the update.

class ETD4RK(equation_nl: Callable[[float, ndarray], ndarray], equation_l: Callable[[float, ndarray], ndarray], save_step: Callable[[float, float, ndarray], bool], init_t: float, y0: ndarray, dt_params: dict[str, float | None])

Bases: _CustIntegrator

Based on Cox & Matthews (2002). Exponential time differencing for stiff systems. Journal of Computational Physics, 176(2), 430-455.

Warning

Only designed for constant-in-time viscosities

Note

See Utils/thresholds.py in the source code for the determination of the thresholds used in computing Taylor series

Parameters:

equation_nl – linear-less update step
equation_l – linear update step
init_t – initial time
y0 – initial 1D array
dt_params – {dt0 = initial time step, dtmin, dtmax}

static dopri_step(fun: Callable[[float, ndarray], ndarray], t: float, y: ndarray, h: float, A: ndarray, B: ndarray, C: ndarray, K: ndarray) → tuple[ndarray, ndarray]

Perform a single Runge-Kutta step. This function computes a prediction of an explicit Runge-Kutta method and also estimates the error of a less accurate method. Notation for Butcher tableau is as in [3].

Parameters:

fun – Right-hand side of the system.
t – Current time.
y – Current state.
h – Step to use.
A – Coefficients for combining previous RK stages to compute the next stage. For explicit methods the coefficients at and above the main diagonal are zeros.
B – Coefficients for combining RK stages for computing the final prediction.
C – Coefficients for incrementing time for consecutive RK stages. The value for the first stage is always zero.
K – Storage array for putting RK stages here. Stages are stored in rows. The last row is a linear combination of the previous rows with coefficients

Returns: tuple (y_new, f_new) where y_new is the solution at t + h computed with a higher accuracy, and f_new is the derivative fun(t + h, y_new). .. rubric:: References

rk() → ndarray: Performs one Rk step

solve() → None: Main solving loop. The solver runs until an interruption is fired by the update.

class ETD35(equation_nl: Callable[[float, ndarray], ndarray], equation_l: Callable[[float, ndarray], ndarray], save_step: Callable[[float, float, ndarray], bool], init_t: float, solver_params: dict, y0: ndarray, dt_params: dict[str, float | None])

Bases: _CustIntegrator

Based on https://github.com/whalenpt/rkstiff

Warning

Only supports diagonal (and constant-in-time) linear terms for now

Parameters:

equation_nl – linear-less update step
equation_l – linear update step
init_t – initial time
solver_params – {"rtol"[relative tolerance]: 1e-2, "adapt_cutoff"[Limits values used in the computation of the suggested step size to those with |u| > adapt_cutoff*max(|u|)]: 1e-2, "minh"[minimum time step]: 1e-16}
y0 – initial 1D array
dt_params – {dt0 = initial time step, dtmin, dtmax}

rk() → tuple[ndarray, float, float]

Performs one Rk step

Returns:: tuple (y_new, h, h_new_suggested)

solve() → None: Main solving loop. The solver runs until an interruption is fired by the update.

exception SolverInterruptedError

Bases: InterruptedError

custom error thrown to change grid size

class Solver(D: int, l_params: dict[str, float | bool], fields_names: list[str], equation_nl: Callable[[float, Grid, dict], dict], equation_l: Callable[[float, Grid, dict], dict[str, ndarray]], k_min: float = None, k0: bool = False, simu_params: dict = None, n_threads: int = None)

Generic top-level object to handle solving log-grid equations

Parameters:

D – dimension of physical space (usually 1, 2 or 3)
l_params – parameters that define the grid parameter ̀`l` as a dict with keys {"a", "b", "plastic"}. See below for more detail.
fields_names – the names of the scalar grid-shaped fields
equation_nl – the function that performs the nonlinear update step. Returns the time derivative of the fields
equation_l – implicit visocsity update step, called after the RK step has ended. Returns the new fields and their new time derivative
k_min – minimum k of the grid. Default is defined by pyloggrid.LogGrid.Grid.Grid.
k0 – whether there’s a k0 mode
simu_params – physical quantities relevant to the simulation, fixed for the whole simulation.
n_threads – Number of threads to use, default is max/2 (not recommended, run benchmarking for better results. If running batch simulations, 1 is optimal.)

In pyloggrid, simulating equations is done through the Solver class.

Initializing the simulation

A Solver object is initialized with parameters of the simulation:

Geometric parameters (D, l_params, k_min, k0)
- l_params is a dict that describes the scaling parameter l of the grid. If the “plastic” key is set to True (l_params={"plastic": True, ...}) then l~=1.3 (the plastic number). If plastic is not True, then keys a, b must be specified and have integer values, in which case l is the solution to l**b-l**a=1. In particular if l_params={"a": 1, "b": 2} then l~=1.6 (the golden number). To get the special case l=2, use l_params={"a":None, "b": None}.
Important

Choosing either "plastic": True, "a": int, "b": int or "a": None, "b": None doesn’t just change the scaling of the grid, but also the number of interactions per point for the convolution (listed here in decreasing order), which affects both the physics and the required computing power.
- k0 is a bool that indicates whether the axes \(k_i=0\) are included. Enabling this comes at a significant performance cost, and greatly modifies the interaction between grid points.
- k_min is the minimum absolute value of a wavevector on each axis, excluding zeros. This is not the minimum magnitude of a wavevector on the grid, which would be \(\sqrt{D}k_{min}\) in \(D\) dimensions.
Equation parameters (field_names, equation_nl, equation_l, simu_params)
- For stability and performance reasons, the equation to solve is divided into two parts: one that is proportional to the scalar input equation_l, and the rest equation_nl. If your equation is \(\partial_t u=A(u)−b\cdot u\), then the linear part is \(−b\), and the nonlinear part is \(A(t, u)\).
- equation_l, equation_nl expect functions that take as arguments t: float, grid: Grid, simu_params: dict and return a dict {"field_name": d/dt_field_name, ...}.
Warning

Note that the linear part must be time-independent.
- simu_params are arbitrary (serializable) constants used in the simulation. The fields named by fields_names, on the other hand, are grid-shaped complex scalars., i.e. “the fields being simulated”.

Running the simulation

Simulations are ran by calling the solve() method.

Warning

Although in theory you can use the solve method several times on a single Solver object, this isn’t a common use case. Do so at your own risks, but please report unexpected behavior to us.

_load_parameters(path: str) → tuple[float, float]

Loads simulation parameters from settings.json. Overwrites instance variables

Parameters:: path – directory to load the settings from
Returns:: (last used time, last used timstep)

_save_step_all(t: float, dt: float, path: str) → None

Save the current step. Saves both fields and settings.

Parameters:

t – current time
dt – timestep
path – folder to save to. Will override any existing file.

solve(save_path: str, initial_conditions: str | Callable[[dict[str, ndarray], Grid, dict], dict[str, ndarray]] | tuple | list, solver_params: dict = None, init_t: float = 0, end_simulation=None, save_one_in: float = 1, update_gridsize_cb: Callable[[Grid], int | None] | None = None, dt_params: dict[str, float | None] = None, solver: Literal[ViscDopri, ETD4RK, ETD35] = 'ETD35') → None

Solve the solver’s equation with the provided numerical parameters.

Parameters:

save_path – where to save the fields and simulation parameters. Will backup if exists, will create otherwise
initial_conditions – if str: either "loadfromsave" to resume simulation, [legacy: or the path of the .npz save file whence to initialize the fields]. If tuple (name, step): load step step from .h5 name. Otherwise: function that returns the initial fields.
solver_params – params forwarded to the solver
init_t – initial simulation time
end_simulation – dict of thresholds to end the simulation: {t: max physical time, elapsed_time: max real time spent computing, step: max save step}
save_one_in – save one step to file every X ode_step
update_gridsize_cb – optional callback to change the grid size after a step. Returns None if no change is to be made, returns the new grid size otherwise.
dt_params – {dt0: initial timestep, dtmin, dtmax}
solver – "ETD35" for ETD35 (default), "ETD4RK" for ETD4RK, "ViscDopri" for ViscDopri

Details on a few parameters

If save_path="<path>/<foldername>" already exists on disk, the existing directory will be renamed to <path>/bk_<foldername>_<TIMESTAMP>. Otherwise, it will be automatically created.
If initial_conditions="loadfromsave" then settings and initial conditions will be loaded from the last step of the existing simulation in save_path. This will override the following arguments: end_simulation, l_params, D, fields_name, simu_params, k_min, k0.
To load from a specific step of a .h5 simulation file, set loadfromsave=(<filepath>:str, <step>:int).
Otherwise, initial_conditions must be a function that takes (fields, grid, simu_params) and returns the initial fields.
update_gridsize_cb is an optional function to dynamically change the grid size. It takes (grid) and returns either None (no change) or a new grid size.

pyloggrid.LogGrid.Framework

Module Contents

Classes

Functions

Initializing the simulation

Running the simulation

Details on a few parameters

`pyloggrid.LogGrid.Framework`