[deal.II] Re: triangulation.add_periodicity takes an enormous amount of time

[Daniel Arndt]

Andreas,

I am trying to read from a mesh file and impose periodic boundary 
> condition. This works out very well for small geometries in three spatial 
> dimensions. However, if the read meshes tend to have more elements, the 
> amount of time needed to create the mesh rises drastically. The time is 
> spent within the routine triangulation.add_periodicity. Thus, I have 
> created a minimal code example where cubes of side lengths 5, 20, and 50 
> are created. For each cube, the periodic face pairs are identified and 
> collected. Afterwards, triangulation.add_periodicity is executed and the 
> needed time is measured using the Timer class. The output of the attached 
> code example looks as follows:
>
> Starting Constructor with side length = 5
> Periodic faces collected: Size of periodicity_vector = 75
> Start triangulation.add_periodicity(periodicity_vector); after 0.018097 
> seconds
>  Done triangulation.add_periodicity(periodicity_vector); after 3.0626 
> seconds
>
> Starting Constructor with side length = 20
> Periodic faces collected: Size of periodicity_vector = 1200
> Start triangulation.add_periodicity(periodicity_vector); after 2.10931 
> seconds
>  Done triangulation.add_periodicity(periodicity_vector); after 15749 
> seconds
>
> Starting Constructor with side length = 50
> Periodic faces collected: Size of periodicity_vector = 7500
> Start triangulation.add_periodicity(periodicity_vector); after 298.734 
> seconds
>

Were you running in Debug or Release mode? I only tried Release mode and 
the timings weren't that terrible. I observed a runtime like three times 
larger in that function than for creating the mesh for up to 100 elements 
per side.
Still, there are some oppurtunities for improvement, see 
https://github.com/dealii/dealii/pull/8665. With these changes, creating 
the mesh and adding periodicity takes roughly about the same time.

Apart from that, it is much better to have a coarse mesh with few cells and 
refine a few times than starting with a rather fine coarse mesh (if you can 
do that).

Best,
Daniel

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2982487e-6e33-4f49-b4e2-a4cc93229d61%40googlegroups.com.

Re: [deal.II] Implement 'w.n = u.n' as a boundary condition, where 'u' is a vector valued FE function/solution from a phase/PDE and 'n' is the unit outward normal in a distributed setting

[Bhanu Teja]

Thanks Daniel.
I have looked at 'VectorTools::compute_nonzero_normal_flux_constraints'. 
But what I have here to compute the 'inhomogeneity' is not a closed form 
function, but a discrete FE solution computed from another PDE(in my case 
Stokes). The 'u' in 'w.n = u.n' is a 'LA::MPI::Vector' FE solution vector. 
Can you please suggest?

Bhanu.

On Wednesday, August 28, 2019 at 9:43:42 PM UTC+5:30, Daniel Arndt wrote:
>
> Bhanu,
>
> Please suggest me in implementing 'w.n = u.n' as a boundary condition 
>> where 'u' is Stokes solution with 'FESystem 
>> fe_stokes(FE_Q(degree+1), dim, FE_Q(degree), 1)' and 'w' for mesh 
>> velocity with 'FESystem fe_move(FE_Q(degree), dim)'. As I have 
>> taken 'degree=1', 'fe_move' has unit support only at geometric vertices. I 
>> want to know if there is a direct way to get outward unit normal at a 
>> boundary vertex which is average of normals of faces that share this 
>> vertex(The faces are flat as I have not assigned any MapingQ). If this 
>> is not possible/straight_forward, the normal at a quadrature point on the 
>> face from which I visit this vertex the first time is good enough. 
>> [...]
>>
>
> Have a look at VectorTools::compute_nonzero_normal_flux_constraints (
> https://www.dealii.org/current/doxygen/deal.II/group__constraints.html#gacb97b6eb5113d649d4daee0867dc82ac
> ).
>
> Best,
> Daniel
>
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/77b9561f-1d89-44d0-9b66-3046b998f30b%40googlegroups.com.

Re: [deal.II] Implement 'w.n = u.n' as a boundary condition, where 'u' is a vector valued FE function/solution from a phase/PDE and 'n' is the unit outward normal in a distributed setting

