function face_component_to_system_index() for Discontinuous Galerkin ?

[Sylvain Mathonnière]

Hello,

I have an equations (named eq1) that I am solving just like in step-12 using the discontinuous Galerkin method but I am solving it using the discrete ordinate method so I actually have N_dir = 40 directions so 40 equations. 

My setup looks like this :

FESystem<dim>        fe_eq1;         

DoFHandler<dim>      dof_handler_eq1;

And the constructor of the class :

template <int dim> DG_FEM<dim>::DG_FEM ()

: fe_eq1 (FE_DGQ<dim>(1),N_dir), dof_handler_eq1 (triangulation)) {}

For the eq1, I am using like in step-12 the MeshWorker::mesh_loop(...) to assemble my matrix :

In the cell_worker function, in order to access the index_dof_i-th component (out of 40) of my shape_function for the direction dir_comp, I have:

...

for (unsigned int dir_comp = 0; dir_comp < N_dir; ++dir_comp)

{

for (unsigned int i = 0; i < n_vertices; ++i) // here I hardcoded n_vertices = 4 for quadrangle

{

unsigned int index_dof_i = fe_eq1.component_to_system_index(dir_comp,i);

for (unsigned int j = 0; j < n_vertices; ++j)

{

unsigned int index_dof_j = fe_eq2.component_to_system_index(dir_comp,j);

  

copy_data_face.cell_matrix_RTE(index_dof_j, index_dof_i) +=

// some calculations involving

fe_face.shape_value_component(index_dof_i, point, dir_comp)

fe_face.shape_value_component(index_dof_j, point, dir_comp)

}

}

}

...

   

I need to loop on the vertices of my cell to do the assembly and I did not find a better way than that. I cannot do like in step-8 looping over all the dof_per_cell, because then

I could not couple the component at one vertice with the same component at another vertice, like I am doing above.

I guess it is fine as long as I only have quadrangle. 

This code works fine it seems and I can apply the same strategy for the boundary_worker. My main issue is that this technique does not work for the face_worker. Here I need to loop on all the vertices twice (because there are two cells sharing the same vertices so two shape function per vertices).

If I write :

for (unsigned int j = 0; j < n_vertices*2; ++j)

{

unsigned int index_dof_jfe_eq2.component_to_system_index(dir_comp,j);

...

this will trigger an error since component_to_system_index(dir_comp,j) will not accept j >= 4

I found in the documentation a fonction face_system_to_component_index() but no function face_component_to_system_index(). What would be the best way to proceed at this point ?

 

Best,

Sylvain M.

PS : As I write this email I realise that I could technically solve 40 times the equation, only changing the direction each time since the components are not coupled to each other (they are howerver coupled with another equation not mentionned here) but that looks cumbersome at best and probably inefficiency.

[Sylvain Mathonnière]

I wanted to add some clarification:

I realised that what I called n_vertices is missleading, it is not the number of vertices I want but the number of "nodes per cell". Numerically, it is however correct if  FE_DGQ<dim>(1).

More generally, it should be something like n_nodes = (fe_eq1.degree+1)^dim, so n_nodes = 4 for  FE_DGQ<dim>(1) but n_nodes = 9 for FE_DGQ<dim>(2).

[Wolfgang Bangerth]

On 9/23/21 5:48 AM, Sylvain Mathonnière wrote:
>
> I found in the documentation a fonction

> *face_system_to_component_index()* but no
> function*face_component_to_system_index()*. What would be the best way

> to proceed at this point ?

It would not be difficult to actually add this function. Would you like
to give it a try? If you see how the non-face functions are implemented,
you'd probably see right away how to add such a function for faces as well.

That said, the way you approach this is a bit unusual. A better way
would be to consider that you have a vector-valued element and, as we
always do, loop over all DoFs i and j that belong to *all* components.
You can then either do something like

for (i=0...dofs_per_cell)
for (j=0...dofs_per_cell)
if (fe.system_to_component_index(i).first ==
fe.system_to_component_index(j).first)
...add matrix term...

or, if you would rather do this, write things as

for (i=0...dofs_per_cell)
for (j=0...dofs_per_cell)
for (c=0...n_components)
copy_data_face.cell_matrix_RTE(i,j) +=
... fe_face[c].shape_value(i,q) ...
... fe_face[c].shape_value(j,q) ...

The first of these ways would correspond to how step-8 assembles its
matrix, whereas the second is used in all other tutorial programs
solving vector-valued problems (e.g., step-22). You might also want to
look at the vector-valued documentation module for discussions around
these sorts of issues.

Best
W.

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

[Sylvain Mathonnière]

About writing the function :

I must confess I did not try to write the function myself but I did check to see how it was implemented and it seems like it is just retrieving the value from a table. I currently only have the user version of deal.II (v 9.2.0),

If I wanted to write this function, the proper way to do it would be to download the developer version 9.3.0, write the function down then recompile the librairies ?

I haven't thought much about contributing before to deal.II because I am fairly new to it (I love it) and not that skilled at c++.

About your first method :

Yes, in the end this is the way I implemented it, but it takes really long time to compute since dof_per_cell = 160 for me (with Q1 mapping and fe_degree =1), so I have 160*160 iterations but the "if condition" is only true 1/40th of the time, so I end up with a lot of unnecessary looping.

My equation system is b.\nabla I + I = RHS where b is a direction vector and I is a scalar function. For the discrete ordinate method, I must solve this equation for 40 values of b.

About your second method :

The problem in doing so, is that for me the 40 components are already included in the dof_per_cell (4 nodes * 40 components = 160 dof_per_cell), so the ith and jth indices of the loops already correspond to a specific component. Hence why I used the 

system_to_component_index() function.

( I have  FESystem<dim>  =  fe_eq1 (FE_DGQ<dim>(1), 40) )

What worked the best :

I mentioned in my previous post that since my equations are not coupled with each other, I could technically loop over the assembly and solving functions 40 times only changing a constant every time (I change the direction vector  for the discrete ordinate method).

I though it was a really unefficient way but I tried it nonetheless and it turns out that it is actually much faster (I gain a speed factor of 16).

In the end, assembling and solving 40 times a system with a sparse matrix of size (n_cells*4)^2 is faster than assembling and solving 1 time a system with a sparse matrix of size (n_cells*4*40)^2.

It sounds correct but does that make sense ?

Caveat :

I use MeshWorker::mesh_loop with queue_length = 1 so I believe I only use one thread, maybe it is a different result with more thread.

Otherwise if I use mesh_loop like in step 12, I get an error

"pure virtual method called
terminate called without an active exception
Abandon (core dumped)"

I believe this error occurs because I create a second iterator in the cell worker to access values bound to another dof_handler like this:

typename DoFHandler<dim>::active_cell_iterator cell_2 (&cell->get_triangulation(), cell->level(), cell->index(), &dof_handler_2);

Best,

Sylvain