Merge pull request #10 from tjhei/lla-sync

pull from upstream

Author: tjhei

.github/workflows/linux.yml

@@ -33,7 +33,7 @@ jobs:
       run: |
         sudo apt-get update
         sudo apt-get install -y software-properties-common
-        sudo add-apt-repository ppa:ginggs/deal.ii-9.5.1-backports
+        sudo add-apt-repository ppa:ginggs/deal.ii-9.6.0-backports
         sudo apt-get update
         sudo apt-get install -yq --no-install-recommends libdeal.ii-dev
     - name: compile
@@ -53,6 +53,9 @@ jobs:
     - uses: actions/checkout@v4
     - name: setup
       run: |
+        sudo apt-get update
+        sudo apt-get install -y software-properties-common
+        sudo add-apt-repository ppa:ginggs/deal.ii-9.6.0-backports
         sudo apt-get update
         sudo apt-get install -yq --no-install-recommends texlive-plain-generic texlive-base texlive-latex-recommended texlive-latex-base texlive-fonts-recommended texlive-bibtex-extra lmodern texlive-latex-extra texlive-science graphviz python3-pip libdeal.ii-dev doxygen latexmk biber inkscape
         doxygen --version

doc/modules/changes/20250618KiralyAgi renamed to doc/modules/changes/20250618_KiralyAgi

@@ -0,0 +1,6 @@
+Fixed: On Mac systems, ASPECT sometimes crashed at the very end of the
+program run with an exception triggered by the CGAL library that was
+caused by the settings for floating point exceptions. This is now
+fixed.
+<br>
+(Timo Heister, Wolfgang Bangerth, 2025/08/04)

doc/modules/changes/20250804_heisterbangerth

@@ -26,6 +26,7 @@
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/multithread_info.h>
 #include <deal.II/base/revision.h>
+#include <deal.II/base/scope_exit.h>
 #include <csignal>
 #include <string>
 #include <thread>
@@ -677,6 +678,16 @@ int main (int argc, char *argv[])
 
   // enable floating point exceptions
   feenableexcept(FE_DIVBYZERO|FE_INVALID);
+
+  // Make sure that at the end of the program run, we
+  // disable floating exceptions again. This fixes an issue
+  // of CGAL checking floating points exceptions during
+  // finalization and throwing an error:
+  ScopeExit fp_except_cleanup(
+    []()
+  {
+    fedisableexcept(FE_DIVBYZERO|FE_INVALID);
+  });
 #endif
 #endif
 
@@ -685,6 +696,12 @@ int main (int argc, char *argv[])
   // required for Numpy interop
   if (_import_array() < 0)
     AssertThrow(false, ExcMessage("Numpy init failed!"));
+
+  ScopeExit python_cleanup(
+    []()
+  {
+    Py_FinalizeEx();
+  });
 #endif
 
   std::string prm_name = "";
@@ -936,10 +953,5 @@ int main (int argc, char *argv[])
       return 1;
     }
 
-
-#ifdef ASPECT_WITH_PYTHON
-  Py_FinalizeEx();
-#endif
-
   return 0;
 }

source/main.cc

@@ -627,14 +627,8 @@ namespace aspect
       // For the moment add constraints from all plugins into one object, then
       // merge that object with the existing constraints (respecting the existing
       // constraints as more important).
-      //
-      // We initialize the constraints object with the locally
-      // relevant DoF indices as the set of constraints to store (also
-      // using the locally relevant indices for the "locally owned"
-      // DoFs). This way, plugins can fill information for all locally
-      // relevant DoFs and we do not need to communicate.
 #if DEAL_II_VERSION_GTE(9,7,0)
