Re: [deal.II] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[shanglong zhang]

Hi Wolfgang,

Thank you for the explanation. Then everything make sense.

Best,
Shanglong

在 2016年7月27日星期三 UTC-4下午6:13:16，bangerth写道：
>
> On 07/27/2016 01:43 PM, Alexander wrote: 
> > 
> > step-17 claims to solve same problem as step-8. If one produces symmetry 
> > matrix and the other does not or both do not -- this is suspicious. 
>
> The reason why the matrix in step-17 is not symmetric is actually 
> explained in the code: 
>
>  // The last argument to the call to 
>  // MatrixTools::apply_boundary_values() below allows for some 
>  // optimizations. It controls whether we should also delete 
>  // entries (i.e., set them to zero) in the matrix columns 
>  // corresponding to boundary nodes, or to keep them (and passing 
>  // true means: yes, do eliminate the columns). If we 
>  // do eliminate columns, then the resulting matrix will be 
>  // symmetric again if it was before; if we don't, then it 
>  // won't. The solution of the resulting system should be the same, 
>  // though. The only reason why we may want to make the system 
>  // symmetric again is that we would like to use the CG method, 
>  // which only works with symmetric matrices. The reason why we may 
>  // not want to make the matrix symmetric is because this 
>  // would require us to write into column entries that actually 
>  // reside on other processes, i.e., it involves communicating 
>  // data. This is always expensive. 
>  // 
>  // Experience tells us that CG also works (and works almost as 
>  // well) if we don't remove the columns associated with boundary 
>  // nodes, which can be explained by the special structure of this 
>  // particular non-symmetry. To avoid the expense of communication, 
>  // we therefore do not eliminate the entries in the affected 
>  // columns. 
>  std::map boundary_values; 
>  VectorTools::interpolate_boundary_values (dof_handler, 
>0, 
>ZeroFunction(dim), 
>boundary_values); 
>  MatrixTools::apply_boundary_values (boundary_values, 
>  system_matrix, 
>  solution, 
>  system_rhs, 
>  false); 
>
> In other words, it is expected that the matrix is not symmetric, and 
> that consequently when you run a direct solver, you can't use the 
> symmetric setting. 
>
> Best 
>   Wolfgang 
>

-- 
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] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[Wolfgang Bangerth]

On 07/27/2016 01:43 PM, Alexander wrote:


step-17 claims to solve same problem as step-8. If one produces symmetry
matrix and the other does not or both do not -- this is suspicious.


The reason why the matrix in step-17 is not symmetric is actually 
explained in the code:


// The last argument to the call to
// MatrixTools::apply_boundary_values() below allows for some
// optimizations. It controls whether we should also delete
// entries (i.e., set them to zero) in the matrix columns
// corresponding to boundary nodes, or to keep them (and passing
// true means: yes, do eliminate the columns). If we
// do eliminate columns, then the resulting matrix will be
// symmetric again if it was before; if we don't, then it
// won't. The solution of the resulting system should be the same,
// though. The only reason why we may want to make the system
// symmetric again is that we would like to use the CG method,
// which only works with symmetric matrices. The reason why we may
// not want to make the matrix symmetric is because this
// would require us to write into column entries that actually
// reside on other processes, i.e., it involves communicating
// data. This is always expensive.
//
// Experience tells us that CG also works (and works almost as
// well) if we don't remove the columns associated with boundary
// nodes, which can be explained by the special structure of this
// particular non-symmetry. To avoid the expense of communication,
// we therefore do not eliminate the entries in the affected
// columns.
std::map boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
  0,
  ZeroFunction(dim),
  boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs,
false);

In other words, it is expected that the matrix is not symmetric, and 
that consequently when you run a direct solver, you can't use the 
symmetric setting.


Best
 Wolfgang

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

[deal.II] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[Alexander]

So why is it not 162x162? 

step-17 claims to solve same problem as step-8. If one produces symmetry 
matrix and the other does not or both do not -- this is suspicious. 
Anyway, I recommend you creating an issue on github and put short 
description with the link to this thread. People who created step-17 will 
more likely see it there than go into details of this thread. 

Best,
Alexander

Am Mittwoch, 27. Juli 2016 20:10:09 UTC+2 schrieb shanglong zhang:
>
> Alexander,
>
> I run step-8 and the system matrix is symmetry at least for the first 
> cycle. I attach the output txt of the system matrix in this reply. Thanks.
>
> Best,
> Shanglong
>
> 在 2016年7月27日星期三 UTC-4上午10:53:46，Alexander写道：
>>
>> Shanglong,
>> I'm not an expert in this field, but I assume it should be.
>>
>> Can you please check if system matrix in step-8 (which was used to 
>> implement step-17) is symmetric?
>>
>> Alexander
>>
>> Am Mittwoch, 27. Juli 2016 15:58:42 UTC+2 schrieb shanglong zhang:
>>>
>>> Hi Alexander,
>>>
>>> Thank you for your explanation. I should have checked the symmetricity 
>>> first.
>>> However, since the step-17 is isotropic linear elasticity problem, 
>>> should the system matrix be symmetry?
>>>
>>> Thank you very much for your help and patience.
>>>
>>> Regards,
>>> Shanglong
>>>
>>> 在 2016年7月27日星期三 UTC-4上午4:49:53，Alexander写道：

 Shanglong,
 I can reproduce this behaviour, but this isn't a bug. I output the 
 162x162 system matrix at first cycle and as far as I can see it fails 
 symmetry check. You can see this already from first few non-zero entries:

 (0,0) 1.3
 (1,1) 1.3
 (2,2) 1.3
 (3,3) 1.3
 (4,4) 1.3
 (5,5) 1.3
 (6,0) -0.67

 See, there is no (0, 6) and symmetry is violated.

 The matrix is still positive-definite, that is probably why CG still 
 gets something out of it. 

 Alexander



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

