Merge pull request #204 from tamiko/poison_invalid_data

StateVector: poison invalid states in debug mode

source/hyperbolic_module.template.h

@@ -1220,6 +1220,10 @@ namespace ryujin
       }
     }
 
+    /* In debug mode poison precomputed values: */
+    Vectors::debug_poison_precomputed_values<Description>(new_state_vector,
+                                                          *offline_data_);
+
     /* Return the time step size tau: */
     return tau;
   }

source/initial_values.template.h

@@ -248,21 +248,6 @@ namespace ryujin
 
     U.update_ghost_values();
 
-#ifdef DEBUG
-    /* Poison constrained degrees of freedom: */
-    {
-      const unsigned int n_owned = offline_data_->n_locally_owned();
-      const auto &partitioner = offline_data_->scalar_partitioner();
-      for (unsigned int i = 0; i < n_owned; ++i) {
-        if (offline_data_->affine_constraints().is_constrained(
-                partitioner->local_to_global(i)))
-          U.write_tensor(dealii::Tensor<1, dim + 2, Number>() *
-                             std::numeric_limits<Number>::signaling_NaN(),
-                         i);
-      }
-    }
-#endif
-
     return U;
   }
 

source/mpi_ensemble.cc

@@ -16,11 +16,15 @@ namespace ryujin
       , n_ensembles_(1)
       , ensemble_rank_(0)
       , n_ensemble_ranks_(1)
+      , initialized_(false)
   {
   }
 
   MPIEnsemble::~MPIEnsemble()
   {
+    if (!initialized_)
+      return;
+
     MPI_Group_free(&world_group_);
     for (auto &it : ensemble_groups_)
       MPI_Group_free(&it);
@@ -34,6 +38,8 @@ namespace ryujin
   void MPIEnsemble::prepare(const int n_ensembles /* = 1 */,
                             const bool global_synchronization /* = true */)
   {
+    Assert(initialized_ == false, dealii::ExcInternalError());
+
     n_ensembles_ = n_ensembles;
     global_synchronization_ = global_synchronization;
 
@@ -113,5 +119,7 @@ namespace ryujin
               << ", e = " << n_ensemble_ranks_ << ", p = " << n_peer_ranks
               << ") -> belongig to ensemble: " << ensemble_ << std::endl;
 #endif
+
+    initialized_ = true;
   }
 } /* namespace ryujin */

source/mpi_ensemble.h

@@ -135,6 +135,8 @@ namespace ryujin
     int ensemble_rank_;
     int n_ensemble_ranks_;
 
+    bool initialized_;
+
     MPI_Group world_group_;
     std::vector<MPI_Group> ensemble_groups_;
     MPI_Group ensemble_leader_group_;

source/solution_transfer.h

@@ -13,8 +13,6 @@
 
 #include <deal.II/base/parameter_acceptor.h>
 
