Unify fracture network generation in 2D and 3D

[jhabriel]

The main idea is to generate a FractureNetwork2D by passing a list of 1D fracture objects.

In this way, we will be unifying the construction of the classes between FractureNetwork2D and FractureNetwork3D.

Furthermore, we introduce a fracture_network.py sub-module that contains a method for creating and returning the appropriate FractureNetwork class, given a domain and list of fractures. The correct class (2D or 3D) is defined implicitly by the arguments.

At last, we need to update the meshing tutorial. This will be postponed, until the unification of the mesh generation and the domain definition are done.

[pschultzendorff]

Hei,

I pushed a first attempt at the unification. My approach so far is to use two functions linefractures_to_pts_edges and pts_edges_to_linefractures that I created in fracs.utils, which convert between LineFractures and tuples of points and edges. These functions are then added to FractureNetwork2d at the necessary places, while the internal functioning is left alone.

However, one loses the functionality to pass tags along the edges when creating the network. One possibility would be to add the option to pass the tags without the edges, but I don't know if this is a sensible approach and would like to hear some expertise from you on this.

Also, I am not quite sure, whether there are some other places in the source code that needs to be adjusted, so far I only found a call to FractureNetwork2d in grids.standard_grids.utils.make_mdg_2d_simplex.

edit: And one more thing: I added the tests for the conversion functions in a seperate file, as the test file for fracs.utils uses unittest instead of pytest. To avoid occlusion of the tests, these should be combined, but I couldn't find any info on whether one can just combine both types of test and haven't tried it out yet.

[keileg]

Thanks for the update, @pschultzendorff. Some thoughts:

• My instinct is that one can expand the __init__ of FractureNetwork2d (or introduce a separate __init__) so that the network can be initialized using LineFractures. This new init would translate to the format internal to FractureNetwork2d, and take it from there. This will amount to a simplification of the user interface, to the price of a slight but localized increase of complexity. Does this make sense, or am I overlooking something.
• Passing along tags will be a problem with the above solution; I think your suggested solution is fair. We could also add a tag attribute to the Fracture classes, indeed this is perhaps the more natural option, but it could have major repercussions in the fracture processing, so I am hesitant to advocate this.
• You are right that there are few other adjustments needed in the source code. For standard usage, the FractureNetwork classes provides the mesh and are then forgotten.

[jhabriel]

Thanks @pschultzendorff. I think this is headed in the right direction.

• I agree with @keileg suggestion. Instantiation of the FractureNetwork2d should be done by passing a list of LineFractures.
• Are tags really necessary for generating FractureNetworks? Arguably, once the MixedDimensionalGrid is created, each subdomain and interface receives a unique identifier.
• This would be a matter of checking how many tests fail after introducing the changes. I'm not too familiar with this part of the code though.

[keileg]

@jhabriel Tags are part of the API (to the degree that exists), so it should not be removed lightheartedly. I for one have used it extensively in some projects.

[pschultzendorff]

Hei, thanks for the comments @keileg and @jhabriel.

• The way you describe the instantiation is exactly the way I implemented it.
• Your latter suggestion sounds more elegant @keileg, what are the repercussions you are talking about?
• I haven't run all tests on the new code so far, but only the ones specific to the fracture networks (and the newly introduced functions). Will do so when I continue working on this.

I think we should just discuss this shortly in this week's usermeeting and go forward from there.

[keileg]

If you start changing the internal data formats for tags, or anything else related to geometry processing, all kind of trouble will break loose. But if you do a translation of the tags, the same way you translate the fractures, that should be fine.

[jhabriel]

Short update:

We agreed on adding a tag: int attribute to the abstract Fracture class. To preserve backward compatibility, the FractureNetwork2d will retrieve said tag from each LineFracture and use it to "translate" back into the current way of generating the network.

[jhabriel]

Hi. I pushed a new module containing the wrapper method create_fracture_network() that unifies the generation of fracture networks in 2d and 3d:

https://github.com/pmgbergen/porepy/blob/780-unify-fracture-network-generation-in-2d-and-3d/src/porepy/fracs/fracture_network.py

Let me know what you think.

[pschultzendorff]

Hei,

I started work on implementing tags into the Fracture class and came around the following problem:

As of now, tags are passed to the FractureNetwork2d class by the edges array. This automatically ensures that each fracture/edge has the same number of tags. However, when each Fracture stores their own tags, the number of tags may vary between fractures, which would create a ragged array when they are internally converted into edges again.

Now my question is, if there is some reserved None tag that can be assigned to edges that have less tags than the edge with the highest number of tags? The easiest way would probably be to assign 0 or -1 to mean "not tagged", but as the documentation mentions that tags can for example be imported from gmsh, I am unsure whether any of these numbers is used in gmsh.

