Bug in Nodal values, Materials and BCs

[WilkAndy]

Nodal values are not working with Materials and BCs in the way I expect. Could someone tell me whether this is a bug, or a deliberate feature, or whether I'm misunderstanding something.

Consider the case with a linear quad element with one of the edges being a boundary. Place an IntegratedBC on this boundary.

(1) Moose loops over _i (0 to 3) and _qp (0 to 1) to evaluate the Residual.

(2) The range of _i presumably comes from the element, and the range of _qp from the boundary.

(3) The QpResidual is multiplied by _test[_i][_qp], and I've verified that _test[_i][_qp] are the quad-element test functions evaluated on the boundary. This is expected.

(4) My Variable for this IntegratedBC, _u[_qp], is correctly sized (_qp = 0, 1).

(5) When i have a coupledNodalValue in my IntegratedBC, it is sized by _i. Viz, it has size 4 in this case. I have verified that the values are the values at the 4 nodes. In my IntegratedBC i just need to index it with [_i], good, and the loop over _i and _qp with _test[_i][_qp] will correctly pick up what i want (the nodal values off the boundary will get killed by _test being zero there).

(6) When I use _mat_prop(getMaterialProperty("blah")), it is sized by _qp. I have verified that the _mat_prop are the material properties correctly projected onto the boundary (except in the case (7), below). So in my IntegratedBC, i'd just use _mat_prop[_qp] and everything would be fine.

(7) HOWEVER, when my MaterialProperty comes from a coupledNodalValue, it is still sized by _qp. It's values are screwy (according to me): in this case they are just the first two values of the 4 values in the element of the coupledNodalValue. That is, they are not the values projected onto the boundary: they are just two random nodal values.

I can provide examples, but since it involves me writing a clean Material and BC and simple test file, which will take a while. So I thought i'd ask first in case there's something wrong in my logic, or someone can see an immediate reason why this is true and can fix it.

Tagging @andrsd because i think he did the coupledNodalValue stuff.

[andrsd]

You can only get variable values at nodes, not material properties. Materials always work with q-rule.

If you compute material property at an element, it is unrelated to the material properties on a side. The material system does not work the same way nodal values do. The 4 internal material property values are stored in a different memory than the 4 * 2 values for edges...

Remember there are 2 (actually 3) instances of each Material, one for the volumetric version, one for the boundary version. When you are doing your volumetric computation you have 4 qps and 4 nodal values. On a side, 2 qps and 4 nodal values. So, in order to store you material properties correctly (so that the BC can pick it up), you will need to related your 2 matprop values to the right nodal values. Look at Material::_bnd variable as a way to determine if your material acts on a volume or boundary.

[WilkAndy]

I think I understand what you're saying, @andrsd , regarding volumetric and boundary Materials. However, how does my IntegratedBC choose to get the boundary version? At present I have:

_pp(getMaterialProperty<std::vector<Real> >("PorousFlow_porepressure_nodal")),

in my IntegratedBC constructor. Here _pp.size() = 2 (correct, i believe), but _pp[0] and _pp[1] are incorrect. They are equal to the _pp[0] and _pp[1] values I would obtain if I had the same line as above in a Kernel instead. (Of course there _pp.size() = 4, but the _pp[0] and _pp[1] values are identical.)

[WilkAndy]

Oh, the std::vector in the above is irrelevant - i just copied-and-pasted from my real example. Just pretend it's simply getMaterialProperty

[WilkAndy]

[WilkAndy]

And in my Material:

_porepressure_nodal(declareProperty<std::vector<Real> >("PorousFlow_porepressure_nodal")),

_porepressure_nodal_var(coupledNodalValue("porepressure")),

[WilkAndy]

Hmmm, i don't think i've woken up yet. One final(?) line from my Material:

_porepressure_nodal[_qp][0] = _porepressure_nodal_var[_qp];

[WilkAndy]

I think you're going to tell me i'm abusing the Material system by storing nodal information in something indexed by _qp (see #6940). But i still don't understand why we can't get this to work.

[andrsd]

how does my IntegratedBC choose to get the boundary version?

It always gets the boundary version, because integrated BCs work only on boundaries. They do not work with volume at all, it does not make sense for them.

I think you're going to tell me i'm abusing the Material system by storing nodal information in something indexed by _qp

You already know that ;-) You need to figure out which _i in the BC corresponds to _qp. You should be actually doing that in your kernel as well, because I think you are just getting lucky on the numbering, see this example (inner numbers are _qp indices, outer nodal _i):

4        3            3        2
 +------+              +------+
 | 4  3 |              | 3  4 |
 | 1  2 |              | 1  2 |
 +------+              +------+
1        2            4        1

  ok                    not ok

Similarly in your BC, you have to figure out that (4, 1) and (3, 1) belong together (assuming your BC is on 3-4 edge)...

[WilkAndy]

Thanks @andrsd . Now i have two questions:

(1) How can i figure out which _i corresponds to which _qp? Actually i have to do this in my Materials, as they are where basically everything happens in porous_flow. All i can think of is to look at which _test is biggest. That seems a bit bodgy - is there a better way?

(2) In my BC i have no choice. It only receives only two things, so somehow i need to set up things beforehand(?) I know my multiple comments above were rather confused, but in my Material i use coupledNodalValue to create a thing called "PorousFlow_porepressure_nodal" which stores nodal porepressure values at the Material's qps, which is then used by the BC. The BC receives two values only.

I think you know all this, but i'm trying to be as clear as possible.

In your "OK" picture, the boundary is on line 1-4, and the BC receives the values porepressure_1 and porepressure_2. Obviously not what i wanted.

Thanks for spending the time guiding me through this, David. I'm looking forward to getting this right, as i suspect some convergence issues in the Richards module are due to this problem too.

a

[andrsd]

If you want to investigate convergence issues and see if you have problems with ordering _i vs. _qp. you can compare the distance of a q-point to a node. Nodes inherit from Points, so (_qpoint[_qp] - node[_i]).size() should give you the distance. The problem here is obvious, you have to compare all _is to all _qps before you do anything. That can happen in Kernel::precalculateResidual, so you are not doing it at each qp. This will sacrifice some performance, of course. We do not have the same function for integrated BCs, but you can mimic what is in Kernel and add that temporarily on some branch (and maybe open a PR for it - it is quite a simple task). You want to build an indirect index from _i to _qp, not necessarily std::map, std::vector should give better performance and it would not waste too much memory, your arrays are small. With the indirect map, your indexing into arrays that take _qp would be _value_expecting_qp_indexing[_my_qp_map[_i]]. Note, that this map is different in Kernels and BCs, so you cannot reuse them...

BTW: The real problem you are fighting is that MOOSE can do only one type of quadrature rule. If you needed just nodal values everywhere, you'd use FIRST TRAP rule and all would work like a charm... But, I think we already talked about it in the past, right?

[WilkAndy]

Thanks @andrsd , i'll try your ideas. Yes, i'm fighting MOOSE's one-q-rule restriction. Maybe one day someone with lots of $$$ will get the libmesh people to generalise.

[andrsd]

We could do multiple q-rules in MOOSE (it is not really a libmesh limitation, I think). The problem is that you need to tell MOOSE which kernels/bcs/materials, etc. are using which q-rule. We could "simply" add another FE object with the other qrule and let MOOSE recompute all of that. There are obvious difficulties with things like stateful material properties. And another whole lot of other small things that we would have to get right. As you can imagine qrule is almost everywhere, so it is quite a big task...