pmgbergen / porepy

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

Limit-cases-of-ad-functions #1169

Closed IvarStefansson closed 5 months ago

IvarStefansson commented 6 months ago

Proposed changes

This PR introduces two minor modifications of the ad.functions:

  1. Uses first argument in the case of equality of the arguments of the max function. This impact the tangential contact mechanics and the aperture calculation.
  2. Assigns unitary jacobian values for l2_norm of an argument with val=0. Thus, derivative information is nonzero and we avoid special cases where a block of the linear system vanishes.

Types of changes

What types of changes does this PR introduce to PorePy? Put an x in the boxes that apply.

Checklist

Put an x in the boxes that apply or explain briefly why the box is not relevant.

keileg commented 6 months ago

@mariusnevland will have a look at the impact of the suggested changes on some of his setups.

keileg commented 6 months ago

@mariusnevland: Will you have time to have a look at this in the nearest future, or should we proceed without your input?

mariusnevland commented 6 months ago

@mariusnevland: Will you have time to have a look at this in the nearest future, or should we proceed without your input?

I have taken a quick look at it already. I ran the same simulations that I presented on the previous MaPSI meeting (where the dilation angle and the angle between two intersecting fractures are varied), and the performance was generally much better with the changes in this PR; the Newton solver converged for all combinations. I have not tested for any other setups though.

I will be quite busy the next couple of weeks, so I might not be able to look further into it for now.

One comment: I am assuming that the change in the max-function is primarily felt by the b_p variable in the tangential equation? I believe that, from a theoretical standpoint, the max-function in b_p is redundant, because the distinction between open and closed cells is already being handled by the characteristic function. In other words, changing the definition to b_p=self.friction_bound(subdomains) should give the same solution. However, I can confirm through some experimentation that this change will affect the behavior of the solver (not sure yet if in a good or bad way).

IvarStefansson commented 6 months ago

@jwboth and I have revisited the usage and confirmed that the changes make sense from a "heuristic theoretical" point of view. We will keep probing the implementation a bit (approx. a week?) before merging. Marking as draft in the meantime.

IvarStefansson commented 5 months ago

@jwboth, have you had a chance to test this?

jwboth commented 5 months ago

@jwboth, have you had a chance to test this?

I integrated the branch in my recent numerical studies. But I am not really able to conclude. I did not really do a proper study, comparing both side by side for many cases. I only quickly compared the variants on this branch and develop for one single test case. And I could not see a real difference in the performance.

It will be a bit difficult to do a proper study this week. If you do not want to wait, I still feel like the proposed changes are meaningful, and would be OK with integrating them. Otherwise, I make a proper study in the beginning of next week.

IvarStefansson commented 5 months ago

Quote reply

@jwboth, have you had a chance to test this?

I integrated the branch in my recent numerical studies. But I am not really able to conclude. I did not really do a proper study, comparing both side by side for many cases. I only quickly compared the variants on this branch and develop for one single test case. And I could not see a real difference in the performance.

It will be a bit difficult to do a proper study this week. If you do not want to wait, I still feel like the proposed changes are meaningful, and would be OK with integrating them. Otherwise, I make a proper study in the beginning of next week.

Any updates?