Sparsity pattern error in the context of diffuse domain with hp-capability and CIP

[Matthias Nicola Henssler]

Dear all,

Background : I am currently trying to implement a diffuse domain formulation to solve the poisson equation whilst penalizing the jumps in gradient in very similar way to what step-85 showcases. The nuance relative to the mentioned tutorial stems the need for AMR close to the boundary.
Having worked through this problem without neglecting the outside of the domain (contrary to what is shown in step-85 with the use of FENothing in the FECollection object for all the elements for which the level set is positive on all vertices), I found myself wanting to implement this feature again.

Issue : Taking good care to assigne FENothing only to the cells that extend beyond an arbitrary threshold to avoid computing the face values on a cell who's neighbor is neglected, I got the following error when computing the flux sparsity pattern. Note that hanging nodes are taken into account. Note that this issue arise only if cells are flagged for ghost penalty and if elements (outside the domain of interest) are attributed the FENothing finite element.



    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
   
    const unsigned int           n_components = fe_collection.n_components();
    Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components);
    Table<2, DoFTools::Coupling> face_coupling(n_components, n_components);
    cell_coupling[0][0] = DoFTools::always;
    face_coupling[0][0] = DoFTools::always;


    const bool                      keep_constrained_dofs = true;
    constraints.clear();
    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
    constraints.close();
    DoFTools::make_flux_sparsity_pattern(dof_handler,
                                          dsp,
                                          constraints,
                                          keep_constrained_dofs,
                                          cell_coupling,
                                          face_coupling,
                                          numbers::invalid_subdomain_id);


-------------------------------------------------------------
The violated condition was:
    matrix_values->column() == column
Additional information:
    You are trying to access the matrix entry with index <14,142>, but
    this entry does not exist in the sparsity pattern of this matrix.
-------------------------------------------------------------

Having tried the following sister make_flux_sparsity_pattern() function call without explicitly passing the cell and face couplings, I get this error:
-------------------------------------------------------------
An error occurred in line <4467> of file </home/mh2294/codes/dealii-9.5.2/include/deal.II/lac/affine_constraints.templates.h> in function
    void dealii::AffineConstraints<number>::add_entries_local_to_global(const std::vector<unsigned int>&, const dealii::AffineConstraints<number>&, const std::vector<unsigned int>&, dealii::SparsityPatternBase&, bool, const dealii::Table<2, bool>&) const [with number = double]
The violated condition was:
    false
Additional information:
    You are trying to use functionality in deal.II that is currently not
    implemented. In many cases, this indicates that there simply didn't
    appear much of a need for it, or that the author of the original code
    did not have the time to implement a particular case. If you hit this
    exception, it is therefore worth the time to look into the code to
    find out whether you may be able to implement the missing
    functionality. If you do, please consider providing a patch to the
    deal.II development sources (see the deal.II website on how to
    contribute).
------------------------------------------------------------

Implying that the both sparsity patterns aren't the same although specifying should result in a full coupling in the case of scalar valued problems.
cell_coupling[0][0] = DoFTools::always;
face_coupling[0][0] = DoFTools::always;

Questions :
1. Is there a way to deal make a correct sparsity pattern accounting for hanging nodes, FENothing elements and coupling of faces that have ghost penalty or is this really a feature that isn't implemented in deal.ii ?

2. Should I assemble the matrix without taking care of the constraints and apply them afterwards ?

3. Why do the sparsity patterns differ from one call to another ?

I have attached the slightly modified step-85 code I am talking about. 

I hope this question isn't redundant,

Thank you for your help,

Matthias

[simon...@gmail.com]

Hi.

1. Is there a way to deal make a correct sparsity pattern accounting for hanging nodes, FENothing elements and coupling of faces that have ghost penalty or is this really a feature that isn't implemented in deal.ii ?

If you do not care so much about performance at this stage, you can just avoid passing the face_has_flux_coupling lambda to make_flux_sparsity_pattern (the version that takes cell and face couplings). You will then get a flux coupling over all faces. This gives a sparsity pattern with more entries than you need, but you will not get an exception.

The alternative is to choose the face_has_flux_coupling lambda more carefully. The problem is that you are trying to assign values to entries in the matrix that does not exist in the sparsity pattern.

2. Should I assemble the matrix without taking care of the constraints and apply them afterwards ?

No, see 1.

3. Why do the sparsity patterns differ from one call to another ?

The implementation of the two functions are different. The simpler function has a bug in it. There is an open issue for this:

https://github.com/dealii/dealii/issues/12037

but since there is a workaround, it has not been prioritized.


Best,
Simon

[Simon Sticko]

Hi.

> 1. Is there a way to deal make a correct sparsity pattern accounting for hanging nodes, FENothing elements and coupling of faces that have ghost penalty or is this really a feature that isn't implemented in deal.ii ?

If you do not care so much about performance at this stage, you can just avoid passing the face_has_flux_coupling lambda to make_flux_sparsity_pattern (the version that takes cell and face couplings). You will then get a flux coupling over all faces. This gives a sparsity pattern with more entries than you need, but you will not get an exception.

