What is the easiest way of topping field values within the framework of the SIR example?

[HannoSpreeuw]

Looking at the SIR model example what would be the easiest way to limit the values of one of the fields?

By means of evolution_rate the three fields s, i and r evolve.
What if one of those three fields cannot exceed some value because of some constraint?

Say that i cannot exceed some value alpha. So if i became larger than alpha it would have to be set equal to alpha.

How would you implement that?

[HannoSpreeuw]

I guess one could use to_scalar for a ScalarField, like this:

j = i.to_scalar(lambda f: numpy.fmin(f,alpha))

[david-zwicker]

I'm not sure what kind of PDE your are solving, but directly limiting the field is not something that is typically possible in PDEs. Nonetheless, there is a workaround in py-pde, which only applies to explicit solvers. Here is a simplified example, which modifies the field after every step to set all values above 1 to 1.

from pde import UnitGrid, MemoryStorage, PDEBase, ScalarField, plot_kymograph


class LinearPDE(PDEBase):

    def evolution_rate(self, state, t=0):
        return pde.ScalarField(state.grid, 1)
    
    def make_modify_after_step(self, state):
        def limit_state(state_data):
            state_data[state_data > 1] = 1
            return 0
        return limit_state

# initialize the equation and the space
grid = UnitGrid([32], periodic=True)
state = ScalarField.random_uniform(grid)

# solve the equation and store the trajectory
storage = MemoryStorage()
eq = LinearPDE()
eq.solve(state, t_range=2, dt=0.01, method='explicit', tracker=storage.tracker(0.1))

# plot the trajectory as a space-time plot
plot_kymograph(storage)

[HannoSpreeuw]

Thanks, I am solving PDEs in the framework of the SIR example, so using a custom class.

They are about the time evolution of compositions and concentrations of molecules in sediments,
governed by chemical self-organisation and stochastic resonance.
Under- and oversaturation factors in these PDEs are some non-integer power of certain variables.
The base of these exponentiations may not become negative. That is why I need to set an upper limit of 1
to product_of_two_scalar_fields that appears as 1 - product_of_two_scalar_fields in the base of the exponentiations.

Is the SIR example an example of an explicit solver? I am not sure how to interpret "explicit".
Anyway, I can still run after adding method="explicit", scheme = "rk" as arguments to my eq.solve(.....).
This seems to give the same solutions, at first sight.

Your fix seems very close to my proposed j = i.to_scalar(lambda f: numpy.fmin(f,alpha)), which passes without error in my code.
I understand that that does not mean that it will work correctly.

[HannoSpreeuw]

Perhaps I should have added that this is not about any field that I am solving for but for the under- and oversaturation factors which cover the same grid. So these are actually "helper fields" and have to be computed for each time step in order to derive the fields I am solving for.

[david-zwicker]

I see, if these are not the dynamic fields, the simplest way might be to use numpy.clip to limit the field. A simple example in py-pde could read

import numpy as np
import pde

grid = pde.UnitGrid([32])
field = pde.ScalarField.random_normal(grid)
np.clip(field.data, 0, 1, out=field.data)
field.plot()

Here, the third line limits the field to the interval [0, 1].