Flow problem with no Dirichlet bc

[alessiofumagalli]

Dear all, I'm planning to implement the condition (zero/or given mean pressure) to fix the uniqueness of the solution in case of only Neumann bc flow problem. Is it already present? If not, is it something you are interested to have in the pp core? Let me know, thanks.

[keileg]

It is not present in the core, but we have talked about it on various occasions, e.g. while working on solvers based on domain decomposition, where fractures not in touch with a Dirichlet boundary form floating subdomains. So in principle, this is definitely of interest; of course it depends a bit on how the implementation ends up. How do you plan to do it - on the matrix level or discretization level?

[alessiofumagalli]

I don't have the implementation yet, but I was thinking to introduce a new discretization that might be added to fix the value with its appropriate Lagrange multiplier. This probably fit well in the framework of pp.

[alessiofumagalli]

after some thoughts, here's my strategy (I assume one single grid so far)

• create one fake 0d grid and a mortar grid that couples all the cells in the real grid with the single dof associated to the fake 0d grid
• make the discretization that fix the pressure mean as an edge discretization
  this should work, all the dof should be set correctly as well as the matrices dimensions. We will be able to access all the variables. I might come with a prettier solution that does not involve the creation of the 0d grid.

[keileg]

I see the reasoning, and it kind of makes sense, but pretty it is not. The underlying issue is that we have assume a 1-1 relation between equations + variables and quantities (cells, faces, nodes) in the entire grid. So the thorough fix is to allow for more general variables and equations, but that is a non-trivial extension which requires some thinking.

My suggestion is that you find an implementation you think will work for your case and push the solution to porepy_labs (which we try to use as a staging ground for sharing ideas and code that may not yet be ready for inclusion in the main repository) - I just gave you write access. Feel free to make a new directory there, with a long and descriptive name. Then we can experiment with how to do this, and move the code to the main repo if that feels natural at some point. Even if we don't, the trick you develop will be available to other interested parties.

Does this sound like reasonable way forward to you?

[alessiofumagalli]

If someone is interested, here's a possible solution

   A, b = assembler.assemble_matrix_rhs()
    A, b = fix_integral_of_variable(A, b, assembler._dof_manager, variable, gb.get_grids())

    solution = sps.linalg.spsolve(A, b)[:-1]

def fix_integral_of_variable(A, b, dof_manager, variables, gs):
    variables = np.atleast_1d(variables).tolist()
    dof_block = dof_manager.grid_and_variable_block_range(gs, variables)

    values = np.zeros((len(variables), dof_manager.num_dofs()))
    for (g, v), dof in dof_block.items():
        values[variables.index(v), np.arange(*dof)] = g.cell_volumes

    new_A = sps.bmat([[A, values.T], [values, None]], format=A.getformat())
    new_b = np.append(b, [0.]*values.shape[0])

    return new_A, new_b

Should we put it in the core of porepy?

[keileg]

Not sure, but this could potentially be useful. @IvarStefansson, what do you think?

[keileg]

Me neither. Possibly the new matirx_operations.py, proposed in #613 (if you approve of it). It is a bit outside the scope of matrix_operations as the module is defined now, which may tell us that either we need another module (but that will give more fragmentation), or we can expand the scope.

[alessiofumagalli]

another possibility is to consider it as somehow related to the bc, otherwise the matrix_operations module could be a possibility

[keileg]

This is numerics, so it should enter in a numerics module.

@alessiofumagalli If you make a PR with the function added in https://github.com/pmgbergen/porepy/blob/9966d696ec41d85c06354ad1c71f93a7c28e7f10/src/porepy/numerics/linalg/matrix_operations.py
I'll take care of it. I can of course copy the code myself, but you made it, so the commit should be associated with you.

[alessiofumagalli]

sure, thanks! I'll add also some integration tests.

[keileg]

Thanks!