make_post_step_hook seems to be not-working

[VarvaraGordeeva]

Hello there,
I tried to use the hook function in my calculations. I found only two exapmles and description that the method "post_step_hook" will start after every time step. In both examples, I see solution which I need: if some values of the solved field are higher that one, they become equal to one (in my case it is mass fraction, which must be less or equal one but becomes slightly higher because of calculation error). My case is a little more complex: my "state" is a FieldCollection, and the mass fraction is one of the fields. I tried few variants (under "reaction" I mean that the mass fraction has no values higher than one):

1. totally copy the example - no reaction
   def make_post_step_hook(self, state):
   """Create a hook function that is called after every time step."""
   def post_step_hook(state_data, t, post_step_data):
   """Limit state 1 and abort when standard deviation exceeds 1."""
   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
   if post_step_data > 400:
   # Abort simulation when correction exceeds 400
   # Note that the post_step_data of the previous step will be returned.
   raise StopIteration
   return post_step_hook, 0.0 # hook function and initial value for data
2. modify the example: state_data -> state (otherwise "state" in "def make_post_step_hook(self, state):" is colored as not used) - no reaction
3. modify the example: state_data -> state[0] (the mass fraction is the first field of the FieldCollection state) - no reaction
   4,5,6) used numpy.clip(state...,0,1) with the same variants: state_data, state, state[0] - no reaction
   I see two variants here: either I use/write the class/method in wrong way, or the method is not working.

[david-zwicker]

What exactly are you asking? It would be helpful to have a minimal self-contained example with a clear description of the expected behavior.

[VarvaraGordeeva]

A plot of mass fraction vs coordinate is attached. As you can see, in some region it slightly higher than 1 (zoom in black rectangular). Obviously, it is result of the calculation error, and this error persists and affect for the next time steps. I would like to eliminate this error at every time step of calculation. The stated method should do it but doesn't. The following is piece of my code:

class evapLyth(PDEBase):
    def __init__(self, tt=1, bc = {"neumann"}):
        super().__init__()
        self.bc = bc
        self.tt = tt

    def make_post_step_hook(self, state):
        def post_step_hook(state, t, post_step_data):
            i = state > 1  # get violating entries
            state[i] = 1  # limit data entries
        return post_step_hook, 0.0  # hook function and initial value for data
    
    def evolution_rate(self, state, t: float = 0):
        Cc,Hhrho = state # mass fraction, film thickness * density
        { a lot of equations }
        return FieldCollection([Cc_t,Hh_t])

## initialize the equation and the space
grid = CartesianGrid([[-Xx/h0,Xx/h0]], [Nn])
Cc1 = ScalarField(grid,data = c0)
Hh1 = ScalarField(grid,data = rho)
state = FieldCollection(fields=[Cc1,Hh1])
eq = evapLyth()

# # solution
storage = MemoryStorage()
result = eq.solve(state, t_range, dt, tracker=storage.tracker(ttrack))
{ outputs and vizualisation }

[david-zwicker]

Thank you for the details! However, the provided code does not yet run since crucial details are missing (particularly the definition of the PDE). Could you please provide a self-contained example, so I can re-create the problem?

[VarvaraGordeeva]

I made a simplification of my problem. It has no physical sence but repeats the thing which I mean: "mass fraction" grows above 1 despite of make_post_step_hook method which should cut the over-values to 1. The plots show two last states of Cc, and I can see that previous state wasn't cut. The following is the code (pde version is 0.21.0 if it is important):

from pde import FieldCollection, PDEBase, ScalarField, CartesianGrid,MemoryStorage
import numpy as np
import matplotlib.pyplot as plt

# coordinates
Xx = 6   # d/less, half of the field
Nn = 600 # coordinate nodes
Nt = 100 # time steps
xx = np.linspace(-Xx,Xx,Nn)

