Use a manifold to generate better meshes for the shell. #242
Conversation
|
(not ready to merge of course) |
| for (std::set<types::boundary_id>::iterator it = ids.begin(); | ||
| it!=ids.end(); ++it) | ||
| if (*it > 1) | ||
| coarse_grid.set_boundary (*it, straight_border); |
There was a problem hiding this comment.
This is the place where I think a comment on why you're doing that (and which parts of the boundary it affects) is in place. I believe the word "border" only applies to the edge of a 2d object. Maybe call it straight_boundary?
|
Looks good with the one added piece of documentation mentioned in my last comment. |
|
The problem was that your code crashed if you were running without a full shell: |
|
I also had to change the manifold id from 2 to something different, because this conflicted with the straight boundaries. Not sure if this is by design. |
|
Sadly we have two tests that now crash: |
|
Do you have any suggestions? |
tjhei
force-pushed
the
pull_229
branch
2 times, most recently
from
3b43213
to
7df516f
Compare
|
it turns out that the two tests don't crash if we use MappingQ also in the interior! |
|
Wow, who knew? How many iterations did you need for the first time step? |
a few more than with the old mesh: |
|
Sadly depth_average_01 (and 02 which is the same computation) are crashing now: What is weird about this is that the mesh is a unit cube refined to be 4x4 square cells so we should have a FlatManifold and it shouldn't make a difference. Switching back to the plain mapping in the interior or reducing the mapping from Q4 to Q1 fixes the problem. Do you have any suggestion, @bangerth ? It might be a bug in deal.II, so maybe I need to sit down and construct a testcase... |
|
I created a deal.II test and I can not find a difference between true and false for the interior mapping: mapping quadrature points, shape_value, shape_gradient of FE_Q, and JxW are identical. Any other ideas? |
|
That's too crazy. Looking into it. |
|
I'm seeing changes on the order of roundoff in the matrix. Need to look further, I assume... |
|
So I now tried to circumvent the problem by using ARGH! |
|
So where do we stand now? I see differences on the order of roundoff in the matrix between using the higher order mapping everywhere or not. I think that's fair. I did figure out which solver actually fails: it's the mass matrix solver (bottom right block) for which we use Trilinos's implementation of CG. I fail to see why this can happen -- the viscosity for this model is constant one. Can you summarize what works and what doesn't? If I understand correctly, then this is the case:
|
Maybe the matrix is not exactly SPD?
A) use |
|
I think we should keep We could try our CG, which is probably more robust, or test with GMRES to see if that is the problem. |
|
Even GMRES fails for the velocity block in periodic-box.prm. The example is periodic in x and has tangential conditions at the top/bottom so this is not too surprising. That the order of the mapping matters is just weird, though. |
|
Yes, I forgot to mention this as well that using GMRES doesn't make a difference. I output the matrix before I left but didn't get to look at it. What could make it singular? It's a mass matrix on the pressure, after all, not a Laplace matrix on the velocities. I would like to understand what exactly is the problem here -- if it worked before, then it seems like that was only by accident. I agree that we should keep |
|
By the way, I use this patch to get the matrix in the step that fails: The matrix is only 25x25, so I plan on just putting it into maple or matlab or something and compute the spectrum. For reference, this is the matrix that fails: Anyone good with matlab and wants to compute whether the matrix is SPD? |
|
Seems to be SPD to me, unless I've made a mistake. Eigenvalues of that matrix: I don't have a good intuition for the matrix, is the degeneracy expected? |
|
No, it's not supposed to be degenerate. In fact, the eigenvalues look just right to me. On the other hand, we work on a square which has four-fold symmetry, so I would expect that there are four-fold degenerate eigenvalues, but I don't see any. Here's another observation: the Trilinos solver only fails in this time step if the |
I already did that yesterday: tjhei@b18093e |
|
My last comment corresponds to the has_curved_cell patch. |
|
okay, periodic_box.prm doesn't crash in the A solve when I change the coarse grid solver of AMG from KLU to "ML symmetric Gauss-Seidel". That is probably a bad idea for performance, but this shows us what the problem is: |
|
Wait, what? I missed the periodic-box thing. What's going wrong there? Is the A matrix singular for this case? I think we really do want the direct KLU solver on the coarse mesh, though, no? |
I think so, yes. I just did this to see what is happening in this test. |
|
@tjhei Regarding the KLU: From the ML guide, I can't see a way to set coarse solver parameters with respect to that (other than defining one's own solver). On the other hand, I can't think of KLU to not do the 'right' thing when it gets a singular matrix and failing only in case the incoming vector is not in the range of the matrix. |
|
@kronbichler : The case @tjhei asks about is for the Laplace matrix which is singular for this case. The mass matrix is a different testcase. |
|
@kronbichler No, if I play around with the MappingQ a different test fails. The velocity-velocity block Trilinos::SolverCG with AMG fails with If I switch to a different coarse solver, it works. |
|
I've isolated the depth_average failure into this deal.II testcase: I think we should merge this PR (after some cleanup of the commit history) to mainline and deal with the remaining periodic-box at a later time. |
tjhei
force-pushed
the
pull_229
branch
2 times, most recently
from
81544cd
to
af3a54b
Compare
This only works with a sufficiently new deal.II, of course.
- set straight boundary for quarter/half shell walls - only set boundary objects as needed for spherical shell - remove deal.II version check - use higher order mapping also in the interior
|
I cleaned up the commits. Are you okay with the code changes? Before we merge I want to a) read through the test changes, and b) wait for the tester. Tomorrow... |
|
Yea, I like it. Thanks for working it over -- merge whenever the tests are to your liking! |
Use a manifold to generate better meshes for the shell.
|
What an ordeal :-) |
This is a rework from #229, rebased to master so I can fix some problems.