Mpsa subface formulation #176
Conversation
Dual_vem fails, but that is from the merge from the merge_develop_multiphysics branch
…f a vector is given
…is and robin was switched. I believe this is a bug that happened when we changed the ordering of the basis and robin weight to place the face index at the last axis compared to the second as was done initially.
…t for the Neumann boundaries.
…ion region is connected to both slave and master
There was a problem hiding this comment.
I have two substantial comments: First, more comments are needed - for instance it should be clear from the documentation of classes and functions what is their intended usage. This should be the case even for private functions. The more technical parts of the code also needs more explanation - it should not be necessary to run a debugger through one line at a time to understand the logic of the code.
Second, I see what you mean to do with the FVEllipticVectorial concept, but I do not think we should do that kind of abstraction yet. If nothing else, it is highly unlikely we will have another FV discretization for elasticity anytime soon. So my preference right now is to move the relevant functions into mpsa, and delete the parent class. However, at this stage it is not clear to me how the corresponding functions will be for the extra terms in poro-elasticity. So perhaps it is better to leave the code in this PR, and revisit it with a critical mind in the biot branch?
| @@ -0,0 +1,54 @@ | |||
| import numpy as np | |||
There was a problem hiding this comment.
We need documentation here, both on the module and class level.
| import porepy as pp | ||
|
|
||
|
|
||
| class FvSubGrid(pp.Grid): |
There was a problem hiding this comment.
Is this a grid that replaces subcells used in mpfa/mpsa with proper cells in the porepy sense, and similar with subfaces and nodes? If so, I am most interested, this is related to another project of mine. But the fact that I can't say for sure shows the documentation need.
There was a problem hiding this comment.
I wrote this function exactly with the proper sub-grid in mind. In 2D it was straight-forward, but in 3D it was a bit more tricky, and I abounded this function and just forgot about it.
I can not decide if it has the right of life or not. It might be a convenient way to set the subface-boundary conditions if you don't know anything about SubfaceTopology, however, then you are not able to interpret the output of Mpsa which is on the subface level.
I think the cleanest would be to actually construct the sub-cell grid, let mpsa/mpfa work on this, and return the discretization on the sub-grid. This would simplify mpsa/mpfa considerably and possible make the code readable. But as you say, this is a project in its own
| @@ -297,7 +297,7 @@ def _check_mappings(self, tol=1e-4): | |||
| # ------------------------------------------------------------------------------# | |||
|
|
|||
|
|
|||
| class BoundaryMortar(object): | |||
| class BoundaryMortar(MortarGrid): | |||
There was a problem hiding this comment.
Do we still need the BoundaryMortar as a separate concept from the standard MortarGrid? It may very well be, I am curious to understand if this implicitly is telling us that the implementation is poor.
There was a problem hiding this comment.
For now, yes, Ideally, no.
It should be quite straight forward to remove the BoundaryMortar, however, MortarGrid was written with only mixed-dim couplings in mind, which means that there still are a few differences that needs to be sorted out first.
| # convert to indices. | ||
| # First check for dimension | ||
| if faces and c.size != g.num_faces: | ||
| raise IndexError("boolean index did not match number of faces") |
There was a problem hiding this comment.
Raising of errors should be documented in the function preamble.
| @@ -344,3 +347,429 @@ def enforce_neumann_int_bound( | |||
| """ | |||
| # Operation is void for finite volume methods | |||
| pass | |||
|
|
|||
|
|
|||
| class FVVectorElliptic(pp.numerics.mixed_dim.solver.Solver): | |||
There was a problem hiding this comment.
Do you actually use this class for anything? If not, it feels like a premature abstraction that may cause more harm than good in the future. As an example, it indicates that mappings between grids in different dimensions is straightforward for elasticity - fact is that nobody knows this.
Unless I am overlooking something here, my preference would be to include what is needed to get contact mechanics working into mpsa, and then remove the rest - but save the commit id.
There was a problem hiding this comment.
I kind off agree, but I don't see the harm in the abstraction. If you want to couple different dimensions you first need to write a mixed dim interface law, and if you do that you know enough to know what you are doing. But I agree most (or all) can be moved to the Mpsa-class. However, I think we should add a VectorEllipticDiscretization to the mixed-dim folder (this folder should not be called mixed-dim, but that is an other issue) as Hau has plans of finite elements for elasticity.
There was a problem hiding this comment.
Also, this specify the contract between the node-discretization and interface law. Keeping this separate from mpsa should be beneficial.
| from porepy.numerics.fv import fvutils | ||
| from porepy.utils import matrix_compression, mcolon, sparse_mat | ||
| from porepy.grids import structured, partition | ||
| from porepy.params import tensor, bc |
There was a problem hiding this comment.
It is good to see those old imports disappearing.
| and give the displacement at the points given by eta, also for the boundary. | ||
| eta (float or ndarray, range=[0,1]): Optional. Parameter determining the point | ||
| at which the displacement is evaluated. If eta is a nd-array it should be on | ||
| the size of subcell_topology.num_subfno. If eta is not given the method will |
There was a problem hiding this comment.
This also means that etas that vary on the sub-face level must be set according to the numbering in that array? Sound like one really has to know what one is doing, but it has to be that way, I suppose. The alternative of leaving it to fvutils is fully acceptable.
There was a problem hiding this comment.
Yes, eta has to be set according to SubcellTopology. I can't really think of any other way of setting this. I guess if you want to set the continuity points on the subcells it can be required that you know what you are doing
| else: | ||
| master_ind = 0 | ||
|
|
||
| self.discr_master.assemble_int_bound_displacement_trace( |
There was a problem hiding this comment.
Now I see where you get the FVEllipticVectorial stuff from, some of those functions are indeed necessary. Nevertheless, I think those should be moved to mpsa - in particular since it seems unlikely we will invent and implement a new fv discretization for elasticity anytime soon.
|
|
||
| """ | ||
| # We don't call the init since we don't have access to the grid. | ||
| # Maybe this is bad style. |
There was a problem hiding this comment.
Quite likely. Let us see if anyone will ever notice.
| @@ -4,144 +4,6 @@ | |||
| import porepy as pp | |||
|
|
|||
|
|
|||
| class TestUpdateDisc(unittest.TestCase): | |||
There was a problem hiding this comment.
Good catch. It must have been deleted by an accident
There was a problem hiding this comment.
Never mind, I now remember. It was deleted on purpose. This test does exactly the same as the other class in the same file. There was a purpose of two tests for the early versions of displacement reconstruction, but at the current version there is no difference between the two tests.
There was a problem hiding this comment.
I agree that it could be useful to state the contract explicitly in an EllipticVectorDiscretization, and that FVVectorElliptic (we should think about the order in which we merge words in order to make this as uniform as possible throughout the code, FVEllipticVector(ial), cf BCs?) is probably superfluous.
Apart from a couple of minor suggestions which I have mentioned directly to Runar, I have no further objections.
|
|
||
| """ | ||
|
|
||
| def __init__(self, keyword, physics=None): |
There was a problem hiding this comment.
The physics keyword was purged in the parameter branch. We only kept the one called "keyword".
…on params to Paramter dictionary. The naming change of eta does not belong here. If the eta name is changed, this has to be done on a module level (not just at the outer most level). The natural choice now is that all parameters for the discretization now belongs in the Paramert dictionary, not the outer dictionary
…ization used slave discretization for mortars
|
I'm ready to merge this if you agree @keileg @IvarStefansson |
There is a new Python 3.5 issue in test_elliptic_interface_laws.py. This just shows the importance of fixing the other one (#177), as it increases the risk of overlooking further 3.5 failures. |
It should work now. It was the recurring error of node_numbers switching randomly in 3.5 making the tests fail at random |
The main part of this pull request is an expansion of Mpsa so that we can retrieve the subcell and subface discretizations. This enables us to specify the boundary conditions for each subface instead of on each face. Further it let us obtain the subface tractions instead of the face traction. This increased accuracy is crucial if you are interested in the boundary (contact mechanics). If you don't care about the subfaces everything should work as before. Mpsa will only return the subfaces if the boundary-condition is given for the subfaces instead of faces.
This pull request also includes some improvements in the subface-displacement reconstruction, notably a chang from 'C' ordering to 'F' and from giving the displacement per gradient on the subfaces, to 1 displacement per subface (the displacement is averaged over the displacement from the gradient on the two sides. At the continuity point this is the same)
We have also included two interface-laws for the elasticity problem, a Robin coupling and a full continuity coupling.
Before accepting this pull request have a look in fv_elliptic and comment on the classes of FVElliptic and FVEllipticVector. Should we have both? Should they be renamed?