[Daniel Arndt]

Bhanu,

Please suggest me in implementing 'w.n = u.n' as a boundary condition where
> 'u' is Stokes solution with 'FESystem fe_stokes(FE_Q(degree+1),
> dim, FE_Q(degree), 1)' and 'w' for mesh velocity with 'FESystem
> fe_move(FE_Q(degree), dim)'. As I have taken 'degree=1', 'fe_move' has
> unit support only at geometric vertices. I want to know if there is a
> direct way to get outward unit normal at a boundary vertex which is average
> of normals of faces that share this vertex(The faces are flat as I have not
> assigned any MapingQ). If this is not possible/straight_forward, the
> normal at a quadrature point on the face from which I visit this vertex the
> first time is good enough.
> [...]
>

Have a look at VectorTools::compute_nonzero_normal_flux_constraints (
https://www.dealii.org/current/doxygen/deal.II/group__constraints.html#gacb97b6eb5113d649d4daee0867dc82ac
).

Best,
Daniel

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAOYDWbK212WjrKR3dL%3Dq-LSP_FjTYt%3DQ9RmJkfPVRqGfykSRFA%40mail.gmail.com.

Re: [deal.II] On using block-diagonal preconditioner for the mixed Poisson problem from step-20

[David Wells]

Hi Charlie,

The issue is that Trilinos really needs an actual Matrix object (not just a
linear operator) to compute the AMG preconditioner: i.e., it can only
construct the preconditioner if it has access to the actual matrix entries,
which are not available unless you explicitly calculate B (DA)^-1 B^T. To
get this to work you will have to multiply those three matrices together
and then pass that result (stored as a Trilinos matrix) to the initialize
function.

It should be possible to do this with TrilinosWrappers::SparseMatrix::mmult
and TrilinosWrappers::SparseMatrix::Tmmult.

Thanks,
David

On Tue, Aug 27, 2019 at 2:39 PM Charlie Z  wrote:

