:py:mod:`pyloggrid.LogGrid.Framework` ===================================== .. py:module:: pyloggrid.LogGrid.Framework .. autoapi-nested-parse:: Solver for the log grids. Module Contents --------------- Classes ~~~~~~~ .. autoapisummary:: pyloggrid.LogGrid.Framework._CustIntegrator pyloggrid.LogGrid.Framework.ViscDopri pyloggrid.LogGrid.Framework.ETD4RK pyloggrid.LogGrid.Framework.ETD35 pyloggrid.LogGrid.Framework.Solver Functions ~~~~~~~~~ .. autoapisummary:: pyloggrid.LogGrid.Framework.fields_to_1D pyloggrid.LogGrid.Framework.D1_to_fields .. py:function:: fields_to_1D(fields: dict[str, numpy.ndarray], field_names: list[str]) -> numpy.ndarray Convert N-dimensional fields to 1D arrays :param field_names: [name1, ...] name of the fields (order is important, as it it used to convert back from 1D) :param fields: N-dimensional named fields, grid-shaped :returns: 1D ndarray containing all the appended 1D fields, as ordered by ``field_names`` .. py:function:: D1_to_fields(linear: numpy.ndarray, field_names: list[str], field_shape: (int, Ellipsis)) -> dict[str, numpy.ndarray] Reverse of :func:`fields_to_1D` Converts one 1D array to a dict of the encoded fields in their original shape :param linear: 1D ndarray of fields :param field_names: [name1, ...] name of the fields (order is important) :param field_shape: shape of the grid / fields :returns: the fields with their full shape .. py:class:: _CustIntegrator(equation_nl: Callable[[float, numpy.ndarray], numpy.ndarray], equation_l: Callable[[float, numpy.ndarray], numpy.ndarray], save_step: Callable[[float, float, numpy.ndarray], bool], init_t: float, solver_params: dict, y0: numpy.ndarray) Abstract integrator class .. py:method:: select_initial_step(t0: float, y0: numpy.ndarray, order: int, rtol: float = 1e-06, atol: float = 1e-06) -> float Empirically select a good initial step. The algorithm is described in [#]_. :param t0: Initial value of the independent variable. :param y0: Initial value of the dependent variable. :param order: Error estimator order. It means that the error controlled by the algorithm is proportional to ``step_size ** (order + 1)``. :param rtol: Desired relative tolerance. :param atol: Desired absolute tolerance. :returns: ``h_abs`` = Absolute value of the suggested initial step. .. rubric:: References .. [#] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential Equations I: Nonstiff Problems", Sec. II.4. .. py:method:: solve() :abstractmethod: .. py:class:: ViscDopri(equation_nl: Callable[[float, numpy.ndarray], numpy.ndarray], equation_l: Callable[[float, numpy.ndarray], numpy.ndarray], save_step: Callable[[float, float, numpy.ndarray], bool], init_t: float, solver_params: dict, y0: numpy.ndarray, dt_params: dict[str, Optional[float]]) Bases: :py:obj:`_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. :param equation_nl: linear-less update step :param equation_l: linear update step :param init_t: initial time :param y0: initial 1D array :param solver_params: = ``{"atol"[absolute tolerance]: 1e-6, "rtol"[relative tolerance]: 1e-4, "exact_viscous_splitting"[if True, use exact viscous splitting]: False}`` :param dt_params: ``{dt0 = initial time step, dtmin, dtmax}`` .. py:method:: rk_step(fun: Callable[[float, numpy.ndarray], numpy.ndarray], t: float, y: numpy.ndarray, h: float, A: numpy.ndarray, B: numpy.ndarray, C: numpy.ndarray, K: numpy.ndarray) -> tuple[numpy.ndarray, numpy.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 [#]_. :param fun: Right-hand side of the system. :param t: Current time. :param y: Current state. :param h: Step to use. :param 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. :param B: Coefficients for combining RK stages for computing the final prediction. :param C: Coefficients for incrementing time for consecutive RK stages. The value for the first stage is always zero. :param 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 .. [#] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential Equations I: Nonstiff Problems", Sec. II.4. .. py:method:: rkdp(force_step: bool = False) -> tuple[numpy.ndarray, float, float] Performs one RKDP step .. py:method:: solve() -> None Main solving loop. The solver runs until an interruption is fired by the update. .. py:class:: ETD4RK(equation_nl: Callable[[float, numpy.ndarray], numpy.ndarray], equation_l: Callable[[float, numpy.ndarray], numpy.ndarray], save_step: Callable[[float, float, numpy.ndarray], bool], init_t: float, y0: numpy.ndarray, dt_params: dict[str, Optional[float]]) Bases: :py:obj:`_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 :param equation_nl: linear-less update step :param equation_l: linear update step :param init_t: initial time :param y0: initial 1D array :param dt_params: ``{dt0 = initial time step, dtmin, dtmax}`` .. py:method:: dopri_step(fun: Callable[[float, numpy.ndarray], numpy.ndarray], t: float, y: numpy.ndarray, h: float, A: numpy.ndarray, B: numpy.ndarray, C: numpy.ndarray, K: numpy.ndarray) -> tuple[numpy.ndarray, numpy.ndarray] :staticmethod: 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 [#]_. :param fun: Right-hand side of the system. :param t: Current time. :param y: Current state. :param h: Step to use. :param 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. :param B: Coefficients for combining RK stages for computing the final prediction. :param C: Coefficients for incrementing time for consecutive RK stages. The value for the first stage is always zero. :param 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 .. [#] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential Equations I: Nonstiff Problems", Sec. II.4. .. py:method:: rk() -> numpy.ndarray Performs one Rk step .. py:method:: solve() -> None Main solving loop. The solver runs until an interruption is fired by the update. .. py:class:: ETD35(equation_nl: Callable[[float, numpy.ndarray], numpy.ndarray], equation_l: Callable[[float, numpy.ndarray], numpy.ndarray], save_step: Callable[[float, float, numpy.ndarray], bool], init_t: float, solver_params: dict, y0: numpy.ndarray, dt_params: dict[str, Optional[float]]) Bases: :py:obj:`_CustIntegrator` Based on https://github.com/whalenpt/rkstiff .. warning:: Only supports diagonal (and constant-in-time) linear terms for now :param equation_nl: linear-less update step :param equation_l: linear update step :param init_t: initial time :param 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}`` :param y0: initial 1D array :param dt_params: ``{dt0 = initial time step, dtmin, dtmax}`` .. py:method:: rk() -> tuple[numpy.ndarray, float, float] Performs one Rk step :returns: tuple ``(y_new, h, h_new_suggested)`` .. py:method:: solve() -> None Main solving loop. The solver runs until an interruption is fired by the update. .. py:exception:: SolverInterruptedError Bases: :py:obj:`InterruptedError` custom error thrown to change grid size .. py:class:: Solver(D: int, l_params: dict[str, Union[float, bool]], fields_names: list[str], equation_nl: Callable[[float, pyloggrid.LogGrid.Grid.Grid, dict], dict], equation_l: Callable[[float, pyloggrid.LogGrid.Grid.Grid, dict], dict[str, numpy.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 :param D: dimension of physical space (usually 1, 2 or 3) :param l_params: parameters that define the grid parameter `̀`l`` as a dict with keys ``{"a", "b", "plastic"}``. See below for more detail. :param fields_names: the names of the scalar grid-shaped fields :param equation_nl: the function that performs the nonlinear update step. Returns the time derivative of the fields :param equation_l: implicit visocsity update step, called after the RK step has ended. Returns the new fields and their new time derivative :param k_min: minimum k of the grid. Default is defined by :class:`pyloggrid.LogGrid.Grid.Grid`. :param k0: whether there's a k0 mode :param simu_params: physical quantities relevant to the simulation, fixed for the whole simulation. :param 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 :class:`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 :math:`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 :math:`\sqrt{D}k_{min}` in :math:`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 :math:`\partial_t u=A(u)−b\cdot u`, then the linear part is :math:`−b`, and the nonlinear part is :math:`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 :func:`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. .. py:method:: _load_parameters(path: str) -> tuple[float, float] Loads simulation parameters from settings.json. Overwrites instance variables :param path: directory to load the settings from :returns: (last used time, last used timstep) .. py:method:: _save_step_all(t: float, dt: float, path: str) -> None Save the current step. Saves both fields and settings. :param t: current time :param dt: timestep :param path: folder to save to. Will override any existing file. .. py:method:: solve(save_path: str, initial_conditions: str | Callable[[dict[str, numpy.ndarray], pyloggrid.LogGrid.Grid.Grid, dict], dict[str, numpy.ndarray]] | tuple | list, solver_params: dict = None, init_t: float = 0, end_simulation=None, save_one_in: float = 1, update_gridsize_cb: Optional[Callable[[pyloggrid.LogGrid.Grid.Grid], Optional[int]]] = None, dt_params: dict[str, Optional[float]] = None, solver: Literal[ViscDopri, ETD4RK, ETD35] = 'ETD35') -> None Solve the solver's equation with the provided numerical parameters. :param save_path: where to save the fields and simulation parameters. Will backup if exists, will create otherwise :param 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. :param solver_params: params forwarded to the solver :param init_t: initial simulation time :param end_simulation: dict of thresholds to end the simulation: ``{t: max physical time, elapsed_time: max real time spent computing, step: max save step}`` :param save_one_in: save one step to file every X ``ode_step`` :param 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. :param dt_params: ``{dt0: initial timestep, dtmin, dtmax}`` :param solver: ``"ETD35"`` for :class:`ETD35` (default), ``"ETD4RK"`` for :class:`ETD4RK`, ``"ViscDopri"`` for :class:`ViscDopri` Details on a few parameters ########################### * If ``save_path="/"`` already exists on disk, the existing directory will be renamed to ``/bk__``. 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=(:str, :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.