Also, am I correct in the assumption that (theoretically) an edge can have an arbitrarily high number of tags?

[jhabriel]

Hi Peter,

Thanks for the update. I understand the problem and I think your suggestion is fair (as long as it is properly documented). I'm sure @keileg has more experience with tags, but it in principle it should be safe to use -1 for untagged fractures.

By reading the documentation of FractureNetwork2d it is indeed possible to assign an arbitrary amount of tags to each fracture. Here's a suggestion that might solve the problem:

1. Define tag as a numpy array with dtype=np.int8.
2. Transverse the list_of_fractures to determine the LineFracture with the highest amount of tags max_num_tag.
3. Pad the bottom rows of edges array of -1 with shape=(max_num_tags, len(list_of_fractures)).
4. Replace the -1 entries where there exists non-empty tags.

[pschultzendorff]

Hei,

thanks for the input! That sounds like a nice way to solve the problem. I will implement it, once @keileg gives the okay to use -1 for untagged fractures.

[keileg]

I think that if you assign -1 as a default tag and document it, this will be both acceptable and an improvement.

[IvarStefansson]

I merged this one into my well branch. There may be commits for you to cherrypick (for now related to test fixes). Let's coordinate!

[pschultzendorff]

The tags work now. I also changed tags that were previously initialized with np.empty in FractureNetwork2d.add_network, to be initialized with -1 for consistency.

There remain 8 test files with errors and two errors during test discovery in the functional tests. Most importantly the index tests for the PlaneFracture class fail and I am not quite sure why. I haven't looked at the other tests, but I assume they are mostly just calls of the FractureNetwork2d class with edges and points instead of fractures.

[jhabriel]

Thanks @pschultzendorff. We're making good progress. I will make a list of the next steps.

[jhabriel]

• [x] Merge develop branch and resolve conflicts.
• [x] Fix failing tests.
• [x] Revisit create_fracture_network() and write tests.
• [x] Revisit pts_edges_to_linefracture() and linefracture_to_pts_edges() and discuss tests.

[OmarDuran]

I merged develop branch and resolved the conflicts. Also, I have included two commits from "wells_in_models" as requested by @IvarStefansson.

[jhabriel]

Hi all. I fixed the failing tests and now all tests pass with the current version of develop.

I fixed a few bugs, which is an indication that we need to write/improve the tests. Let's discuss specific actions tomorrow.

[jhabriel]

@keileg : Design choice question. We are unifying the creation of fracture networks by introducing the function pp.create_fracture_network(). Now, do we want to still allow for the construction of fracture networks 2d via pp.FractureNetwork2d (and pp.FractureNetwork3d)? Or do we make these classes "private" sort to speak?

Lingering a bit on the reason of the question, what we want to avoid is a scenario where we introduce even more diversification (i.e., providing 3 ways of constructing fracture networks) instead of unification. One possibility that we considered is to remove the classes FractureNetwork2d and FractureNetwork3d from the porepy init. But we are not sure if this would be enough.

[keileg]

Good question. My instinct is to lead the user towards the new function, and I think that removing the classes from the top-level init could be a good idea; but that depends on possible changes to the general approach to inits - I can explain a bit more in a more suitable channel.

Two questions:

1. Does the function returns a FractureNetwork,?
2. If yes on 1., will Cartesian meshing also be based on a FractureNetwork?

[jhabriel]

Answering the questions:

1. The function will return FractureNetwork2d or FractureNetwork3d, depending on the dimension of the fractures (or alternatively, the domain when there are no fractures).
2. Indeed, that's the plan (via #795). Basically, create_mdg() will take a fracture_network, a grid_type and a mesh_size, and it will generate the mdg. Cartesian grids are included in grid_type, of course. Could you confirm that this is the case @OmarDuran?

[keileg]

Sounds good. If possible, I think you should then make it unnatural (whatever that means) to generate a network directly through the class inits - the init classes can be considered private to that degree that this makes sense. However, the FractureNetwork classes should still be easily accessible in documentation, since users will interact wiht this.

I don't know how to best combine these goals, but I think we will go with whatever solution you find acceptable.

[jhabriel]

Yes, I agree. I will go for removing FractureNetwork2d and FractureNetwork3d from the top init. This will force anyone that wants to use these classes in a runscript to import them directly (which IMO is tedious enough to avoid such practice). Thanks for the feedback :)

[OmarDuran]

Answering the questions:

