[DEV] UFEM: Constant particle coagulation term

plugins/UFEM/src/UFEM/particles/ParticleConcentration.cpp

@@ -61,6 +61,16 @@ ParticleConcentration::ParticleConcentration(const std::string& name) :
     .pretty_name("Concentration Variable")
     .description("Variable for the particle concentration")
     .attach_trigger(boost::bind(&ParticleConcentration::trigger_set_expression, this));
+
+  options().add("source_term_tag", "concentration_source_terms")
+    .pretty_name("Source Term Tag")
+    .description("Tag for the field containing a source term")
+    .attach_trigger(boost::bind(&ParticleConcentration::trigger_set_expression, this));
+
+  options().add("source_term_variable", "c_src")
+    .pretty_name("Source Term Variable")
+    .description("Variable name for the source term")
+    .attach_trigger(boost::bind(&ParticleConcentration::trigger_set_expression, this));
 
   options().add("theta", m_theta)
     .pretty_name("Theta")
@@ -109,6 +119,11 @@ void ParticleConcentration::trigger_set_expression()
 
   // Particle concentration
   FieldVariable<1, ScalarField> c(options().value<std::string>("concentration_variable"), solution_tag());
+
+  // Source term
+  FieldVariable<2, ScalarField> s(options().value<std::string>("source_term_variable"), options().value<std::string>("source_term_tag"));
+
+//   std::cout << "options summary: " << options().value<std::string>("velocity_variable") << ", " <<  options().value<std::string>("velocity_tag") << ", " << options().value<std::string>("concentration_variable") << ", " << solution_tag() << ", " << options().value<std::string>("source_term_variable") << ", " << options().value<std::string>("source_term_tag") << std::endl;
 
   const Real theta = options().value<Real>("theta");
 
@@ -128,15 +143,16 @@ void ParticleConcentration::trigger_set_expression()
     AllowedElementTypesT(),
     group
     (
-      _A = _0, _T = _0,
+      _A = _0, _T = _0, _a = _0,
       compute_tau.apply(v, 0., lit(dt()), lit(tau_su)),
       element_quadrature
       (
         _A(c,c) +=  transpose(N(c) + lit(tau_su)*(v*nabla(c))) * (v*nabla(c) + divergence(v)*N(c)),
-        _T(c,c) +=  transpose(N(c) + lit(tau_su)*(v*nabla(c))) * N(c)
+        _T(c,c) +=  transpose(N(c) + lit(tau_su)*(v*nabla(c))) * N(c),
+        _a[c]   +=  transpose(N(c) + lit(tau_su)*(v*nabla(c))) * s
       ),
       system_matrix += invdt()*_T + theta*_A,
-      system_rhs += -_A*_x
+      system_rhs += -_A*_x + _a
     )
   ));
 

plugins/UFEM/src/UFEM/particles/Polydisperse.cpp

@@ -18,6 +18,8 @@
 
 #include "math/Consts.hpp"
 
+#include "mesh/Field.hpp"
+
 #include "physics/PhysModel.hpp"
 
 #include "solver/actions/Proto/ProtoAction.hpp"
@@ -39,6 +41,87 @@ common::ComponentBuilder < Polydisperse, common::Action, LibUFEMParticles > Poly
 
 ////////////////////////////////////////////////////////////////////////////////////////////
 
