Adding effective tangential/normal permeability to fractures

[Cirdans-Home]

Hi,
I've been trying to follow the tutorials and the examples to reproduce the Validation 5.1 test from the paper "PorePy: an open-source software for simulation of multiphysics processes in fractured porous media".

I think I added both the 2D network fractures and the boundary conditions, but I do not understand how to set the different permeability to the fractures, i.e., what is said about "the fractures identified in red (1, 2, 3, 6, 7, 8, 9, 10) as high permeable with K_i = 1 and κ_i = 2 × 10^8, whereas the fractures identified in blue (4, 5) act as low permeable barriers, with K_i = 1 × 10^−8 and κ_i = 2. The 0d-intersections between fractures are divided into two classes, if both fractures are high permeable we assign κi = 2 × 10^4, otherwise, we set κ_i = 4/(10^4 + 10^−4)."
manyfractures2D.txt

Any help in this regard would be greatly appreciated.

[IvarStefansson]

Hi,
I admit this may not be the most intuitive part of the code. First of all, the tangential (in-plane) permeability is simply referred to as "permeability" in PorePy, whereas the normal (out-of-plane) permeability is "normal_permeability". To set different values for the former for different subdomains, you'll have to override the "permeability" method of your model (the default is defined inside constitutive_laws.py). The latter can be modified by passing a "normal_permeability" to a SolidConstants object. Note that internally, (meth. interface_darcy_flux_equation) the passed value will be divided by the distance between fracture centre and fracture wall, i.e., \kappa_i=normal_permeability/(a/2), with a denoting aperture.
Hope this helps, and don't hesitate to ask if further clarification is needed.

[Cirdans-Home]

Dear Ivar,

Thank you very much for the answer. I went and took a look at the code in the constituve_laws.py file. The ConstantPermeability class has two inside methods "permeability" and "normal_permeability", they seem to loop between the different domains to select the constants, To perform the overload you are mentioning do I have to pass the constants as;

fluid_constants = pp.FluidConstants({"viscosity": 1,
	"permeability" : np.array([1,1,1,1e-8,1e-8,1,1,1,1,1]),
	"normal_permeability" : np.array([2e+8,2e+8,2e+8,2,2,2e+8,2e+8,2e+8,2e+8,2e+8])})
material_constants = {"fluid": fluid_constants}
params = {"material_constants": material_constants}
model = SinglePhaseFlowFractureSource(params)

