Inactive cells in Porepy?

[rafaelmarch3]

Hello all,

I have one quick question, does PorePy support "active cells", as usual in typical reservoir simulators? I have cartesian grids that are relatively large for the typical scipy sparse solvers, but in some cases many cells do not contribute to fluid flow (k~=0). I would like to assemble a smaller system that leave out these cells. In reservoir simulators you usually do this with active/inactive cells. Any way to achieve this in PorePy? Also, any suggestions of sparse solvers that go a bit beyond spsolve? It's in my todo list to try pyamg and PETSC for Python...

Thank you!
Rafael March.

[keileg]

We have not implemented an active cell concept; the kind of problems we have been interested in has not called for such tricks. However, If you have the indices of the cells to be eliminated, making a restriction matrix is not difficult; the code will roughly be

A = matrix with too many rows and columns
inds = indices_of_active_cells  # assumed to be a numpy array
sz = inds.size
R = sps.coo_matrix((np.ones(sz), (np.arange(sz), inds), shape=(sz, A.shape[1]))
A_smaller = R.T * A * R
b_smaller = R * b

x_smaller = spsolve(A_smaller, b_smaller)
x_full = R.T * x_smaller  # This will have zeros for the inactive cells

There surely is a typo or two in the snippet, but it should not be far from working.

As for the solvers, some suggestions:

• spsolve uses superlu as its backend (direct) solver. If you use conda for installation, you can try pardiso instead, it should be much faster if it is applicable. See this project for install instructions.
  Side note, if conda is not an option: Consider umfpack, accessible through scikit-umfpack. This is also much faster than superlu.
• I have played with both pyamg and petsc4py; getting started with these should not be difficult. How well it will work depends on grid quality, discretization etc. If you have fractures present, I am not sure if using amg as a black box is the best approach, but I can't say which adjustments (transformations etc.) of the linear systems are necessary for better behavior. This is something I have wanted to look into for a long time now, but there is always something else.

I hope this is of some help

Cheers,
Eirik

[rafaelmarch3]

Very helpful Eirik, thanks for that! Indeed, I have access to the assembled matrix, so this code will be useful to restrict the cells.

As for the solvers, I'll give it a go. I'm not playing with fractures at the moment, just testing some concepts around upscaling for digital rocks with strong diagenetic features, like microcracks, vugs...but it all translates into a segmented image which you can assign porosity/permeability trends.

Thanks again!
Rafael.