4.4.4 pde.pdes.diffusion module

A simple diffusion equation.

class DiffusionPDE(diffusivity: float = 1, *, bc: BoundariesData | None = None, noise: float = 0, rng: np.random.Generator | None = None)[source]

Bases: SDEBase

A simple diffusion equation.

The mathematical definition is

\[\partial_t c = D \nabla^2 c\]

where \(c\) is a scalar field and \(D\) denotes the diffusivity.

Parameters:

  • diffusivity (float) – The diffusivity of the described species

  • bc – The boundary conditions applied to the field. Boundary conditions are generally given as a dictionary with one condition for each axis side. For periodic axes, only periodic boundary conditions are allowed (indicated by ‘periodic’ and ‘anti-periodic’). For non-periodic axes, different boundary conditions can be specified for the lower and upper end (using specific identifiers, like x- and y+). For instance, Dirichlet conditions enforcing a value NUM (specified by {‘value’: NUM}) and Neumann conditions enforcing the value DERIV for the derivative in the normal direction (specified by {‘derivative’: DERIV}) are supported. Note that the special value ‘auto_periodic_neumann’ imposes periodic boundary conditions for periodic axis and a vanishing derivative otherwise. More information can be found in the boundaries documentation.

  • noise (float) – Variance of the (additive) noise term

  • rng (Generator) – Random number generator (default: default_rng()) used for stochastic simulations. Note that this random number generator is only used for numpy functions, while compiled numba code uses the random number generator of numba. Moreover, in simulations using multiprocessing, setting the same generator in all processes might yield unintended correlations in the simulation results.

check_rhs_consistency(rhs_implementation: Callable, state: TField, t: float = 0, *, tol: float = 1e-07) → None

Checks an implementation of the right hand side versus the numpy variant.

Parameters:

  • rhs_implementation (callable) – The implementation of the numba variant that is to be checked.

  • state (FieldBase) – The state for which the evolution rates should be compared

  • t (float) – The associated time point

  • tol (float) – Acceptance tolerance. The check passes if the evolution rates differ by less then this value

determine_backend(state: TField, backend: str | pde.backends.base.BackendBase = 'auto', *, use_mpi: bool = False) → pde.backends.base.BackendBase

Returns backend that will be chosen automatically for this PDE.

Parameters:

  • state (FieldBase) – An example for the state from which the grid and other information can be extracted.

  • backend (str) – Information about which backend to choose. The special value auto tries various backends and returns one for which the evolution rate is implemented for this PDE.

  • use_mpi (bool) – Request backend with MPI support for parallel simulation.

Returns:

The backend used automatically

Return type:

str

evolution_rate(state: ScalarField, t: float = 0) → ScalarField[source]

Evaluate the right hand side of the PDE.

Parameters:

  • state (ScalarField) – The scalar field describing the concentration distribution

  • t (float) – The current time point

Returns:

Scalar field describing the evolution rate of the PDE

Return type:

ScalarField

property expression: str

the expression of the right hand side of this PDE

Type:

str

property is_sde: bool

flag indicating whether this is a stochastic differential equation

The SDEBase class supports additive Gaussian white noise, whose magnitude is controlled by the noise property. In this case, is_sde is True if self.noise != 0. Sub-classes might need to set is_sde explicitly to signal that they define a stochastic partial differential equation.

Type:

bool

make_evolution_rate(state: pde.fields.scalar.ScalarField, backend: pde.backends.base.BackendBase) → Callable[[TNativeArray, float], TNativeArray][source]

Create a compiled function evaluating the right hand side of the PDE.

Parameters:

  • state (ScalarField) – An example for the state defining the grid and data types

  • backend (str or BackendBase) – The backend used for numerical operations

Returns:

A function with signature (state_data, t), which can be called with an instance of the state data and time to obtain the associated evolution rate.

make_noise_variance(state: TField, *, backend: pde.backends.base.BackendBase[TNativeArray], ret_diff: bool = False) → Callable[[TNativeArray, float], TNativeArray | tuple[TNativeArray, TNativeArray]]

Make function that calculates noise variance.

Parameters:

  • state (FieldBase) – An example for the state from which the grid and other information can be extracted.

  • backend (str) – Determines the backend.

  • ret_diff (bool) – Determines whether only the noise variance or also its derivative with respect to the field at this position is returned.

Returns:

A function with signature (state_data, t) that either returns just the noise variance or also its derivative, depending on ret_diff.

make_pde_rhs(state: TField, backend: str | pde.backends.base.BackendBase = 'auto') → Callable[[TNativeArray, float], TNativeArray]

Return a function for evaluating the right hand side of the PDE.

