Fix issue with constraints in the VANS solver and update boycott (cha…

…os-polymtl#895)

Description of the problem
There was a major issue with the void fraction constraints in the VANS and CFD-DEM solver which led to an assertion being thrown in debug mode. The constraints were not reinited with the right size

Description of the solution
Fix this issue by reiniting the constraints correctly
Also fixed an issue with an example while I was there

How Has This Been Tested?
All tests should now pass in debug

Comments
Some part of the VANS solver and the unresolved CFD-DEM solver are still a bit of a mess. I will revisit this slowly but surely within the coming months.

CHANGELOG.md

@@ -3,6 +3,12 @@ All notable changes to the Lethe project will be documented in this file.
 
 The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
+## [Master] - 2023-09-30
+
+### Fixed
+
+- MINOR Affine constraints used to bound the void fraction and used as boundary conditions within the heat transfer and the block Navier-stokes solver were corrected [#885](https://github.com/lethe-cfd/lethe/pull/895): There was an error in the application of the affine constraints used to clip the void fraction in the lethe-fluid-vans and lethe-fluid-particles solvers. This led to assertions being thrown in debug. This has been corrected by reiniting the constraints with the appropriate size. For the heat transfer and the block Navier-stokes, the issue was that the constraints were never reinited with the correct size, so they contained ghosted elements. This was caught by a new assert introduced in 2023-09 within deal.II master.
+
 ## [Master] - 2023-09-20
 
 ### Deprecated
@@ -36,4 +42,4 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
 ### Fixed
 
-- MAJOR/MINOR/PATCH Description (#PR).
+- MAJOR/MINOR/PATCH Description (#PR).

examples/unresolved-cfd-dem/boycott-effect/boycott-effect.prm

@@ -106,7 +106,6 @@ subsection cfd-dem
   set shear force                   = true
   set pressure force                = true
   set drag model                    = difelice
-  set post processing               = true
   set coupling frequency            = 250
   set grad-div length scale         = 0.005
   set vans model                    = modelA

source/fem-dem/gls_vans.cc

@@ -325,8 +325,14 @@ GLSVANSSolver<dim>::update_solution_and_constraints()
                                                 system_rhs_void_fraction);
 
   void_fraction_constraints.clear();
+
+  // reinitialize affine constraints
+  void_fraction_constraints.reinit(locally_relevant_dofs_voidfraction);
+
+  // Remake hanging node constraints
+  DoFTools::make_hanging_node_constraints(void_fraction_dof_handler,
+                                          void_fraction_constraints);
   active_set.clear();
-  std::vector<bool> dof_touched(void_fraction_dof_handler.n_dofs(), false);
 
   for (const auto &cell : void_fraction_dof_handler.active_cell_iterators())
     {

source/solvers/gd_navier_stokes.cc

@@ -559,6 +559,8 @@ GDNavierStokesSolver<dim>::setup_dofs_fd()
   auto &nonzero_constraints = this->nonzero_constraints;
   {
     nonzero_constraints.clear();
+    nonzero_constraints.reinit(locally_relevant_dofs_acquisition);
+
 
     DoFTools::make_hanging_node_constraints(this->dof_handler,
                                             nonzero_constraints);
@@ -627,6 +629,8 @@ GDNavierStokesSolver<dim>::setup_dofs_fd()
 
   {
     this->zero_constraints.clear();
+    this->zero_constraints.reinit(locally_relevant_dofs_acquisition);
+
     DoFTools::make_hanging_node_constraints(this->dof_handler,
                                             this->zero_constraints);
 

source/solvers/heat_transfer.cc

@@ -1018,6 +1018,7 @@ HeatTransfer<dim>::setup_dofs()
 
   {
     nonzero_constraints.clear();
+    nonzero_constraints.reinit(this->locally_relevant_dofs);
     DoFTools::make_hanging_node_constraints(this->dof_handler,
                                             nonzero_constraints);
 
@@ -1043,6 +1044,7 @@ HeatTransfer<dim>::setup_dofs()
   // Boundary conditions for Newton correction
   {
     zero_constraints.clear();
+    zero_constraints.reinit(this->locally_relevant_dofs);
     DoFTools::make_hanging_node_constraints(this->dof_handler,
                                             zero_constraints);