Dynamic time step with scheduling causes dividing by zero

[Yuriyzabegaev]

Dynamic time stepping is partly implemented out of the box, but can trigger a hard-to-catch error related to the time scheduling. The time step is adjusted in TimeManager.compute_time_step, which is not called anywhere, so to make it work we must do something like:

from porepy.models.poromechanics import Poromechanics
import porepy as pp
import numpy as np

class MyModel(Poromechanics):

    def after_nonlinear_convergence(
        self, solution: np.ndarray, errors: float, iteration_counter: int
    ) -> None:
        super().after_nonlinear_convergence(solution, errors, iteration_counter)
        self.time_manager.compute_time_step(iteration_counter, recompute_solution=False)

    def after_nonlinear_failure(
        self, solution: np.ndarray, errors: float, iteration_counter: int
    ) -> None:
        self.time_manager.compute_time_step(recompute_solution=True)

Now, we can set some schedule we want to match exactly and the limits for the time step, e.g., time_manager = pp.TimeManager(dt_init=1, dt_min_max=(0.1, 1), schedule=[0, 5, 10], constant_dt=False). However, at the time step 5, the scheduler will notice that we're approaching the scheduled time=5, and decrease the time step to match it exactly, see this. In our case, it will set dt to 0, which will lead to hardly traceable errors when parsing the AD equations. The full example to reproduce this is:

from porepy.models.poromechanics import Poromechanics
import porepy as pp
import numpy as np


class MyModel(Poromechanics):

    def after_nonlinear_convergence(
        self, solution: np.ndarray, errors: float, iteration_counter: int
    ) -> None:
        super().after_nonlinear_convergence(solution, errors, iteration_counter)
        self.time_manager.compute_time_step(iteration_counter, recompute_solution=False)

    def after_nonlinear_failure(
        self, solution: np.ndarray, errors: float, iteration_counter: int
    ) -> None:
        self.time_manager.compute_time_step(recompute_solution=True)


time_manager = pp.TimeManager(
    dt_init=1, dt_min_max=(0.1, 1), schedule=[0, 5, 10], constant_dt=False
)
model = MyModel(
    {
        "time_manager": time_manager,
    }
)

pp.run_time_dependent_model(
    model,
    {
        "progressbars": True,
    },
)

I realize that using the dynamic time stepping is something not fully supported now, but if you're interested, I can implement a fix for it.

[IvarStefansson]

Let's discuss this issue in person, @Yuriyzabegaev.