Linear Raviart-Thomas Finite Elements#3638
Conversation
|
Job Coverage on 8bd7b33 wanted to post the following: Coverage
Warnings
This comment will be updated on new commits. |
||||||||||||||||||||||||||
|
I'm in between chores and errands on my day off, and I probably won't get to a review here until Monday, but let me at least immediately say: Thank you!!! Adding H(div) and Raviart-Thomas have been something in the "on my todo list but getting buried lower faster than it bubbles higher" category for ... ouch, literally a decade now. This is wonderful to see. |
|
You copied an 11-year-old typo from one of @pbauman's Nedelec code comments. I'd have preferred it if you'd kept re-bootstrap results in separate commits. Unit test coverage would have been nice (just like it would have been for the other vector elements I decided not to cover when I was writing fe_test.h...) But I'm not seeing anything here yet to really criticize. I'd like to play around a bit with this tomorrow afternoon, but if anyone else hits approve before then there's no need to wait on me. |
I caught a few, let me know if I tackled the one you spotted. :)
Done.
Please do! :) |
|
Not sure whether he has the time, but it would be valuable to get @pbauman's input on this PR as well given all the prior work that he did in adding vector-valued FE to libmesh. |
jwpeterson
left a comment
jwpeterson
left a comment
There was a problem hiding this comment.
I'm not really qualified to comment on the FEM aspects here, but I definitely appreciate that you've added an example (including the associated html documentation!) and this doesn't appear to break any of our existing tests, so definitely 💯% approval there. Thanks for this very nice contribution!
|
Didn't really have time to check the math in detail, but the structure all seems to jive with what I remember. The main testing I'd done with Nedelec was to make sure I achieved the correct convergence rates in the different norms, so I would suggest the same be done with the example in the PR, namely for the example problem check the convergence rates of L2 and HDiv norms to the exact solution and make sure the expected rates are achieved. Definitely mimic @roystgnr's sentiment. Was always on my TODO list and never got there. :( So happy to someone pick that up and add it. |
Would love to AMR and Dirichlet BCs get properly tackled for these and Nedelec. It was indeed a lot of work, esp. handling 2D vs. 3D, was going to be a real PITA IIRC. |
Both the one I spotted, and some more egregious ones I missed. :-P Thanks! |
|
The It was also a bit annoying to play with that; rather than forcing some parameters to be loaded only from input file and others only from command line, why not just have the same GetPot object read the input file and then override on the command line? I've pushed some commits for those things to my fork of your branch, if you want to look at them. I'm happy with Paul's "check the convergence rates" test plan for now. Results: L2 errors
HDiv semi-norm errors
The H(div) errors look roughly like what I'd have expected, but not quite. Shouldn't we be leveling off closer to C*h by the time we have h=1/32? 50% of the error after each subdivision, not 55%? And that 55-ish% seems pretty consistent through the table, there's not some obvious higher-order terms that are disappearing. I'm afraid I'm ignorant of what I should expect in L2. Do we only expect O(h) here, not O(h^2)? I know the Nitsche lift "drop a derivative and add an h" trick doesn't always work even for scalar-valued elements ... thanks to stumbling upon an exception (Powell-Sabin-Heindl in H2) in grad school and then wasting weeks looking for a software bug where no software bug existed ... but I don't recall the rules for H(div) or H(curl) elements. |
Yeah, I noticed the same, we can lower the grid size as you suggest.
Well, yeah... I must admit I was bit confused by the situation w/ the Nedelec examples, so I just went with the same approach. I didn't realise the GetPot object allowed overrides via the command line.
I'll apply your two commits at the tip of my branch and push in a sec. Thank you.
Ah! You were quicker, I was arriving at exactly the same conclusions... This is for TRI6, the black circles are C*h, which is what I'd expect for both L2 and H(div).... But you seem to have shown the same deviation happens across all element types :/
According to this it should be O(h) for both (scroll to the pure Neumann on the scalar field, which equates to pure Dirichlet on the vector field, which is what I used). |
|
The natural thing to check first is always "where are we trying to cheat" and then "are we getting away with it" ... with that in mind, let's try modifying With triangles, |
I pushed this to my branch, but I'd say don't bother cherry-picking it unless we find something actually useful to push alongside. |
|
It works but I don't particularly like it... I'm now solving a problem defined only up to a constant... The solver somehow finds the solution we seek but I'm a bit uncomfortable... Maybe I should use the penalty method to fix one of the scalar field's dofs.... |
|
We were misconstraining p? I admit that boundary term wasn't what I'd have chosen, and in hindsight "evaluate a discontinuous variable at element boundaries" should have been an obvious code smell ... but surely we're now underconstrained instead? Is this example now going to scream and die the first time someone tries LU? (or even ILU? I often leave |
Yeah... I thought I was being smart by doing that...
Indeed.
Yes... |
|
I pushed a fix for the unused variables issues to my fork. And I've just noticed that |
|
I guess it wouldn't kill us to not be compatible with LU in a simple example; we're already hard-coding |
|
We should probably put some comment about the problem in the example code, though. |
nmnobre
force-pushed
the
rt
branch
3 times, most recently
from
5e71742 to
8cf1525
Compare
|
Looking better now. To fix it, I'm instead adding the constraint that the integral over the volume (not the surface) must match that of the known solution. :) I've updated the html docs to match. The plot for the tets is shorter just because I didn't want to wait, no other reason. It'd be nice if you could independently regenerate the table above. :)
|
|
In case the exact solution isn't known, would you then just constrain the scalar variable average value to be 0? This is sometimes done for the pressure in the case of enclosed incompressible flows. |
It is a problem if the exact solution is not known. I tried setting the average value to zero using the scalar Lagrange constraint in MOOSE. I don't believe it did anything to fix the issue. |
This is precisely the net effect of the constraint I'm putting in, as the scalar field integrates to zero over the domain. |
|
Regarding "setting the average value to zero": you can use a SCALAR variable for that kind of integral constraint, e.g. see Another example is |
Thanks @dknez, I hadn't looked into that, no. I think I spotted a couple of typos in that example, when repositioning the submatrices (I didn't know you could do that either, that's something I could've used here too!) between lines 241 and 249, whenever you multiply, it should be |
Nah; I'm happy enough as-is. I'm going to need to think harder about the problem before writing a RT-based mixed method myself, but this certainly looks adequate for demonstrating and testing the RT elements. Back in the day I would always just constrain away pressure kernels by adding a constraint pinning a single arbitrary point. :-D |
While I was waiting for you I actually implemented the example using the Lagrange multiplier. 😅 |
|
Either way is fine; whichever you prefer.
Or I need to just think harder about it because it's bugging the hell out of me. I can understand how trying to impose a Dirichlet BC on the boundaries of the discontinuous p variable could hurt the L2 convergence (constant, not rate) of p, adding a constant error term to it, because that's not the right way to get a high-order reconstruction there ... but I don't actually understand how that was hurting the H(div) convergence of u. We're solving u=-∇p, right? And ∇(constant) is still 0? None of this should block merging, my ignorance here is a me problem not a fe_raviart.C problem, but it's like an itch... |
I'd be in favor of having the Lagrange multiplier approach in the example (especially if you already coded it), since it seems like something that would be needed by anyone trying to extend the example. |
Right.
Right.
I think it's simply because I was doing this in a non-standard way. I'll try to explain my understanding. |
But the catch is that you're not just adding that one new equation, you're really adding one equation per boundary scalar DoF? And so then the solver doesn't just cancel the one-dimensional kernel, it's being forced to tweak every boundary value a bit? The error in p then isn't just a constant, it's actually got a (weak) gradient, and that leaks into u. That would explain it, if I've got it right. |
Yeah, that's it. It's not exactly one equation per boundary scalar dof, but I'm amending the existing equations with an extra term.
Exactly. (I think I've learnt my lesson, Lagrange multipliers it is from now on 😥) |
|
Lastly, there's something else I wanted to check with you. |
Nope. All you need is 5 hexes sharing an edge, and there's no way to get that by deforming any subset of a cartesian Hex mesh. If you need the deformation to also have a continuous inverse then the situation is even worse; just look at our And there's no way to set up a perfectly consistent single orientation rule in general, even if you can handle the high and low valence edge cases somehow. Imagine starting with a Nx1x1 row of hexes, giving it a 90 degree twist from one end to the other, then transforming from x,y,z to theta,r,z and stitching where the last hex meets the first to make it a twisted torus. There's no set of orientations you can give those hexes without there being some point at which the orientation has to locally twist. I was just discussing this on another ticket the other day - it's easy to split a hex into 6 tets, but the hex-to-tet code in libMesh splits them into 24 tets, because that way every face is symmetric and we don't have to worry about orientation, whereas the otherwise-more-reasonable 6-tet splits have to worry about whether there's any twists in the mesh that will lead to selection of two different diagonal edges from different sides of the same face. |
|
I don't know whether you want to do anything to modify the preconditioning. With the global saddle point, domain decomposition is not great; bjacobi/asm take > 400 linear iterations to converge. Interestingly mumps does not converge even with some shift on the diagonal. superLU and strumpack work however. The defaults do still work (albeit slowly they do converge) for the problem size, so I think what you have is fine |
I'd say leave it be; it's fast enough for a small scale example. I'd be curious as to what you'd recommend, though, for cases like this in general. On 16 cores, with 3375 hexes, I see 683 Krylov iterations with default block Jacobi + ILU0; |
|
I was just mentioning this PR to other coworkers and realized I skipped over some important bits:
It ought to be trivial. Add entries to a couple
You're in the UK, so 6 hours ahead of me in Austin? Early morning this Friday (anywhere 8-10am my time, so 2-4pm yours) would work for me for a telecon; otherwise I'll have more time most mornings next week. |
|
I've never done a Schur complement for a single global constraint. It would be interesting to try at some point. Since the off-diagonal blocks have the relation B = C^T probably common Schur complement preconditioning strategies would be applicable |
I think I can see this. Please correct me if I'm wrong, but you start by building what looks like the Schlegel diagram of a cube (2d) or a tesseract (3d) and, yes, for the latter, now seeing it as a pure 3d object, all 20 "internal" edges are valence 3, which I believe is what you are referring to?
I understand the construction of the twisted torus (the 90 degree rotation being along an axis parallel to the x axis). But I don't get the conclusion, can you please expand?
This is easy to understand I think, it could be the case that the diagonal edge on the common face of two neighbouring hexes didn't match in the 6-tet split case, but that won't happen with the 24-tet split. Got it. |
|
Correct about the sphere. For the torus, picture the eta-zeta axes on each cube (choosing xi to be the direction pointing from one cube to the next). For most of the cubes, eta and zeta on one cube will (on their shared face) match those on the neighbor, but you have to have one "twist" at some point, where e.g. eta on one cube is zeta on its neighbor and zeta on the cube is negative eta on its neighbor. Whether that matters for you depends on what you're doing, of course. Our C0 elements don't care; this new 6-tet split will scream and die. |
I think I was complicating. I misread/misunderstood what you said as meaning each cube would be rotated in relation to its neighbours, whereas you meant a single 90 degree rotation as seen from the end faces... My bad. All good! |
|
BTW I see O(h^2) convergence for the L2 norm in 'p'. I believe this is superconvergence |
Is this for the current mixed formulation or your new hybrid formulation? |
|
The former. Yes, seeing it for both BC types |
|
Interesting, do you have a branch I can play with? |
|
It turned out to be under-integration during the error calculation. Adding some extra integration order removed the superconvergence. I will be pushing some more commits to #3668 tomorrow that adds the calculation of L2 norms for 'p' |
8 participants
Hi @roystgnr, @jwpeterson, @lindsayad,
Christmas came earlier this year! :)
Hopefully it's not too hard to review.
Includes Raviart-Thomas basis functions for 2d and 3d, Piola mappings for H(div), an example that tests the functionality in both 2d and 3d for all supported element types (TRI6/7, QUAD8/9, and TET14 and HEX27), and a dump of the example output for the html docs.
We do have the required additions to MOOSE ready to go as well as a standalone app which tests the functionality. We're just working on tidying things up and in bringing the example under one of the MOOSE modules, most likely Navier-Stokes or Electromagnetics. We also noticed TET14s are missing from MOOSE but, at least at first sight, it doesn't look like that would be a massive amount of work?
Lastly, we do have plans (and funding) to do a bit more with libMesh and MOOSE beyond Raviart-Thomas. It'd be great if we could have a chat about these... it seems it'd be quite a lot of work, so we'd like to better understand the required effort, but also to get your opinion on how relevant these new features would be, how these align with your vision for the future of libMesh/MOOSE, etc....
Cheers,
Nuno and Karthik (cc @karthichockalingam)