pmgbergen / porepy

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

Negative defined mechanics submatrix #849

Closed Yuriyzabegaev closed 1 year ago

Yuriyzabegaev commented 1 year ago

Hei,

I noticed that the matrix arising from the mechanics problem in PorePy is negatively defined. It has a dominant negative diagonal and gives a negative scalar when multiplied by a random vector from left and right. That's alright if we solve the mechanics block separately, but it might lead to suboptimal solver performance if we couple mechanics with something else. For example, the single-phase flow problem matrix is positively defined.

Providing the code to reproduce:

import matplotlib.pyplot as plt
import porepy as pp

model = pp.poromechanics.Poromechanics()
pp.run_stationary_model(model, {})
mat, rhs = model.linear_system
plt.matshow(mat.A)
plt.colorbar()

Result: image

jhabriel commented 1 year ago

Interesting observation @Yuriyzabegaev. Do you get the same results if you use pp.run_time_dependent_model() instead?

Yuriyzabegaev commented 1 year ago

@jhabriel yes the result is the same.

Yuriyzabegaev commented 1 year ago

@IvarStefansson if we test the incompressible single-phase flow setup (elliptic problem), the diagonal values remain positive. As you suggested, it must be a scalar and vector divergence inconsistency.

IvarStefansson commented 1 year ago

Thanks for checking! @keileg, could it be that the sign convention is different between mpfa.flux and mpsa.stress?