+namespace detail
+{
+
+Real pow_int(const Real in, const int pow)
+{
+  Real result = 1.;
+  for(int i = 0; i < pow; ++i)
+    result *= in;
+  return result;
+}
+
+/// Functor to compute the moment source terms
+struct MomentSourceFunctor : FunctionBase
+{
+  typedef void result_type;
+
+  MomentSourceFunctor(mesh::Region& region, const std::vector<std::string>& concentration_names, const std::vector<std::string>& weighted_volume_names) :
+    m_nb_phases(concentration_names.size()),
+    m_mat(2*m_nb_phases, 2*m_nb_phases),
+    m_rhs(2*m_nb_phases)
+  {
+    mesh::Mesh& mesh = common::find_parent_component<mesh::Mesh>(region);
+    mesh::Dictionary& dict = mesh.geometry_fields(); // TODO: Generalize to other dicts
+    m_source_field = &(dict.field("moment_source_terms"));
+
+    m_concentration_fields.reserve(m_nb_phases);
+    m_weighted_volume_fields.reserve(m_nb_phases);
+    for(Uint i = 0; i != m_nb_phases; ++i)
+    {
+      m_concentration_fields.push_back(&(dict.field(concentration_names[i])));
+      m_weighted_volume_fields.push_back(&(dict.field(weighted_volume_names[i])));
+    }
+  }
+
+  void operator()(const Uint node_idx)
+  {
+    const Uint nb_moments = 2*m_nb_phases;
+    m_rhs.setZero();
+    for(Uint alpha = 0; alpha != m_nb_phases; ++alpha) // For all phases
+    {
+      const Real n_alpha = (*m_concentration_fields[alpha])[node_idx][0];
+      const Real vol_alpha = (*m_weighted_volume_fields[alpha])[node_idx][0] / n_alpha;
+      for(int ki = 0; ki != nb_moments; ++ki) // For all moments
+      {
+        const Real k = static_cast<Real>(ki);
+        m_mat(ki, alpha) = (1.-k) * pow_int(vol_alpha, ki);
+        m_mat(ki, alpha+m_nb_phases) = k * pow_int(vol_alpha, ki-1);
+      }
+      for(Uint gamma = 0; gamma != m_nb_phases; ++gamma)
+      {
+        const Real n_gamma = (*m_concentration_fields[gamma])[node_idx][0];
+        const Real vol_gamma = (*m_weighted_volume_fields[gamma])[node_idx][0] / n_gamma;
+        for(int k = 0; k != nb_moments; ++k ) // For all moments
+        {
+          m_rhs[k] += (0.5*pow_int(vol_alpha + vol_gamma, k ) - pow_int(vol_alpha, k )) * (1.*n_alpha*n_gamma);
+        }
+      }
+    }
+
+//     std::cout << "source mat:\n" << m_mat << std::endl;
+//     std::cout << "source rhs: " << m_rhs.transpose() << std::endl;
+
+    Eigen::Map<RealVector> x(&(*m_source_field)[node_idx][0], nb_moments);
+    Eigen::FullPivLU<RealMatrix> lu(m_mat);
+    if(!lu.isInvertible())
+      throw common::FloatingPointError(FromHere(), "Matrix is not invertible.");
+    x = lu.solve(m_rhs);
+//     std::cout << "source sol: " << x.transpose() << std::endl;
+  }
+
+  mesh::Field* m_source_field;
+  std::vector<mesh::Field*> m_concentration_fields;
+  std::vector<mesh::Field*> m_weighted_volume_fields;
+
+  const Uint m_nb_phases;
+  RealMatrix m_mat;
+  RealVector m_rhs;
+};
+
+}
+
 Polydisperse::Polydisperse(const std::string& name) :
   solver::Action(name),
   m_nb_phases(0)
@@ -79,6 +162,8 @@ Polydisperse::Polydisperse(const std::string& name) :
   m_concentration_solver->set_solution_tag("particle_concentration_0");
   m_concentration_solver->options().set("concentration_variable", std::string("c_0"));
   m_concentration_solver->options().set("velocity_variable", std::string("vp_0"));
+  m_concentration_solver->options().set("source_term_tag", std::string("moment_source_terms"));
+  m_concentration_solver->options().set("source_term_variable", std::string("c_src_0"));
 }
 
 Polydisperse::~Polydisperse()
@@ -113,6 +198,8 @@ void Polydisperse::trigger_nb_phases()
   m_weighted_volume_variables.clear();
   m_tau_variables.clear();
   m_velocity_variables.clear();
+  m_concentration_src_variables.clear();
+  m_weighted_volume_src_variables.clear();
 
   m_ic_actions->clear();
 
@@ -125,6 +212,8 @@ void Polydisperse::trigger_nb_phases()
     m_weighted_volume_variables.push_back("zeta_"+phase_label);
     m_tau_variables.push_back("tau_"+phase_label);
     m_velocity_variables.push_back("vp_"+phase_label);
+    m_concentration_src_variables.push_back("c_src_"+phase_label);
+    m_weighted_volume_src_variables.push_back("zeta_src_"+phase_label);
 
     const Real init_c = initial_concentrations[i];
     const Real dp = initial_diameters[i];
