Merge pull request #76 from tobiasleibner/threadsafe_expression

Thread-safe expression function

dune/stuff/functions/expression.hh

@@ -16,6 +16,7 @@
 
 #include <dune/stuff/common/configuration.hh>
 #include <dune/stuff/common/exceptions.hh>
+#include <dune/stuff/common/parallel/threadstorage.hh>
 
 #include "expression/base.hh"
 #include "interfaces.hh"
@@ -215,15 +216,15 @@ public:
     bool failure = false;
     std::string error_type;
     for (size_t rr = 0; rr < dimRange; ++rr) {
-      tmp_row_ = ret[rr];
+      *tmp_row_ = ret[rr];
       for (size_t cc = 0; cc < dimRangeCols; ++cc) {
-        if (DSC::isnan(tmp_row_[cc])) {
+        if (DSC::isnan(tmp_row_->operator[](cc))) {
           failure    = true;
           error_type = "NaN";
-        } else if (DSC::isinf(tmp_row_[cc])) {
+        } else if (DSC::isinf(tmp_row_->operator[](cc))) {
           failure    = true;
           error_type = "inf";
-        } else if (std::abs(tmp_row_[cc]) > (0.9 * std::numeric_limits<double>::max())) {
+        } else if (std::abs(tmp_row_->operator[](cc)) > (0.9 * std::numeric_limits<double>::max())) {
           failure    = true;
           error_type = "an unlikely value";
         }
@@ -236,13 +237,13 @@ public:
                          << function_->variable()
                          << "\n"
                          << "The expression of this functional is: "
-                         << function_->expression().at(0)
+                         << function_->expression().at(rr * dimRangeCols + cc)
                          << "\n"
                          << "You tried to evaluate it with:   xx = "
                          << xx
                          << "\n"
                          << "The result was:                       "
-                         << ret
+                         << tmp_row_->operator[](cc)
                          << "\n\n"
                          << "You can disable this check by defining DUNE_STUFF_FUNCTIONS_EXPRESSION_DISABLE_CHECKS\n");
       }
@@ -296,11 +297,11 @@ private:
   template <size_t rC>
   void evaluate_helper(const DomainType& xx, RangeType& ret, internal::ChooseVariant<rC>) const
   {
-    function_->evaluate(xx, tmp_vector_);
+    function_->evaluate(xx, *tmp_vector_);
     for (size_t rr = 0; rr < dimRange; ++rr) {
       auto& retRow = ret[rr];
       for (size_t cc = 0; cc < dimRangeCols; ++cc)
-        retRow[cc] = tmp_vector_[rr * dimRangeCols + cc];
+        retRow[cc] = (*tmp_vector_)[rr * dimRangeCols + cc];
     }
   } // ... evaluate_helper(...)
 
@@ -381,9 +382,8 @@ private:
   std::shared_ptr<const MathExpressionFunctionType> function_;
   size_t order_;
   std::string name_;
-  mutable FieldVector<RangeFieldType, dimRange * dimRangeCols> tmp_vector_;
-  mutable FieldVector<RangeFieldType, dimRangeCols> tmp_row_;
-  mutable FieldVector<RangeFieldType, dimDomain> tmp_gradient_row_;
+  mutable typename DS::PerThreadValue<FieldVector<RangeFieldType, dimRange * dimRangeCols>> tmp_vector_;
+  mutable typename DS::PerThreadValue<FieldVector<RangeFieldType, dimRangeCols>> tmp_row_;
   std::vector<std::vector<std::shared_ptr<const MathExpressionGradientType>>> gradients_;
 }; // class Expression
 

dune/stuff/functions/expression/base.hh

@@ -86,6 +86,7 @@ public:
   void evaluate(const Dune::FieldVector<DomainFieldType, dimDomain>& arg,
                 Dune::FieldVector<RangeFieldType, dimRange>& ret) const
   {
+    std::lock_guard<std::mutex> guard(mutex_);
     // copy arg
     for (typename Dune::FieldVector<DomainFieldType, dimDomain>::size_type ii = 0; ii < dimDomain; ++ii)
       *(arg_[ii]) = arg[ii];
@@ -99,6 +100,7 @@ public:
    */
   void evaluate(const Dune::DynamicVector<DomainFieldType>& arg, Dune::DynamicVector<RangeFieldType>& ret) const
   {
+    std::lock_guard<std::mutex> guard(mutex_);
     // check for sizes
     assert(arg.size() > 0);
     if (ret.size() != dimRange)
@@ -114,6 +116,7 @@ public:
   void evaluate(const Dune::FieldVector<DomainFieldType, dimDomain>& arg,
                 Dune::DynamicVector<RangeFieldType>& ret) const
   {
+    std::lock_guard<std::mutex> guard(mutex_);
     // check for sizes
     if (ret.size() != dimRange)
       ret = Dune::DynamicVector<RangeFieldType>(dimRange);
@@ -130,6 +133,7 @@ public:
    */
   void evaluate(const Dune::DynamicVector<DomainFieldType>& arg, Dune::FieldVector<RangeFieldType, dimRange>& ret) const
   {
+    std::lock_guard<std::mutex> guard(mutex_);
     assert(arg.size() > 0);
     // copy arg
     for (size_t ii = 0; ii < std::min(dimDomain, arg.size()); ++ii)
@@ -208,6 +212,7 @@ private:
   RVar* var_arg_[dimDomain];
   RVar* vararray_[dimDomain];
   ROperation* op_[dimRange];
+  mutable std::mutex mutex_;
 }; // class MathExpressionBase
 
 } // namespace Functions