The alternative is to choose the face_has_flux_coupling lambda more carefully. The problem is that you are trying to assign values to entries in the matrix that does not exist in the sparsity pattern.

> 2. Should I assemble the matrix without taking care of the constraints and apply them afterwards ?

No, see 1.

> 3. Why do the sparsity patterns differ from one call to another ?

The implementation of the two functions are different. The simpler function has a bug in it. There is an open issue for this:

https://github.com/dealii/dealii/issues/12037

but since there is a workaround, it has not been prioritized.


Best,
Simon

On 01/05/2024 23:39, 'Matthias Nicola Henssler' via deal.II User Group wrote:
> Dear all,
>
> *Background :* I am currently trying to implement a diffuse domain formulation to solve the poisson equation whilst penalizing the jumps in gradient in very similar way to what step-85 showcases. The nuance relative to the mentioned tutorial stems the need for AMR close to the boundary.

> Having worked through this problem without neglecting the outside of the domain (contrary to what is shown in step-85 with the use of FENothing in the FECollection object for all the elements for which the level set is positive on all vertices), I found myself wanting to implement this feature again.
>

> *Issue :* Taking good care to assigne FENothing only to the cells that extend beyond an arbitrary threshold to avoid computing the face values on a cell who's neighbor is neglected, I got the following error when computing the flux sparsity pattern. Note that hanging nodes are taken into account. Note that this issue arise only if cells are flagged for ghost penalty and if elements (outside the domain of interest) are attributed the FENothing finite element.

> *Questions :*

> 1. Is there a way to deal make a correct sparsity pattern accounting for hanging nodes, FENothing elements and coupling of faces that have ghost penalty or is this really a feature that isn't implemented in deal.ii ?
> 2. Should I assemble the matrix without taking care of the constraints and apply them afterwards ?
> 3. Why do the sparsity patterns differ from one call to another ?
>
> I have attached the slightly modified step-85 code I am talking about.
>
> I hope this question isn't redundant,
>
> Thank you for your help,
>
> Matthias
>
>

> --
> The deal.II project is located at http://www.dealii.org/ <http://www.dealii.org/>
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en <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 <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/fb5f6087-b96c-47f1-b1d2-36444f3192b3n%40googlegroups.com <https://groups.google.com/d/msgid/dealii/fb5f6087-b96c-47f1-b1d2-36444f3192b3n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Simon Sticko]

Hi.

> 1. Is there a way to deal make a correct sparsity pattern accounting for hanging nodes, FENothing elements and coupling of faces that have ghost penalty or is this really a feature that isn't implemented in deal.ii ?

If you do not care so much about performance at this stage, you can just avoid passing the face_has_flux_coupling lambda to make_flux_sparsity_pattern (the version that takes cell and face couplings). You will then get a flux coupling over all faces. This gives a sparsity pattern with more entries than you need, but you will not get an exception.

The alternative is to choose the face_has_flux_coupling lambda more carefully. The problem is that you are trying to assign values to entries in the matrix that does not exist in the sparsity pattern.

> 2. Should I assemble the matrix without taking care of the constraints and apply them afterwards ?

No, see 1.

> 3. Why do the sparsity patterns differ from one call to another ?

The implementation of the two functions are different. The simpler function has a bug in it. There is an open issue for this:

https://github.com/dealii/dealii/issues/12037

but since there is a workaround, it has not been prioritized.


Best,
Simon


On 01/05/2024 23:39, 'Matthias Nicola Henssler' via deal.II User Group wrote:
> Dear all,
>

> *Background :* I am currently trying to implement a diffuse domain formulation to solve the poisson equation whilst penalizing the jumps in gradient in very similar way to what step-85 showcases. The nuance relative to the mentioned tutorial stems the need for AMR close to the boundary.

> Having worked through this problem without neglecting the outside of the domain (contrary to what is shown in step-85 with the use of FENothing in the FECollection object for all the elements for which the level set is positive on all vertices), I found myself wanting to implement this feature again.
>

> *Issue :* Taking good care to assigne FENothing only to the cells that extend beyond an arbitrary threshold to avoid computing the face values on a cell who's neighbor is neglected, I got the following error when computing the flux sparsity pattern. Note that hanging nodes are taken into account. Note that this issue arise only if cells are flagged for ghost penalty and if elements (outside the domain of interest) are attributed the FENothing finite element.

> *Questions :*

> 1. Is there a way to deal make a correct sparsity pattern accounting for hanging nodes, FENothing elements and coupling of faces that have ghost penalty or is this really a feature that isn't implemented in deal.ii ?
> 2. Should I assemble the matrix without taking care of the constraints and apply them afterwards ?
> 3. Why do the sparsity patterns differ from one call to another ?
>
> I have attached the slightly modified step-85 code I am talking about.
>
> I hope this question isn't redundant,
>
> Thank you for your help,
>
> Matthias
>
>