1. The function will return FractureNetwork2d or FractureNetwork3d, depending on the dimension of the fractures (or alternatively, the domain when there are no fractures).
2. Indeed, that's the plan (via Unify mixed-dimensional grid generation #795). Basically, create_mdg() will take a fracture_network, a grid_type and a mesh_size, and it will generate the mdg. Cartesian grids are included in grid_type, of course. Could you confirm that this is the case @OmarDuran?

2. Yes, I confirm that is the plan.

[jhabriel]

Second part

For this second part of the project, we will be working on the branch unify_fracture_network_part_2.

The roadmap is the following one:

• [x] Remove FractureNetwork2d and FractureNetwork3d from the top level __init__.
• [x] Replace instances pp.FractureNetwork2d and pp.FractureNetwork3d by pp.create_fracture_network() in src.
• [x] Fix affected tests.
• [x] Fix affected tutorials.
• [ ] Update documentation of the fracs subpackage to be in line with current standards.

[jhabriel]

Documentation roadmap

Files and assignees for documentation update of the fracs subpckage:

• [x] fracture (@jhabriel)
• [x] fracture_importer (@jhabriel)
• [x] fracture_network (@jhabriel)
• [x] fracture_network_2d (@jhabriel)
• [x] fracture_network_3d (@jhabriel)
• [x] line_fracture (@pschultzendorff)
• [x] meshing (@pschultzendorff)
• [x] mesh_2_grid (@pschultzendorff)
• [x] non_conforming (@pschultzendorff)
• [x] plane_fracture (@pschultzendorff)
• [x] gmsh_interface (@vlipovac )
• [x] simplex (@vlipovac)
• [x] split_grid (@vlipovac)
• [x] structured (@vlipovac)
• [x] tools (@vlipovac)
• [x] utils (@vlipovac)
• [x] wells_3d (@jhabriel)

Reminder: Changes should be incorporated in the unify_fracture_network_mirror branch.

[pschultzendorff]

Hei,

I uploaded my part documentation. The non_conforming module is not finished, however, as I struggled quite a lot with the missing documentation. Haven't managed to get all the types right so far and there are also more underlying problems, as, e.g., update_nodes has an argument in the signature that is not used and update_face_tags never accesses the values of the inner list of its argument new_faces (this is kind of weird).

To finish documentation of this file, I will need some input. Won't be able to continue work before Friday though.

Best, Peter

[vlipovac]

Hi,

I have merged the docs-rework I was responsible for. Some things remain open:

1. The function split_grids.split_specific_faces is unused, and most likely unfinished. The docs are accordingly more or less guessing work
2. I am unsure about the helper function split_grids._sort_sub_list. This is some indexing voodoo related to sparse matrix and I am unsure if I captured the docs well. Can someone take a look?
3. Mypy was complaining on several occasions with regards to other files. pp.matrix_operators.slice_indices f.e. needs an overload documentation. Also, the return type Union[FractureNetwork2D, FractureNetwork3D] causes mypy errors, because these two classes have apparently different members and it cannot figure out in the logic which one is addressed.

Also, before the PR we should compile the whole documentation of the subpackage fracs and have a final look at it and polish the results.

[jhabriel]

Hi @pschultzendorff and @vlipovac, thanks a lot for the effort.

Peter:

I have a quick look at the non_conforming module and is indeed heavily undocumented. Let me know when you're done with the module and will have a go.

Veljko:

1. The function split_grid.split_specific_faces is used in _split_fracture_extension() from the propagate_fracture module. I've added a warning that the function should be used with care.
2. I had a look but I couldn't figure out what is going on there. The code is too obscure and the names of the variables don't help much either. I think we're good with the Todo you wrote. I've complemented it a bit to alert the reader of the low level of certainty we have regarding the docstrings.
3. Indeed, the union causes mypy some headaches. The workaround when the fracture network object demands to be dimension-specific is to assert whether this is a FractureNetwork2d or a FractureNetwork3d. This is the price we sometimes have to pay after unifying the fracture network generation.

P.D.: The branch is now up to date with develop. Moreover, I fixed all the static typing errors.

Best,

[pschultzendorff]

Everything should now be documented. As discussed, non_conforming.py is left with deprecation warnings. @vlipovac: I did not add any overload signatures, as mypy would always complain and I couldn't figure out the problem.

[pschultzendorff]

Pushed the file now. I get a lot of pytest errors on the branch, but this could also be a local problem with my python env.

[vlipovac]

I took a look at the docs in the subpackage and reworked quite a bit. The guidelines for documentation and the linking should now be enforced to a satisfying level.

There are a lot of open Todo notes though, possibly critical for the functionality of porepy. Respective parts have Todo: directives in their docstrings.

I pushed the changes to the mirror branch @jhabriel

[jhabriel]

Thanks a lot @vlipovac. I will update the branch and open the PR.

[jhabriel]

Closed after merging of #816