Hi everyone,
I am back with good & bad news!
- Good news first, I implemented a parallel direct solver (mumps via petsc)
and going from a 2x2 blocked system to
a regular one (actually 1x1 still blocked system)
one just needs to adapt the setup phase, not reordering per component & not
use the "get_view"
but rather use all the same vector and matrix (types) with 1 block instead
of 2. Once you realize that, it takes ~15 mins for a
non-experienced deal.ii user like me.
- Bad news is, that I still do not obtain the optimal convergence rate in
the pressure, i.e., 2 for Q1Q1 PSPG-stabilized Stokes.
The obtained rates in the Poisuille benchmark (quadratic velocity profile &
linear pressure) using the direct solver are:
L2-errors in v & p: Q1Q1 + PSPG
--------------------------------------------------------------------------------------------------
lvl 0: | e_v = 1 | eoc_v = 0 | e_p = 1 |
eoc_p = 0 |
--------------------------------------------------------------------------------------------------
lvl 1: | e_v = 0.531463 | eoc_v = 0.91196 | e_p = 0.479585 |
eoc_p = 1.06014 |
--------------------------------------------------------------------------------------------------
lvl 2: | e_v = 0.218196 | eoc_v = 1.28434 | e_p = 0.208925 |
eoc_p = 1.1988 |
--------------------------------------------------------------------------------------------------
lvl 3: | e_v = 0.0661014 | eoc_v = 1.72287 | e_p = 0.0674415 |
eoc_p = 1.63128 |
--------------------------------------------------------------------------------------------------
lvl 4: | e_v = 0.0176683 | eoc_v = 1.90352 | e_p = 0.0194532 |
eoc_p = 1.79363 |
--------------------------------------------------------------------------------------------------
lvl 5: | e_v = 0.00452588 | eoc_v = 1.96489 | e_p = 0.00549756 |
eoc_p = 1.82315 |
--------------------------------------------------------------------------------------------------
lvl 6: | e_v = 0.00114238 | eoc_v = 1.98616 | e_p = 0.00158448 |
eoc_p = 1.79478 |
--------------------------------------------------------------------------------------------------
lvl 7: | e_v = 0.000286765 | eoc_v = 1.9941 | e_p = 0.000475072 |
eoc_p = 1.73779 |
--------------------------------------------------------------------------------------------------
lvl 8: | e_v = 7.1824e-05 | eoc_v = 1.99733 | e_p = 0.000149148 |
eoc_p = 1.67141 |
--------------------------------------------------------------------------------------------------
lvl 9: | e_v = 1.79717e-05 | eoc_v = 1.99874 | e_p = 4.88044e-05 |
eoc_p = 1.61166 |
--------------------------------------------------------------------------------------------------
lvl 10: | e_v = 4.49483e-06 | eoc_v = 1.99939 | e_p = 1.64706e-05 |
eoc_p = 1.56712 |
--------------------------------------------------------------------------------------------------
where e_v & e_p are relative errors, i.e., int(p-ph)dx/int(p)dx
The obtained rates using stable Taylor-Hood Q2Q1 are:
L2-errors in v & p: Q2Q1
--------------------------------------------------------------------------------------------------
lvl 0: | e_v = 2.20205e-16 | eoc_v = 0 | e_p = 3.58011e-16 |
eoc_p = 0 |
--------------------------------------------------------------------------------------------------
lvl 1: | e_v = 2.1471e-16 | eoc_v = 0.0364598 | e_p = 1.02455e-15 |
eoc_p = -1.51692 |
--------------------------------------------------------------------------------------------------
lvl 2: | e_v = 4.23849e-16 | eoc_v = -0.98116 | e_p = 4.1472e-16 |
eoc_p = 1.30479 |
--------------------------------------------------------------------------------------------------
lvl 3: | e_v = 6.9843e-16 | eoc_v = -0.720565 | e_p = 4.49311e-15 |
eoc_p = -3.4375 |
--------------------------------------------------------------------------------------------------
lvl 4: | e_v = 1.53993e-15 | eoc_v = -1.14068 | e_p = 9.82242e-15 |
eoc_p = -1.12837 |
--------------------------------------------------------------------------------------------------
lvl 5: | e_v = 7.60182e-15 | eoc_v = -2.30348 | e_p = 5.44953e-14 |
eoc_p = -2.47198 |
--------------------------------------------------------------------------------------------------
lvl 6: | e_v = 1.60118e-14 | eoc_v = -1.07472 | e_p = 1.86328e-13 |
eoc_p = -1.77364 |
--------------------------------------------------------------------------------------------------
... which shows, that I get the exact solution, which is expected.
Using Q2Q2 + PSPG i get:
L2-errors in v & p: Q2Q2 + PSPG
--------------------------------------------------------------------------------------------------
lvl 0: | e_v = 3.5608e-16 | eoc_v = 0 | e_p = 6.5321e-16 |
eoc_p = 0 |
--------------------------------------------------------------------------------------------------
lvl 1: | e_v = 7.51502e-16 | eoc_v = -1.07758 | e_p = 1.59993e-15 |
eoc_p = -1.29239 |
--------------------------------------------------------------------------------------------------
lvl 2: | e_v = 7.47141e-16 | eoc_v = 0.00839586 | e_p = 7.92418e-16 |
eoc_p = 1.01367 |
--------------------------------------------------------------------------------------------------
lvl 3: | e_v = 1.36773e-15 | eoc_v = -0.872327 | e_p = 4.87009e-15 |
eoc_p = -2.61961 |
--------------------------------------------------------------------------------------------------
lvl 4: | e_v = 2.16714e-15 | eoc_v = -0.664013 | e_p = 2.01012e-14 |
eoc_p = -2.04526 |
--------------------------------------------------------------------------------------------------
lvl 5: | e_v = 1.14488e-14 | eoc_v = -2.40134 | e_p = 2.04823e-14 |
eoc_p = -0.0270959 |
--------------------------------------------------------------------------------------------------
lvl 6: | e_v = 1.15543e-14 | eoc_v = -0.0132335 | e_p = 1.7867e-13 |
eoc_p = -3.12484 |
--------------------------------------------------------------------------------------------------
also (!) giving me the exact solution for all meshes used, but actually one
has to use the PSPG here.
The only suspicious thing I found, is that when I use in the preparation of
the testfunctions:
if (get_laplace_residual)
{
Tensor<3,dim> phi_hessian_v;
phi_hessian_v =
fe_hessians[velocities].hessian(k,q);
for (int d = 0; d < dim; ++d)
{
phi_laplacians_v[k][d] =
trace(phi_hessian_v[d]);
}
if (phi_hessian_v.norm() > 1e-14)
std::cout << "|H| = " << std::scientific <<
phi_hessian_v.norm() << std::endl;
}
There is actually a non-zero norm printed during execution. Linear shape
functions on a rectangular(!) grid should
despite being subject to a bilinear mapping still give zero second
derivatives, right?
Since my understanding of the hessian() function is somewhat limited, I
also tried with
if (get_laplace_residual)
{
phi_hessian_v = 0;
phi_hessian_v =
fe_hessians[velocities].hessian(k,q);
for (int d = 0; d < dim; ++d)
{
phi_laplacians_v[k][d] =
trace(phi_hessian_v[d]); // ###
if (fe.system_to_component_index(k).first < dim)
phi_laplacians_v[k][d] =
trace(fe_hessians.shape_hessian_component(k,q,d));
else
phi_laplacians_v[k][d] = 0.0;
Tensor<1,dim> diff;
diff = 0;
diff = trace(phi_hessian_v[d]) -
trace(fe_hessians.shape_hessian_component(k,q,d));
if (diff.norm () > 1e-14)
std::cout << "|d| = " << std::scientific <<
diff.norm() << std::endl;
}
if (phi_hessian_v.norm() > 1e-14)
std::cout << "|H| = " << std::scientific <<
phi_hessian_v.norm() << std::endl;
}
Showing, that I get same (correct?) results using any function to compute
the laplacian, since nothing is printed to screen when executing.
Can someone confirm, that the use of these 2 functions is correct?
I (almost) use the same call as the proposed code 'lethe'.
Does the Q2Q2+PSPG not show, that the implementation is correct? Is the
(internal) treatment of linear shape functions different?
The non-zero second derivatives are still somewhat suspicious to me.
If all else fails, I am gonna strip off all the extra stuff I have in my
flow solver and post a simple Q(p)xQ(q) + pspg.
Kind regards & greetings from Austria,
Richard
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/9ad2e869-da5a-4de6-be12-55d19aa92d5b%40googlegroups.com.
