Wrong indices assumed from scipy (after pypardiso install)

[rbe051]

I came over a weird bug. Here is a minimal example:

import numpy as np
import porepy as pp

# Define a rectangular domain in terms of range in the two dimensions
bounding_box = {'xmin':0, 'xmax':1, 'ymin':0, 'ymax': 1}
domain = pp.Domain(bounding_box=bounding_box)

# Define each individual fracture, collect into a list.
frac1 = pp.LineFracture(np.array([[0.5, 0.5],
                                  [0, 1]]))
fractures = [frac1]

# Define a fracture network in 2d
network_2d = pp.create_fracture_network(fractures, domain)

# Set overall target cell size and target cell size close to the fracture.
mesh_args: dict[str, float] = {"cell_size": 0.5}

# Generate a mixed-dimensional grid
mdg = pp.create_mdg("cartesian", mesh_args, network_2d)
pp.set_local_coordinate_projections(mdg)

Which may give out the error (I will get back to may):

  File c:\users\runar\repositories\porepy\src\porepy\utils\tangential_normal_projection.py:284 in _invert_3d_matrix
    M_inv[:, :, i] = np.linalg.inv(M[:, :, i])

  File ~\anaconda3\envs\pp\Lib\site-packages\numpy\linalg\linalg.py:561 in inv
    ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)

  File ~\anaconda3\envs\pp\Lib\site-packages\numpy\linalg\linalg.py:112 in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix

I think I have traced down the error to line 857 in mortar_grid.py:

secondary_f, primary_f, data = sps.find(primary_secondary)

which gives out:
secondary_f = [0, 0, 1, 1]
primary_f = [ 4, 13, 1, 12]

while the code following the line assumes the ordering:
secondary_f = [1, 0, 1, 0]
primary_f = [ 1, 4, 12, 13]

This started to happen after I install pypardiso. I then created an identical conda enviorment, except for installing pypardiso and everything works just fine. I have no idea why pypardiso affects anything in this case. The version number of both scipy and numpy is the same in both environments.

[keileg]

Thanks for the report, and not least for digging this deep into things! I will look into it, but in the meantime, could you check the versions of scipy in the two environments? Specifically, do you have version 1.11 in the environment where you have trouble?

[rbe051]

I have scipy 1.10 in both environments.

[OmarDuran]

If possible could test your example on the following branch and If the issue has been resolved, could you please let us know?

[IvarStefansson]

Can this be closed, @OmarDuran and @rbe051?

[rbe051]

Yes, it think so