[FIX] UFEM: NaN problems

coolfluid · Apr 23, 2014 · b2d2bea · b2d2bea
1 parent 8b7f678
commit b2d2bea
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 7 deletions.
diff --git a/plugins/UFEM/src/UFEM/particles/Polydisperse.cpp b/plugins/UFEM/src/UFEM/particles/Polydisperse.cpp
@@ -184,7 +184,8 @@ struct DNSCollisionKernel
     {
       std::stringstream msg;
       msg << "Unknown combination of eigenvalues: " << ev[0] << ", " << ev[1] << ", " << ev[2];
-      throw common::ShouldNotBeHere(FromHere(), msg.str());
+      CFdebug << msg.str() << CFendl;
+      beta2 = 0.;
     }
 
     return beta1 + beta2;

diff --git a/plugins/UFEM/src/UFEM/particles/RelaxationTime.cpp b/plugins/UFEM/src/UFEM/particles/RelaxationTime.cpp
@@ -15,7 +15,6 @@
 #include "common/OptionArray.hpp"
 #include "common/PropertyList.hpp"
 
-#include "math/Consts.hpp"
 #include "math/LSS/SolveLSS.hpp"
 #include "math/LSS/ZeroLSS.hpp"
 #include "mesh/LagrangeP1/ElementTypes.hpp"
@@ -36,7 +35,6 @@ namespace particles {
 
 using namespace solver::actions::Proto;
 using boost::proto::lit;
-using math::Consts::pi;
 
 ////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -67,10 +65,12 @@ RelaxationTime::RelaxationTime(const std::string& name) :
     .pretty_name("Tau Variable")
     .description("Variable for the relaxation time");
 
-  options().add("reference_volume", m_reference_volume)
+  r2.m_reference_volume = 1.;
+
+  options().add("reference_volume", r2.m_reference_volume)
     .pretty_name("Reference Volume")
     .description("Reference volume, all particle volumes are divided by this")
-    .link_to(&m_reference_volume);
+    .link_to(&(r2.m_reference_volume));
 }
 
 RelaxationTime::~RelaxationTime()
@@ -92,7 +92,7 @@ void RelaxationTime::on_regions_set()
   (
     group
     (
-      tau = lit(2./9.)*rho_p/mu*_std_pow(lit(3./(4.*pi()))*zeta/c*lit(m_reference_volume), lit(2./3.))
+      tau = lit(2./9.)*rho_p/mu*lit(r2)(c, zeta)
     )
   ));
 }

diff --git a/plugins/UFEM/src/UFEM/particles/RelaxationTime.hpp b/plugins/UFEM/src/UFEM/particles/RelaxationTime.hpp
@@ -7,6 +7,8 @@
 #ifndef cf3_UFEM_RelaxationTime_hpp
 #define cf3_UFEM_RelaxationTime_hpp
 
+#include "math/Consts.hpp"
+
 #include "solver/actions/Proto/ProtoAction.hpp"
 #include "solver/actions/Proto/Expression.hpp"
 
@@ -16,6 +18,23 @@ namespace cf3 {
 namespace UFEM {
 namespace particles {
 
+// Functor to compute the squared particle radius
+struct SquaredParticleRadius : solver::actions::Proto::FunctionBase
+{
+  typedef Real result_type;
+
+  Real operator()(const Real c, const Real zeta)
+  {
+    const Real r2 = ::pow(3./(4.*math::Consts::pi())*zeta/c*m_reference_volume, 2./3.);
+    if(!std::isfinite(r2))
+      return 0.;
+
+    return r2;
+  }
+
+  Real m_reference_volume;
+};
+
 /// Compute the particle relaxation time
 class RelaxationTime : public solver::actions::Proto::ProtoAction
 {
@@ -33,7 +52,7 @@ class RelaxationTime : public solver::actions::Proto::ProtoAction
 private:
   virtual void on_regions_set();
 
-  Real m_reference_volume;
+  SquaredParticleRadius r2;
 };
 
 } // particles