with the numpy array ordered as I have ordered the fractures? (I'll have to multiply the entries of the normal_permeability array by the correct value to have the value used in the benchmark). What about the 0th-order intersections? Do I add them in the "normal_permeability" entry in between the corresponding fractures?

If I run the code in this way, I also get the warning:

UserWarning: Is it really meaningful to ask for signs of a one sided mortar grid?

so I'm probably doing something silly somewhere.

[Cirdans-Home]

Thank you for your patience. I redid the calculation considering what you told me and arranging the order of the points in the fractures. Now I feel like I have the right number of pieces. The thing I can't get to work though is passing normal permeabilities, the SolidConstants object only accepts a scalar as a dictionary argument associated with the "normal_permeability" keyword, that is,

solid_constants = pp.SolidConstants({"normal_permeability" : np.array([2e+8,2e+8,2e+8,2,2,2e+8,2e+8,2e+8,2e+8,2e+8])})

returns an error of

TypeError: only length-1 arrays can be converted to Python scalars

The normal_permeability method seems to be called just once in the code, how do I make it return to normal permeability for each different interface? I guess I also need to overload the normal_permeability method somehow.

[keileg]

Ivar is rather busy today and tomorrow, so although he has the best feeling for how to set up such models, I will do my best to give a fair answer: You need first to feed the values of the normal permeability to the model, and then assign the right value to the different interfaces. The SolidConstant class is not really made to assign different values to different domains, so I believe the best approach is to override the normal_permeability method.

As you have noted, normal_permeability is called once, and then with all interfaces (MortarGrids) as an argument. You need a way to identify individual interfaces, and this is most readily done via the subdomain grids of the neighboring fractures: These have an attribute frac_num which reflects the order in which you feed them to the mesh generation (note: you need to check if this is 0- of 1-offset, and I also recommend that you verify that the order is indeed preserved, if not, we will consider it a bug and will appreciate a report).

From what I recall, the example you are looking at has about 10 fractures, some of them intersecting, and some of the fractures having low normal permeability. An overridden method to set the normal permeability would look something like:

def normal_permeability(self, interfaces: list[pp.MortarGrid]) -> pp.ad.Operator:
    for intf in interfaces:
         # Get hold of the fracture subdomain
         sd_high, sd_low = self.mdg.interface_to_subdomain_pair(intf)
        if intf.dim == 1:
           if sd_low in list_of_fracture_indexes_with_high_permeability:
                 return pp.ad.Scalar(HIGH_PERM_VALUE)
           else:
                 return pp.ad.Scalar(LOW_PERM_VALUE)
       else:  # 0d-1d coupling
            # Not sure exactly how the problem is defined here, but you can identify the intersection by the 
            # fracture, which is now the higher-dimensional neighbor, like
            if sd_high.frac_num == SOME_VALUE:
                return pp.ad.Scalar(??)

The numerical values can be hardcoded in the implementation, as can the list of fracture indices that should be assigned high values. While simple, this is perhaps not best practice (though while debugging, everything is permissible), but you can also pass arbitrary arguments to the model parameter dictionary (the same that you use to send in your SolidConstants object), and access it inside the model by self.params['high_normal_permeability_value'] or similar.

Let me know if this makes sense. I will try to follow this up, but I may not be able to respond until the middle of next week.

[Cirdans-Home]

Dear Keileg,

Thank you very much for the answer, the idea seems clear. What I find strange is that the normal_permeability routine is only called once in the code. So if I do a for loop with returns in it don't I get an operator that contains the first normal permeability value I end up in due to the if? Shouldn't the result instead be a sort of piecewise constant operator on the various interfaces?

Namely, it seems I'm doing something like:

def function():
    for intf in range(1,10):
        return intf

and then calling function(), which returns only a 1, and that's it. It could well be a misunderstanding due to my poor knowledge of Python.

Regarding the sd.frac_num the indexing starts from 0 and the order is the one in which I passed them to the mesher, I printed some information from inside the permeability function and things seem fine.

[keileg]

Sorry for the late answer, I was out of office for a few days. You are right, the solution I sketched out is not complete, in that the returned quantity should be a single array valid for all interfaces, and actually for all cells on each of the interfacses
This must be built up as you iterate over the interfaces. Trying again, it would look like this:

def normal_permeability(self, interfaces: list[pp.MortarGrid]) -> pp.ad.Operator:

    # temporary storage for the values
    tmp_array = []

    for intf in interfaces:
        # Make a numpy array of ones, one element per cell in the interface grid
        e = np.ones(intf.num_cells) 
        # Get hold of the fracture subdomain
        sd_high, sd_low = self.mdg.interface_to_subdomain_pair(intf)
        if intf.dim == 1:
           if sd_low in list_of_fracture_indexes_with_high_permeability:
                 # Scale the array of ones with the right permeability, append it to the list
                 tmp_array.append(HIGH_PERM_VALUE * e)
           else:
                 tmp_array.append(LOW_PERM_VALUE * e)
       else:  # 0d-1d coupling
            pass  # Insert similar code here

    # tmp_array is a list of numpy arrays, but we need to convert it to a single array of type pp.ad.Operator:
    # First convert to a single numpy array, then wrap it into the right format.
    single_array = pp.wrap_as_ad_array(np.hstack(tmp_array), name='normal_permeability_array')
    return single_array

If you want to check that the array is indeed of the right size or similar, you can inspect the result of np.hstack(tmp_array).

Does this make more sense?

Cheers,
Eirik

[Cirdans-Home]

Dear Eirik,

I am sorry for the very long delay in answering. First of all, I want to thank you for the help, I think I managed to use all your suggestions to obtain the test problem:
manyfractures2D.txt

When generating the matrix I observed that there were few zeros on the main diagonal. From the discussion in the article, it is not clear to me whether this is the expected behaviour, or if I introduced some errors into the classes I overloaded. Again, every suggestion is really helpful.

All the best,
Fabio

[keileg]

Dear Fabio,

I'm glad to hear that you made it work, from a quick glance, the code you attached looked reasonable.. This is definitely not an easy problem to set up, due to the geometric complexity, parameter heterogeneity etc., so it is impressive to see what you have done. We have discussed this case among the core developers, and agreed it would be very interesting to have this as an application example that ships with the code. Would you be interested in contributing your script to that end, with the understanding that we take the responsibility for cleaning up the code and keeping it up to date (at least for the foreseeable future)? I should mention two caveats:
• We will likely do substantial changes to your script to bring it closer to the standard style we use, but potentially also to exploit synergies with related setups, or for other reasons. You will of course be fully credited for your work.
• It may take us a while before we get around to updating the script; in addition to the usual activities, we are in the middle of a release cycle and have limited capacity for the moment.

If this is of interest to you, I suggest you simply make a pull request of the script into the folder src/porepy/applications, then we will move it to the appropriate subfolder and start modifying.

Regarding the script, you should be aware that we have changed how boundary conditions are specified, so if you pull the latest version of the develop branch, your script will likely stop working as intended. The new setup is outlined in the tutorial for single phase flow. To be clear, if you make the PR, we will take care of this change (when we have time).

As for the zero diagonal elements, we expect some of these, but not too many: In the 0d grids, there is no internal flow, thus mass balance only considers fluxes over interfaces to the neighboring 1d grids, with no contribution from the 0d pressure, hence a zero diagonal (this is only true for incompressible fluids, if there is compressibility, the pressure in the 0d cell will enter the mass balance through a dp/dt term). The 0d pressure will contribute to the interface fluxes, though, so its full column will not be zero, but should contain 4 elements. Summarized, there should be as many 0 elements on the diagonal as there are 0d cells in the grid. If there are more than that, my hunch is that something is wrong.

Best Regards,
Eirik

[Cirdans-Home]

Dear Erik,

I finally managed to prepare and send the pull request, I have updated the code to work with the new boundary condition definition and substituted the old array wrap command with wrap_as_dense_ad_array. I have done some testing, and it seems to work.

Thank you for the further information on the matrix, I'm trying to do some work on preconditioners/linear solvers for this problem with a library of which I am one of the developers (PSCToolkit). This kind of problem seems to me to be a very interesting test bed.

Best regards,
Fabio

[keileg]

Dear Fabio,

Thanks for the effort (see also comments in the PR #1081). Linear solvers for these problems are indeed interesting, we dabble in it ourselves from time to time, though not on the scales you are working on.

Best,
Eirik