-      AffineConstraints<double> plugin_constraints(mesh_vertex_constraints.get_local_lines(),
+      AffineConstraints<double> plugin_constraints(mesh_deformation_dof_handler.locally_owned_dofs(),
                                                    mesh_vertex_constraints.get_local_lines());
 #else
       AffineConstraints<double> plugin_constraints(mesh_vertex_constraints.get_local_lines());
@@ -647,7 +641,7 @@ namespace aspect
           for (const auto &model : boundary_id.second)
             {
 #if DEAL_II_VERSION_GTE(9,7,0)
-              AffineConstraints<double> current_plugin_constraints(mesh_vertex_constraints.get_local_lines(),
+              AffineConstraints<double> current_plugin_constraints(mesh_deformation_dof_handler.locally_owned_dofs(),
                                                                    mesh_vertex_constraints.get_local_lines());
 #else
               AffineConstraints<double> current_plugin_constraints(mesh_vertex_constraints.get_local_lines());
@@ -656,7 +650,7 @@ namespace aspect
               model->compute_velocity_constraints_on_boundary(mesh_deformation_dof_handler,
                                                               current_plugin_constraints,
                                                               boundary_id_set);
-              if ((this->is_stokes_matrix_free()))
+              if (this->is_stokes_matrix_free())
                 {
                   mg_constrained_dofs.make_zero_boundary_constraints(mesh_deformation_dof_handler,
                                                                      boundary_id_set);

source/mesh_deformation/interface.cc

@@ -190,9 +190,8 @@ namespace aspect
 
 
     template <int dim>
-    // We need to pass the arguments by value, as this function can be called on a separate thread:
-    void Particles<dim>::writer (const std::string &filename, //NOLINT(performance-unnecessary-value-param)
-                                 const std::string &temporary_output_location, //NOLINT(performance-unnecessary-value-param)
+    void Particles<dim>::writer (const std::string &filename,
+                                 const std::string &temporary_output_location,
                                  const std::string &file_contents)
     {
       std::string tmp_filename = filename;

source/postprocess/particles.cc

@@ -52,7 +52,7 @@ subsection Material model
     set Reference temperature         = 1613
     set Thermal conductivity          = 1e-6
     set Thermal expansion coefficient = 2e-5
-    set Viscosity                     = 1
+    set Viscosity                     = 1e21
   end
 end
 

tests/adiabatic_conditions.prm

@@ -18,14 +18,14 @@ Number of degrees of freedom: 40,836 (25,090+3,201+12,545)
      Reference thermal expansion coefficient (1/K): 2e-05
      Reference specific heat capacity (J/(K*kg)):   1250
      Reference thermal conductivity (W/(m*K)):      1e-06
-     Reference viscosity (Pa*s):                    1
+     Reference viscosity (Pa*s):                    1e+21
      Reference thermal diffusivity (m^2/s):         2.42424e-13
      Rayleigh number:                               0
 
      Writing graphical output:           output-adiabatic_conditions/solution/solution-00000
-     RMS, max velocity:                  2.81e+04 m/s, 1.12e+05 m/s
+     RMS, max velocity:                  2.69e-17 m/s, 1.11e-16 m/s
      Temperature min/avg/max:            1613 K, 1613 K, 1613 K
-     Heat fluxes through boundary parts: 6721 W, -99.27 W, 0 W, 0 W
+     Heat fluxes through boundary parts: 9.544e-18 W, 8.936e-19 W, 0 W, 0 W
      Writing depth average:              output-adiabatic_conditions/depth_average
 
 Termination requested by criterion: end time

tests/adiabatic_conditions/depth_average.gnuplot

@@ -18,4 +18,4 @@
 # 18: Outward heat flux through boundary with indicator 1 ("top") (W)
 # 19: Outward heat flux through boundary with indicator 2 ("left") (W)
 # 20: Outward heat flux through boundary with indicator 3 ("right") (W)
-0 0.000000000000e+00 0.000000000000e+00 3072 28291 12545 0 14 15 15 output-adiabatic_conditions/solution/solution-00000 2.81189741e+04 1.12336096e+05 1.61300000e+03 1.61300000e+03 1.61300000e+03 6.72093540e+03 -9.92725149e+01 0.00000000e+00 0.00000000e+00 
+0 0.000000000000e+00 0.000000000000e+00 3072 28291 12545 0 14 15 15 output-adiabatic_conditions/solution/solution-00000 2.68839569e-17 1.10748720e-16 1.61300000e+03 1.61300000e+03 1.61300000e+03 9.54446020e-18 8.93580077e-19 0.00000000e+00 0.00000000e+00