Refactor void fraction out of solver (#1392)

Description
The facilities required to calculate the void fraction were all contained within the unresolved CFD-DEM solver and the VANS solver. This is very problematic, because it makes it difficult to develop new void fraction schemes, but also prevents the development of new VANS schemes (e.g. matrix free). Consequently, I Have used the new interface developed by @AmishgaAlphonius to refactor the void fraction calculation completely out of the VANS and unresolved CFD-DEM solver. The result is very lightweight.

I also took the opportunity to deprecate the bounding of the void fraction. This feature was never used and was actually dangerous for convergence. Consequently, it has now been fully deprecated and the documentation has been adapted.

Testing
None of the tests results were changed, this is great since nothing should have changed.

Documentation
I have removed the documentation related to the bounding of the void fraction and all of its mention within the examples.

Co-authored-by: Amishga Alphonius <[email protected]>

CHANGELOG.md

@@ -4,6 +4,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] - 2024-12-03
+
+### Removed
+
+- MAJOR The ability to bound the void fraction from below and above (using l2 lower bound and l2 upper bound parameters) has been removed. This bounding of the void fraction was highly problematic, since it could create discontinuities in the time derivative of the void fraction and contaminate the solution. In reality, this feature was never used. Importantly, this PR refactors the calculation of the void fraction outside of the VANS and CFD-DEM solvers of Lethe so that it exists only independently in a seperate subequation solver. This change is in preperation for the matrix free implementation of the VANS equations. [#1392](https://github.com/chaos-polymtl/lethe/pull/1392)
+
 ## Release of Lethe v1.0 - 2024-11-30
 
 The lethe v1.0 release marks the transition of lethe to a numbered released format. Following numbered releases will be generated more often and we aim to release a numbered version 4 to 6 times per year or when sufficient changes have been realized.

applications_tests/lethe-fluid-particles/adaptive_sparse_contacts.prm

@@ -164,7 +164,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.0001
-  set l2 lower bound      = 0
 end
 
 # --------------------------------------------------

applications_tests/lethe-fluid-particles/conserve_phase_volumes.prm

@@ -74,9 +74,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.00000125
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
-  set bound void fraction = false
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/dynamic_contact_search.prm

@@ -73,9 +73,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.00001
-  set bound void fraction = false
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/liquid_fluidized_bed.prm

@@ -75,9 +75,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.00000125
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
-  set bound void fraction = false
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/particle_sedimentation.prm

@@ -73,9 +73,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.00001
-  set bound void fraction = false
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/periodic_particles_qcm.prm

@@ -146,7 +146,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.0005
-  set bound void fraction = false
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/periodic_particles_qcm_opposite_flow.prm

@@ -146,7 +146,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.0005
-  set bound void fraction = false
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/restart-gas-solid-fluidized-bed.prm

@@ -78,9 +78,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.000005
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
-  set bound void fraction = false
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/restart_particle_sedimentation.prm

@@ -73,9 +73,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.00001
-  set bound void fraction = false
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
 end
 
 #---------------------------------------------------

applications_tests/lethe-fluid-particles/spouted_bed_load_balancing.prm

@@ -59,7 +59,6 @@ subsection void fraction
   set qcm sphere equal cell volume = true
   set read dem                     = true
   set dem file name                = dem
-  set bound void fraction          = false
 end
 
 #---------------------------------------------------

doc/source/examples/unresolved-cfd-dem/boycott-effect/boycott-effect.rst

@@ -274,7 +274,7 @@ Void Fraction
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Since we are calculating the void fraction using the particle insertion of the DEM simulation, we set the ``mode`` to ``dem``. For this, we need to read the dem files which we already wrote using check-pointing. We, therefore, set the ``read dem`` to ``true`` and specify the prefix of the dem files to be dem.
 We choose to use the quadrature centered method (`QCM <../../../theory/unresolved_cfd-dem/unresolved_cfd-dem.html>`_) to calculate the void fraction. For this, we specify the ``mode`` to be ``qcm``. We want the radius of our volume averaging sphere to be equal to the length of the element where the void fraction is being calculated. We don't want the volume of the sphere to be equal to the volume of the element.
-For this, we set the ``qcm sphere equal cell volume`` equals to ``false``. Since we want to keep the mass conservative properties of the :math:`L^2` projection, we do not bound the void fraction and as such we set ``bound void fraction`` to ``false``. Unlike the other schemes, we do not smooth the void fraction as we usually do using the PCM and SPM void fraction schemes since QCM is continuous in time and space.
+For this, we set the ``qcm sphere equal cell volume`` equals to ``false``. Unlike the other schemes, we do not smooth the void fraction as we usually do using the PCM and SPM void fraction schemes since QCM is continuous in time and space.
 
 .. code-block:: text
 
@@ -283,7 +283,6 @@ For this, we set the ``qcm sphere equal cell volume`` equals to ``false``. Since
       set qcm sphere equal cell volume = false
       set read dem                     = true
       set dem file name                = dem
-      set bound void fraction          = false
     end
 
 CFD-DEM

doc/source/examples/unresolved-cfd-dem/cylindrical-packed-bed/cylindrical-packed-bed.rst

@@ -296,7 +296,7 @@ The additional sections that define the VANS solver are the void fraction subsec
 Void Fraction
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Since we are calculating the void fraction using the packed bed of the DEM simulation, we set the mode to ``dem``. For this, we need to read the dem files which we already wrote using check-pointing. We therefore set the read dem to ``true`` and specify the prefix of the ``dem`` files to be read. In order to ensure that our void fraction projection is bounded, we choose an upper bound limit of 1. We decide not to lower bound the void fraction and thus attributed a value of 0 to the L2 lower bound parameter. We now choose a smoothing factor for the void fraction as to reduce discontinuity which can lead to oscillations in the velocity. The factor we choose is around the square of twice the particle's diameter. 
+Since we are calculating the void fraction using the packed bed of the DEM simulation, we set the mode to ``dem``. For this, we need to read the dem files which we already wrote using check-pointing. We therefore set the read dem to ``true`` and specify the prefix of the ``dem`` files to be read. We now choose a smoothing factor for the void fraction as to reduce discontinuity which can lead to oscillations in the velocity. The factor we choose is around the square of twice the particle's diameter. 
 
 .. code-block:: text
 

doc/source/examples/unresolved-cfd-dem/gas-solid-fluidized-bed/gas-solid-fluidized-bed.rst

@@ -308,7 +308,7 @@ The additional sections for the CFD-DEM simulations are the void fraction subsec
 Void Fraction
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Since we are calculating the void fraction using the packed bed of the DEM simulation, we set the mode to "dem". For this, we need to read the dem files which we already wrote using check-pointing. We, therefore, set the read dem to "true" and specify the prefix of the dem files to be dem. In order to ensure that our void fraction projection is bounded, we choose an upper bound limit of 1. We decide not to lower bound the void fraction and thus attributed a value of 0 to the L2 lower bound parameter. We now choose a smoothing factor for the void fraction to reduce discontinuity which can lead to oscillations in the velocity. The factor we choose is around the square of twice the particle's diameter. 
+Since we are calculating the void fraction using the packed bed of the DEM simulation, we set the mode to "dem". For this, we need to read the dem files which we already wrote using check-pointing. We, therefore, set the read dem to "true" and specify the prefix of the dem files to be dem. We now choose a smoothing factor for the void fraction to reduce discontinuity which can lead to oscillations in the velocity. The factor we choose is around the square of twice the particle's diameter. 
 
 .. code-block:: text
 
@@ -317,9 +317,6 @@ Since we are calculating the void fraction using the packed bed of the DEM simul
         set read dem            = true
         set dem file name       = dem
         set l2 smoothing factor = 0.000005
-        set l2 lower bound      = 0
-        set l2 upper bound      = 1
-        set bound void fraction = false
     end
 
 CFD-DEM

doc/source/examples/unresolved-cfd-dem/gas-solid-spouted-bed/gas-solid-spouted-bed.rst

@@ -307,7 +307,7 @@ Void Fraction
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Since we are calculating the void fraction using the packed bed of the DEM simulation, we set the ``mode`` to ``dem``. For this, we need to read the dem files which we already wrote using check-pointing. We, therefore, set the ``read dem`` to ``true`` and specify the prefix of the dem files to be dem. We choose to use the quadrature centered method (QCM) to calculate the void fraction. This method does not require smoothing the void fraction as it is space and time continuous. For this simulation, we use a reference sphere having the same volume as the mesh elements as the averaging volume to calculate the void fraction.
-For this, we specify the ``mode`` to be ``qcm``. We want the volume of the volume averaging sphere to be equal to the volume of the element. For this, we set the ``qcm sphere equal cell volume`` equals to ``true``. Since we want to keep the mass conservative properties of the  :math:`L^2` projection, we do not bound the void fraction and as such we set ``bound void fraction`` to ``false``.
+For this, we specify the ``mode`` to be ``qcm``. We want the volume of the volume averaging sphere to be equal to the volume of the element. For this, we set the ``qcm sphere equal cell volume`` equals to ``true``.
 
 .. code-block:: text
 
@@ -316,7 +316,6 @@ For this, we specify the ``mode`` to be ``qcm``. We want the volume of the volum
       set qcm sphere equal cell volume = true
       set read dem                     = true
       set dem file name                = dem
-      set bound void fraction          = false
     end
 
 CFD-DEM

doc/source/examples/unresolved-cfd-dem/liquid-solid-fluidized-bed/liquid-solid-fluidized-bed.rst

@@ -332,8 +332,6 @@ We choose the `particle centroid method (PCM) <../../../parameters/unresolved-cf
       set l2 smoothing factor = 2.8387584e-5
     end
 
-.. note::
-    Note that void fraction is not bound in this case. The size of the particles used in this example forces us to use a very coarse mesh. Bounding void fraction would lead to instability in the present case.
 
 CFD-DEM
 ~~~~~~~~~~

doc/source/parameters/unresolved-cfd-dem/void-fraction.rst

@@ -10,9 +10,6 @@ In this subsection, all parameters required for the calculation of the void frac
     set read dem                   = true
     set dem file name              = dem
     set l2 smoothing factor        = 0
-    set l2 lower bound             = 0
-    set l2 upper bound             = 1
-    set bound void fraction        = false
     set particle refinement factor = 0
   end
 
@@ -38,8 +35,6 @@ If the ``mode`` chosen is ``pcm``, then the void fraction is calculated using th
 * The ``read dem`` allows us to read an already existing dem simulation result which can be obtained from checkpointing the Lethe-DEM simulation. This is important as the `lethe-fluid-vans` solver requires reading an initial dem triangulation and particle information to simulate flows in the presence of particles. 
 * The ``dem_file_name`` parameter specifies the prefix of the dem files that must be read.
 * The ``l2 smoothing factor`` is a smoothing length used for smoothing the L2 projection of the void fraction to avoid sharp discontinuities which can lead to instabilities in the simulation.
-* The ``l2 lower bound`` and ``l2 upper bound`` are the minimum and maximum values around which the void fraction is bounded. This is important especially for upper bounds as the void fraction can sometimes slightly exceed a value of 1 when projected.
-* The ``bound void fraction`` parameter determines whether or not to bound the void fraction between the lower and upper bounds specified in the previous two parameters. As the void fraction is calculated and then projected using L2 projection, it can sometimes exceeds the maximum value of 1. In order to prevent this, we use an active set method to bound the void fraction.
 * The ``qcm sphere diameter`` allows us to fix the diameter of all reference spheres in the simulation to a given value. If this option is used (a value other than 0 is specified), it overrides the default calculation of the size of the sphere and sets its diameter to the value specified.
 * The ``qcm sphere equal cell volume`` determines whether or not we want to use a reference sphere with the same volume as the element in which it is located. If it is disabled, then each sphere will have a radius equal to the size of the element in which it is located. This parameter is important only when the ``qcm sphere diameter`` is not used or is set to 0.
 * The ``particle refinement factor`` is only required for the ``spm``. It allows to determine the number of pseudo-particles that we want to divide our particle into. By default, it is set to 0 refinements, and results in no refinement of the original meshed particle (division into 7 particles in 3D). Every additional refinement results in a :math:`2^{dim}` times more particles. The figure below shows how the number of pseudo-particles change with every refinement. Every cell in the particle mesh represents a pseudo-particle in the satellite point method.

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

@@ -94,7 +94,6 @@ subsection void fraction
   set qcm sphere equal cell volume = false
   set read dem                     = true
   set dem file name                = dem
-  set bound void fraction          = false
 end
 
 #---------------------------------------------------

examples/unresolved-cfd-dem/cylindrical-packed-bed/flow-in-bed.prm

@@ -55,8 +55,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.00001
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
 end
 
 #---------------------------------------------------

examples/unresolved-cfd-dem/gas-solid-fluidized-bed/gas-solid-fluidized-bed.prm

@@ -58,9 +58,6 @@ subsection void fraction
   set read dem            = true
   set dem file name       = dem
   set l2 smoothing factor = 0.000005
-  set l2 lower bound      = 0
-  set l2 upper bound      = 1
-  set bound void fraction = false
 end
 
 #---------------------------------------------------

examples/unresolved-cfd-dem/gas-solid-spouted-bed/gas-solid-spouted-bed.prm

@@ -70,7 +70,6 @@ subsection void fraction
   set qcm sphere equal cell volume = true
   set read dem                     = true
   set dem file name                = dem
-  set bound void fraction          = false
 end
 
 #---------------------------------------------------

examples/unresolved-cfd-dem/gas-solid-spouted-cylinder-bed/gas-solid-spouted-cylinder-bed.prm

@@ -77,7 +77,6 @@ subsection void fraction
   set qcm sphere equal cell volume = true
   set read dem                     = true
   set dem file name                = dem
-  set bound void fraction          = false
 end
 
 #---------------------------------------------------

include/fem-dem/cfd_dem_simulation_parameters.h

@@ -24,8 +24,8 @@ class CFDDEMSimulationParameters
   SimulationParameters<dim> cfd_parameters;
   DEMSolverParameters<dim>  dem_parameters;
 
-  std::shared_ptr<Parameters::VoidFraction<dim>> void_fraction;
-  Parameters::CFDDEM                             cfd_dem;
+  std::shared_ptr<Parameters::VoidFractionParameters<dim>> void_fraction;
+  Parameters::CFDDEM                                       cfd_dem;
 
   void
   declare(ParameterHandler             &prm,
@@ -34,7 +34,7 @@ class CFDDEMSimulationParameters
     cfd_parameters.declare(prm, size_of_subsections);
     dem_parameters.declare(prm);
 
-    void_fraction = std::make_shared<Parameters::VoidFraction<dim>>();
+    void_fraction = std::make_shared<Parameters::VoidFractionParameters<dim>>();
     void_fraction->declare_parameters(prm);
     Parameters::CFDDEM::declare_parameters(prm);
   }