Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-30 Thread Michał Wichrowski 
Dear Martin,
I finally managed to make it work by following Your idea. Thanks!

Michał 

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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-28 Thread 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 
  void LaplaceProblem::solve ()
  {
  SystemMatrixType system_matrix_no_dirichlet;
{
  typename MatrixFree::AdditionalData additional_data;
  additional_data.tasks_parallel_scheme =
MatrixFree::AdditionalData::none;
  additional_data.mapping_update_flags = (update_gradients | 
update_JxW_values |
  update_quadrature_points);
  std::shared_ptr >
  system_mf_storage(new MatrixFree());
  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());

//
constraints.distribute(solution);
LinearAlgebra::distributed::Vector residual=system_rhs;
system_matrix_no_dirichlet.vmult(residual, solution);
system_rhs -= residual;
LinearAlgebra::distributed::Vector 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.
 

> 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 
>

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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-28 Thread 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

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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-28 Thread 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
0I
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
  
(subblock(src,j).local_element(edge_constrained_indices[j][i]),
  
 subblock(dst,j).local_element(edge_constrained_indices[j][i]));
subblock(const_cast(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.


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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-28 Thread 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? 

W dniu piątek, 27 października 2017 15:02:29 UTC+2 użytkownik Wolfgang 
Bangerth napisał:
>
> 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/ 
>
>

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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-27 Thread 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: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-27 Thread 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


On 27.10.2017 14:44, Bruno Turcksin wrote:
> 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
>
> On Friday, October 27, 2017 at 3:57:09 AM UTC-4, Martin Kronbichler
> wrote:
>
> 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
>
> On 26.10.2017 22:07, Michał Wichrowski wrote:
>> 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(),
>> constraints);
>> to
>> VectorTools::interpolate_boundary_values (dof_handler,
>> 0,
>> Functions::ConstantFunction(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.  
>> -- 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
>> . For more options,
>> visit https://groups.google.com/d/optout
>> . 
>
> -- 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
> . For more options, visit
> https://groups.google.com/d/optout. 

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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-27 Thread 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

On Friday, October 27, 2017 at 3:57:09 AM UTC-4, Martin Kronbichler wrote:
>
> 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
> On 26.10.2017 22:07, Michał Wichrowski wrote:
>
> 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(),
> constraints);
> to
> VectorTools::interpolate_boundary_values (dof_handler,
> 0,
> Functions::ConstantFunction(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.  
> -- 
> 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.
> For more options, visit https://groups.google.com/d/optout.
>
>
>

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

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

2017-10-27 Thread 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

On 26.10.2017 22:07, Michał Wichrowski wrote:
> 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(),
> constraints);
> to
> VectorTools::interpolate_boundary_values (dof_handler,
> 0,
> Functions::ConstantFunction(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.  
> -- 
> 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
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2009 - 2017 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Authors: Katharina Kormann, Martin Kronbichler, Uppsala University,
 * 2009-2012, updated to MPI version with parallel vectors in 2016
 */


#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 


namespace Step37
{
  using namespace dealii;


  const unsigned int degree_finite_element = 3;
  const unsigned int dimension = 3;

  typedef double number;
  typedef float level_number;


  // ==
  // Functions
  // ==

  template 
  class Solution : public Function
  {
  public:
Solution ()  : Function() {}

virtual double value (const Point   ,
  const unsigned int  component = 0) const;

virtual Tensor<1,dim> gradient (const Point   ,
const unsigned int  component = 0) const;


virtual Tensor<1,dim,VectorizedArray> gradient (const 
Point  ,
const unsigned int  
component = 0) const;

virtual double laplacian (const Point   ,
  const unsigned int  component = 0) const;

virtual VectorizedArray laplacian (const 
Point   ,
   const unsigned int  component = 
0)