Hi everyone!
I am currently implementing some stationary Stokes solvers based on step-55.
Therein, Taylor-Hood elements are being used. One can check the optimal
order of
convergence easily comparing to the Kovasznay or Poisuille flow solution,
the first one being already implemented in this step.
I went ahead and added PSPG and Bochev-Dohrmann stabilizations to check the
errors using those.
For the latter, i get the expected suboptimal (!) convergence rates (vel:2
and p~1.5-2).
Using PSPG, one adds
tau * residual momentum dot gradient(q)
(q being the pressure test function)
per element like:
local_matrix (i,j) += fe_values.JxW(q) * tau_sups * phi_grads_p [j] *
phi_grads_p [i];
if (ADD_THE_LAPLACIAN)
{
local_matrix (i,j) += fe_values.JxW(q) * tau * (
- viscosity * phi_laplacians_v [j]
) * phi_grads_p [i];
}
Using the laplacian of the velocity u defined as below:
if (get_laplace_residual)
{
phi_hessians_v[k] = fe_hessians[velocities].hessian(k,q);
for (int d = 0; d < dim; ++d)
{ phi_laplacians_v[k][d] = trace(phi_hessians_v[k][d]);
}
// Check laplacian, both give zero for rectangular grid and identical
values for distorted grids.
unsigned int comp_k = fe.system_to_component_index(k).first;
Tensor<1,dim> vector_laplace_v;
vector_laplace_v = 0;
if (comp_k<dim)
{
Tensor<2,dim> nonzero_hessian_k =
fe_hessians.shape_hessian_component(k,q,comp_k);
vector_laplace_v [comp_k] += trace(nonzero_hessian_k);
}
Tensor<1,dim> diff;
diff = phi_laplacians_v[k] - vector_laplace_v;
if (diff.norm() > 1e-6 || vector_laplace_v.norm() > 1e-16)
std::cout << "|divgrad_v| = " << std::setw(15) << vector_laplace_v.norm()
<< " , diff = " << diff.norm() << std::endl;
}
Which also includes a second option to compute the wanted laplacian.
Using a perfectly rectangular grid, the laplacian is actually zero for all
elements, thus everything reduces to just adding a pressure stiffness in
the pressure-pressure block.
Nonetheless, i observe convergence rates of vel:2 and p:1.5-1.8 which is
suboptimal and not(!) expected.
I implemented the same methods in a matlab inhouse code & observed perfect
convergence rates of 2&2 for v&p using a direct solver.
The stabilization parameter is in both cases (matlab or deal.ii) chosen as
tau = 0.1*h^2/viscosity
with
double h = std::pow(cell->measure(), 1.0/ static_cast<double>(dim))
or h = cell->diameter() giving similar results.
Could the observed behaviour be caused by applying an iterative solver*?
(using ReductionControl and a reduction factor of 1e-13 which is really low)
* Block-triangular Schur-complement-based approach, similar as in step-55
suggested in the further possibilities for extension, basically using
S~-(1/viscosity)*Mass_pressure giving 20-40 iterations up to 3mio dofs
tested.
I have been checking the code for a week now, and i really doubt, that the
integration of just tau*grad(p)*grad(p) is wrong (again, the laplace drops
out for a rectangular grid).
Just to check, i used a manufactured solution from Poisuille flow and added
the laplacian from PSPG to the rhs, giving me the exact (linear) solution
in the pressure and quadratic convergence in the velocity, thus I assume,
that the grad(p)grad(q) term added is correctly implemeted. (exact solution
meaning, that i get an L2 error of err<1e-9 for any number of elements used)
Anyone has ever experienced this or has anyone some tipps for further
debugging?
Thanks for taking the time to read the post, any help or comment would be
greatly appreciated!
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/2ded867b-a2e7-4041-b90c-2dd774f1f5dd%40googlegroups.com.