@@ -138,6 +227,10 @@ void Polydisperse::trigger_nb_phases()
     Handle<ProtoAction> wv_action(m_ic_actions->create_initial_condition(m_weighted_volume_tags.back(), "cf3.solver.ProtoAction"));
     wv_action->set_expression(nodes_expression(wv = lit(init_wv)));
 
+    FieldVariable<2, ScalarField> c_src(m_concentration_src_variables.back(), "moment_source_terms");
+    Handle<ProtoAction> c_src_action(m_ic_actions->create_initial_condition(m_concentration_tags.back(), "cf3.solver.ProtoAction"));
+    c_src_action->set_expression(nodes_expression(c_src = lit(0.)));
+
     Handle<common::Component> tau_calc = m_compute_velocities->create_component("ComputeTau"+phase_label, "cf3.UFEM.particles.RelaxationTime");
     tau_calc->options().set("concentration_tag", m_concentration_tags.back());
     tau_calc->options().set("concentration_variable", m_concentration_variables.back());
@@ -157,6 +250,13 @@ void Polydisperse::trigger_nb_phases()
       eq_euler->options().set("compute_gradient", false);
     }
   }
+  // Outside the loop for per-block ordering
+  for(Uint i = 0; i != m_nb_phases; ++i)
+  {
+    FieldVariable<3, ScalarField> wv_src(m_weighted_volume_src_variables[i], "moment_source_terms");
+    Handle<ProtoAction> wv_src_action(m_ic_actions->create_initial_condition(m_weighted_volume_tags[i], "cf3.solver.ProtoAction"));
+    wv_src_action->set_expression(nodes_expression(wv_src = lit(0.)));
+  }
 }
 
 void Polydisperse::trigger_initial_conditions()
@@ -168,16 +268,29 @@ void Polydisperse::trigger_initial_conditions()
 void Polydisperse::execute()
 {
   m_compute_velocities->execute();
+
+  BOOST_FOREACH(const Handle<mesh::Region>& region, m_loop_regions)
+  {
+    detail::MomentSourceFunctor compute_source_terms(*region, m_concentration_tags, m_weighted_volume_tags);
+    for_each_node(*region, lit(compute_source_terms)(node_index));
+  }
+
   for(Uint i = 0; i != m_nb_phases; ++i)
   {
     m_concentration_solver->set_solution_tag(m_concentration_tags[i]);
     m_concentration_solver->options().set("concentration_variable", m_concentration_variables[i]);
     m_concentration_solver->options().set("velocity_variable", m_velocity_variables[i]);
+    m_concentration_solver->options().set("source_term_variable", m_concentration_src_variables[i]);
     m_concentration_solver->execute();
+//     std::cout << "concentration system for phase " << i << ":" << std::endl;
+//     Handle<math::LSS::System>(m_concentration_solver->get_child("LSS"))->print(std::cout);
 
     m_concentration_solver->set_solution_tag(m_weighted_volume_tags[i]);
     m_concentration_solver->options().set("concentration_variable", m_weighted_volume_variables[i]);
+    m_concentration_solver->options().set("source_term_variable", m_weighted_volume_src_variables[i]);
     m_concentration_solver->execute();
+//     std::cout << "weighted volume system for phase " << i << ":" << std::endl;
+//     Handle<math::LSS::System>(m_concentration_solver->get_child("LSS"))->print(std::cout);
   }
 }
 

plugins/UFEM/src/UFEM/particles/Polydisperse.hpp

@@ -53,6 +53,8 @@ class Polydisperse : public solver::Action
   std::vector<std::string> m_weighted_volume_variables;
   std::vector<std::string> m_tau_variables;
   std::vector<std::string> m_velocity_variables;
+  std::vector<std::string> m_concentration_src_variables;
+  std::vector<std::string> m_weighted_volume_src_variables;
 
   Handle<LSSAction> m_concentration_solver;
 };

plugins/UFEM/test/particles/CMakeLists.txt