[deal.II] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[shanglong zhang]

Hi Alexander,

I switch to :
dealii 8.4.1
petsc 3.6.4
mumps 5.0.1

The result is still not correct.
The modification I did for the code is replacing the following code in 
step-17.cc:

*SolverControl   solver_control (solution.size(),*
*1e-8*system_rhs.l2_norm());*
*PETScWrappers::SolverCG cg (solver_control,*
*mpi_communicator);*
*PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);*
*cg.solve (system_matrix, solution, system_rhs,*
*  preconditioner);*
with
*  SolverControl solver_control (**solution.size()**, 1e-12);*
*  PETScWrappers::SparseDirectMUMPS solver(solver_control, 
mpi_communicator);*
*  solver.set_symmetric_mode(true);*
*  solver.solve(system_matrix, solution, system_rhs);*

Thank you very much.
Regards,
Shanglong

在 2016年7月26日星期二 UTC-4上午3:40:09，Alexander写道：
>
>
> deal.II 8.4.0
>> petsc 3.3.0.3
>> MUMPS 4.10.0-p3
>>
>
> This PETSc and MUMPS > 4 years old. Update to the latest versions and try 
> again. If problem persists, provide the relevant snippet of the code.
>  
>
>> Sorry, I am not familiar with the petsc library. You meant if I call 
>> *solver.solve(system_matrix,solution,system_rhs)* twice, internally, it 
>> would use the factorization from the first time for the second solving 
>> process? (of course the system matrix is the same.) Do i understand this 
>> correctly? Thanks.
>>
>
> Yes, correct.
>
> @Denis: I agree, the logic is somewhat deficient. I also think 
> documentation and interface itself (e.g. I have no idea why does method 
> set_solver_type 
> 
>  
> exist and exposed to the public) could be improved. 
>

-- 
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] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[Denis Davydov]

Hi Alexander,


> On 26 Jul 2016, at 09:40, Alexander  wrote:
> 
> 
> @Denis: I agree, the logic is somewhat deficient. I also think documentation 
> and interface itself (e.g. I have no idea why does method set_solver_type 
> 
>  exist and exposed to the public) could be improved. 

I created a GitHub issue to track this 
https://github.com/dealii/dealii/issues/2869

As per set_solver_type, it is protected so it is not exposed to users. Should 
be ok, AFAICT.

Regards,
Denis.

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

[deal.II] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[Alexander]

> deal.II 8.4.0
> petsc 3.3.0.3
> MUMPS 4.10.0-p3
>

This PETSc and MUMPS > 4 years old. Update to the latest versions and try 
again. If problem persists, provide the relevant snippet of the code.
 

> Sorry, I am not familiar with the petsc library. You meant if I call 
> *solver.solve(system_matrix,solution,system_rhs)* twice, internally, it 
> would use the factorization from the first time for the second solving 
> process? (of course the system matrix is the same.) Do i understand this 
> correctly? Thanks.
>

Yes, correct.

@Denis: I agree, the logic is somewhat deficient. I also think 
documentation and interface itself (e.g. I have no idea why does method 
set_solver_type 

 
exist and exposed to the public) could be improved. 

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

[deal.II] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[shanglong zhang]

Hi Denis,

Thank you for your information.

Best,
Shanglong

在 2016年7月25日星期一 UTC-4下午12:02:09，Denis Davydov写道：
>
> Hi,
>
> I can't comment on the first part of your question, but as for the second 
> part -- you are right, 
> currently you can't reuse factorization. I you want it, it should be 
> relatively easy to add.
> See Direct solver in TrilinosWrappers which was recently modified to 
> support this: 
>
> https://www.dealii.org/developer/doxygen/deal.II/classTrilinosWrappers_1_1SolverDirect.html#a6c0fecc921bc93f2c005e265cfbdd882
>
> Regards,
> Denis.
>
> On Monday, July 25, 2016 at 5:52:06 PM UTC+2, shanglong zhang wrote:
>>
>> Dear all,
>>
>> I modify the step-17 so that it uses PETScWrappers::SparseDirectMUMPS as 
>> solver. However, when I active the option of set_symmetric_mode, the 
>> solution has discrepancy from CG or MUMPS without the symmetry option. (The 
>> solutions from CG and MUMPS without symmetry option are identical though.) 
>> Any suggestion?
>>
>> And it seems to me that  PETScWrappers::SparseDirectMUMPS doesn't have 
>> the functionality of reusing the factorization of system matrix for 
>> different right-hand side vector. Do I miss anything?
>>
>> Thank you very much.
>>
>> Best,
>> Shanglong
>>
>

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