Various issues

[IngeborgGjerde]

Hi,

In checking out this repo I made note of various issues/concerns/suggestions, hope this is of help 🙂

Mathematical issues

• The variational form in the demos use a finite element pair that is not stable. For the hydraulic network model (i.e. a dual mixed Poisson equation on a network) you should use P1-P0 elements. P2-P1 elements are appropriate for Stokes type network models, as they contain a viscous term. See the original fenics implementation for more details.

• The variational forms you consider are straightforward implement if you use the primal form of the model, which reads:

Find $p \in H^1(\Lambda)$ such that

$$(\mathcal{R} \partial_s p, \partial_s v) _\Lambda= (f,v)_\Lambda$$

for all $v \in H^1(\Lambda)$, where $\Lambda$ is the network domain. This can be discretized with P1-elements on the global mesh.

For more information you should check out the graphnics introduction demo. This approach works really well in practice, see for example this repo containing numerical comparisons against analytical solutions.

With all of this in mind, I recommend not spending too much time on optimizing the dual mixed assembly time, given that there isn't strictly a need for it.

Implementation issues

• I can't get the following demos to run using the provided docker image:
  • demo_honeycomb.py, demo_double_Y_bifurcation.py and demo_arterial_tree.py fail with

  A = _cpp.fem.petsc.create_matrix_block(a)
  ValueError: vector::_M_range_insert
  

  • demo_tree.py fails with

  A = _cpp.fem.petsc.create_matrix_block(a)
  RuntimeError: Cannot insert rows that do not exist in the IndexMap.
  

• There is an issue with running out of MPI communicators as soon as there ~100 edges in the network. This causes the code to crash.

Copyright issues

Currently, the code is licensed under an MIT license stating Copyright 2022 Cécile Daversin-Catty. However, as the code contained in this repo is directly based on graphnics / fenics-networks, you need to retain the original GNU General Public License v3.0 along with original copyright notices.

I suggest adding in the copyright notice in graphnics where relevant, along with information on what was changed and when. At a minimum, this needs to be done for the following files:

• networks_fenicsx/mesh_generation.py, being a copy of mesh generators in graphnics / fenics-networks, see for example here,
• networks_fenicsx/mesh.py, being a copy of the FenicsGraph class in graphnics / fenics-networks adapted to fenicsx,
• solvers/assembly.py, being a copy of the flow model implemented e.g. here. Moreover,
  • the jump_vector function was copied verbatim from here,
  • the method you use for implementing the bifurcation conditions, by inserting the jump vector into the global matrix, was taken from here,
• utils/petsc_utils.py, being verbatim copied from here.