pmgbergen / porepy

Python Simulation Tool for Fractured and Deformable Porous Media
GNU General Public License v3.0
251 stars 88 forks source link

Dynamic time step with scheduling causes dividing by zero #1152

Closed Yuriyzabegaev closed 5 months ago

Yuriyzabegaev commented 7 months ago

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 commented 6 months ago

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