> Dear all,
>
> I'm trying to use a block-diagonal preconditioner for the mixed Poisson
> problem from step-20, but haven't made it yet.
>
> The diagonal preconditioner of interest is from Powell & Silvester,2003
> . Specifically,
> for the matrix system after discretization with the RT0-DG0 element pair
>   [A  BT
>B  0],
> they suggested a block diagonal preconditioner
>   P =
>   [DA  0
>0 PS],
> where DA is diag(A), and PS is AMG applied to the approximated Schur
> complement SD = B(DA^-1)BT.
>
> As the above idea is quite similar to that in step-55, so I borrowed quite
> a few code blocks from there, e.g. InverseMatrix and
> BlockDiagonalPreconditioner. The resulting solve() with a SolverMinRes solver
> is as follows. (The rest of the code is mainly from step-20)
>
>   template 
>  void MixedLaplaceProblem::solve()
>  {
>TrilinosWrappers::PreconditionJacobi diagA;
>diagA.initialize(system_matrix.block(0, 0));//DA
>
> ApproximateSchurComplement approximate_schur(system_matrix, solution);
> //SD
>
> InverseMatrix approximate_inverse(
> approximate_schur);//PS
>
> const BlockDiagonalPreconditioner ,
>  InverseMatrix >
>  preconditioner(diagA, approximate_inverse);
>
> SolverControl solver_control(system_matrix.m(),
> 1e-6 * system_rhs.l2_norm());
>
> SolverMinRes solver(solver_control
> );
>solver.solve(system_matrix, solution, system_rhs,
>   preconditioner
>  );
>
> std::cout<<"last step = "<  }
>
> However, the code can't compile, and gives the following error message:
>
>   /home/dealii/test.cc:336:5: error: no matching function for call to ‘
> dealii::TrilinosWrappers::PreconditionAMG::initialize(const Step20::
> ApproximateSchurComplement&, dealii::TrilinosWrappers::PreconditionAMG::
> AdditionalData&)’
>   preconditioner.initialize(*matrix, data);
>   ^
>  In file included from /home/dealii/dealii-v9.0.0/include/deal.II/lac/
> generic_linear_algebra.h:141:0,
>   from /home/dealii/test.cc:31:
>  /home/dealii/dealii-v9.0.0/include/deal.II/lac/trilinos_precondition.h:
> 1516:10: note: candidate: void dealii::TrilinosWrappers::PreconditionAMG::
> initialize(const dealii::TrilinosWrappers::SparseMatrix&, const dealii::
> TrilinosWrappers::PreconditionAMG::AdditionalData&)
>   void initialize (const SparseMatrix   ,
>
> The error is related to the use of PreconditionAMG in my implementation of
> InverseMatrix::vmult:
>
>   template 
>  void InverseMatrix::vmult (LinearAlgebraTrilinos::MPI::Vector
>   ,
> const LinearAlgebraTrilinos::MPI::
> Vector ) const
>  {
>SolverControl solver_control(src.size(),
> 1e-6 * src.l2_norm());
>SolverCGcg (solver_control);
>dst = 0;
>
> LinearAlgebraTrilinos::MPI::PreconditionAMG preconditioner;
>LinearAlgebraTrilinos::MPI::PreconditionAMG::AdditionalData data;
>preconditioner.initialize(*matrix, data);
>
> // PreconditionIdentity preconditioner;//A Identity preconditioner
> works fine here
>
> cg.solve (*matrix, dst, src, preconditioner);
>
> std::cout<<"inner last step = "< ;
>  }
>
> It seems that a TrilinosWrappers::PreconditionAMG can't be initialized
> with a user-defined ApproximateSchurComplement object. Instead, a matrix
> form of the approximated Schur complement SD=B(DA^-1)BT is expected.
>
> Could you please give a hint on how to form this matrix?
>
> Attached is a minimal (but not working) code.
>
> Best regards,
> Charlie
>
> --
> 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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/f45e2027-3a16-4424-8cfd-956381725502%40googlegroups.com
> 
> .
>

-- 
The deal.II project is located at 

Re: [deal.II] On using block diagonal preconditioner for the mixed Poisson problem

[Wolfgang Bangerth]

On 8/28/19 6:23 AM, Charlie Z wrote:
> 
> The code works fine now with the help of 
> TrilinosWrappers::SparseMatrix::mmult(), though, I understand, mmult is 
> expensive.

Glad to hear! Yes, mmult is expensive, but then you're only doing it once. 
Unless you have concrete evidence to the contrary, performance problems are 
not worth chasing.

Best
  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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/15388c89-f99c-da0a-5a18-b3e0bebf2948%40colostate.edu.

Re: [deal.II] On using block diagonal preconditioner for the mixed Poisson problem

[Charlie Z]

Hi Wolfgang,

Thanks very much for your clarification.

The code works fine now with the help of 
TrilinosWrappers::SparseMatrix::mmult(), though, I understand, mmult is 
expensive. 

Regards,
Charlie

On Wednesday, August 28, 2019 at 12:43:04 AM UTC+2, Wolfgang Bangerth wrote:
>
> On 8/27/19 2:23 PM, Charlie Z wrote: 
> > 
> > I find a conflict here is that:: 
> > 1) TrilinosWrappers::PreconditionAMG can only be initialized with a 
> > matrix, but not a ApproximateSchurComplement; 
> > 2) step-20 uses ApproximateSchurComplement because we actually don't 
> > want to/can't assemble the Schur matrix or its inverse. 
>
> Yes, you correctly identified the issue. step-20 gets away with using 
> ApproximateSchurComplement as a preconditioner because, from the 
> perspective of GMRES, a preconditioner only has be able to provide a 
> matrix-vector product (i.e., apply itself to a vector). 
>
> But if you want to use an object to initialize a algebraic multigrid, 
> the AMG implementation needs to be able to do much more than just apply 
> the object to a vector -- it needs to build a multilevel hierarchy, and 
> for that it needs access to the elements that make up that matrix. The 
> problem is of course that ApproximateSchurComplement does not offer 
> access to the elements of the matrix it represents (because it never 
> explicitly computes them). 
>
>
> > Could you please give me a hint on how to solve this conflict in 
> > deal.II? Is it actually possible to assemble this B(DA^-1)BT directly? 
>
> Yes, of course. You can just multiply it out as the matrix-matrix 
> product of three matrices. The result is a square matrix of the right 
> size, with a rather large number of nonzero entries per row. 
>
> Best 
>   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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4fadbe93-ff6c-4653-8b22-180bf71f095f%40googlegroups.com.