-#include <optional>
-
 namespace ryujin
 {
   /**

source/state_vector.h

@@ -51,6 +51,67 @@ namespace ryujin
         BlockVector<Number> /*parabolic state vector*/>;
 
 
+    template <
+        typename Description,
+        int dim,
+        typename Number,
+        typename View =
+            typename Description::template HyperbolicSystemView<dim, Number>,
+        int problem_dimension = View::problem_dimension,
+        int prec_dimension = View::n_precomputed_values>
+    void debug_poison_constrained_dofs(
+        StateVector<Number, problem_dimension, prec_dimension> &state_vector
+        [[maybe_unused]],
+        const OfflineData<dim, Number> &offline_data [[maybe_unused]])
+    {
+#ifdef DEBUG
+      auto &[U, precomputed, V] = state_vector;
+
+      const unsigned int n_owned = offline_data.n_locally_owned();
+      const auto &partitioner = offline_data.scalar_partitioner();
+
+      for (unsigned int i = 0; i < n_owned; ++i) {
+        if (!offline_data.affine_constraints().is_constrained(
+                partitioner->local_to_global(i)))
+          continue;
+        constexpr auto nan = std::numeric_limits<Number>::signaling_NaN();
+        U.write_tensor(dealii::Tensor<1, problem_dimension, Number>() * nan, i);
+      }
+#endif
+    }
+
+
+    template <
+        typename Description,
+        int dim,
+        typename Number,
+        typename View =
+            typename Description::template HyperbolicSystemView<dim, Number>,
+        int problem_dimension = View::problem_dimension,
+        int prec_dimension = View::n_precomputed_values>
+    void debug_poison_precomputed_values(
+        StateVector<Number, problem_dimension, prec_dimension> &state_vector
+        [[maybe_unused]],
+        const OfflineData<dim, Number> &offline_data [[maybe_unused]])
+    {
+#ifdef DEBUG
+      auto &[U, precomputed, V] = state_vector;
+
+      constexpr auto nan = std::numeric_limits<Number>::signaling_NaN();
+      const unsigned int n_owned = offline_data.n_locally_owned();
+      const auto block_size = offline_data.n_parabolic_state_vectors();
+
+      for (unsigned int i = 0; i < n_owned; ++i) {
+        precomputed.write_tensor(
+            dealii::Tensor<1, prec_dimension, Number>() * nan, i);
+        for (unsigned int b = 0; b < block_size; ++b) {
+          V.block(b).local_element(i) = nan;
+        }
+      }
+#endif
+    }
+
+
     /**
      * Helper function that (re)initializes all components of a StateVector
      * to proper sizes.
@@ -61,10 +122,10 @@ namespace ryujin
         typename Number,
         typename View =
             typename Description::template HyperbolicSystemView<dim, Number>,
-        int problem_dim = View::problem_dimension,
-        int prec_dim = View::n_precomputed_values>
+        int problem_dimension = View::problem_dimension,
+        int prec_dimension = View::n_precomputed_values>
     void reinit_state_vector(
-        StateVector<Number, problem_dim, prec_dim> &state_vector,
+        StateVector<Number, problem_dimension, prec_dimension> &state_vector,
         const OfflineData<dim, Number> &offline_data)
     {
       auto &[U, precomputed, V] = state_vector;
@@ -76,6 +137,23 @@ namespace ryujin
       for (unsigned int i = 0; i < block_size; ++i) {
         V.block(i).reinit(offline_data.scalar_partitioner());
       }
+
+#ifdef DEBUG
+      /* Poison all vectors: */
+      using state_type = typename View::state_type;
+
+      constexpr auto nan = std::numeric_limits<Number>::signaling_NaN();
+
+      const unsigned int n_owned = offline_data.n_locally_owned();
+      for (unsigned int i = 0; i < n_owned; ++i) {
+        U.write_tensor(state_type{} * nan, i);
+        precomputed.write_tensor(
+            dealii::Tensor<1, prec_dimension, Number>() * nan, i);
+        for (unsigned int b = 0; b < block_size; ++b) {
+          V.block(b).local_element(i) = nan;
+        }
+      }
+#endif
     }
   } // namespace Vectors
 

source/time_loop.template.h

@@ -7,6 +7,7 @@
 
 #include "scope.h"
 #include "solution_transfer.h"
+#include "state_vector.h"
 #include "time_loop.h"
 #include "version_info.h"
 
@@ -306,6 +307,15 @@ namespace ryujin
       }
     }
 
+    /*
+     * In debug mode poison constrained degrees of freedom and precomputed
+     * values:
+     */
+    Vectors::debug_poison_constrained_dofs<Description>(state_vector,
+                                                        offline_data_);
+    Vectors::debug_poison_precomputed_values<Description>(state_vector,
+                                                          offline_data_);
+
     unsigned int cycle = 1;
     Number last_terminal_output = (terminal_update_interval_ == Number(0.)
                                        ? std::numeric_limits<Number>::max()
@@ -338,8 +348,6 @@ namespace ryujin
       /* Perform various tasks whenever we reach a timer tick: */
 
       if (t >= relax * timer_cycle * timer_granularity_) {
-        output(state_vector, base_name_ensemble_ + "-solution", t, timer_cycle);
-
         if (enable_compute_error_) {
           StateVector analytic;
           {
@@ -354,12 +362,20 @@ namespace ryujin
             std::get<0>(analytic) =
                 initial_values_.interpolate_hyperbolic_vector(t);
           }
+
+          /*
+           * FIXME: a call to output() will also write a checkpoint (if
+           * enabled). So as a workaround we simply call the output()
+           * function for the analytic solution first...
+           */
           output(analytic,
                  base_name_ensemble_ + "-analytic_solution",
                  t,
                  timer_cycle);
         }
 
+        output(state_vector, base_name_ensemble_ + "-solution", t, timer_cycle);
+
         if (enable_compute_quantities_ &&
             (timer_cycle % timer_compute_quantities_multiplier_ == 0)) {
           Scope scope(computing_timer_,
@@ -541,7 +557,7 @@ namespace ryujin
     Vectors::reinit_state_vector<Description>(state_vector, offline_data_);
 
     SolutionTransfer<Description, dim, Number> solution_transfer(
-        mpi_ensemble_.ensemble_communicator(),
+        mpi_ensemble_,
         /* we need write access: */ discretization_.triangulation(),
         offline_data_,
         hyperbolic_system_,