residual for a static condensation problem

[Yann Jobic]

Dear Dealii group,

I am currently learning Dealii (and i really like this lib). I would
like to bench it in order to compare its performance with an in house
FEM code at my lab.

I want to solve the stokes equation with a static condensation of the
pressure over the velocity dof. I started from the step-22.cc example.
It's working fine. I also computed the "pressure", to within one
constant. I put the code in attachment.

However, i don't know how to compute the residual ! It should be easy,
but i still haven't succeeded. If I extract the different submatrices of
the system matrix at the element level, i don't have the right type of
{U} to compute {P}:
{P} = -[Kpp]^{-1}*[Kpu]*{U}
The matrices are in FullMatrix format, and {u} is of type
std::vector<Tensor<1, dim>> if I'm using "get_function_values". Thus i
can't compute [Kpu]*{U}. I may have the nodal velocity values, but they
have to be ordered as [Kup], [Kpu], ...

Can you give me some advice on the correct way to solve this problem?

Thanks,

Yann

PS :
1) {U} : column vector U
2) I'm using a penalisation method for the pressure, by defining div u =
\lambda p, \lambda the penalisation parameter. Thus the weak form looks
like (v shape function of u, and phi shape function of p),
2*a(grad v, grad u)-a(div v,p)-a(phi,div u)+a(phi,\lambda p) = l(v*f)
then
{P}=-[Kpp]^{-1}[Kpu]*{U}
and finally
([Kuu]-[Kup][Kpp]^-1[Kpu]){U}={f}

[Wolfgang Bangerth]

On 10/14/24 14:03, Yann Jobic wrote:
>
> However, i don't know how to compute the residual ! It should be easy,
> but i still haven't succeeded. If I extract the different submatrices of
> the system matrix at the element level, i don't have the right type of
> {U} to compute {P}:
> {P} = -[Kpp]^{-1}*[Kpu]*{U}
> The matrices are in FullMatrix format, and {u} is of type
> std::vector<Tensor<1, dim>> if I'm using "get_function_values". Thus i
> can't compute [Kpu]*{U}. I may have the nodal velocity values, but they
> have to be ordered as [Kup], [Kpu], ...

Yann -- the term "residual" is used in different meanings in numerical
analysis. Can you define precisely what it is you want to compute?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Yann Jobic]

My next step is to compute that way (static condensation of the
pressure) the navier-stokes equations with a Newton non linear solver,
as in step-57, but i do it on a simpler problem first (stokes).
Thus if i define
{R(u,p)} = {-2*div u + \nabla p -f},
I want to compute ||R(u,p)||.

Many thanks,

Yann

[Wolfgang Bangerth]

On 10/15/24 01:19, Yann Jobic wrote:
>
> My next step is to compute that way (static condensation of the
> pressure) the navier-stokes equations with a Newton non linear solver,
> as in step-57, but i do it on a simpler problem first (stokes).
> Thus if i define
> {R(u,p)} = {-2*div u + \nabla p -f},
> I want to compute ||R(u,p)||.

Yann,
so the residual for you is a function of x, and what you're looking to compute
is an integral (for the norm). Integrals in deal.II are always computed using
FEValues, and what you want here is going to look something like the following:

double integral = 0;
for (cell=...)
{
fe_values.reinit(cell);
fe_values[velocities].get_function_divergences(solution, div_u);
fe_values[pressure].get_function_gradients(solution, grad_p);
rhs.value_list(fe_values.quadrature_points(), f_values);

for(q=...)
integral += std::pow(-2*div_u[q] + grad_p[q] - f_values[q], 2) *
fe_values.JxW(q);
}
double norm = std::sqrt(integral);

(Your formula for the residual isn't correct because div u is a scalar and
you're adding grad p to it, which is a vector. This can't be right. I'll let
you fix this, and just wrote the code matching your formula.)

You might want to take a look at step-15 to see some of the parts I left out
above.

[Yann Jobic]

Many thanks for your answer. I may not have all the terms of the static condensation in your solution, as i made some mistakes. Sorry.

My email is way too long! So the short question is:

In order to have nodal values of the solution vector, i've got to work with the type "std::vector<Tensor<1,dim>>" which is difficult to handle with the "FullMatrix<double>" of the extracted submatrix of the system matrix at the element level, as i have to compute [Kpu]{u}. Is it possible to transform the std::vector<Tensor<1,dim>> to the std::vector<double> type ?

The detailed but i hope comprehensive  explanation of why i want that:

The way i compute the static condensation of the pressure at the element level is as follow:
*) the weak form is 2*a(grad v, grad u)-a(div v,p)-a(phi,div u)+a(phi,\lambda p) = a(v,f), lambda being small
*) I create the system matrix [K] at the element level.
*) I extract [Kuu], [Kpu], [Kup], and [Kpp] as FullMatrix, and condense the pressure (I create A=-[Kup][Kpp]^{^-1}[Kpu]) and scatter the values to the velocity dof)

Then, i would like to compute the elementary residual at each cell <R(u),v>, that is to say compute
For each cell, {R(u)} = 2*a(grad v, grad u) - a(v,f) - [Kup][Kpp]^{^-1}[Kpu]{u}
So i've got to, {u} given:
*) compute {R(u)}=2*a(grad v, grad u) - a(v,f) at gauss points
*) compute [Kup][Kpp]^{^-1}[Kpu]{u} at solution node values
*) subtract the latter to take into account the condensed terms, then distribute_local_to_global
*) compute the L2 error of the resulting vector.

The problem is that i can have the nodal values of {u} with the following code:

const std::vector<Point<dim>> support_points = fe.base_element(0).get_unit_support_points();
std::vector<double>        weights(support_points.size(),1);
const Quadrature<dim> node_quadrature_formula(support_points,weights);
FEValues<dim> fe_node_values (fe, node_quadrature_formula,
        date_values   | update_gradients |
        date_quadrature_points | update_JxW_values);

But then i will work with the type "std::vector<Tensor<1,dim>>" which is difficult to handle with the "FullMatrix" of the extracted submatrix, as i have to compute [Kpu]{u}. Then how can i compute [Kpu]{u} ? Is it possible to transform the std::vector<Tensor<1,dim>> to the std::vector<double> type ? Is it the way to do the static condensation at the element level ? (step-44 does it but not at the element level)

Do that make more sense ? It could still be confusing.

Many thanks,

Best regards,

Yann

[Wolfgang Bangerth]

Yann,
you did not read my email. The approach to compute the norm of the
residual goes through *integration with FEValues*, not with computing
nodal values. You are stuck thinking about how the approach you are
trying to use can be made to work; but the approach you are using is not
the right approach.

> Then, i would like to compute the elementary residual at each cell
> <R(u),v>, that is to say compute
> For each cell, {R(u)} = 2*a(grad v, grad u) - a(v,f) - [Kup][Kpp]^-1[Kpu]{u}

This formula does not make sense. The terms
a(grad v, grad u)
and
a(v,f)
are numbers (they are the result of integration) whereas

[Kup][Kpp]^-1[Kpu]{u}

is a vector. You cannot add numbers and vectors.

> std::vector<double>        weights(support_points.size(),1);
> const Quadrature<dim> node_quadrature_formula(support_points,weights);
> FEValues<dim> fe_node_values (fe, node_quadrature_formula,
>         date_values   | update_gradients |
>         date_quadrature_points | update_JxW_values);
>
> But then i will work with the type "std::vector<Tensor<1,dim>>" which is
> difficult to handle with the "FullMatrix" of the extracted submatrix, as
> i have to compute [Kpu]{u}.

You have a conceptual misunderstanding. The std::vector<Tensor<1,dim>>
are the vector-valued values of the velocity at the evaluation points.
But when you want to compute...

> Then how can i compute [Kpu]{u} ?

...the {u} are the *nodal values* of the velocity. These are not the
same. It is not a problem of how these two things are represented in
data structures, but that they are conceptually different. You simply
cannot multiply [Kpu] by a collection of velocity vectors; [Kpu] needs
to be multiplied by a vector of nodal values.

Best
W.