Non-homogeneus boundary conditions with-matrix-free not working

[Michał Wichrowski]

Dear deal.II developers,

I think I found problem with non-homogeneus boundary conditions.  

By changing:

VectorTools::interpolate_boundary_values (dof_handler,

0,

Functions::ZeroFunction<dim>(),

constraints);

to

VectorTools::interpolate_boundary_values (dof_handler,

0,

Functions::ConstantFunction<dim>(5),

constraints);

The results shown on included picture were outputted. I've tried to figure out what is going on, and it looks like matrix-free framework is solving problem with homogeneous constrains and then Dirichlet nodes are overwritten.  

[Martin Kronbichler]

Dear Michal,

This is expected: the matrix-free operator evaluation cannot apply non-homogeneous boundary conditions while solving (at least I have never figured out how to do that). To solve such a problem, you need to bring the non-homogeneous part on the right hand side first. I often solve this by first creating a vector that is filled with the Dirichlet values and apply the operator on that (without setting the zero constraints) and then solve for an increment that has zero boundary conditions. I attach an example which is an extension of step-37 (and uses coefficients described here, if you want some more background: http://www.it.uu.se/research/publications/reports/2017-006/2017-006-nc.pdf

In the program, you will find that LaplaceOperator not only contains a vmult() function, but also a compute_residual() function that reads the data in from an FEEvaluation object without Dirichlet conditions, due the operation, and assembles the right hand side into an FEEvaluation object with Dirichlet conditions (-> local_compute_residual()).

Please let me know in case you have questions.

Best,
Martin

--
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+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bruno Turcksin]

Martin,

Could you add this explanation to the documentation and/or to a tutorial. I think mentioning it step-37 would be great.

Best,

Bruno

To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscribe@googlegroups.com.

[Martin Kronbichler]

Hi Bruno,

I definitely agree - we should mention this in step-37. I will open an issue since I will not get to do it today.

Best,
Martin

-- 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+un...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

On 10/27/2017 06:44 AM, Bruno Turcksin wrote:
>
> Could you add this explanation to the documentation and/or to a tutorial. I
> think mentioning it step-37 would be great.

Yes, I was actually thinking the same. Maybe in "Possibilities for extensions"
in the results section of step-37, like we have for other programs?

Cheers
W.

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

[Michał Wichrowski]

That's very bad news for me. This workaround will work for Laplace problem but not for my case. I'm dealing with Stokes problem and my method requires that second block of right hand side is zero, otherwise MinRes crashes and GMRES have to be used (that is far less effective). 

Can You give me some details what is the problem? 

[Michał Wichrowski]

I think that the problem can be avoided by splitting operator A into two parts: dirichlet nodes and other. The structure will be:

A1 0

0    I

and the operator A1 applies matrix-vector product using fixed values on Dirichlet nodes.

If I got it right, the problem is with this part of code:

Line:  1186 of  operatrors.h:

    // set zero Dirichlet values on the input vector (and remember the src and

    // dst values because we need to reset them at the end)

    for (unsigned int j=0; j<n_blocks(dst); ++j)

      {

        for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)

          {

            edge_constrained_values[j][i] =

              std::pair<Number,Number>

              (subblock(src,j).local_element(edge_constrained_indices[j][i]),

               subblock(dst,j).local_element(edge_constrained_indices[j][i]));

            subblock(const_cast<VectorType &>(src),j).local_element(edge_constrained_indices[j][i]) = 0.;

          }

      }

By replacing 0 with Dirichlet boundary value and then applying the operator without Dirichlet constrains the multiplication by A1 can be done. Multiplication by identity for Dirichlet nodes is already implemented. 

This will require some changes in matrix-free framework. MatrixFree object will handle only homogeneous constrains and Operator will handle the non-homogeneous one.  The same change have to be done in vmult_interface_down  and  vmult_interface_up.

[Martin Kronbichler]

Dear Michal,

You mean because once you apply an inhomogeneous Dirichlet condition on
the velocity you also get a contribution because the divergence is not
zero? You could work around that by applying a function in the whole
domain that satisfies the Dirichlet constraints and is divergence free,
in which case you should not get a contribution to right right hand side
of the divergence equation but only the momentum equation. On the other
hand, I'm not sure how easy it would be to construct such a function...