@@ -2,4 +2,7 @@ coolfluid_add_test( ATEST atest-ufem-particles-taylor-green
                     PYTHON atest-ufem-particles-taylor-green.py)
 
 coolfluid_add_test( ATEST atest-ufem-particles-taylor-green-polydisperse
-                    PYTHON atest-ufem-particles-taylor-green-polydisperse.py)
+                    PYTHON atest-ufem-particles-taylor-green-polydisperse.py)
+
+coolfluid_add_test( ATEST atest-ufem-particles-polydisperse-uniform
+                    PYTHON atest-ufem-particles-polydisperse-uniform.py)

plugins/UFEM/test/particles/atest-ufem-particles-polydisperse-uniform.py

@@ -0,0 +1,92 @@
+import sys
+import coolfluid as cf
+import numpy as np
+import os
+from optparse import OptionParser
+import pylab as pl
+
+env = cf.Core.environment()
+
+# Global configuration
+env.assertion_throws = False
+env.assertion_backtrace = False
+env.exception_backtrace = False
+env.regist_signal_handlers = False
+env.log_level = 3
+
+nu = 1./5000.
+
+segs = 32
+D = 0.5
+Vs = 1./(4.*np.pi)
+Ua = 0.
+Va = 0.
+
+tau = 0.25
+beta = 3.
+
+dt = 0.001
+
+numsteps = 200
+write_interval = 50
+
+
+# Create the model and solvers
+model = cf.Core.root().create_component('NavierStokes', 'cf3.solver.ModelUnsteady')
+domain = model.create_domain()
+physics = model.create_physics('cf3.UFEM.NavierStokesPhysics')
+physics.options.dimension = 2
+solver = model.create_solver('cf3.UFEM.Solver')
+
+polydisp = solver.add_unsteady_solver('cf3.UFEM.particles.Polydisperse')
+polydisp.options.velocity_variable = 'Velocity'
+polydisp.options.velocity_tag = 'navier_stokes_u_solution'
+
+# Set up the physical constants
+physics.density = 1.
+physics.dynamic_viscosity = nu
+
+polydisp.initial_diameters = [(6./np.pi)**(1./3.), (12./np.pi)**(1./3.)]
+polydisp.initial_concentrations = [1., 0.0001]
+polydisp.nb_phases = 2
+
+# Create the mesh
+mesh = domain.create_component('Mesh', 'cf3.mesh.Mesh')
+mesh_generator = domain.create_component("mesh_generator","cf3.mesh.SimpleMeshGenerator")
+mesh_generator.options().set("mesh",mesh.uri())
+mesh_generator.options().set("nb_cells",[1, 1])
+mesh_generator.options().set("lengths",[1.,1.])
+mesh_generator.options().set("offsets",[0.,0.])
+mesh_generator.execute()
+
+polydisp.regions = [mesh.topology.interior.uri()]
+
+polydisp.children.ConcentrationSolver.LSS.SolutionStrategy.Parameters.linear_solver_type = 'Amesos'
+
+series_writer = solver.TimeLoop.create_component('TimeWriter', 'cf3.solver.actions.TimeSeriesWriter')
+writer = series_writer.create_component('Writer', 'cf3.mesh.VTKXML.Writer')
+writer.file = cf.URI('polydisperse-uniform-{iteration}.pvtu')
+writer.mesh = mesh
+series_writer.interval = write_interval
+
+solver.create_fields()
+writer.fields = [mesh.geometry.particle_concentration_1.uri(), mesh.geometry.weighted_particle_volume_1.uri()]
+
+# Time setup
+time = model.create_time()
+time.time_step = dt
+time.end_time = numsteps*dt
+
+model.simulate()
+model.print_timing_tree()
+
+c0 = np.array(mesh.geometry.particle_concentration_0)
+c1 = np.array(mesh.geometry.particle_concentration_1)
+zeta0 = np.array(mesh.geometry.weighted_particle_volume_0)
+zeta1 = np.array(mesh.geometry.weighted_particle_volume_1)
+
+print 'zeta:', zeta0[0,0], zeta1[0,0]
+print 'c:', c0[0,0], c1[0,0]
+print 'v:', zeta0[0,0]/c0[0,0], zeta1[0,0]/c1[0,0]
+
+domain.write_mesh(cf.URI('polydisperse-uniform-end.pvtu'))