## evaporation
Aa = 0.01
Bb = 1
Jj = -Aa*xx**2+Bb
# plt.plot(xx,Jj)
# plt.show()

## consts
cg = 0.7
c0 = 0.4
h0 = 1

class evapLyth(PDEBase):

    def __init__(self, tt=1, bc = {"neumann"}):
        super().__init__()
        self.bc = bc
        self.tt = tt

    def make_post_step_hook(self, state):
        """Create a hook function that is called after every time step."""

        def post_step_hook(state_data, t, post_step_data):
            """Limit state 1 and abort when standard deviation exceeds 1."""
            # i = state > 1  # get violating entries
            # state[i] = 1  # limit data entries
            np.clip(state_data,0,1)

            '''
            # overshoot = (state_data[i] - 1).sum()  # get total correction
            # post_step_data += overshoot  # accumulate total correction
            # if post_step_data > 400:
            #     # Abort simulation when correction exceeds 400
            #     # Note that the `post_step_data` of the previous step will be returned.
            #     raise StopIteration
            '''

        return post_step_hook, 0.0  # hook function and initial value for data
    
    def evolution_rate(self, state, t: float = 0):
        
        assert state.grid.dim == 1
        Cc,Hhrho = state # mass fraction, film thickness * density
        c_x = Cc.gradient("neumann")[0]
        c1_x = (Hhrho*c_x).gradient("neumann")[0]
        hvc = np.heaviside(cg-Cc,0)
        Cc_t = 1/Hhrho*c1_x*hvc+Jj*Cc/Hhrho
        Hh_t = -Jj-(Hhrho).gradient("neumann")[0]
        
        return FieldCollection([Cc_t,Hh_t])

## Initialization
grid = CartesianGrid([[-Xx,Xx]], [Nn])
# print("grid size ",grid.shape)
Cc1 = ScalarField(grid,data = c0)
Hh1 = ScalarField(grid,data = h0)
state = FieldCollection(fields=[Cc1,Hh1])
assert state.grid.dim ==1
eq = evapLyth()

# # solution
storage = MemoryStorage()
t_range = 0.11
dt = t_range/Nt
ttrack = t_range/5
result = eq.solve(state, t_range, dt, tracker=storage.tracker(ttrack))
storLen = len(storage.data)

## visualization
plt.figure(figsize=(5,6))

# # Mass part
plt.subplot(211)
plt.ylabel("$c$")
plt.grid()
cVorLetz = storage.data[-2][0]
cLetz = storage.data[-1][0]
plt.plot(xx,cVorLetz)
plt.plot(xx,cLetz,':')

# Thickness
plt.subplot(212)
plt.ylabel("$h$, $\mu$m")
plt.grid()
plt.xlabel("$x$, mm")
for i in range(len(storage.data)):
    h1 = storage.data[i][1]
    plt.plot(xx,h1)

plt.tight_layout()
plt.show()

[david-zwicker]

I had a brief look and can reproduce the problem. The problem could be fixed on my side by replacing np.clip(state_data,0,1) by np.clip(state_data,0,1,out=state_data) so state_data is replaced in-place. However, note that py-pde version 0.21.0 is severely outdated, so if you run into issues, please first update to a more up-to-date version.

[VarvaraGordeeva]

Thank you very much for your help! пт, 14 мар. 2025 г. в 18:06, David Zwicker ***@***.***>:

…

I had a brief look and can reproduce the problem. The problem could be fixed on my side by replacing np.clip(state_data,0,1) by np.clip(state_data,0,1,out=state_data) so state_data is replaced in-place. However, note that py-pde version 0.21.0 is severely outdated, so if you run into issues, please first update to a more up-to-date version. — Reply to this email directly, view it on GitHub <#675 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BNILDLQPTRQOQO3Z2YHMPZT2ULH6HAVCNFSM6AAAAABY5EWJB2VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTENJQGA4DGOA> . You are receiving this because you authored the thread.Message ID: ***@***.***>