Parameters:

  • state (FieldBase) – An example for the state from which the grid and other information can be extracted.

  • backend (str) – The backend that is used to create the function. The special value numpy uses the method evaluation_rate. Other backends are only available if make_evolution_rate is defined for the PDE. If this is the case, the special value auto selects the numba backend, otherwise it defaults to numpy.

Returns:

Function determining the right hand side of the PDE

Return type:

callable

make_post_step_hook(state: FieldBase, backend: str | pde.backends.base.BackendBase[TNativeArray] = 'numpy') → tuple[StepperHook, Any]

Returns a function that is called after each step.

This function receives three arguments: the current state as a numpy array, the current time point, and a numpy array that can store data for the hook function. The function must return the state data and the hook data, which it can both modify in place.

The hook can also be used to abort the simulation when a user-defined condition is met by raising StopIteration. Note that this interrupts the inner-most loop, so that some final information might be still reflect the values they assumed at the last tracker interrupt. Additional information (beside the current state) should be returned by the post_step_data. Note that raising StopIteration only works for some backends.

Example

The following code provides an example that creates a hook function that limits the state to a maximal value of 1 and keeps track of the total correction that is applied. This is achieved using post_step_data, which is initialized with the second value (0) returned by the method and incremented each time the hook is called.

def make_post_step_hook(self, state, backend):
    def post_step_hook(state_data, t, post_step_data):
        i = state_data > 1  # get violating entries
        overshoot = (state_data[i] - 1).sum()  # get total correction
        state_data[i] = 1  # limit data entries
        post_step_data += overshoot  # accumulate total correction
        return state_data, post_step_data

    return post_step_hook, 0.0  # hook function and initial value
Parameters:

  • state (FieldBase) – An example for the state from which the grid and other information can be extracted

  • backend (str) – Determines how the function is created (like ‘numpy’ and ‘numba’)

Returns:

The first entry is the function that implements the hook. The second

entry gives the initial data that is used as auxiliary data in the hook. This can be None if no data is used.

Return type:

tuple

solve(state: TField, t_range: TRangeType, dt: float | None = None, tracker: TrackerCollectionDataType = 'auto', *, backend: str | pde.backends.base.BackendBase = 'auto', solver: str | SolverBase = 'euler', ret_info: bool = False, **kwargs) → TField | tuple[TField | None, dict[str, Any]] | None

Solves the partial differential equation.

The method constructs a suitable solver (SolverBase) and controller (Controller) to advance the state over the temporal range specified by t_range. This method only exposes the most common functions, so explicit construction of these classes might offer more flexibility.

Parameters:

  • state (FieldBase) – The initial state (which also defines the spatial grid).

  • t_range (float or tuple) – Sets the time range for which the PDE is solved. This should typically be a tuple of two numbers, (t_start, t_end), specifying the initial and final time of the simulation. If only a single value is given, it is interpreted as t_end and the time range is (0, t_end).

  • dt (float) – Time step of the chosen stepping scheme. If None, the solver chooses a default value when constructing the stepping function. If adaptive stepping is enabled (e.g., supported by EulerSolver), dt sets the initial time step.

  • tracker – Defines trackers that process the state of the simulation at specified times. A tracker is either an instance of TrackerBase or a string identifying a tracker (possible identifiers can be obtained by calling registered_trackers()). Multiple trackers can be specified as a list. The default value auto checks the state for consistency (tracker ‘consistency’) and displays a progress bar (tracker ‘progress’) when tqdm is installed. More general trackers are defined in trackers, where all options are explained in detail. In particular, the time points where the tracker analyzes data can be chosen when creating a tracker object explicitly.

  • backend (str) – Determines how the function is created. Accepted values are ‘numpy’ and ‘numba’. Alternatively, ‘auto’ lets the code pick the optimal backend.

  • solver (SolverBase or str) – Specifies the persistent numerical strategy used for solving the differential equation. This can either be a solver factory or a descriptive name like ‘explicit’ or ‘scipy’. The valid names are given by pde.solvers.registered_solvers(). Details of the solver classes and additional features (like adaptive stepping) are explained in solvers.

  • ret_info (bool) – Flag determining whether diagnostic information about the solver and stepping process should be returned. Note that the same information is also available as the diagnostics attribute.

  • **kwargs – Additional keyword arguments are forwarded to the solver chosen with the solver argument. In particular, adaptive stepping can often be enabled using adaptive=True.

Returns:

The state at the final time point. If ret_info == True, a tuple with the final state and a dictionary with additional information is returned. Note that None instead of a field is returned in multiprocessing simulations if the current node is not the main MPI node.

Return type:

FieldBase