[deal.II] Implement 'w.n = u.n' as a boundary condition, where 'u' is a vector valued FE function/solution from a phase/PDE and 'n' is the unit outward normal in a distributed setting

[Bhanu Teja]

Dear deal.ii,
Please suggest me in implementing 'w.n = u.n' as a boundary condition where 
'u' is Stokes solution with 'FESystem fe_stokes(FE_Q(degree+1), 
dim, FE_Q(degree), 1)' and 'w' for mesh velocity with 'FESystem 
fe_move(FE_Q(degree), dim)'. As I have taken 'degree=1', 'fe_move' has 
unit support only at geometric vertices. I want to know if there is a 
direct way to get outward unit normal at a boundary vertex which is average 
of normals of faces that share this vertex(The faces are flat as I have not 
assigned any MapingQ). If this is not possible/straight_forward, the 
normal at a quadrature point on the face from which I visit this vertex the 
first time is good enough. 

I attempted implementing in the following way
MappingQGeneric mapping(fe.degree); 
Quadrature face_quadrature(femove.
get_unit_face_support_points());
FEFaceValues fe_face_values_stokes(mapping, fe, 
face_quadrature, update_values | update_normal_vectors);
const unsigned int n_face_dofs = face_quadrature.size();   
 
std::vector> stokes_values(n_face_dofs,  Vector
(dim + 1));
std::vector local_dof_indices(femove.
dofs_per_face);

Tensor<1, dim> normal_vector1;
Tensor<1, dim> col_indices;
Tensor<1, dim> stokes_value_atvertex;

auto cell_move = move_dof_handler.begin_active();
auto endcell_move = move_dof_handler.end();
auto cell_stokes = dof_handler.begin_active();

std::vector vertex_touched(triangulation.n_vertices(), 
false);
for (; cell_move != endcell_move; ++cell_stokes, ++cell_move)
{
if (!cell_move->is_artificial())
{
for (unsigned int f = 0; f < GeometryInfo::
faces_per_cell; ++f) 
{
if (cell_move->face(f)->boundary_id() == 1006)
{
fe_face_values_stokes.reinit(cell_stokes, f);
normal_vector1 = fe_face_values_stokes.
normal_vector(1);
if(fabs(normal_vector1[0]) < 1e-10)
normal_vector1[0] = 1e-10;
fe_face_values_stokes.get_function_values(
lr_solution, stokes_values);
cell_move->face(f)->get_dof_indices(
local_dof_indices);

for (unsigned int v=0; v < GeometryInfo::
vertices_per_face; ++v)
{
if(vertex_touched[cell_move->face(f)->
vertex_index(v)] == false)
{
vertex_touched[cell_move->face(f)->
vertex_index(v)] = true;
for(int j=0; j < dim; ++j)
{
col_indices[j] = cell_move->face(f
)->vertex_dof_index(v, j);
if(j==0)
moveconstraints.add_line(
col_indices[0]);
const unsigned int component = 
femove.system_to_component_index(v*dim+j).first;
stokes_value_atvertex[j] = 
stokes_values[v*dim+j](component);
if(j!=0)
moveconstraints.add_entry(
col_indices[0], col_indices[j], normal_vector1[j]/normal_vector1[0]);
}
moveconstraints.set_inhomogeneity(
col_indices[0], -stokes_value_atvertex*normal_vector1/normal_vector1[0]); 
}
}
}
}
}
}

but I get the following error.
void dealii::AffineConstraints::add_entry(dealii::AffineConstraints
::size_type, dealii::AffineConstraints::size_type, number) [
with number = double; dealii::AffineConstraints::size_type = 
unsigned int]
The violated condition was: 
std::abs(p.second - value) < 1.e-14
Additional information: 
The entry for the indices 1143 and 1144 already exists, but the values 
0.387236 (old) and 0 (new) differ by -0.387236.

Please correct any gaps in my understanding and suggest if there is a 
general way which works with any order finite element.

Thanks,
Bhanu.


-- 
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.
To view this discussion on the web visit