Regarding the workaround: It could work by the modification you
suggested, but we cannot really fix that in the library because the
operators are meant for linear operators and not affine ones. The
problematic issue is that there are some cases where you want to apply
an operator with a non-zero constraint and others where you don't. For
example, all Krylov solvers expect you to provide a matrix-vector for
the linear operator, not the affine one as you would get when you have
inhomogeneous constraints inside the mat-vec. In particular, I do not
really understand the part in your second post where you talk about the
vmult_interface_{up,down} functions because those are only used for
multigrid where the library always assumes to work on homogeneous problems.

I have not gone down this route and as far as I know, nobody else has
previously - including the matrix-based solvers we have available in
deal.II whose way of operation resembles the setup I described in my
previous post, but do it inside the
ConstraintMatrix::distribute_local_to_global call - so I cannot say how
much work that would be. We would be happy to provide a solution that is
generic if you find one, but I'm not sure I will be able to help much
along this direction.

Best,
Martin

[Michał Wichrowski]

W dniu sobota, 28 października 2017 18:52:57 UTC+2 użytkownik Martin Kronbichler napisał:

Dear Michal,

You mean because once you apply an inhomogeneous Dirichlet condition on
the velocity you also get a contribution because the divergence is not
zero?

 Exactly.

You could work around that by applying a function in the whole
domain that satisfies the Dirichlet constraints and is divergence free,
in which case you should not get a contribution to right right hand side
of the divergence equation but only the momentum equation. On the other
hand, I'm not sure how easy it would be to construct such a function...

Constructing such a function is not an option (for now), it is too expensive. 

One question about that way of imposing Dirichlet boundary conditions. I've modified step-37 following the same idea, but in different way:

  template <int dim>

  void LaplaceProblem<dim>::solve ()

  {

  SystemMatrixType system_matrix_no_dirichlet;

    {

      typename MatrixFree<dim,double>::AdditionalData additional_data;

      additional_data.tasks_parallel_scheme =

        MatrixFree<dim,double>::AdditionalData::none;

      additional_data.mapping_update_flags = (update_gradients | update_JxW_values |

                                              update_quadrature_points);

      std::shared_ptr<MatrixFree<dim,double> >

      system_mf_storage(new MatrixFree<dim,double>());

      system_mf_storage->reinit (dof_handler, constraints_without_dirichlet, QGauss<1>(fe.degree+1),

                                 additional_data);

      system_matrix_no_dirichlet.initialize (system_mf_storage);

    }

    system_matrix_no_dirichlet.evaluate_coefficient(Coefficient<dim>());

//

    constraints.distribute(solution);

    LinearAlgebra::distributed::Vector<double> residual=system_rhs;

    system_matrix_no_dirichlet.vmult(residual, solution);

    system_rhs -= residual;

    LinearAlgebra::distributed::Vector<double> solution_update=system_rhs;

    solution_update=0;

......

    cg.solve (system_matrix, solution_update, system_rhs,

              preconditioner);

    solution+= solution_update;

    constraints.distribute(solution);

In this way I don't have to modify my operator, bum I'm not sure if it will work for all cases. For now it works and the results looks OK.

 

Regarding the workaround: It could work by the modification you
suggested, but we cannot really fix that in the library because the
operators are meant for linear operators and not affine ones. The
problematic issue is that there are some cases where you want to apply
an operator with a non-zero constraint and others where you don't. For
example, all Krylov solvers expect you to provide a matrix-vector for
the linear operator, not the affine one as you would get when you have
inhomogeneous constraints inside the mat-vec. In particular, I do not
really understand the part in your second post where you talk about the
vmult_interface_{up,down} functions because those are only used for
multigrid where the library always assumes to work on homogeneous problems.

You're right,  forget about what I written previously. I've totally misunderstood a few things.

[Michał Wichrowski]

Dear Martin,

I finally managed to make it work by following Your idea. Thanks!

Michał