Singular matrix error in THM simulations with multiple intersecting fractures

[adyngo]

Hi maintainers, and thank you for your work on this project.

I'm working on a THM coupling in a 3D DFM model, where I simulate water injection from a well to study the resulting rock responses. I'm encountering convergence failures with the Newton solver, and I'd like to ask if this is expected behavior for this type of problem, or if it suggests an issue with my configuration or implementation.

Problem Description

I'm encountering singular matrix errors when running THM simulations with multiple intersecting fractures. The simulation fails during the first or earlier nonlinear iteration with the error:

Matrix A is singular, because it contains empty row(s)

My simulations converge with 3 or fewer fractures, but consistently fail with 4 or more. I've confirmed that during the nonlinear iterations, the residual norm grows significantly, and the matrix becomes singular.I would like to understand if encountering a singular matrix in my situation is expected behavior, or if this points to a configuration issue or code error that I can resolve.

Are there specific settings that must be adjusted when working with such complex fracture networks? I would greatly appreciate any guidance, best practices, or workarounds for handling multiple intersecting fractures in these THM simulations.

If you need any information, please let me know.
Thank you for your help!

Below, I have attached an example of my situation:

Simulation Setup

(I've tried multiple similar simulations with varying configurations, but the following represents a typical example that demonstrates the issue.)

Domain

• Size: 300 m × 300 m × 300 m
• Material : Granite

Fracture Network

• Number of fractures: 4-20 (tested with various subsets)
• Fracture geometry: Hexagonal polygons (6 vertices each)
• Fracture sizes:
  • Fracture area: 6,500-60,000 m² per fracture
  • Edge lengths: 50-150 m (mean: ~100 m)
• Fracture intersections:
  • Fractures intersect both each other and some of them are intersecting domain boundaries
  • Injection wells intersect fractures at discrete points

Mesh Parameters

• Matrix cell size: 30-50 m
• Fracture cell size: 20-35 m
• Boundary cell size: 30-50 m
• Minimum cell size: 6-10 m

Boundary Conditions

• Fluid (Darcy flux):
  • Matrix (3D): Dirichlet (pore pressure) on all external boundaries
  • Fractures (2D): Dirichlet (pore pressure) at fracture-domain intersections
• Thermal: Dirichlet (reservoir temperature) on all boundaries(20℃)
• Mechanical: Dirichlet (zero displacement) at bottom, Neumann (far-field stress traction) on other boundaries

Initial Conditions

• Fluid: Hydrostatic pore pressure with overpressure
• Thermal: Uniform reservoir temperature (20℃)
• Mechanical: Initial stress state set to far-field stress values (anisotropic, constant per axis, 15 ~ 25 MPa)

Sources

• Injection well: 1 well intersecting 1-3 fractures
• Injection rate: 0.1-1.0 kg/s (water at 15°C)

Solver Configuration

• Nonlinear solver: Newton solver (also tried with linear search)
• Convergence tolerance: 1e-3 (increment), 1e-3 (residual)

[IvarStefansson]

Hi,
First of all: I'm impressed by your brave attempt at using our tool for what are really quite complicated setups! And thanks for your thorough description of your issues. I will not pretend to be highly surprised you run into convergence issues, we have encountered them aplenty ourselves! Here are a few things you might want to look into in order to diagnose the problem further.

1. Can you identify some qualitative change between your 3 and 4 fracture setups? E.g. fractures at external boundaries, different kinds of intersections.
2. Try to run with a poromechanics model instead of thermoporomechanics. If that works, it would indicate that the problem is in some sense "nonlinearly difficult".
3. Did you modify your model apart from initial conditions, boundary conditions and geometry?
   I'll be happy to provide more reflections and suggestions once you have considered these three.
   Ivar

[adyngo]

Hi Ivar,

Thank you for the quick response and for the helpful suggestions.

Here's what I found based on your questions:

1. I re-ran the program with several network models, and it seems the key factor is the nature of the intersections, not just the total number of fractures. (The injection well isn't displayed, but it runs through the center of the domain and intersects the horizontal fractures in all cases.)

When I use up to 5 layered, non-intersecting fractures (like the first image), the simulation converges.

When I increase this to 6 or more layered fractures (second image), the calculation did not converge within 20 iterations. (I also tried small time steps, down to 1, with no success).

When I use networks with multiple intersecting fractures (like the last two images), I consistently get the singular matrix error.

2 Poromechanics Model: I haven't had time to test the pure poromechanics model yet, but I will try to this week and report back. However, I did check the MassAndEnergyBalance model, and it converged successfully, even with the complex, intersecting networks.

3 Model Modifications: No, the only modification I made was adding pp.constitutive_laws.CubicLawPermeability to the basic model.

Thanks again for your guidance!