Merge branch 'develop' into feature/adaptive_resolution_universe

# Conflicts:
#	src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C

.gitignore

@@ -5,6 +5,8 @@
 log_*
 log.*
 *~
+*.swp
+*.swo
 
 **/linux*Gcc*/
 **/.vscode

README.md

@@ -1,6 +1,6 @@
 # CFDEMcoupling
 
-CFDEM®coupling stands for Computational Fluid Dynamics (CFD) - Discrete Element Method (DEM) coupling. It combines the open source packages OpenFOAM® (CFD) and LIGGGHTS® (DEM) to simulate particle-laden flows. CFDEM®coupling is part of the [CFDEM®project](www.cfdem.com).
+CFDEM®coupling stands for Computational Fluid Dynamics (CFD) - Discrete Element Method (DEM) coupling. It combines the open source packages OpenFOAM® (CFD) and LIGGGHTS® (DEM) to simulate particle-laden flows. CFDEM®coupling is part of the [CFDEM®project](https://www.cfdem.com).
 
 [![CircleCI](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling.svg?style=shield&circle-token=e4b6af30d3aa7aee109d206116f01600bf9ee9c6)](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling)
 [![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.html)
@@ -8,7 +8,7 @@ CFDEM®coupling stands for Computational Fluid Dynamics (CFD) - Discrete Element
 ## Disclaimer
 
 > This is an academic adaptation of the CFDEM®coupling software package, released by the
-[Department of Particulate Flow Modelling at Johannes Kepler University in Linz, Austria.](www.jku.at/pfm)
+[Department of Particulate Flow Modelling at Johannes Kepler University in Linz, Austria.](https://www.jku.at/pfm)
 > LIGGGHTS® and CFDEM® are registered trademarks, and this offering is not approved or
 endorsed by DCS Computing GmbH, the official producer of the LIGGGHTS® and CFDEM®coupling software.
 > This offering is not approved or endorsed by OpenCFD Limited, producer and distributor of the OpenFOAM software via www.openfoam.com, and owner of the OPENFOAM®  and OpenCFD®  trade marks.
@@ -29,5 +29,5 @@ endorsed by DCS Computing GmbH, the official producer of the LIGGGHTS® and CFDE
 - Some parts of CFDEM®coupling are based on OpenFOAM® and Copyright on these
   parts is held by the OpenFOAM® Foundation (www.openfoam.org)
   and potentially other parties.
-- Some parts of CFDEM®coupling are contributied by other parties, which are
+- Some parts of CFDEM®coupling are contributed by other parties, which are
   holding the Copyright. This is listed in each file of the distribution.

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -54,6 +54,9 @@
 
     thermo.correct();
 
+    Info << "Qsource" << max(Qsource).value() << " " << min(Qsource).value() << endl;
+    Info << "QCoeff" << max(QCoeff).value() << " " << min(QCoeff).value() << endl;
+    Info << "Cpv" << max(Cpv).value() << " " << min(Cpv).value() << endl;
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
 
     particleCloud.clockM().start(31,"energySolve");

applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H

@@ -10,8 +10,8 @@ thCond=particleCloud.thermCondM().thermCond();
 Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
 // correct source for the thermodynamic reference temperature
-dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+// dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
-Qsource += QCoeff*Tref;
+// Qsource += QCoeff*Tref;
 
 fvScalarMatrix EEqn
 (
@@ -49,7 +49,13 @@ fvScalarMatrix EEqn
 
     thermo.correct();
 
-    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+    Info << "Qsource :" << max(Qsource).value() << " " << min(Qsource).value() << endl;
-    Info << "he min/max : " << max(he).value() << "  " << min(he).value() << endl;
+    Info << "QCoeff :" << max(QCoeff).value() << " " << min(QCoeff).value() << endl;
+    Info << "Cpv :" << max(Cpv).value() << " " << min(Cpv).value() << endl;
+    Info<<  "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+    Info << "he max/min : " << max(he).value() << "  " << min(he).value() << endl;
 
+    particleCloud.clockM().start(31,"energySolve");
+    particleCloud.solve();
+    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C

@@ -52,7 +52,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-//    #include "postProcess.H"
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -61,7 +60,6 @@ int main(int argc, char *argv[])
     #include "createRDeltaT.H"
 
     #include "createFields.H"
-    #include "createFieldRefs.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/cfdemSolverRhoPimpleChem/createFields.H

@@ -8,7 +8,6 @@
     );
 
     rhoReactionThermo& thermo = combustion->thermo();
-
     thermo.validate(args.executable(), "h", "e");
 
     basicSpecieMixture& composition = thermo.composition();
@@ -34,6 +33,8 @@
     }
 
     volScalarField& p = thermo.p();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -57,9 +58,6 @@
     );
 
     //  kinematic fields
-
-
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -88,18 +86,8 @@
         mesh
     );
 
-    volScalarField rhoeps
+    volScalarField rhoeps ("rhoeps", rho*voidfraction);
-    (
+
-        IOobject
-        (
-            "rhoeps",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        rho*voidfraction
-    );
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Qsource
@@ -282,9 +270,18 @@
         mesh,
         dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
     );
+
+    volScalarField dSauter
+    (
+        IOobject
+        (
+            "dSauter",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 1, 0, 0, 0,0,0),0)
+    );
 //===============================
-
-
-
- //   singlePhaseTransportModel laminarTransport(U, phi);
-

doc/_themes/lammps_theme/breadcrumbs.html

@@ -0,0 +1,34 @@
+<div role="navigation" aria-label="breadcrumbs navigation">
+  <ul class="wy-breadcrumbs">
+    <li><a href="{{ pathto(master_doc) }}">Docs</a> &raquo;</li>
+      {% for doc in parents %}
+          <li><a href="{{ doc.link|e }}">{{ doc.title }}</a> &raquo;</li>
+      {% endfor %}
+    <li>{{ title }}</li>
+      <li class="wy-breadcrumbs-aside">
+        {% if pagename != "search" %}
+          {% if display_github %}
+            <a href="https://{{ github_host|default("github.com") }}/{{ github_user }}/{{ github_repo }}/blob/{{ github_version }}{{ conf_py_path }}{{ pagename }}{{ source_suffix }}" class="fa fa-github"> Edit on GitHub</a>
+          {% elif display_bitbucket %}
+            <a href="https://bitbucket.org/{{ bitbucket_user }}/{{ bitbucket_repo }}/src/{{ bitbucket_version}}{{ conf_py_path }}{{ pagename }}{{ source_suffix }}" class="fa fa-bitbucket"> Edit on Bitbucket</a>
+          {% elif show_source and source_url_prefix %}
+            <a href="{{ source_url_prefix }}{{ pagename }}{{ source_suffix }}">View page source</a>
+          {% elif show_source and has_source and sourcename %}
+            <a href="{{ pathto('_sources/' + sourcename, true)|e }}" rel="nofollow"> View page source</a>
+          {% endif %}
+            <a href="https://www.cfdem.com">Website</a>
+        {% endif %}
+      </li>
+  </ul>
+  <hr/>
+  {% if next or prev %}
+    <div class="rst-footer-buttons" style="margin-bottom: 1em" role="navigation" aria-label="footer navigation">
+      {% if next %}
+        <a href="{{ next.link|e }}" class="btn btn-neutral float-right" title="{{ next.title|striptags|e }}" accesskey="n">Next <span class="fa fa-arrow-circle-right"></span></a>
+      {% endif %}
+      {% if prev %}
+        <a href="{{ prev.link|e }}" class="btn btn-neutral" title="{{ prev.title|striptags|e }}" accesskey="p"><span class="fa fa-arrow-circle-left"></span> Previous</a>
+      {% endif %}
+    </div>
+  {% endif %}
+</div>

doc/_themes/lammps_theme/footer.html

@@ -0,0 +1,36 @@
+<footer>
+  {% if next or prev %}
+    <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
+      {% if next %}
+        <a href="{{ next.link|e }}" class="btn btn-neutral float-right" title="{{ next.title|striptags|e }}" accesskey="n">Next <span class="fa fa-arrow-circle-right"></span></a>
+      {% endif %}
+      {% if prev %}
+        <a href="{{ prev.link|e }}" class="btn btn-neutral" title="{{ prev.title|striptags|e }}" accesskey="p"><span class="fa fa-arrow-circle-left"></span> Previous</a>
+      {% endif %}
+    </div>
+  {% endif %}
+
+  <hr/>
+
+  <div role="contentinfo">
+    <p>
+    {%- if show_copyright %}
+      {%- if hasdoc('copyright') %}
+        {% trans path=pathto('copyright'), copyright=copyright|e %}&copy; <a href="{{ path }}">Copyright</a> {{ copyright }}.{% endtrans %}
+      {%- else %}
+        {% trans copyright=copyright|e %}&copy; Copyright {{ copyright }}.{% endtrans %}
+      {%- endif %}
+    {%- endif %}
+
+    {%- if last_updated %}
+      {% trans last_updated=last_updated|e %}Last updated on {{ last_updated }}.{% endtrans %}
+    {%- endif %}
+    </p>
+  </div>
+
+  {%- if show_sphinx %}
+  {% trans %}Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/snide/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>{% endtrans %}.
+  {%- endif %}
+
+</footer>
+

doc/_themes/lammps_theme/layout.html

@@ -0,0 +1,195 @@
+{# TEMPLATE VAR SETTINGS #}
+{%- set url_root = pathto('', 1) %}
+{%- if url_root == '#' %}{% set url_root = '' %}{% endif %}
+{%- if not embedded and docstitle %}
+  {%- set titlesuffix = " — "|safe + docstitle|e %}
+{%- else %}
+  {%- set titlesuffix = "" %}
+{%- endif %}
+
+<!DOCTYPE html>
+<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
+<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
+<head>
+  <meta charset="utf-8">
+  {{ metatags }}
+  <meta name="viewport" content="width=device-width, initial-scale=1.0">
+  {% block htmltitle %}
+  <title>{{ title|striptags|e }}{{ titlesuffix }}</title>
+  {% endblock %}
+
+  {# FAVICON #}
+  {% if favicon %}
+    <link rel="shortcut icon" href="{{ pathto('_static/' + favicon, 1) }}"/>
+  {% endif %}
+
+  {# CSS #}
+
+  {# OPENSEARCH #}
+  {% if not embedded %}
+    {% if use_opensearch %}
+      <link rel="search" type="application/opensearchdescription+xml" title="{% trans docstitle=docstitle|e %}Search within {{ docstitle }}{% endtrans %}" href="{{ pathto('_static/opensearch.xml', 1) }}"/>
+    {% endif %}
+
+  {% endif %}
+
+  {# RTD hosts this file, so just load on non RTD builds #}
+  {% if not READTHEDOCS %}
+    <link rel="stylesheet" href="{{ pathto('_static/' + style, 1) }}" type="text/css" />
+  {% endif %}
+
+  {% for cssfile in css_files %}
+    <link rel="stylesheet" href="{{ pathto(cssfile, 1) }}" type="text/css" />
+  {% endfor %}
+
+  {% for cssfile in extra_css_files %}
+    <link rel="stylesheet" href="{{ pathto(cssfile, 1) }}" type="text/css" />
+  {% endfor %}
+
+  {%- block linktags %}
+    {%- if hasdoc('about') %}
+        <link rel="author" title="{{ _('About these documents') }}"
+              href="{{ pathto('about') }}"/>
+    {%- endif %}
+    {%- if hasdoc('genindex') %}
+        <link rel="index" title="{{ _('Index') }}"
+              href="{{ pathto('genindex') }}"/>
+    {%- endif %}
+    {%- if hasdoc('search') %}
+        <link rel="search" title="{{ _('Search') }}" href="{{ pathto('search') }}"/>
+    {%- endif %}
+    {%- if hasdoc('copyright') %}
+        <link rel="copyright" title="{{ _('Copyright') }}" href="{{ pathto('copyright') }}"/>
+    {%- endif %}
+    <link rel="top" title="{{ docstitle|e }}" href="{{ pathto('index') }}"/>
+    {%- if parents %}
+        <link rel="up" title="{{ parents[-1].title|striptags|e }}" href="{{ parents[-1].link|e }}"/>
+    {%- endif %}
+    {%- if next %}
+        <link rel="next" title="{{ next.title|striptags|e }}" href="{{ next.link|e }}"/>
+    {%- endif %}
+    {%- if prev %}
+        <link rel="prev" title="{{ prev.title|striptags|e }}" href="{{ prev.link|e }}"/>
+    {%- endif %}
+  {%- endblock %}
+  {%- block extrahead %} {% endblock %}
+
+  {# Keep modernizr in head - http://modernizr.com/docs/#installing #}
+  <script src="{{ pathto('_static/js/modernizr.min.js', 1) }}"></script>
+
+</head>
+
+<body class="wy-body-for-nav" role="document">
+
+  <div class="wy-grid-for-nav">
+
+    {# SIDE NAV, TOGGLES ON MOBILE #}
+    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
+      <div class="wy-side-nav-search">
+        {% block sidebartitle %}
+
+        {% if logo and theme_logo_only %}
+          <a href="{{ pathto(master_doc) }}">
+        {% else %}
+          <a href="{{ pathto(master_doc) }}" class="icon icon-home"> {{ project }}
+        {% endif %}
+
+        {% if logo %}
+          {# Not strictly valid HTML, but it's the only way to display/scale it properly, without weird scripting or heaps of work #}
+          <img src="{{ pathto('_static/' + logo, 1) }}" class="logo" />
+        {% endif %}
+        </a>
+
+        {% if theme_display_version %}
+          {%- set nav_version = version %}
+          {% if READTHEDOCS and current_version %}
+            {%- set nav_version = current_version %}
+          {% endif %}
+          {% if nav_version %}
+            <div class="version">
+              {{ nav_version }}
+            </div>
+          {% endif %}
+        {% endif %}
+
+        {% include "searchbox.html" %}
+
+        {% endblock %}
+      </div>
+
+      <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
+        {% block menu %}
+          {% set toctree = toctree(maxdepth=4, collapse=theme_collapse_navigation, includehidden=True) %}
+          {% if toctree %}
+              {{ toctree }}
+          {% else %}
+              <!-- Local TOC -->
+              <div class="local-toc">{{ toc }}</div>
+          {% endif %}
+        {% endblock %}
+      </div>
+      &nbsp;
+    </nav>
+
+    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">
+
+      {# MOBILE NAV, TRIGGLES SIDE NAV ON TOGGLE #}
+      <nav class="wy-nav-top" role="navigation" aria-label="top navigation">
+        <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
+        <a href="{{ pathto(master_doc) }}">{{ project }}</a>
+      </nav>
+
+
+      {# PAGE CONTENT #}
+      <div class="wy-nav-content">
+        <div class="rst-content">
+          {% include "breadcrumbs.html" %}
+          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
+           <div itemprop="articleBody">
+            {% block body %}{% endblock %}
+           </div>
+          </div>
+          {% include "footer.html" %}
+        </div>
+      </div>
+
+    </section>
+
+  </div>
+  {% include "versions.html" %}
+
+  {% if not embedded %}
+
+    <script type="text/javascript">
+        var DOCUMENTATION_OPTIONS = {
+            URL_ROOT:'{{ url_root }}',
+            VERSION:'{{ release|e }}',
+            COLLAPSE_INDEX:false,
+            FILE_SUFFIX:'{{ '' if no_search_suffix else file_suffix }}',
+            HAS_SOURCE:  {{ has_source|lower }}
+        };
+    </script>
+    {%- for scriptfile in script_files %}
+      <script type="text/javascript" src="{{ pathto(scriptfile, 1) }}"></script>
+    {%- endfor %}
+
+  {% endif %}
+
+  {# RTD hosts this file, so just load on non RTD builds #}
+  {% if not READTHEDOCS %}
+    <script type="text/javascript" src="{{ pathto('_static/js/theme.js', 1) }}"></script>
+  {% endif %}
+
+  {# STICKY NAVIGATION #}
+  {% if theme_sticky_navigation %}
+  <script type="text/javascript">
+      jQuery(function () {
+          SphinxRtdTheme.StickyNav.enable();
+      });
+  </script>
+  {% endif %}
+
+  {%- block footer %} {% endblock %}
+
+</body>
+</html>

doc/_themes/lammps_theme/layout_old.html

@@ -0,0 +1,205 @@
+{#
+    basic/layout.html
+    ~~~~~~~~~~~~~~~~~
+
+    Master layout template for Sphinx themes.
+
+    :copyright: Copyright 2007-2013 by the Sphinx team, see AUTHORS.
+    :license: BSD, see LICENSE for details.
+#}
+{%- block doctype -%}
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
+  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
+{%- endblock %}
+{%- set reldelim1 = reldelim1 is not defined and ' &raquo;' or reldelim1 %}
+{%- set reldelim2 = reldelim2 is not defined and ' |' or reldelim2 %}
+{%- set render_sidebar = (not embedded) and (not theme_nosidebar|tobool) and
+                         (sidebars != []) %}
+{%- set url_root = pathto('', 1) %}
+{# XXX necessary? #}
+{%- if url_root == '#' %}{% set url_root = '' %}{% endif %}
+{%- if not embedded and docstitle %}
+  {%- set titlesuffix = " &mdash; "|safe + docstitle|e %}
+{%- else %}
+  {%- set titlesuffix = "" %}
+{%- endif %}
+
+{%- macro relbar() %}
+    <div class="related">
+      <h3>{{ _('Navigation') }}</h3>
+      <ul>
+        {%- for rellink in rellinks %}
+        <li class="right" {% if loop.first %}style="margin-right: 10px"{% endif %}>
+          <a href="{{ pathto(rellink[0]) }}" title="{{ rellink[1]|striptags|e }}"
+             {{ accesskey(rellink[2]) }}>{{ rellink[3] }}</a>
+          {%- if not loop.first %}{{ reldelim2 }}{% endif %}</li>
+        {%- endfor %}
+        {%- block rootrellink %}
+        <li><a href="{{ pathto(master_doc) }}">{{ shorttitle|e }}</a>{{ reldelim1 }}</li>
+        {%- endblock %}
+        {%- for parent in parents %}
+          <li><a href="{{ parent.link|e }}" {% if loop.last %}{{ accesskey("U") }}{% endif %}>{{ parent.title }}</a>{{ reldelim1 }}</li>
+        {%- endfor %}
+        {%- block relbaritems %} {% endblock %}
+      </ul>
+    </div>
+{%- endmacro %}
+
+{%- macro sidebar() %}
+      {%- if render_sidebar %}
+      <div class="sphinxsidebar">
+        <div class="sphinxsidebarwrapper">
+          {%- block sidebarlogo %}
+          {%- if logo %}
+            <p class="logo"><a href="{{ pathto(master_doc) }}">
+              <img class="logo" src="{{ pathto('_static/' + logo, 1) }}" alt="Logo"/>
+            </a></p>
+          {%- endif %}
+          {%- endblock %}
+          {%- if sidebars != None %}
+            {#- new style sidebar: explicitly include/exclude templates #}
+            {%- for sidebartemplate in sidebars %}
+            {%- include sidebartemplate %}
+            {%- endfor %}
+          {%- else %}
+            {#- old style sidebars: using blocks -- should be deprecated #}
+            {%- block sidebartoc %}
+            {%- include "localtoc.html" %}
+            {%- endblock %}
+            {%- block sidebarrel %}
+            {%- include "relations.html" %}
+            {%- endblock %}
+            {%- block sidebarsourcelink %}
+            {%- include "sourcelink.html" %}
+            {%- endblock %}
+            {%- if customsidebar %}
+            {%- include customsidebar %}
+            {%- endif %}
+            {%- block sidebarsearch %}
+            {%- include "searchbox.html" %}
+            {%- endblock %}
+          {%- endif %}
+        </div>
+      </div>
+      {%- endif %}
+{%- endmacro %}
+
+{%- macro script() %}
+    <script type="text/javascript">
+      var DOCUMENTATION_OPTIONS = {
+        URL_ROOT:    '{{ url_root }}',
+        VERSION:     '{{ release|e }}',
+        COLLAPSE_INDEX: false,
+        FILE_SUFFIX: '{{ '' if no_search_suffix else file_suffix }}',
+        HAS_SOURCE:  {{ has_source|lower }}
+      };
+    </script>
+    {%- for scriptfile in script_files %}
+    <script type="text/javascript" src="{{ pathto(scriptfile, 1) }}"></script>
+    {%- endfor %}
+{%- endmacro %}
+
+{%- macro css() %}
+    <link rel="stylesheet" href="{{ pathto('_static/' + style, 1) }}" type="text/css" />
+    <link rel="stylesheet" href="{{ pathto('_static/pygments.css', 1) }}" type="text/css" />
+    {%- for cssfile in css_files %}
+    <link rel="stylesheet" href="{{ pathto(cssfile, 1) }}" type="text/css" />
+    {%- endfor %}
+{%- endmacro %}
+
+<html xmlns="http://www.w3.org/1999/xhtml">
+  <head>
+    <meta http-equiv="Content-Type" content="text/html; charset={{ encoding }}" />
+    {{ metatags }}
+    {%- block htmltitle %}
+    <title>{{ title|striptags|e }}{{ titlesuffix }}</title>
+    {%- endblock %}
+    {{ css() }}
+    {%- if not embedded %}
+    {{ script() }}
+    {%- if use_opensearch %}
+    <link rel="search" type="application/opensearchdescription+xml"
+          title="{% trans docstitle=docstitle|e %}Search within {{ docstitle }}{% endtrans %}"
+          href="{{ pathto('_static/opensearch.xml', 1) }}"/>
+    {%- endif %}
+    {%- if favicon %}
+    <link rel="shortcut icon" href="{{ pathto('_static/' + favicon, 1) }}"/>
+    {%- endif %}
+    {%- endif %}
+{%- block linktags %}
+    {%- if hasdoc('about') %}
+    <link rel="author" title="{{ _('About these documents') }}" href="{{ pathto('about') }}" />
+    {%- endif %}
+    {%- if hasdoc('genindex') %}
+    <link rel="index" title="{{ _('Index') }}" href="{{ pathto('genindex') }}" />
+    {%- endif %}
+    {%- if hasdoc('search') %}
+    <link rel="search" title="{{ _('Search') }}" href="{{ pathto('search') }}" />
+    {%- endif %}
+    {%- if hasdoc('copyright') %}
+    <link rel="copyright" title="{{ _('Copyright') }}" href="{{ pathto('copyright') }}" />
+    {%- endif %}
+    <link rel="top" title="{{ docstitle|e }}" href="{{ pathto('index') }}" />
+    {%- if parents %}
+    <link rel="up" title="{{ parents[-1].title|striptags|e }}" href="{{ parents[-1].link|e }}" />
+    {%- endif %}
+    {%- if next %}
+    <link rel="next" title="{{ next.title|striptags|e }}" href="{{ next.link|e }}" />
+    {%- endif %}
+    {%- if prev %}
+    <link rel="prev" title="{{ prev.title|striptags|e }}" href="{{ prev.link|e }}" />
+    {%- endif %}
+{%- endblock %}
+{%- block extrahead %} {% endblock %}
+  </head>
+  <body>
+{%- block header %}{% endblock %}
+
+{%- block relbar1 %}{{ relbar() }}{% endblock %}
+
+{%- block content %}
+  {%- block sidebar1 %} {# possible location for sidebar #} {% endblock %}
+
+    <div class="document">
+  {%- block document %}
+      <div class="documentwrapper">
+      {%- if render_sidebar %}
+        <div class="bodywrapper">
+      {%- endif %}
+          <div class="body">
+            {% block body %} {% endblock %}
+          </div>
+      {%- if render_sidebar %}
+        </div>
+      {%- endif %}
+      </div>
+  {%- endblock %}
+
+  {%- block sidebar2 %}{{ sidebar() }}{% endblock %}
+      <div class="clearer"></div>
+    </div>
+{%- endblock %}
+
+{%- block relbar2 %}{{ relbar() }}{% endblock %}
+
+{%- block footer %}
+    <div class="footer">
+    {%- if show_copyright %}
+      {%- if hasdoc('copyright') %}
+        {% trans path=pathto('copyright'), copyright=copyright|e %}&copy; <a href="{{ path }}">Copyright</a> {{ copyright }}.{% endtrans %}
+      {%- else %}
+        {% trans copyright=copyright|e %}&copy; Copyright {{ copyright }}.{% endtrans %}
+      {%- endif %}
+    {%- endif %}
+    {%- if last_updated %}
+      {% trans last_updated=last_updated|e %}Last updated on {{ last_updated }}.{% endtrans %}
+    {%- endif %}
+    {%- if show_sphinx %}
+      {% trans sphinx_version=sphinx_version|e %}Created using <a href="http://sphinx-doc.org/">Sphinx</a> {{ sphinx_version }}.{% endtrans %}
+    {%- endif %}
+    </div>
+    <p>asdf asdf asdf asdf 22</p>
+{%- endblock %}
+  </body>
+</html>
+

doc/_themes/lammps_theme/search.html

@@ -0,0 +1,50 @@
+{#
+    basic/search.html
+    ~~~~~~~~~~~~~~~~~
+
+    Template for the search page.
+
+    :copyright: Copyright 2007-2013 by the Sphinx team, see AUTHORS.
+    :license: BSD, see LICENSE for details.
+#}
+{%- extends "layout.html" %}
+{% set title = _('Search') %}
+{% set script_files = script_files + ['_static/searchtools.js'] %}
+{% block footer %}
+  <script type="text/javascript">
+    jQuery(function() { Search.loadIndex("{{ pathto('searchindex.js', 1) }}"); });
+  </script>
+  {# this is used when loading the search index using $.ajax fails,
+     such as on Chrome for documents on localhost #}
+  <script type="text/javascript" id="searchindexloader"></script>
+  {{ super() }}
+{% endblock %}
+{% block body %}
+  <noscript>
+  <div id="fallback" class="admonition warning">
+    <p class="last">
+      {% trans %}Please activate JavaScript to enable the search
+      functionality.{% endtrans %}
+    </p>
+  </div>
+  </noscript>
+
+  {% if search_performed %}
+    <h2>{{ _('Search Results') }}</h2>
+    {% if not search_results %}
+      <p>{{ _('Your search did not match any documents. Please make sure that all words are spelled correctly and that you\'ve selected enough categories.') }}</p>
+    {% endif %}
+  {% endif %}
+  <div id="search-results">
+  {% if search_results %}
+    <ul>
+    {% for href, caption, context in search_results %}
+      <li>
+        <a href="{{ pathto(item.href) }}">{{ caption }}</a>
+        <p class="context">{{ context|e }}</p>
+      </li>
+    {% endfor %}
+    </ul>
+  {% endif %}
+  </div>
+{% endblock %}

doc/_themes/lammps_theme/searchbox.html

@@ -0,0 +1,9 @@
+{%- if builder != 'singlehtml' %}
+<div role="search">
+  <form id="rtd-search-form" class="wy-form" action="{{ pathto('search') }}" method="get">
+    <input type="text" name="q" placeholder="Search docs" />
+    <input type="hidden" name="check_keywords" value="yes" />
+    <input type="hidden" name="area" value="default" />
+  </form>
+</div>
+{%- endif %}

doc/_themes/lammps_theme/versions.html

@@ -0,0 +1,37 @@
+{% if READTHEDOCS %}
+{# Add rst-badge after rst-versions for small badge style. #}
+  <div class="rst-versions" data-toggle="rst-versions" role="note" aria-label="versions">
+    <span class="rst-current-version" data-toggle="rst-current-version">
+      <span class="fa fa-book"> Read the Docs</span>
+      v: {{ current_version }}
+      <span class="fa fa-caret-down"></span>
+    </span>
+    <div class="rst-other-versions">
+      <dl>
+        <dt>Versions</dt>
+        {% for slug, url in versions %}
+          <dd><a href="{{ url }}">{{ slug }}</a></dd>
+        {% endfor %}
+      </dl>
+      <dl>
+        <dt>Downloads</dt>
+        {% for type, url in downloads %}
+          <dd><a href="{{ url }}">{{ type }}</a></dd>
+        {% endfor %}
+      </dl>
+      <dl>
+        <dt>On Read the Docs</dt>
+          <dd>
+            <a href="//{{ PRODUCTION_DOMAIN }}/projects/{{ slug }}/?fromdocs={{ slug }}">Project Home</a>
+          </dd>
+          <dd>
+            <a href="//{{ PRODUCTION_DOMAIN }}/builds/{{ slug }}/?fromdocs={{ slug }}">Builds</a>
+          </dd>
+      </dl>
+      <hr/>
+      Free document hosting provided by <a href="http://www.readthedocs.org">Read the Docs</a>.
+
+    </div>
+  </div>
+{% endif %}
+

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=19.02
+export CFDEM_VERSION=19.09
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade

etc/cshrc

@@ -15,7 +15,7 @@
 #------------------------------------------------------------------------------
 
 setenv CFDEM_PROJECT CFDEM
-setenv CFDEM_VERSION 19.02
+setenv CFDEM_VERSION 19.09
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -34,12 +34,12 @@ Description
 #ifndef versionInfo_H
 #define versionInfo_H
 
-word CFDEMversion="PFM 19.02";
+word CFDEMversion="PFM 19.09";
-word compatibleLIGGGHTSversion="PFM 19.02";
+word compatibleLIGGGHTSversion="PFM 19.09";
 word OFversion="4.x";
 
 Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;
-Info << ", compatible to LIGGGHTS version: " << compatibleLIGGGHTSversion <<  endl;
+Info << "compatible to LIGGGHTS version: " << compatibleLIGGGHTSversion <<  endl;
-Info << ", compatible to OpenFOAM version: " << OFversion <<  endl;
+Info << "compatible to OpenFOAM version: " << OFversion <<  endl;
 
 #endif

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -45,6 +45,7 @@ Description
 #include "smoothingModel.H"
 #include "liggghtsCommandModel.H"
 #include "otherForceModel.H"
+#include "IOmanip.H"
 
 namespace Foam
 {
@@ -660,16 +661,14 @@ bool cfdemCloud::evolve
         }
 
         //============================================
-        //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD
+        //CHECK JUST TIME-INTERPOLATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD
         //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
         //      QUANTITIES AT THE GRID!
-        scalar timeStepFrac = dataExchangeM().timeStepFraction();
+        const scalar timeStepFrac = dataExchangeM().timeStepFraction();
+        int old_precision = Info().precision(10);
         Info << "\n timeStepFraction() = " << timeStepFrac << endl;
-        if(timeStepFrac > 1.0000001)
+        Info().precision(old_precision);
-        {
+
-   //         FatalError << "cfdemCloud::dataExchangeM().timeStepFraction()>1: Do not do this, since dangerous. This might be due to the fact that you used a adjustable CFD time step. Please use a fixed CFD time step." << abort(FatalError);
-              Warning << "cfdemCloud::dataExchangeM().timeStepFraction() = " << timeStepFrac << endl;
-        }
         clockM().start(24,"interpolateEulerFields");
 
         // update voidFractionField

src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C

@@ -72,24 +72,12 @@ diffusionCoefficient::diffusionCoefficient
     P_(sm.mesh().lookupObject<volScalarField>(pressureFieldName_)),
     partPressureName_(propsDict_.lookupOrDefault<word>("partPressureName","partP")),
     partPressure_(NULL),
-    densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
-    rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
-    molarConcFieldName_(propsDict_.lookupOrDefault<word>("totalMoleFieldName","molarConc")),
-    molarConc_(sm.mesh().lookupObject<volScalarField>(molarConcFieldName_)),
     X_(speciesNames_.size()),
     diffusantGasNames_(propsDict_.lookup("diffusantGasNames")),
     diffusionCoefficients_(diffusantGasNames_.size(),NULL),
+    Xdiffusant_(diffusantGasNames_.size()),
     initialized_(false)
 {
-    if(verbose_)
-    {
-        Info << " Reading diffusionCoefficient list: " << diffusantGasNames_ << endl;
-        for (int i = 0; i < diffusantGasNames_.size(); i++)
-        {
-            Info << "Diffusant names: " << diffusantGasNames_[i] << endl;
-        }
-    }
-
     particleCloud_.checkCG(false);
     allocateMyArrays();
     createCoeffs();
@@ -108,7 +96,8 @@ diffusionCoefficient::~diffusionCoefficient()
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
- void diffusionCoefficient::allocateMyArrays() const
+
+void diffusionCoefficient::allocateMyArrays() const
 {
     double initVal=0.0;
     if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
@@ -123,14 +112,11 @@ diffusionCoefficient::~diffusionCoefficient()
 
 void diffusionCoefficient::reAllocMyArrays() const
 {
-    if (particleCloud_.numberOfParticlesChanged())
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partPressure_,initVal,1,"nparticles");
+    for (int i=0; i<diffusantGasNames_.size(); i++)
     {
-        double initVal=0.0;
+        particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partPressure_,initVal,1,"nparticles");
-        for (int i=0; i<diffusantGasNames_.size(); i++)
-        {
-            particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1);
-        }
     }
 }
 
@@ -142,15 +128,22 @@ void diffusionCoefficient::init()
         volScalarField& X = const_cast<volScalarField&>
                 (mesh_.lookupObject<volScalarField>("X_"+speciesNames_[i]));
         X_.set(i, &X);
-
-         if(verbose_)
-         {
-             Info << " Read species list from: " << specDict_.name() << endl;
-             Info << " Reading species list: " << speciesNames_ << endl;
-             Info << " Looking up species fields: " << "X_"+speciesNames_[i] << endl;
-             Info << "The molar fraction fields (X_i): " << X_[i].name() << nl << endl;
-         }
     }
+
+    for (int j = 0; j < diffusantGasNames_.size(); j++)
+    {
+        volScalarField& Xdiffusant = const_cast<volScalarField&>
+                (mesh_.lookupObject<volScalarField>("X_"+diffusantGasNames_[j]));
+        Xdiffusant_.set(j, &Xdiffusant);
+
+        if (verbose_)
+        {
+            Info << " Reading diffusing gas species list: " << diffusantGasNames_ << endl;
+            Info << " Looking up diffusin gas species fields: " << "X_"+diffusantGasNames_[j] << endl;
+            Info << " The molar fraction fields (Xdiffusant_i): " << Xdiffusant_[j].name() << nl << endl;
+        }
+    }
+
     initialized_ = true;
 }
 
@@ -168,196 +161,115 @@ void diffusionCoefficient::execute()
 
     label  cellI=0;
     scalar Tfluid(0);
-    scalar rhofluid(0);
     scalar Pfluid(0);
-    scalar molarConcfluid(0);
     scalar Texp(0);
     scalar dBinary_(0);
     scalar Xnegative(0);
 
-    List<scalar> Xfluid_(0);
+    List<scalar> TotalFraction_(diffusantGasNames_.size(),Zero);
-    Xfluid_.setSize(speciesNames_.size());
-    List<scalar> XfluidDiffusant_(0);
-    XfluidDiffusant_.setSize(diffusantGasNames_.size());
-    List<scalar> MixtureBinaryDiffusion_;
-    MixtureBinaryDiffusion_.setSize(diffusantGasNames_.size());
-    List<scalar> TotalFraction_;
-    TotalFraction_.setSize(diffusantGasNames_.size());
 
-    // defining interpolators for T, rho, voidfraction, N
+    // defining interpolators for T and Pressure
     interpolationCellPoint <scalar> TInterpolator_(tempField_);
-    interpolationCellPoint <scalar> rhoInterpolator_(rho_);
     interpolationCellPoint <scalar> PInterpolator_(P_);
-    interpolationCellPoint <scalar> molarConcInterpolator_(molarConc_);
 
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
-        cellI=particleCloud_.cellIDs()[index][0];
+        cellI = particleCloud_.cellIDs()[index][0];
         if (cellI >=0)
         {
             if(interpolation_)
             {
-                vector position     =   particleCloud_.position(index);
+                vector position = particleCloud_.position(index);
-                Tfluid              =   TInterpolator_.interpolate(position,cellI);
+                Tfluid          = TInterpolator_.interpolate(position,cellI);
-                rhofluid            =   rhoInterpolator_.interpolate(position,cellI);
+                Pfluid          = PInterpolator_.interpolate(position,cellI);
-                Pfluid              =   PInterpolator_.interpolate(position,cellI);
-                molarConcfluid      =   molarConcInterpolator_.interpolate(position,cellI);
             }
             else
             {
-                Tfluid          =   tempField_[cellI];
+                Tfluid = tempField_[cellI];
-                rhofluid        =   rho_[cellI];
+                Pfluid = P_[cellI];
-                Pfluid          =   P_[cellI];
-                molarConcfluid  =   molarConc_[cellI];
 
-                for (int i = 0; i<speciesNames_.size();i++)
+                for (int i = 0; i<speciesNames_.size(); i++)
                 {
-                    Xfluid_[i] = X_[i][cellI];
+                    // total amount of negative molar fractions in the domain
-
+                    // check it and then delete it
-                    for (int j=0; j<diffusantGasNames_.size();j++)
+                    if (X_[i][cellI] < 0.0)
                     {
-                        if (diffusantGasNames_[j] == speciesNames_[i])
+                        Xnegative += X_[i][cellI] * mesh_.time().deltaTValue();
-                        {
+                        Info << "total negative molar fractions = " << Xnegative << endl;
-                            XfluidDiffusant_[j] = X_[i][cellI];
-                        }
-
-                        // total amount of negative molar fractions in the domain
-                        // check it and then delete it
-                        scalar timestep = mesh_.time().deltaTValue();
-                        if (Xfluid_[i] < 0.0)
-                        {
-                            Xnegative += Xfluid_[i]*timestep;
-                            Info << "total negative molar fractions =" << Xnegative << endl;
-                        }
-
-
-                        /*if (Xfluid_[i] <= 0.) Xfluid_[i] = 0.0;
-                        if (XfluidDiffusant_[j] <= 0.) XfluidDiffusant_[j] = 0.0;*/
-
-                        if(verbose_)
-                        {
-                            Info << "X fluid for species " << speciesNames_[i] << " : " << Xfluid_[i] << nl << endl;
-                            Info << "X fluid for diffusant species " << diffusantGasNames_[j] << " : " << XfluidDiffusant_[j] << nl << endl;
-                        }
                     }
                 }
             }
 
             partPressure_[index][0] = Pfluid;
             // change fluid pressure to 1 bar instead of Pa
-            Pfluid  =  Pfluid/100000.0;
+            Pfluid = Pfluid / 100000.0;
-            Texp  = Tfluid*sqrt(sqrt(Tfluid*Tfluid*Tfluid));
+            Texp   = Tfluid * sqrt(sqrt(Tfluid*Tfluid*Tfluid));
 
-            if(verbose_)
+            for (int j=0; j<diffusantGasNames_.size(); j++)
             {
-                Info << "partPressure_[index][0] = " << partPressure_[index][0] << endl;
+                TotalFraction_[j] = 0.0;
-                Info << "pressure field" << Pfluid << nl << endl;
+                dBinary_ = 0.0;
-                Info << "T - exponent calculation" << Texp << nl << endl;
-            }
 
-            for (int i=0; i<diffusantGasNames_.size();i++)
+                for (int i=0; i < speciesNames_.size(); i++)
-            {
-                MixtureBinaryDiffusion_[i]  =   0.0;
-                TotalFraction_[i]   =   0.0;
-                for (int j=0; j < speciesNames_.size();j++)
                 {
                     // get molecular diffusion coefficients if diffusant gas and reactant gas are not equal
-                    if (diffusantGasNames_[i] != speciesNames_[j])
+                    if (diffusantGasNames_[j] != speciesNames_[i])
                     {
-                        if(verbose_)
+                        if (verbose_)
                         {
-                            Info << "molar weights diffuser gases: " << molWeight(speciesNames_[j]) << nl << endl;
+                            Info << "molar weights diffuser gases: "  << molWeight(speciesNames_[i])      << nl << endl;
-                            Info << "molarConc fluid: " << molarConcfluid << nl << endl;
+                            Info << "molar weights diffusant gases: " << molWeight(diffusantGasNames_[j]) << nl << endl;
-                            Info << "rho fluid:  " << rhofluid << nl << endl;
-                            Info << "molar weights diffusant gases: " << molWeight(diffusantGasNames_[i]) << nl << endl;
                             Info << "Pressure: " << Pfluid << nl << endl;
-                            Info << "Temperature: " << Tfluid << nl << endl;
                         }
 
-                        if(coeffs.found(diffusantGasNames_[i]) && coeffs.found(speciesNames_[j]))
+                        if (coeffs.found(diffusantGasNames_[j]) && coeffs.found(speciesNames_[i]))
                         {
                             //  Fuller-Schettler-Giddings Equation
                             //  Unit of dBinary is [m^2/s]
                             //  INFO:: Normally unit of dBinary is cm^2/s, but the 1st term in RHS is 10^-3 instead
                             //  So here it is already converted
-                            dBinary_ = 1e-7*Texp*calcMolNum(i,j)/(Pfluid*calcDiffVol(i,j));
+                            dBinary_ = 1e-7 * Texp * calcMolNum(j,i) / (Pfluid * calcDiffVol(j,i));
 
-                            if(verbose_)
+                            if (verbose_)
                             {
-                                Info << "Molecular diffusion for species " << diffusantGasNames_[i] << " in "
+                                Info << "Xfluid for "     << speciesNames_[i]      << " : " << X_[i][cellI]          << nl << endl;
-                                     << speciesNames_[j] << " is : " << dBinary_ << nl << endl;
+                                Info << "Xdiffusant for " << diffusantGasNames_[j] << " : " << Xdiffusant_[j][cellI] << nl << endl;
-
-                                Info << "Molar fraction of species (speciesNames)" << speciesNames_[j] << " : " << Xfluid_[j] << nl << endl;
-                                Info << "Molar fraction of species (diffusantNames)" << diffusantGasNames_[i] << " : " << XfluidDiffusant_[i] << nl << endl;
                             }
 
-                            TotalFraction_[i] += Xfluid_[j]/dBinary_;
+                            TotalFraction_[j] += X_[i][cellI] / dBinary_;
-                            if (TotalFraction_[i] < VSMALL)
-                                MixtureBinaryDiffusion_[i] = VSMALL;
-                            else
-                                MixtureBinaryDiffusion_[i] = (1.0-XfluidDiffusant_[i])/TotalFraction_[i];
 
-                            if(verbose_)
+                            if (verbose_)
-                            {
+                                Info << "Total Fraction = " << TotalFraction_[j] << nl << endl;
-                                Info << "Total fraction calculated (ratio of stag. gas to binary diffusion : " << TotalFraction_[i] << nl << endl;
-                                Info << "Molar fraction of species diffusant gas " << diffusantGasNames_[i] << " : " << XfluidDiffusant_[i] << nl << endl;
-                                Info << "Multicomp. mix diffusion for species " << diffusantGasNames_[i]
-                                     << " is: " << MixtureBinaryDiffusion_[i] << nl << endl;
-                            }
-
-
-                            /*if (!(Xfluid_[j] <= 0.0))
-                            {
-                                // sum of all stagnant gases to sum of binary diffusion
-                                TotalFraction_[i]   +=  Xfluid_[j]/dBinary_;
-
-                                // dCoeff -- diffusion component of diffusant gas
-                                MixtureBinaryDiffusion_[i]  =   (1.0-XfluidDiffusant_[i])/TotalFraction_[i];
-
-                                if(verbose_)
-                                {
-                                    Info << "Total fraction calculated (ratio of stag. gas to binary diffusion : " << TotalFraction_[i] << nl << endl;
-                                    Info << "Molar fraction of species diffusant gas " << diffusantGasNames_[i] << " : " << XfluidDiffusant_[i] << nl << endl;
-                                }
-                            }
-
-                            if(verbose_)
-                            {
-                                Info << "Multicomp. mix diffusion for species " << diffusantGasNames_[i]
-                                     << " is: " << MixtureBinaryDiffusion_[i] << nl << endl;
-                            }*/
 
                             // pass on dCoeff values to array
-                            diffusionCoefficients_[i][index][0]= MixtureBinaryDiffusion_[i];
+                            if (TotalFraction_[j] < VSMALL)
-
+                                diffusionCoefficients_[j][index][0] = VSMALL;
-                        }else
+                            else
+                                diffusionCoefficients_[j][index][0] = (1.0 - Xdiffusant_[j][cellI]) / TotalFraction_[j];
+                        }
+                        else
                         {
-                        FatalError
+                            FatalError
                                 << "check tables for species diffusion volume"
                                 << endl
                                 << abort(FatalError);
                         }
                     }
                 }
-            }
-        }
 
-        if(verbose_)
+                if(verbose_)
-        {
+                    Info << "diffusionCoefficient of species " << diffusantGasNames_[j] << " = " << diffusionCoefficients_[j][index][0] << endl;
-            for(int i =0; i<diffusantGasNames_.size();i++)
-            {
-                Info << "diffusionCoefficient of species " << diffusantGasNames_[i] << " = " << diffusionCoefficients_[i][index][0] << endl;
             }
         }
     }
 
-    particleCloud_.dataExchangeM().giveData(partPressureName_, "scalar-atom",partPressure_);
+    particleCloud_.dataExchangeM().giveData(partPressureName_,"scalar-atom",partPressure_);
 
-    for (int i=0; i<diffusantGasNames_.size();i++)
+    for (int j=0; j<diffusantGasNames_.size(); j++)
     {
-        word pushName = diffusantGasNames_[i] + "_diffCoeff";
+        word pushName = diffusantGasNames_[j] + "_diffCoeff";
-        particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_[i]);
+        particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_[j]);
-    };
+    }
 
     Info << "give data done" << endl;
 }
@@ -385,10 +297,10 @@ void diffusionCoefficient::createCoeffs()
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-double diffusionCoefficient::calcDiffVol(int i, int j)
+double diffusionCoefficient::calcDiffVol(int j, int i)
 {
-    double sqrtvolDiff  =   coeffs(diffusantGasNames_[i])+coeffs(speciesNames_[j]);
+    double sqrtvolDiff = coeffs(diffusantGasNames_[j])+coeffs(speciesNames_[i]);
-    return sqrtvolDiff*sqrtvolDiff;
+    return sqrtvolDiff * sqrtvolDiff;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -408,15 +320,15 @@ void diffusionCoefficient::molWeightTable()
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // either calculate Molecular Weight addition (eq. D_ij) or consturct hashtable with diffusant and diffuser species
-double diffusionCoefficient::calcMolNum(int i, int j)
+double diffusionCoefficient::calcMolNum(int j, int i)
 {
     double molNum_ = 0.0;
 
-    double W1 = molWeight(diffusantGasNames_[i]);
+    double W1 = molWeight(diffusantGasNames_[j]);
-    double W2 = molWeight(speciesNames_[j]);
+    double W2 = molWeight(speciesNames_[i]);
 
-    molNum_   =   (W1 + W2) / (W1 * W2);
+    molNum_ = (W1 + W2) / (W1 * W2);
-    molNum_   =   sqrt(molNum_);
+    molNum_ = sqrt(molNum_);
 
     return molNum_;
 }

src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H

@@ -78,21 +78,14 @@ private:
 
     mutable double **partPressure_;
 
-    word densityFieldName_;
-
-    const volScalarField& rho_;
-
-    // total mole field
-    word molarConcFieldName_;
-
-    const volScalarField& molarConc_;
-
     UPtrList<volScalarField> X_;
 
     wordList diffusantGasNames_;
 
     mutable List<double**> diffusionCoefficients_;
 
+    UPtrList<volScalarField> Xdiffusant_;
+
     HashTable<scalar, word> coeffs;
 
     HashTable<scalar, word> molWeight;

src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C

@@ -69,6 +69,7 @@ massTransferCoeff::massTransferCoeff
     scaleDia_(1)
 {
     particleCloud_.checkCG(true);
+    allocateMyArrays();
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -95,12 +96,9 @@ void massTransferCoeff::allocateMyArrays() const
 
 void massTransferCoeff::reAllocMyArrays() const
 {
-    if (particleCloud_.numberOfParticlesChanged())
+    double initVal=0.0;
-    {
+    particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-        double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-    }
 }
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
@@ -188,10 +186,10 @@ void massTransferCoeff::execute()
             partRe_[index][0]   =   Rep;
         }
 
-        if (verbose_ && index >=0 && index < 2)
+        if (verbose_)
         {
-            Info << "Nufield = " << nuf << endl;
+            Pout << "Nufield = " << nuf << endl;
-            Info << "Rep = "     << Rep << endl;
+            Pout << "Rep = "     << Rep << endl;
         }
     }
 

src/lagrangian/cfdemParticle/subModels/chemistryModel/noChemistry/noChemistry.C

@@ -57,8 +57,64 @@ noChemistry::~noChemistry()
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
+tmp<volScalarField> noChemistry::Smi(label i) const
+{
+    tmp<volScalarField> Smi_
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "smi_",
+                particleCloud_.mesh().time().timeName(),
+                particleCloud_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+             ),
+             particleCloud_.mesh(),
+             dimensionedScalar
+             (
+              "zero",
+              dimMass/(dimVol*dimTime),
+              0.0
+             )
+        )
+    );
+
+    return Smi_;
+}
+
+tmp<volScalarField> noChemistry::Sm() const
+{
+    tmp<volScalarField> Sm_
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Sm_",
+                particleCloud_.mesh().time().timeName(),
+                particleCloud_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+             ),
+             particleCloud_.mesh(),
+             dimensionedScalar
+             (
+              "zero",
+              dimMass/(dimVol*dimTime),
+              0.0
+             )
+        )
+    );
+
+    return Sm_;
+}
+
 void noChemistry::execute()
-{}
+{
+    //do nothing
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/cfdemParticle/subModels/chemistryModel/noChemistry/noChemistry.H

@@ -66,6 +66,10 @@ public:
 
         void execute();
 
+        tmp <volScalarField> Smi(const label i) const;
+
+        tmp <volScalarField> Sm() const;
+
 };
 
 

src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C

@@ -56,6 +56,7 @@ reactantPerParticle::reactantPerParticle
     chemistryModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
     mesh_(sm.mesh()),
+    verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
     reactantPerParticle_(NULL),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_)),
@@ -98,11 +99,8 @@ void reactantPerParticle::allocateMyArrays() const
 
 void reactantPerParticle::reAllocMyArrays() const
 {
-    if (particleCloud_.numberOfParticlesChanged())
+    double initVal=0.0;
-    {
+    particleCloud_.dataExchangeM().allocateArray(reactantPerParticle_,initVal,1);
-        double initVal=0.0;
-        particleCloud_.dataExchangeM().allocateArray(reactantPerParticle_,initVal,1);
-    }
 }
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
@@ -145,6 +143,8 @@ void reactantPerParticle::execute()
             particlesPerCell=   particlesPerCell_[cellI];
             reactantPerParticle_[index][0] = voidfraction * cellvolume / particlesPerCell;
         }
+
+        if (verbose_) Info << "reactantPerParticle_" << reactantPerParticle_[index][0] << endl;
     }
 
         // give DEM data

src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.H

@@ -54,6 +54,8 @@ private:
 
     const fvMesh& mesh_;
 
+    bool verbose_;
+
     mutable double **reactantPerParticle_;
 
     word voidfractionFieldName_;

src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C

@@ -69,7 +69,7 @@ species::species
     // create a list from the Species table in the specified species dictionary
     speciesNames_(specDict_.lookup("species")),
     mod_spec_names_(speciesNames_.size()),
-    X_(speciesNames_.size()),                           //volumeScalarFields
+    X_(speciesNames_.size()),                           //volumeScalarFields of molarFractions
     molarFractions_(speciesNames_.size(),NULL),         //the value of molar fractions for every species
     changeOfSpeciesMass_(speciesNames_.size(),NULL),    //the values that are received from DEM with the name of Modified_+species name
     changeOfSpeciesMassFields_(speciesNames_.size()),   //the scalar fields generated with the values from Modified_+species names
@@ -125,6 +125,7 @@ species::~species()
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
 void species::allocateMyArrays() const
 {
     double initVal=0.0;
@@ -139,25 +140,21 @@ void species::allocateMyArrays() const
         {
             particleCloud_.dataExchangeM().allocateArray(molarFractions_[i],initVal,1,"nparticles");
             particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1,"nparticles");
-
         }
     }
 }
 
 void species::reAllocMyArrays() const
 {
-    if (particleCloud_.numberOfParticlesChanged())
+    double initVal=0.0;
-    {
+    particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1);
-        double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1);
+    particleCloud_.dataExchangeM().allocateArray(partMolarConc_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partMolarConc_,initVal,1);
 
-        for (int i=0; i<speciesNames_.size(); i++)
+    for (int i=0; i<speciesNames_.size(); i++)
-        {
+    {
-            particleCloud_.dataExchangeM().allocateArray(molarFractions_[i],initVal,1);
+        particleCloud_.dataExchangeM().allocateArray(molarFractions_[i],initVal,1);
-            particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1);
+        particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1);
-        }
     }
 }
 
@@ -169,7 +166,6 @@ void species::init()
         Info << " Reading species list: " << speciesNames_ << endl;
     }
 
-
     for (int i=0; i<speciesNames_.size(); i++)
     {
         // Define the Species volume scalar fields
@@ -182,13 +178,12 @@ void species::init()
 
         if(verbose_)
         {
-            Info << " Looking up species fields \n " << "X_"+speciesNames_[i] << endl;
             Info << "The molar fraction fields (X_i): \n" << X_[i].name() << endl;
             // Check if mod species are correct
             Info << "Modified species names are: \n" << mod_spec_names_[i] << endl;
         }
 
-        // Create new volScalarFields for the changed values of the species mass fields
+        // Create new volScalarFields for the changed values of the species mass fields -- gas species source term
         changeOfSpeciesMassFields_.set
         (
             i,
@@ -196,7 +191,7 @@ void species::init()
             (
                 IOobject
                 (
-                "ModSpeciesMassField_"+X_[i].name(),
+                "TotalChangeOfMassField_"+speciesNames_[i],
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::NO_READ,
@@ -205,7 +200,7 @@ void species::init()
                 mesh_,
                 dimensionedScalar("zero",dimMass/(dimVol*dimTime), 0.0)
             )
-         );
+        );
     }
     initialized_ = true;
 }
@@ -228,12 +223,10 @@ void species::execute()
     reAllocMyArrays();
 
     // get X_i, T, rho at particle positions
-    label  cellI=0;
+    label  cellI = 0;
     scalar Tfluid(0);
     scalar rhofluid(0);
-    List<scalar> Xfluid_;
     scalar voidfraction(1);
-    Xfluid_.setSize(speciesNames_.size());
     scalar molarConcfluid(0);
 
     // defining interpolators for T, rho, voidfraction, molarConc
@@ -250,53 +243,49 @@ void species::execute()
         {
             if(interpolation_)
             {
-                vector position     =   particleCloud_.position(index);
+                vector position = particleCloud_.position(index);
-                Tfluid              =   TInterpolator_.interpolate(position,cellI);
+                Tfluid          = TInterpolator_.interpolate(position,cellI);
-                rhofluid            =   rhoInterpolator_.interpolate(position,cellI);
+                rhofluid        = rhoInterpolator_.interpolate(position,cellI);
-                voidfraction        =   voidfractionInterpolator_.interpolate(position,cellI);
+                voidfraction    = voidfractionInterpolator_.interpolate(position,cellI);
-                molarConcfluid      =   molarConcInterpolator_.interpolate(position,cellI);
+                molarConcfluid  = molarConcInterpolator_.interpolate(position,cellI);
             }
             else
             {
-                Tfluid          =   tempField_[cellI];
+                Tfluid          = tempField_[cellI];
-                rhofluid        =   rho_[cellI];
+                rhofluid        = rho_[cellI];
-                voidfraction    =   voidfraction_[cellI];
+                voidfraction    = voidfraction_[cellI];
-                molarConcfluid	=   molarConc_[cellI];
+                molarConcfluid  = molarConc_[cellI];
-                for (int i = 0; i<speciesNames_.size();i++)
-                {
-                    Xfluid_[i] = X_[i][cellI];
-                   // if(Xfluid_[i] < 0.0) Xfluid_[i] = 0.0;
-                }
             }
 
-            partTemp_[index][0] =   Tfluid;
+            partTemp_[index][0] = Tfluid;
-            partRho_[index][0]  =   rhofluid;
+            partRho_[index][0]  = rhofluid;
-            partMolarConc_[index][0]    =   molarConcfluid;
+            partMolarConc_[index][0] = molarConcfluid;
 
             for (int i=0; i<speciesNames_.size();i++)
             {
                 // attention for indices when not communicating all species
-                molarFractions_[i][index][0]=Xfluid_[i];
+                molarFractions_[i][index][0] = X_[i][cellI];
             }
         }
     }
 
+
     if(verbose_)
     {
-        for(int i =0; i<speciesNames_.size();i++)
+        for(int i =0; i<speciesNames_.size(); i++)
         {
             Info << "X_i = " << X_[i].name() << endl;
             Info << "molarFractions_= " << molarFractions_[i][0][0] << endl;
             Info << "partRho_[index][0] = " << partRho_[0][0] << endl;
-            Info << "rhofluid =" << rhofluid << endl;
+            Info << "rhofluid = " << rhofluid << endl;
-            Info << "Xfluid = " << Xfluid_[i] << endl;
             Info << "partTemp_[index][0] = " << partTemp_[0][0] << endl;
-            Info << "Tfluid = " << Tfluid << endl  ;
+            Info << "Tfluid = " << Tfluid << endl;
-            Info << "voidfraction =" << voidfraction << endl;
+            Info << "voidfraction = " << voidfraction << endl;
         }
     }
 
-        // give DEM data
+    // give DEM data
+    {
         particleCloud_.dataExchangeM().giveData(partTempName_, "scalar-atom", partTemp_);
         particleCloud_.dataExchangeM().giveData(partRhoName_,  "scalar-atom", partRho_);
         particleCloud_.dataExchangeM().giveData(partMolarConcName_, "scalar-atom", partMolarConc_);
@@ -304,19 +293,25 @@ void species::execute()
         for (int i=0; i<speciesNames_.size();i++)
         {
             particleCloud_.dataExchangeM().giveData("X_"+speciesNames_[i],"scalar-atom",molarFractions_[i]);
-        };
+        }
 
-        Info << "give data done" << endl;
+        if (verbose_) Info << "give data done" << endl;
+    }
 
-        // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
+    // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
+    {
         scalar timestep = mesh_.time().deltaTValue();
         changeOfGasMassField_.primitiveFieldRef() = 0.0;
         changeOfGasMassField_.boundaryFieldRef() = 0.0;
         for (int i=0; i<speciesNames_.size();i++)
         {
-            particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_[i]);
             changeOfSpeciesMassFields_[i].primitiveFieldRef() = 0.0;
             changeOfSpeciesMassFields_[i].boundaryFieldRef() = 0.0;
+
+            particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_[i],particleCloud_.dataExchangeM().couplingInterval());
+
+            if (verbose_) Info << "changeOfSpeciesMass received from DEM = " << changeOfSpeciesMass_[i][0][0] << endl;
+
             particleCloud_.averagingM().setScalarSumCentre
             (
                 changeOfSpeciesMassFields_[i],
@@ -330,18 +325,22 @@ void species::execute()
             changeOfSpeciesMassFields_[i].primitiveFieldRef() /= (changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep);
             changeOfSpeciesMassFields_[i].correctBoundaryConditions();
             changeOfGasMassField_ += changeOfSpeciesMassFields_[i];
+
             if (verbose_)
             {
-                Info << "total conversion of species" << speciesNames_[i] << " = " << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep) << endl;
+                Info << "total conversion of species" << speciesNames_[i] << " = "
+                     << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep) << endl;
             }
         }
         massSourceCurr_ = gSum(changeOfGasMassField_*1.0*changeOfGasMassField_.mesh().V() * Nevery_ * timestep);
         massSourceTot_ += massSourceCurr_;
+
         if (verbose_)
         {
             Info << "total conversion of mass:\n\tcurrent source = " << massSourceCurr_ << "\n\ttotal source = " << massSourceTot_ << "\n" << endl;
             Info << "get data done" << endl;
         }
+    }
 }
 
 tmp<volScalarField> species::Smi (const label i) const

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C

@@ -202,7 +202,7 @@ bool dataExchangeModel::couple(int i)
 scalar dataExchangeModel::timeStepFraction() const
 {
     //return fraction between previous coupling TS and actual TS
-    return ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() - (couplingStep_-1) * couplingTime() ) / couplingTime();
+    return clamp( (particleCloud_.mesh().time().timeIndex() - timeIndexOffset_) * particleCloud_.mesh().time().deltaTValue() / couplingTime() - (couplingStep_ - 1) );
 }
 
 int dataExchangeModel::getNumberOfParticles() const
@@ -242,7 +242,8 @@ dataExchangeModel::dataExchangeModel
     maxNumberOfParticles_(0),
     couplingStep_(0),
     DEMts_(-1.),
-    couplingInterval_(readScalar(dict_.lookup("couplingInterval")))
+    couplingInterval_(readScalar(dict_.lookup("couplingInterval"))),
+    timeIndexOffset_(particleCloud_.mesh().time().timeIndex())
 {}
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H

@@ -68,8 +68,16 @@ protected:
 
         int couplingInterval_;
 
+        const label timeIndexOffset_;
+
     // Protected member functions
 
+        // note: C++17 has std::clamp
+        scalar clamp(scalar v, scalar lo=0.0, scalar hi=1.0) const
+        {
+            return ((v < lo) ? lo : (hi < v) ? hi : v);
+        }
+
 public:
 
     //- Runtime type information
@@ -185,38 +193,40 @@ public:
 
         inline int couplingInterval() const { return couplingInterval_; }
 
-        inline scalar couplingTime() const { return couplingInterval_*DEMts_; }
+        inline scalar couplingTime() const { return couplingInterval_ * DEMts_; }
 
         inline scalar TSstart() const { return particleCloud_.mesh().time().startTime().value()
                                         + (couplingStep_-1) * couplingTime(); }
 
-        inline scalar TSend() const { return particleCloud_.mesh().time().startTime().value()
+        inline scalar TSend() const   { return particleCloud_.mesh().time().startTime().value()
-                                        + (couplingStep_) * couplingTime(); }
+                                        + (couplingStep_  ) * couplingTime(); }
 
-        inline int DEMstepsTillT(scalar t) const { return (t - (particleCloud_.mesh().time().value() - couplingTime())  + SMALL) / DEMts_; }
+        inline int DEMstepsTillT(scalar t) const { return (t - (particleCloud_.mesh().time().value() - couplingTime()) + SMALL) / DEMts_; }
 
         inline void checkTSsize() const
         {
-            if (particleCloud_.mesh().time().deltaT().value() >  couplingInterval_ * DEMts_ + SMALL)
+            if (particleCloud_.mesh().time().deltaTValue() >  couplingInterval_ * DEMts_ + SMALL)
             {
-                Info << "particleCloud_.mesh().time().deltaT().value() = " << particleCloud_.mesh().time().deltaT().value() << endl;
+                Info << "particleCloud_.mesh().time().deltaTValue() = " << particleCloud_.mesh().time().deltaTValue() << endl;
                 Info << "couplingInterval_ = " << couplingInterval_ << endl;
                 Info << "DEMts_ = " << DEMts_ << endl;
                 FatalError << "\nError - CFD time-step bigger than coupling time (= DEM time step * coupling interval)!\n" << abort(FatalError);
             }
-            if (std::fabs((round(couplingTime()/particleCloud_.mesh().time().deltaT().value())*particleCloud_.mesh().time().deltaT().value())-couplingTime()) > SMALL)
+
+            if (std::fabs((round(couplingTime()/particleCloud_.mesh().time().deltaTValue()) * particleCloud_.mesh().time().deltaTValue()) - couplingTime()) > SMALL)
             {
-                Info << "particleCloud_.mesh().time().deltaT().value() = " << particleCloud_.mesh().time().deltaT().value() << endl;
+                Info << "particleCloud_.mesh().time().deltaTValue() = " << particleCloud_.mesh().time().deltaTValue() << endl;
                 Info << "couplingInterval_ = " << couplingInterval_ << endl;
                 Info << "DEMts_ = " << DEMts_ << endl;
-                Warning << "\nWarning - Coupling time (= DEM time step * coupling interval) is not a multiple of  CFD time-step!\n" << endl;
+                Warning << "\nWarning - Coupling time (= DEM time step * coupling interval) is not a multiple of CFD time-step!\n" << endl;
             }
+
             if (!particleCloud_.allowCFDsubTimestep())
-                if (particleCloud_.mesh().time().deltaT().value() < couplingInterval_ * DEMts_ + SMALL)
+                if (particleCloud_.mesh().time().deltaTValue() < couplingInterval_ * DEMts_ + SMALL)
                     FatalError << "\nYour models require: CFD time-step = coupling interval (= DEM time step * coupling interval)! \n" << abort(FatalError);
 
             // warn if sub-TS
-            if (particleCloud_.mesh().time().deltaT().value() < couplingTime() - SMALL)
+            if (particleCloud_.mesh().time().deltaTValue() < couplingTime() - SMALL)
                 Warning << "You are using sub-time-steps (i.e. CFD TS < coupling time)! Check your settings properly." << endl;
         }
 
@@ -228,15 +238,10 @@ public:
 
         inline bool doCoupleNow() const
         {
-            if (particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value()-SMALL
+            return ( (particleCloud_.mesh().time().timeIndex() - timeIndexOffset_) * particleCloud_.mesh().time().deltaTValue() - SMALL
-                > couplingStep_*DEMts_*couplingInterval_)
+                     >
-            {
+                     couplingStep_ * DEMts_ * couplingInterval_
-                return true;
+                   );
-            }
-            else
-            {
-                return false;
-            }
         }
 
         virtual int getNumberOfParticles() const;

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C

@@ -29,13 +29,17 @@ Description
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "twoWayMPI.H"
 #include "addToRunTimeSelectionTable.H"
+#include "OFstream.H"
 #include "clockModel.H"
-#include "pair.h"
+#include "liggghtsCommandModel.H"
-#include "force.h"
+
-#include "forceModel.H"
+//LAMMPS/LIGGGHTS
+#include <input.h>
+#include <library_cfd_coupling.h>
+#include <update.h>
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -66,7 +70,7 @@ twoWayMPI::twoWayMPI
     propsDict_(dict.subDict(typeName + "Props")),
     lmp(NULL)
 {
-    Info<<"Starting up LIGGGHTS for first time execution"<<endl;
+    Info << "Starting up LIGGGHTS for first time execution" << endl;
 
     MPI_Comm_dup(MPI_COMM_WORLD, &comm_liggghts);
 
@@ -74,8 +78,8 @@ twoWayMPI::twoWayMPI
     const fileName liggghtsPath(propsDict_.lookup("liggghtsPath"));
 
     // open LIGGGHTS input script
-    Info<<"Executing input script '"<< liggghtsPath.c_str() <<"'"<<endl;
+    Info << "Executing input script '" << liggghtsPath.c_str() << "'" << endl;
-    lmp = new LAMMPS_NS::LAMMPS(0,NULL,comm_liggghts);
+    lmp = new LAMMPS_NS::LAMMPS(0, NULL, comm_liggghts);
     lmp->input->file(liggghtsPath.c_str());
 
     // get DEM time step size

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H

@@ -41,27 +41,15 @@ SourceFiles
 #ifndef twoWayMPI_H
 #define twoWayMPI_H
 
+#include <mpi.h>
+
 #include "dataExchangeModel.H"
-#include "liggghtsCommandModel.H"
-#include "OFstream.H"
-#include <sys/stat.h>
-#include "pair.h"
-#include "force.h"
-#include "forceModel.H"
 
 //=================================//
 //LAMMPS/LIGGGHTS
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <mpi.h>
 #include <lammps.h>         // these are LAMMPS include files
-#include <input.h>
+#include <pair.h>
-#include <atom.h>
+#include <force.h>
-#include <library.h>
-#include <error.h>
-#include <library_cfd_coupling.h>
-#include <update.h>
 //=================================//
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C

@@ -33,13 +33,19 @@ Contributing authors
 
 \*---------------------------------------------------------------------------*/
 
-
 #include "twoWayOne2One.H"
 #include "addToRunTimeSelectionTable.H"
+#include "OFstream.H"
 #include "clockModel.H"
-#include "pair.h"
+#include "liggghtsCommandModel.H"
-#include "force.h"
+
-#include "forceModel.H"
+//LAMMPS/LIGGGHTS
+#include <atom.h>
+#include <input.h>
+#include <library.h>
+#include <library_cfd_coupling.h>
+#include <memory_ns.h>
+#include <update.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,7 +88,10 @@ twoWayOne2One::twoWayOne2One
     foam2lig_scl_tmp_(nullptr),
     staticProcMap_(propsDict_.lookupOrDefault<Switch>("useStaticProcMap", false)),
     cellIdComm_(propsDict_.lookupOrDefault<Switch>("useCellIdComm", false)),
+    prev_cell_ids_(nullptr),
+    dbl_cell_ids_(nullptr),
     my_prev_cell_ids_fix_(nullptr),
+    boundBoxMargin_(propsDict_.lookupOrDefault("boundingBoxMargin", 0.)),
     verbose_(propsDict_.lookupOrDefault("verbose", false)),
     lmp(nullptr)
 {
@@ -120,6 +129,22 @@ twoWayOne2One::twoWayOne2One
         tmpBoundaryFaces.localPoints()
     );
     thisFoamBox_ = treeBoundBox(boundaryFaces.localPoints());
+    if (boundBoxMargin_ > 0.)
+    {
+        thisFoamBox_.max() = thisFoamBox_.max() + boundBoxMargin_ * vector::one;
+        thisFoamBox_.min() = thisFoamBox_.min() - boundBoxMargin_ * vector::one;
+
+        Info<< "twoWayOne2One: Inflating bounding boxes by "
+            << boundBoxMargin_ << "."
+            << endl;
+    }
+    else if (boundBoxMargin_ < 0.)
+    {
+        FatalError
+        << "twoWayOne2One: Bound box margin cannot be smaller than 0."
+        << abort(FatalError);
+    }
+
     if (staticProcMap_)
     {
         createProcMap();
@@ -158,6 +183,11 @@ void twoWayOne2One::createProcMap()
         point(ligbb[0][0], ligbb[0][1], ligbb[0][2]),
         point(ligbb[1][0], ligbb[1][1], ligbb[1][2])
     );
+    if (boundBoxMargin_ > 0.)
+    {
+        thisLigBox.max() = thisLigBox.max() + boundBoxMargin_ * vector::one;
+        thisLigBox.min() = thisLigBox.min() - boundBoxMargin_ * vector::one;
+    }
     ligBoxes[Pstream::myProcNo()] = thisLigBox;
     Pstream::gatherList(ligBoxes);
     Pstream::scatterList(ligBoxes);
@@ -209,6 +239,9 @@ twoWayOne2One::~twoWayOne2One()
     destroy(foam2lig_vec_tmp_);
     destroy(foam2lig_scl_tmp_);
 
+    destroy(prev_cell_ids_);
+    destroy(dbl_cell_ids_);
+
     delete lmp;
 }
 
@@ -408,7 +441,7 @@ void twoWayOne2One::allocateArray
 ) const
 {
     int len = max(length,1);
-    lmp->memory->grow(array, len, width, "o2o:dbl**");
+    LAMMPS_MEMORY_NS::grow(array, len, width);
     for (int i = 0; i < len; i++)
         for (int j = 0; j < width; j++)
             array[i][j] = initVal;
@@ -423,7 +456,7 @@ void twoWayOne2One::allocateArray
 ) const
 {
     int len = max(particleCloud_.numberOfParticles(),1);
-    lmp->memory->grow(array, len, width, "o2o:dbl**:autolen");
+    LAMMPS_MEMORY_NS::grow(array, len, width);
     for (int i = 0; i < len; i++)
         for (int j = 0; j < width; j++)
             array[i][j] = initVal;
@@ -431,7 +464,7 @@ void twoWayOne2One::allocateArray
 
 void inline twoWayOne2One::destroy(double** array,int len) const
 {
-    lmp->memory->destroy(array);
+    LAMMPS_MEMORY_NS::destroy(array);
 }
 
 //============
@@ -445,7 +478,7 @@ void twoWayOne2One::allocateArray
 ) const
 {
     int len = max(length,1);
-    lmp->memory->grow(array, len, width, "o2o:int**");
+    LAMMPS_MEMORY_NS::grow(array, len, width);
     for (int i = 0; i < len; i++)
         for (int j = 0; j < width; j++)
             array[i][j] = initVal;
@@ -460,7 +493,7 @@ void twoWayOne2One::allocateArray
 ) const
 {
     int len = max(particleCloud_.numberOfParticles(),1);
-    lmp->memory->grow(array, len, width, "o2o:int**:autolen");
+    LAMMPS_MEMORY_NS::grow(array, len, width);
     for (int i = 0; i < len; i++)
         for (int j = 0; j < width; j++)
             array[i][j] = initVal;
@@ -468,7 +501,7 @@ void twoWayOne2One::allocateArray
 
 void inline twoWayOne2One::destroy(int** array,int len) const
 {
-    lmp->memory->destroy(array);
+    LAMMPS_MEMORY_NS::destroy(array);
 }
 
 //============
@@ -476,14 +509,14 @@ void inline twoWayOne2One::destroy(int** array,int len) const
 void twoWayOne2One::allocateArray(double*& array, double initVal, int length) const
 {
     int len = max(length,1);
-    lmp->memory->grow(array, len, "o2o:dbl*");
+    LAMMPS_MEMORY_NS::grow(array, len);
     for (int i = 0; i < len; i++)
         array[i] = initVal;
 }
 
 void inline twoWayOne2One::destroy(double* array) const
 {
-    lmp->memory->destroy(array);
+    LAMMPS_MEMORY_NS::destroy(array);
 }
 
 //==============
@@ -491,14 +524,14 @@ void inline twoWayOne2One::destroy(double* array) const
 void twoWayOne2One::allocateArray(int*& array, int initVal, int length) const
 {
     int len = max(length,1);
-    lmp->memory->grow(array, len, "o2o:int*");
+    LAMMPS_MEMORY_NS::grow(array, len);
     for (int i = 0; i < len; i++)
         array[i] = initVal;
 }
 
 void inline twoWayOne2One::destroy(int* array) const
 {
-    lmp->memory->destroy(array);
+    LAMMPS_MEMORY_NS::destroy(array);
 }
 //==============
 
@@ -646,6 +679,7 @@ bool twoWayOne2One::couple(int i)
         particleCloud_.clockM().stop("LIGGGHTS");
         Info<< "LIGGGHTS finished" << endl;
 
+        particleCloud_.clockM().start(4,"One2One_setup");
         if (!staticProcMap_)
         {
             createProcMap();
@@ -656,6 +690,7 @@ bool twoWayOne2One::couple(int i)
         locateParticles();
 
         setupFoam2LigCommunication();
+        particleCloud_.clockM().stop("One2One_setup");
 
         if (verbose_)
         {
@@ -730,12 +765,11 @@ void twoWayOne2One::locateParticles()
     destroy(my_flattened_positions);
 
     double* my_prev_cell_ids = nullptr;
-    double* prev_cell_ids = nullptr;
     if (cellIdComm_)
     {
         my_prev_cell_ids = my_prev_cell_ids_fix_->vector_atom;
-        allocateArray(prev_cell_ids, -1, lig2foam_->ncollected_);
+        allocateArray(prev_cell_ids_, -1, lig2foam_->ncollected_);
-        lig2foam_->exchange(my_prev_cell_ids, prev_cell_ids);
+        lig2foam_->exchange(my_prev_cell_ids, prev_cell_ids_);
     }
 
     if (lig2foam_mask_)
@@ -767,7 +801,7 @@ void twoWayOne2One::locateParticles()
             position,
             cellIdComm_
             ?   // don't know whether using round is efficient
-                (roundedCelli = round(prev_cell_ids[atomi])) < nCells
+                (roundedCelli = round(prev_cell_ids_[atomi])) < nCells
                 ?
                 roundedCelli
                 :
@@ -784,10 +818,6 @@ void twoWayOne2One::locateParticles()
             cellIds.append(cellI);
         }
     }
-    if (cellIdComm_)
-    {
-        destroy(prev_cell_ids);
-    }
 
     particleCloud_.setNumberOfParticles(n_located);
     particleCloud_.reAllocArrays();
@@ -817,7 +847,7 @@ void twoWayOne2One::locateParticles()
         3
     );
     particleCloud_.setPositions(getNumberOfParticles(), extracted_flattened_positions);
-    destroy(extracted_flattened_positions);
+    delete [] extracted_flattened_positions;
     destroy(collected_flattened_positions);
 
     particleCloud_.setCellIDs(cellIds);
@@ -871,14 +901,12 @@ void twoWayOne2One::setupFoam2LigCommunication()
 
     if (cellIdComm_)
     {
-        double** dbl_cell_ids = new double*[getNumberOfParticles()];
+        allocateArray(dbl_cell_ids_, -1, 1);
         for (int atomi = 0; atomi < getNumberOfParticles(); atomi++)
-        {   // TEMPORARY: if this persists after 19.07.2018, call me.
+        {
-            dbl_cell_ids[atomi] = new double[1];
+            dbl_cell_ids_[atomi][0] = particleCloud_.cellIDs()[atomi][0];
-            dbl_cell_ids[atomi][0] = particleCloud_.cellIDs()[atomi][0];
         }
-        giveData("prev_cell_ids", "scalar-atom", dbl_cell_ids, "double");
+        giveData("prev_cell_ids", "scalar-atom", dbl_cell_ids_, "double");
-        delete [] dbl_cell_ids;
     }
 }
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H

@@ -45,33 +45,19 @@ Contributing authors
 #ifndef twoWayOne2One_H
 #define twoWayOne2One_H
 
+#include <mpi.h>
+
 #include "dataExchangeModel.H"
-#include "liggghtsCommandModel.H"
-#include "OFstream.H"
-#include <sys/stat.h>
-#include "pair.h"
-#include "force.h"
-#include "forceModel.H"
 #include "one2one.H"
 
-#include "meshSearch.H"
+#include <treeBoundBox.H>
 
 //=================================//
 //LAMMPS/LIGGGHTS
-#include <stdio.h>
+#include <lammps.h>          // these are LAMMPS include files
-#include <stdlib.h>
+#include <pair.h>
-#include <string.h>
+#include <force.h>
-#include <mpi.h>
+#include <fix_property_atom.h> /// needed
-#include <lammps.h>         // these are LAMMPS include files
-#include <input.h>
-#include <atom.h>
-#include <memory.h>
-#include <library.h>
-#include <library_cfd_coupling.h>
-#include <update.h>
-#include <comm.h>
-#include <fix.h>
-#include <fix_property_atom.h>
 //=================================//
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -115,9 +101,13 @@ private:
 
     Switch    staticProcMap_;
     Switch    cellIdComm_;
+    double*   prev_cell_ids_;
+    double**  dbl_cell_ids_;
     LAMMPS_NS::FixPropertyAtom* my_prev_cell_ids_fix_;
 
-    treeBoundBox       thisFoamBox_;
+    const scalar      boundBoxMargin_;
+
+    treeBoundBox      thisFoamBox_;
 
     Switch            verbose_;
 

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C

@@ -242,6 +242,12 @@ void heatTransferGunn::calcEnergyContribution()
     // get DEM data
     particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
 
+    if(particleCloud_.cg() > 1.)
+    {
+        scaleDia_ = particleCloud_.cg();
+        Info << "Heat Transfer Gunn is using scale from liggghts cg = " << scaleDia_ << endl;
+    }
+
     if(calcPartTempField_)
     {
         partTempField_.primitiveFieldRef() = 0.0;
@@ -353,18 +359,24 @@ void heatTransferGunn::calcEnergyContribution()
                     partNu_[index][0] = Nup;
                 }
 
-                if(particleCloud_.verbose() && index >=0 && index <2)
+                if(verbose_ && index >=0 && index <2)
                 {
-                    Info << "partHeatFlux = " << partHeatFlux_[index][0] << endl;
+                    Pout << "partHeatFlux = " << partHeatFlux_[index][0] << endl;
-                    Info << "magUr = " << magUr << endl;
+                    Pout << "magUr = " << magUr << endl;
-                    Info << "As = " << As << endl;
+                    Pout << "kf0 = " << kf0_ << endl;
-                    Info << "muf = " << muf << endl;
+                    Pout << "Cp = " << Cp_ << endl;
-                    Info << "Rep = " << Rep << endl;
+                    Pout << "rho = " << rho_[cellI] << endl;
-                    Info << "Pr = " << Pr << endl;
+                    Pout << "h = " << h << endl;
-                    Info << "Nup = " << Nup << endl;
+                    Pout << "ds = " << ds << endl;
-                    Info << "voidfraction = " << voidfraction << endl;
+                    Pout << "ds_scaled = " << ds_scaled << endl;
-                    Info << "partTemp_[index][0] = " << partTemp_[index][0] << endl;
+                    Pout << "As = " << As << endl;
-                    Info << "Tfluid = " << Tfluid << endl  ;
+                    Pout << "muf = " << muf << endl;
+                    Pout << "Rep = " << Rep << endl;
+                    Pout << "Pr = " << Pr << endl;
+                    Pout << "Nup = " << Nup << endl;
+                    Pout << "voidfraction = " << voidfraction << endl;
+                    Pout << "partTemp_[index][0] = " << partTemp_[index][0] << endl;
+                    Pout << "Tfluid = " << Tfluid << endl  ;
                 }
             }
     }
@@ -444,6 +456,13 @@ void heatTransferGunn::calcEnergyContribution()
     }
 
     QPartFluid_.correctBoundaryConditions();
+
+    volScalarField minParticleWeights = particleCloud_.averagingM().UsWeightField();
+    Info << "Minimum Particle Weight " << gMin(minParticleWeights) << endl;
+    Info << "Minimum Particle Temperature: " << gMin(partTempField_) << endl;
+    Info << "Maximum Particle Temperature: " << gMax(partTempField_) << endl;
+    Info << "Minimum Fluid Temperature: " << gMin(tempField_) << endl;
+    Info << "Maximum Fluid Temperature: " << gMax(tempField_) << endl;
 }
 
 void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
@@ -451,11 +470,11 @@ void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
     Qsource += QPartFluid_;
 }
 
-void heatTransferGunn::addEnergyCoefficient(volScalarField& Qsource) const
+void heatTransferGunn::addEnergyCoefficient(volScalarField& Qcoeff) const
 {
     if(implicit_)
     {
-        Qsource += QPartFluidCoeff_;
+        Qcoeff += QPartFluidCoeff_;
     }
 }
 

src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C

@@ -46,6 +46,7 @@ reactionHeat::reactionHeat
     energyModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
+    verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
     mesh_(sm.mesh()),
     maxSource_(1e30),
     reactionHeatName_(propsDict_.lookupOrDefault<word>("reactionHeatName","reactionHeat")),
@@ -98,6 +99,14 @@ void reactionHeat::calcEnergyContribution()
 
     particleCloud_.dataExchangeM().getData(reactionHeatName_,"scalar-atom",reactionHeat_);
 
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+        if (verbose_ && index>=0 && index < 2)
+        {
+            Pout << "reactionHeat = " << reactionHeat_[index][0] << endl;
+        }
+    }
+
     reactionHeatField_.primitiveFieldRef() = 0.0;
     reactionHeatField_.boundaryFieldRef() = 0.0;
 

src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H

@@ -48,6 +48,8 @@ protected:
 
     bool interpolation_;
 
+    bool verbose_;
+
     const fvMesh& mesh_;
 
     scalar maxSource_;

src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C

@@ -148,6 +148,7 @@ void BeetstraDrag::setForce() const
     else if (particleCloud_.cg() > 1)
     {
         scaleDia_=particleCloud_.cg();
+        typeCG_[0] = scaleDia_;
         Info << "Beetstra using scale from liggghts cg = " << scaleDia_ << endl;
     }
 

src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H

@@ -17,6 +17,11 @@ Description
     drag law for monodisperse systems according to
     Beetstra et al. AIChE J 53.2 (2007)
 
+Reference
+    R. Beetstra, M.A. van der Hoef, J.A.M. Kuipers
+    "Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres"
+    AIChE Journal, 53 (2007)
+
 SourceFiles
     BeetstraDrag.C
 \*---------------------------------------------------------------------------*/
@@ -61,7 +66,7 @@ protected:
 
     mutable scalar scaleDia_;
 
-    scalarList typeCG_;
+    mutable scalarList typeCG_;
 
     mutable scalar scaleDrag_;
 
@@ -86,12 +91,12 @@ protected:
     virtual void adaptVoidfraction(double&, label) const {}
 
     virtual scalar effDiameter(double d, label cellI, label index) const {return d;}
-    
+
     virtual scalar meanSauterDiameter(double d, label cellI) const {return d;}
 
     double F(double, double) const;
 
-    double terminalVelocity(double, double, double, double, double, double) const; 
+    double terminalVelocity(double, double, double, double, double, double) const;
 
     double a(double) const;
 

src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H

@@ -34,6 +34,11 @@ Description
     - splits off explicit drag component due to fluctuation in fluid and particle
       velocity (optional via forceSubModel "ImExCorr")
 
+Reference
+    D. Gidaspow
+    "Multiphase flow and fluidization: Continuum and kinetic theory description",
+    1994
+
 Class
     GidaspowDrag
 

src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.H

@@ -82,6 +82,7 @@ protected:
 
         List<Switch> switchesList_;         // switches which are requested from dict
 
+        // a Switch is set to false in the default Switch constructor
         List<Switch> switches_;             // switch status
 
         mutable volScalarField nu_;

src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.H

@@ -28,6 +28,12 @@ Description
     This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 
+Reference
+    G. Wierink, C. Goniva, B. Niceno, K. Heiskanen,
+    "Mechanistic Modelling of Particle-Interface Interaction in Three-Phase Flows",
+    Proc. of the 8th Int. Conf. on CFD in Oil and Gas, Metallurgical and Process Industries,
+    Trondheim, 2011
+
 Class
     interface
 

src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H

@@ -28,6 +28,10 @@ Description
     This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 
+Reference
+    Crowe, Clayton T.; Sommerfeld, Martin; Tsuji, Yutaka (1998).
+    "Multiphase flows with droplets and particles"
+
 Class
     virtualMassForce
 
@@ -69,11 +73,10 @@ private:
 
     mutable double **UrelOld_;
 
-    const bool splitUrelCalculation_;    //indicator to split calculation of Urel between CFDEM and LIGGGHTS
+    const bool splitUrelCalculation_;      //indicator to split calculation of Urel between CFDEM and LIGGGHTS
                                            //requires the integration fix to take dv/dt into account!
 
-    mutable double Cadd_;                  //indicator to split calculation of Urel between CFDEM and LIGGGHTS
+    mutable double Cadd_;                  //virtual mass constant
-                                           //requires the integration fix to take dv/dt into account!
 
 public:
 

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.C

@@ -100,7 +100,7 @@ execute::execute
     parseCommandList(commandList_, labelList_, scalarList_, command_, propsDict_, timeStamp_);
     Info << "liggghtsCommand " << command_ << endl;
 
-    strCommand_ = string(command_);
+    strCommand_ = command_;
 
     checkTimeMode(propsDict_);
 

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/newLiggghtsCommandModel.C

@@ -50,7 +50,7 @@ autoPtr<liggghtsCommandModel> liggghtsCommandModel::New
 )
 {
     Info<< "Selecting liggghtsCommandModel "
-         << liggghtsCommandType << " ,provide dicts, numbered from 0 to n" << endl;
+         << liggghtsCommandType << ", provide dicts, numbered from 0 to n" << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(liggghtsCommandType);

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.C

@@ -94,7 +94,7 @@ writeLiggghts::writeLiggghts
     {
         overwrite_ = Switch(propsDict_.lookup("overwrite"));
     }
-    
+
     if (writeLast_)
     {
         runLast_ = true;
@@ -103,13 +103,13 @@ writeLiggghts::writeLiggghts
     {
         //Warning << "Using invalid options of writeLiggghts, please use 'writeLast' option." << endl;
         runEveryWriteStep_ = true;
-    }    
+    }
 
 
     command_ += " " + path_ + "/" + writeName_;
     if (overwrite_)
     {
-        strCommand_=string(command_);
+        strCommand_ = command_;
     }
     else
     {

tutorials/cfdemSolverMultiphase/damBreak/run.config

@@ -12,8 +12,7 @@
       "solver" : "cfdemSolverMultiphase",
       "type" : "CFDEMcoupling/mpi",
       "nprocs" : 4,
-      "pre_scripts" : ["prerun.sh"],
+      "pre_scripts" : ["prerun.sh"]
-      "post_scripts" : ["postrun.sh"]
     }
   ]
 }

tutorials/cfdemSolverMultiphase/particlesWaterEntry/run.config

@@ -12,8 +12,7 @@
       "solver" : "cfdemSolverMultiphase",
       "type" : "CFDEMcoupling/mpi",
       "nprocs" : 4,
-      "pre_scripts" : ["prerun.sh"],
+      "pre_scripts" : ["prerun.sh"]
-      "post_scripts" : ["postrun.sh"]
     }
   ]
 }

tutorials/cfdemSolverMultiphase/tankDrainage/run.config

@@ -12,8 +12,7 @@
       "solver" : "cfdemSolverMultiphase",
       "type" : "CFDEMcoupling/mpi",
       "nprocs" : 4,
-      "pre_scripts" : ["prerun.sh"],
+      "pre_scripts" : ["prerun.sh"]
-      "post_scripts" : ["postrun.sh"]
     }
   ]
 }

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/Allclean.sh

@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cd CFD
+cleanCase
+cd - 
+
+rm -f log*
+
+rm ./DEM/post/dump*
+
+# ----------------------------------------------------------------- end-of-file

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/Allrun.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+#===================================================================#
+# allrun script for testcase
+# M. Efe Kinaci - Sep 2018
+#===================================================================#
+
+#- source CFDEM env vars
+. ~/.bashrc
+
+#- include functions
+source $CFDEM_PROJECT_DIR/etc/functions.sh
+
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+export casePath
+
+cd $casePath/CFD/
+blockMesh
+
+if [ -f "$casePath/DEM/post/restart/liggghts.restart" ]; then
+    echo "LIGGGHTS init was run before - using existing restart file"
+else
+    #- run DEM in new terminal
+    $casePath/parDEMrun.sh
+fi
+
+echo "Run Simulation"
+bash $casePath/parCFDDEMrun.sh

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/CO

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      CO;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.3801847328;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+	type            fixedValue;
+	value		uniform	0.3801847328;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+/*        type            inletOutlet;  
+	value		$internalField;
+        inletValue	uniform 0.;*/
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/CO2

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      CO2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.4376348444;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+	type            fixedValue;
+	value		uniform	0.4376348444;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+        /*type            inletOutlet;  
+	value		$internalField;
+        inletValue	uniform 0.;*/
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/H2

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      H2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.0116562393;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+	type            fixedValue;
+	value		uniform 0.0116562393;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+/*        type            inletOutlet;  
+	value		$internalField;
+        inletValue	uniform 0.; */
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/H2O

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      H2O;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.0444583296;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+	type            fixedValue;
+	value		uniform 0.0444583296;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+/*        type            inletOutlet;  
+	value		$internalField;
+        inletValue	uniform 0.; */
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/Ksl

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      Ksl;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -3 -1 0 0 0 0];
+
+internalField   uniform 0.0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/N2

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      N2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.1260658538;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value		uniform 0.1260658538;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/Qsource

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      Qsource;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/T

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 1023.15;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value		uniform 1023.15;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/Tsource

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Tsource;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 1 0 0 0];
+
+// Tsource[K/s] = q[W/m3]/(rho[kg/m3]*cp[W*s/(kg*K)])
+// Q=10 W
+// V=0.01 m3
+// q=1000 W/m3
+// cp=1000
+// rho=1
+// --> Tsource = 1
+
+internalField   uniform 0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/U

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0.0 0.0 0.0);
+
+boundaryField
+{
+    side-walls
+    {
+        type            noSlip;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0.0 0.25 0.0);
+/*        type            interstitialInletVelocity;
+        inletVelocity      uniform (0 0.25 0);
+        alpha              voidfraction;
+        value              $internalField;  */
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/Us

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      Us;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+    
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/addSource

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location	"0";
+    object      addSource;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            zeroGradient;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/alphat

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/dSauter

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      dSauter;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform 0.00022;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/dSmoothing

@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      dSmoothing;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+	type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;  
+    }
+}
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/epsilon

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.003;
+
+boundaryField
+{
+    side-walls
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.003;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        /*type            inletOutlet;
+        inletValue      uniform 0.003;
+        value           uniform 0.003;*/
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/gradT

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      gradT;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 -1 0 1 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/k

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.001;
+
+boundaryField
+{
+    side-walls
+    {
+        type            kqRWallFunction;
+        value           uniform 0.001;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        /*type            inletOutlet;
+        inletValue      uniform 0.001;
+        value           uniform 0.001;*/
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/mut

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/nut

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    side-walls
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/p

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];	
+
+internalField   uniform 140000; 
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 140000;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/p.org

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 140000; 
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 140000;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/rho

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      rho;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -3 0 0 0 0 0];
+
+internalField   uniform 1.58;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/voidfraction

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      voidfraction;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    side-walls
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        //type            zeroGradient;
+        type            fixedValue;
+        value           uniform 1;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/wallQFactor

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      wallQFactor;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    side-walls
+    {
+        type            fixedValue;
+        value           uniform 1;
+    }
+   
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 1;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 1;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/chemistryProperties

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      chemistryProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+chemistry       off; //on;
+
+/*
+chemistryType
+{
+    chemistrySolver   noChemistrySolver;
+    chemistryThermo   rho;
+}
+
+initialChemicalTimeStep 1e-07;
+
+/*
+EulerImplicitCoeffs
+{
+    cTauChem        0.05;
+    equilibriumRateLimiter off;
+}
+
+odeCoeffs
+{
+    solver          seulex;
+    eps             0.05;
+} */
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/combustionProperties

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      combustionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+combustionModel   noCombustion<rhoThermoCombustion>; //PaSR<rhoChemistryCombustion>;  //PaSR<psiChemistryCombustion>;
+
+active  false; //true;
+
+noCombustionCoeffs
+{
+}
+
+
+PaSRCoeffs
+{
+    Cmix                1.0;
+    turbulentReaction   off;
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/couplingProperties

@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          couplingProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//===========================================================================//
+// sub-models & settings
+
+modelType "A"; // A or B
+
+couplingInterval 50;
+
+voidFractionModel divided;
+
+locateModel engine;	//turboEngineM2M;//
+
+meshMotionModel noMeshMotion;
+
+regionModel allRegion;
+
+IOModel basicIO;
+
+probeModel off;
+
+dataExchangeModel twoWayMPI;//twoWayM2M;//twoWayFiles;//oneWayVTK;//
+
+averagingModel dense;	
+
+clockModel off; //standardClock;//
+
+smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; //
+
+forceModels
+(
+    GidaspowDrag
+    gradPForce
+    viscForce
+);
+
+energyModels
+(
+    heatTransferGunn
+    reactionHeat
+);
+
+thermCondModel SyamlalThermCond;
+
+chemistryModels
+(
+    species
+    diffusionCoefficients
+    massTransferCoeff
+    //off
+);
+
+momCoupleModels
+(
+    implicitCouple
+);
+
+turbulenceModelType "turbulenceProperties";//"LESProperties";//
+
+//===========================================================================//
+// sub-model properties
+reactionHeatProps
+{
+    reactionHeatName            "reactionHeat";
+    verbose                     false;
+}
+
+heatTransferGunnProps
+{
+    partTempName "Temp";
+    partHeatFluxName "convectiveHeatFlux";
+    calcPartTempField true;
+    partRefTemp 1023;
+    implicit true;
+    verbose false;
+    interpolation false;
+}
+
+speciesProps
+{
+    ChemistryFile   "$FOAM_CASE/constant/foam.inp";
+    Nevery          1;
+    verbose         false;
+    interpolation   false;
+}
+
+diffusionCoefficientsProps
+{
+    verbose         false;
+    interpolation   false;
+    ChemistryFile   "$FOAM_CASE/constant/foam.inp";
+    diffusantGasNames       (    CO
+                                 H2
+                            ); 
+}
+
+massTransferCoeffProps
+{
+    verbose         false;
+    interpolation   false;
+}
+
+reactantPerParticleProps
+{
+}
+
+
+SyamlalThermCondProps
+{
+    voidfractionFieldName   "voidfraction";
+    rhoFieldName            "rho";
+}
+
+LaEuScalarTempProps
+{
+    velFieldName "U";
+    tempFieldName "T";
+    voidfractionFieldName "voidfraction";
+    partTempName "Temp";
+    partHeatFluxName "convectiveHeatFlux";
+    lambda 0.0256;
+    Cp 1007;
+}
+
+localPSizeDiffSmoothingProps
+{
+    lowerLimit 0.1;
+    upperLimit 1e10;
+    dSmoothingLength 1.5e-3;
+    Csmoothing    1.0; 
+}
+
+constDiffSmoothingProps
+{
+    lowerLimit 0.1;
+    upperLimit 1e10;
+    smoothingLength 1.5e-3;
+}
+
+implicitCoupleProps
+{
+    velFieldName "U";
+    granVelFieldName "Us";
+    voidfractionFieldName "voidfraction";
+}
+
+ArchimedesProps
+{
+    gravityFieldName "g";
+}
+
+gradPForceProps
+{
+    pFieldName "p";
+    voidfractionFieldName "voidfraction";
+    velocityFieldName "U";
+    interpolation true;
+}
+
+viscForceProps
+{
+    velocityFieldName "U";
+    interpolation true;
+}
+
+volWeightedAverageProps
+{
+    scalarFieldNames
+    (
+        voidfraction
+    );
+    vectorFieldNames
+    (
+    );
+    upperThreshold 0.999;
+    lowerThreshold 0;
+    verbose true;
+}
+
+totalMomentumExchangeProps
+{
+    implicitMomExFieldName "Ksl";
+    explicitMomExFieldName "none";
+    fluidVelFieldName "U";
+    granVelFieldName "Us";
+}
+
+GidaspowDragProps
+{
+    verbose true;
+    velFieldName "U";
+    granVelFieldName "Us";
+    voidfractionFieldName "voidfraction";
+    interpolation true;
+    phi 1;
+}
+
+BeetstraDragProps
+{
+    velFieldName "U";
+    granVelFieldName "Us";
+    gravityFieldName "g";
+    voidfractionFieldName "voidfraction";
+    interpolation false;
+    useFilteredDragModel ;
+    useParcelSizeDependentFilteredDrag ;
+    rhoP     4630.;
+    dPrim   0.00022;
+    rho     1.58;
+    nuf     1.126e-05;
+    g       9.81;
+    k       0.05;
+    aLimit  0.0;
+    verbose false;
+}
+
+DEMbasedDragProps
+{
+    velFieldName "U";
+    voidfractionFieldName "voidfraction";
+}
+
+DiFeliceDragProps
+{
+    //verbose true;
+    velFieldName "U";
+    voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
+    interpolation false;
+}
+
+KochHillDragProps
+{
+    verbose true;
+    velFieldName "U";
+    voidfractionFieldName "voidfraction";
+    interpolation false;
+    //forceSubModels
+    //(
+    //    ImExCorr
+    //);
+    implForceDEM true;
+    //implForceDEMaccumulated true;
+    //explicitCorr true;
+}
+
+RongDragProps
+{
+    verbose true;
+    velFieldName "U";
+    voidfractionFieldName "voidfraction";
+    interpolation false;
+    implForceDEM true;
+    implForceDEMaccumulated true;
+    granVelFieldName "Us";
+}
+
+virtualMassForceProps
+{
+    velFieldName "U";
+}
+
+particleCellVolumeProps
+{
+    upperThreshold 0.999;
+    lowerThreshold 0.;
+    verbose true;
+}
+
+fieldStoreProps
+{
+    scalarFieldNames
+    (
+    );
+
+    vectorFieldNames
+    (
+        "U"
+    );
+}
+
+oneWayVTKProps
+{
+    couplingFilename "vtk_out%4.4d.vtk";
+    maxNumberOfParticles 30000;
+}
+
+twoWayFilesProps
+{
+    maxNumberOfParticles 10100;
+}
+
+centreProps
+{
+    alphaMin 0.1;
+}
+
+engineProps
+{
+    treeSearch true;
+}
+
+turboEngineProps
+{
+    treeSearch true;
+}
+
+turboEngineM2MProps
+{
+    turboEngineProps
+    {
+        treeSearch true;
+    }
+}
+
+dividedProps
+{
+    alphaMin 0.01;
+    scaleUpVol 1.0;
+}
+
+twoWayMPIProps
+{
+    liggghtsPath "../DEM/in.liggghts_run";
+}
+
+twoWayM2MProps
+{
+    maxNumberOfParticles 10100;
+    liggghtsPath "../DEM/in.liggghts_run";
+}
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/foam.dat

@@ -0,0 +1,111 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      foam.dat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+H2
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       2.02;
+    }
+    thermodynamics
+    {
+        Cp              15016.8;
+        Hf              2.544e+06;
+    }
+    transport
+    {
+        mu              8.42e-05;
+        Pr              0.76;
+    }
+}
+H2O
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       18.015;
+    }
+    thermodynamics
+    {
+        Cp              2304.1;
+        Hf              2.544e+06;
+    }
+    transport
+    {
+        mu              1.7e-05;
+        Pr              0.99;
+    }
+}
+
+CO
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       28.01;
+    }
+    thermodynamics
+    {
+        Cp              1189.14;
+        Hf              2.544e+06;
+    }
+    transport
+    {
+        mu              1.66e-05;
+        Pr              0.792;
+    }
+}
+
+CO2
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       44.01;
+    }
+    thermodynamics
+    {
+        Cp              1239.98;
+        Hf              2.544e+06;
+    }
+    transport
+    {
+        mu              1.37e-05;
+        Pr              0.69;
+    }
+}
+
+N2
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       28.0134;
+    }
+    thermodynamics
+    {
+        Cp              1171.6;
+        Hf              2.544e+06;
+    }
+    transport
+    {
+        mu              1.66e-05;
+        Pr              0.69;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/foam.inp

@@ -0,0 +1,14 @@
+species
+(
+        CO
+        CO2
+        H2
+        H2O
+        N2
+);
+
+reactions
+{
+}
+
+

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/liggghtsCommands

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          liggghtsCommands;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+liggghtsCommandModels
+(
+    runLiggghts
+    writeLiggghts
+);
+// ************************************************************************* //
+
+/*runLiggghtsProps
+{
+    preNo false;
+}*/
+
+writeLiggghtsProps
+{
+    writeLast on;
+    writeName "post/restart/liggghts.restartCFDEM";
+    overwrite on;
+}

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/polyMesh/blockMeshDict

@@ -0,0 +1,144 @@
+// blockMesh :  Block mesh description file
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant/polyMesh";
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+   meshGenApp blockMesh;
+   convertToMeters 1.0;
+
+    //68 mm column diameter
+    //150 mm length
+
+
+    //Width of middle square section
+
+    //how many cells in the square section
+    //how many cells from square section to perimeter
+    // how many cells from top to bottom
+
+   vertices
+   (
+    ( 0.017 0.0  0.017) // Vertex fiveoclocksqb = 0
+    (-0.017 0.0  0.017) // Vertex sevenoclocksqb = 1
+    (-0.017 0.0 -0.017) // Vertex elevenoclocksqb = 2
+    ( 0.017 0.0 -0.017) // Vertex oneoclocksqb = 3
+
+    ( 0.0240416305819187 0.0  0.0240416305387665) // Vertex fiveoclockcb = 4
+    (-0.0240416305819187 0.0  0.0240416305387665) // Vertex sevenoclockcb = 5
+    (-0.0240416305819187 0.0 -0.0240416305387665) // Vertex elevenoclockcb = 6
+    ( 0.0240416305819187 0.0 -0.0240416305387665) // Vertex oneoclockcb = 7
+
+    ( 0.017 0.150  0.017) // Vertex fiveoclocksqt = 8
+    (-0.017 0.150  0.017) // Vertex sevenoclocksqt = 9
+    (-0.017 0.150 -0.017) // Vertex elevenoclocksqt = 10
+    ( 0.017 0.150 -0.017) // Vertex oneoclocksqt = 11
+
+    ( 0.0240416305819187 0.150  0.0240416305387665) // Vertex fiveoclockct = 12
+    (-0.0240416305819187 0.150  0.0240416305387665) // Vertex sevenoclockct = 13
+    (-0.0240416305819187 0.150 -0.0240416305387665) // Vertex elevenoclockct = 14
+    ( 0.0240416305819187 0.150 -0.0240416305387665) // Vertex oneoclockct = 15
+   );
+
+   blocks
+   (
+    //square block
+    hex (
+       1 0 3 2
+       9 8 11 10
+       )
+    (7 7 24)
+    simpleGrading (1 1 1)
+
+    //slice1
+    hex (
+       5 4 0 1
+       13 12 8 9
+       )
+    (7 3 24)
+    simpleGrading (1 1 1)
+
+    //slice2
+    hex (
+       1 2 6 5
+       9 10 14 13
+       )
+   (7 3 24)
+   simpleGrading (1 1 1)
+
+   //slice3
+   hex (
+         2 3 7 6
+         10 11 15 14
+       )
+   (7 3 24)
+   simpleGrading (1 1 1)
+
+   //slice4
+   hex (
+         3 0 4 7
+         11 8 12 15
+       )
+   (7 3  24)
+   simpleGrading (1 1 1)
+
+   );
+
+
+   //create the quarter circles
+   edges
+   (
+    arc 4 5 (0.0 0.0 0.034)
+    arc 5 6 (-0.034 0.0 0.0)
+    arc 6 7 (0.0 0.0 -0.034)
+    arc 7 4 (0.034 0.0 0.0)
+
+    arc 12 13 (0.0 0.150 0.034)
+    arc 13 14 (-0.034 0.150 0.0)
+    arc 14 15 (0.0 0.150 -0.034)
+    arc 15 12 (0.034 0.150 0.0)
+
+   );
+
+   patches
+   (
+    patch inlet
+    (
+    (0 3 2 1)
+     (0 4 7 3)
+     (4 0 1 5)
+     (1 2 6 5)
+     (3 7 6 2)
+    )
+
+    patch outlet
+    (
+     (8 11 10 9)
+     (8 12 15 11)
+     (12 8 9 13)
+     (9 10 14 13)
+     (11 15 14 10)
+    )
+
+    wall side-walls
+    (
+     (5 4 12 13)
+     (5 13 14 6)
+     (6 14 15 7)
+     (7 15 12 4)
+    )
+
+);
+

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/polyMesh/cylinderMesh.m4

@@ -0,0 +1,153 @@
+// blockMesh :  Block mesh description file
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant/polyMesh";
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+changecom(//)changequote([,])
+define(calc, [esyscmd(perl -e 'printf ($1)')])
+define(VCOUNT, 0)
+define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])
+
+
+   meshGenApp blockMesh;
+   convertToMeters 1.0;
+
+   define(D, 0.068) //68 mm column diameter
+   define(L, 0.150)  //150 mm length
+   define(PI, 3.14159265)
+
+   define(R, calc(D/2))
+   define(CW, calc(D/4)) //Width of middle square section
+
+   define(CX, calc(R*cos((PI/180)*45)))
+   define(CZ, calc(R*sin((PI/180)*45)))
+
+   define(NPS, 7) //how many cells in the square section
+   define(NPD, 3) //how many cells from square section to perimeter
+   define(NPY, 24) // how many cells from top to bottom
+
+   vertices
+   (
+    ( CW 0.0  CW) vlabel(fiveoclocksqb)
+    (-CW 0.0  CW) vlabel(sevenoclocksqb)
+    (-CW 0.0 -CW) vlabel(elevenoclocksqb)
+    ( CW 0.0 -CW) vlabel(oneoclocksqb)
+
+    ( CX 0.0  CZ) vlabel(fiveoclockcb)
+    (-CX 0.0  CZ) vlabel(sevenoclockcb)
+    (-CX 0.0 -CZ) vlabel(elevenoclockcb)
+    ( CX 0.0 -CZ) vlabel(oneoclockcb)
+
+    ( CW L  CW) vlabel(fiveoclocksqt)
+    (-CW L  CW) vlabel(sevenoclocksqt)
+    (-CW L -CW) vlabel(elevenoclocksqt)
+    ( CW L -CW) vlabel(oneoclocksqt)
+
+    ( CX L  CZ) vlabel(fiveoclockct)
+    (-CX L  CZ) vlabel(sevenoclockct)
+    (-CX L -CZ) vlabel(elevenoclockct)
+    ( CX L -CZ) vlabel(oneoclockct)
+   );
+
+   blocks
+   (
+    //square block
+    hex (
+       sevenoclocksqb fiveoclocksqb oneoclocksqb elevenoclocksqb
+       sevenoclocksqt fiveoclocksqt oneoclocksqt elevenoclocksqt
+       )
+    (NPS NPS NPY)
+    simpleGrading (1 1 1)
+
+    //slice1
+    hex (
+       sevenoclockcb fiveoclockcb fiveoclocksqb sevenoclocksqb
+       sevenoclockct fiveoclockct fiveoclocksqt sevenoclocksqt
+       )
+    (NPS NPD NPY)
+    simpleGrading (1 1 1)
+
+    //slice2
+    hex (
+       sevenoclocksqb elevenoclocksqb elevenoclockcb sevenoclockcb
+       sevenoclocksqt elevenoclocksqt elevenoclockct sevenoclockct
+       )
+   (NPS NPD NPY)
+   simpleGrading (1 1 1)
+
+   //slice3
+   hex (
+         elevenoclocksqb oneoclocksqb oneoclockcb elevenoclockcb
+         elevenoclocksqt oneoclocksqt oneoclockct elevenoclockct
+       )
+   (NPS NPD NPY)
+   simpleGrading (1 1 1)
+
+   //slice4
+   hex (
+         oneoclocksqb fiveoclocksqb fiveoclockcb oneoclockcb
+         oneoclocksqt fiveoclocksqt fiveoclockct oneoclockct
+       )
+   (NPS NPD  NPY)
+   simpleGrading (1 1 1)
+
+   );
+
+
+   //create the quarter circles
+   edges
+   (
+    arc fiveoclockcb sevenoclockcb (0.0 0.0 R)
+    arc sevenoclockcb elevenoclockcb (-R 0.0 0.0)
+    arc elevenoclockcb oneoclockcb (0.0 0.0 -R)
+    arc oneoclockcb fiveoclockcb (R 0.0 0.0)
+
+    arc fiveoclockct sevenoclockct (0.0 L R)
+    arc sevenoclockct elevenoclockct (-R L 0.0)
+    arc elevenoclockct oneoclockct (0.0 L -R)
+    arc oneoclockct fiveoclockct (R L 0.0)
+
+   );
+
+   patches
+   (
+    patch inlet
+    (
+    (fiveoclocksqb oneoclocksqb elevenoclocksqb sevenoclocksqb)
+     (fiveoclocksqb fiveoclockcb oneoclockcb oneoclocksqb)
+     (fiveoclockcb fiveoclocksqb sevenoclocksqb sevenoclockcb)
+     (sevenoclocksqb elevenoclocksqb elevenoclockcb sevenoclockcb)
+     (oneoclocksqb oneoclockcb elevenoclockcb elevenoclocksqb)
+    )
+
+    patch outlet
+    (
+     (fiveoclocksqt oneoclocksqt elevenoclocksqt sevenoclocksqt)
+     (fiveoclocksqt fiveoclockct oneoclockct oneoclocksqt)
+     (fiveoclockct fiveoclocksqt sevenoclocksqt sevenoclockct)
+     (sevenoclocksqt elevenoclocksqt elevenoclockct sevenoclockct)
+     (oneoclocksqt oneoclockct elevenoclockct elevenoclocksqt)
+    )
+
+    wall side-walls
+    (
+     (sevenoclockcb fiveoclockcb fiveoclockct sevenoclockct)
+     (sevenoclockcb sevenoclockct elevenoclockct elevenoclockcb)
+     (elevenoclockcb elevenoclockct oneoclockct oneoclockcb)
+     (oneoclockcb oneoclockct fiveoclockct fiveoclockcb)
+    )
+
+);
+

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/thermophysicalProperties

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo; //hePsiThermo;
+    mixture         reactingMixture;    //species and reactions are listed in chemistry file
+    transport       const; //sutherland;         //calcualtes viscosity as a function of temperature from Stuherland coefficient As and Ts
+    thermo          hConst; //janaf;
+    equationOfState perfectGas;
+    energy          sensibleInternalEnergy; //sensibleEnthalpy;
+    specie          specie;
+}
+
+chemistryReader foamChemistryReader;
+
+foamChemistryThermoFile     "$FOAM_CASE/constant/foam.dat";
+
+foamChemistryFile           "$FOAM_CASE/constant/foam.inp";
+
+inertSpecie     N2;
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/transportProperties

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu      nu [ 0 2 -1 0 0 0 0 ]  1.126e-05;
+
+kf      kf [ 1 1 -3 -1 0 0 0 ] 0.0507;          // comes from energy model
+
+Cp      Cp [ 0 2 -2 -1 0 0 0 ] 1118.15;         // comes from energy model
+
+
+// ******* Non-Newtonian transport properties ************************ //
+/*CrossPowerLawCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    m               m [ 0 0 1 0 0 0 0 ] 1;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+
+BirdCarreauCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    k               k [ 0 0 1 0 0 0 0 ] 0;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+*/
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/constant/turbulenceProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/system/controlDict

@@ -0,0 +1,266 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     cfdemSolverRhoPimpleChem;
+
+startFrom       startTime;
+
+startTime       0.0;
+
+stopAt          endTime;
+
+endTime         120.0;
+
+deltaT          0.00025;
+
+writeControl    runTime; //timeStep;
+
+writeInterval   0.25;    //5000;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+adjustTimeStep  no;
+
+maxCo           0.5;
+
+maxDeltaT       0.1;
+
+
+// ************************************************************************* //
+libs
+(
+    "libfieldFunctionObjects.so"
+);
+
+functions
+{
+    moleFrac
+    {
+        type    rhoReactionThermoMoleFractions;
+    }
+
+    probes
+    {
+        type probes;
+
+        functionObjectLibs ("libsampling.so");
+
+        // Name of the directory for probe data
+        name probes;
+
+        // Locations to be probed.
+        probeLocations
+        (
+            (0.03 0.025 0.0)
+            (0.03 0.050 0.0)
+            (0.03 0.075 0.0)
+            (0.03 0.090 0.0)
+            (0.03 0.120 0.0)
+        );
+
+        fields
+        (
+            rho
+            p
+            T
+            Cpv
+            N2
+            CO2
+            CO
+            H2
+            H2O
+            ModSpeciesMassField_N2
+            ModSpeciesMassField_CO2
+            ModSpeciesMassField_CO
+            ModSpeciesMassField_H2
+            ModSpeciesMassField_H2O
+            X_CO
+            X_CO2
+            X_H2
+            X_N2
+            X_H2O
+        );
+
+        writeControl   timeStep;
+        writeInterval  10000;
+    }
+
+
+    fieldMinMax_T
+    {
+        type            fieldMinMax;
+        libs            ("libfieldFunctionObjects.so");
+        fields          (T);
+        // Report the location of the field extrema
+        location        yes;
+        // Type of extrema for rank > 0 primitives
+        mode            magnitude; // magnitude | component
+    }
+
+    fieldMinMax_Cpv
+    {
+        type            fieldMinMax;
+        libs            ("libfieldFunctionObjects.so");
+        fields          (Cpv);	
+        // Report the location of the field extrema
+        location   yes;
+        // Type of extrema for rank > 0 primitives
+        mode            magnitude; // magnitude | component
+    }
+
+    fieldMinMax_partTemp
+    {
+        type       fieldMinMax;
+        libs       ("libfieldFunctionObjects.so");
+        fields     (partTemp);
+        location   yes;
+        mode       magnitude;
+    }
+
+    globalMassFrac
+    {
+        type            volRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl    timeStep;//outputTime;
+        writeInterval   1;
+        log             true;
+        writeFields     false;
+        regionType      all;
+        name            c0;
+        operation       weightedVolAverage;
+        weightField     rhoeps;
+        fields
+        (
+            H2
+            H2O
+            CO2
+            CO
+            N2
+        );
+    }
+
+    globalMass
+    {
+        type            volRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl   timeStep;//outputTime;
+        writeInterval   1;
+        log             true;
+        writeFields     false;
+        regionType      all;
+        name            c1;
+        operation       volIntegrate;
+        fields
+        (
+            rhoeps
+            rho
+        );
+    }
+
+    inflow_GasSpecies
+    {
+        type            surfaceRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl    timeStep;
+        writeInterval   100;
+        log             true;
+        // Output field values as well
+        writeFields     false;
+        regionType      patch;
+        name            inlet;
+        operation       sum;//areaIntegrate;
+
+        fields
+        (
+            CO
+            CO2
+            H2
+            H2O
+            N2
+        );
+    }
+
+    inflow_phi
+    {
+        type            surfaceRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl   timeStep;
+        writeInterval   100;
+        log             true;
+        // Output field values as well
+        writeFields     false;
+        regionType      patch;
+        name            inlet;
+        operation       sum;//areaIntegrate;
+
+        fields
+        (
+            phi
+        );
+    }
+
+    inflow_U
+    {
+        type            surfaceRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl    timeStep;
+        writeInterval   100;
+        log             true;
+        // Output field values as well
+        writeFields     false;
+        regionType      patch;
+        name            inlet;
+        operation       sum;//areaIntegrate;
+
+        fields
+        (
+            U
+        );
+    }
+
+    outflow_GasSpecies
+    {
+        $inflow_GasSpecies;
+        name    outlet;
+    }
+
+    outflow_phi
+    {
+        $inflow_phi;
+        name    outlet;
+    }
+
+    outflow_U
+    {
+        $inflow_U;
+        name    outlet;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/system/controlDict_test

@@ -0,0 +1,266 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     cfdemSolverRhoPimpleChem;
+
+startFrom       startTime;
+
+startTime       0.0;
+
+stopAt          endTime;
+
+endTime         0.2;
+
+deltaT          0.00025;
+
+writeControl    runTime; //timeStep;
+
+writeInterval   0.02;    //5000;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+adjustTimeStep  no;
+
+maxCo           0.5;
+
+maxDeltaT       0.1;
+
+
+// ************************************************************************* //
+libs
+(
+    "libfieldFunctionObjects.so"
+);
+
+functions
+{
+    moleFrac
+    {
+        type    rhoReactionThermoMoleFractions;
+    }
+
+    probes
+    {
+        type probes;
+
+        functionObjectLibs ("libsampling.so");
+
+        // Name of the directory for probe data
+        name probes;
+
+        // Locations to be probed.
+        probeLocations
+        (
+            (0.03 0.025 0.0)
+            (0.03 0.050 0.0)
+            (0.03 0.075 0.0)
+            (0.03 0.090 0.0)
+            (0.03 0.120 0.0)
+        );
+
+        fields
+        (
+            rho
+            p
+            T
+            Cpv
+            N2
+            CO2
+            CO
+            H2
+            H2O
+            ModSpeciesMassField_N2
+            ModSpeciesMassField_CO2
+            ModSpeciesMassField_CO
+            ModSpeciesMassField_H2
+            ModSpeciesMassField_H2O
+            X_CO
+            X_CO2
+            X_H2
+            X_N2
+            X_H2O
+        );
+
+        writeControl   timeStep;
+        writeInterval  50;
+    }
+
+
+    fieldMinMax_T
+    {
+        type            fieldMinMax;
+        libs            ("libfieldFunctionObjects.so");
+        fields          (T);
+        // Report the location of the field extrema
+        location        yes;
+        // Type of extrema for rank > 0 primitives
+        mode            magnitude; // magnitude | component
+    }
+
+    fieldMinMax_Cpv
+    {
+        type            fieldMinMax;
+        libs            ("libfieldFunctionObjects.so");
+        fields          (Cpv);	
+        // Report the location of the field extrema
+        location   yes;
+        // Type of extrema for rank > 0 primitives
+        mode            magnitude; // magnitude | component
+    }
+
+    fieldMinMax_partTemp
+    {
+        type       fieldMinMax;
+        libs       ("libfieldFunctionObjects.so");
+        fields     (partTemp);
+        location   yes;
+        mode       magnitude;
+    }
+
+    globalMassFrac
+    {
+        type            volRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl    timeStep;//outputTime;
+        writeInterval   1;
+        log             true;
+        writeFields     false;
+        regionType      all;
+        name            c0;
+        operation       weightedVolAverage;
+        weightField     rhoeps;
+        fields
+        (
+            H2
+            H2O
+            CO2
+            CO
+            N2
+        );
+    }
+
+    globalMass
+    {
+        type            volRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl   timeStep;//outputTime;
+        writeInterval   1;
+        log             true;
+        writeFields     false;
+        regionType      all;
+        name            c1;
+        operation       volIntegrate;
+        fields
+        (
+            rhoeps
+            rho
+        );
+    }
+
+    inflow_GasSpecies
+    {
+        type            surfaceRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl    timeStep;
+        writeInterval   100;
+        log             true;
+        // Output field values as well
+        writeFields     false;
+        regionType      patch;
+        name            inlet;
+        operation       sum;//areaIntegrate;
+
+        fields
+        (
+            CO
+            CO2
+            H2
+            H2O
+            N2
+        );
+    }
+
+    inflow_phi
+    {
+        type            surfaceRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl   timeStep;
+        writeInterval   100;
+        log             true;
+        // Output field values as well
+        writeFields     false;
+        regionType      patch;
+        name            inlet;
+        operation       sum;//areaIntegrate;
+
+        fields
+        (
+            phi
+        );
+    }
+
+    inflow_U
+    {
+        type            surfaceRegion;
+        libs ("libfieldFunctionObjects.so");
+        writeControl    timeStep;
+        writeInterval   100;
+        log             true;
+        // Output field values as well
+        writeFields     false;
+        regionType      patch;
+        name            inlet;
+        operation       sum;//areaIntegrate;
+
+        fields
+        (
+            U
+        );
+    }
+
+    outflow_GasSpecies
+    {
+        $inflow_GasSpecies;
+        name    outlet;
+    }
+
+    outflow_phi
+    {
+        $inflow_phi;
+        name    outlet;
+    }
+
+    outflow_U
+    {
+        $inflow_U;
+        name    outlet;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/system/decomposeParDict

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+//method          scotch;
+method simple;
+
+simpleCoeffs
+{
+    n               (2 1 2);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/system/fvOptions

@@ -0,0 +1,28 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvOptions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+limitedTemperature
+{
+    type        limitTemperature;
+    active      yes;
+    limitTemperatureCoeffs
+    {
+        active          yes;
+        selectionMode   all;
+        Tmin     1020;
+        Tmax     1500;
+    }
+} 

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/system/fvSchemes

@@ -0,0 +1,83 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         cellMDLimited leastSquares 0.5;
+    grad(U)         cellMDLimited leastSquares 0.5; //cellLimited Gauss linear 1; //
+    grad(h)         cellMDLimited leastSquares 0.5;
+    grad(e)         cellMDLimited leastSquares 0.5;
+}
+
+divSchemes
+{
+    default         Gauss linear;
+
+    div(phi,U)         Gauss limitedLinear 1; //Gauss linear; //Gauss limitedLinear 1; // //Gauss limitedLinearV 1; //
+    div(phid,p)        Gauss limitedLinear 1; //Gauss upwind; //Gauss limitedLinearV 1;
+    div(phi,K)         Gauss limitedLinear 1;
+    div(phi,h)         Gauss limitedLinear 1; //Gauss upwind;
+    div(phi,k)         Gauss limitedLinear 1; //Gauss upwind;
+    div(phi,epsilon)   Gauss limitedLinear 1; //Gauss upwind;
+    div(U)             Gauss limitedLinear 1;
+    div(phi,Yi_h)   Gauss multivariateSelection
+    {
+        CO limitedLinear01 1;
+        H2 limitedLinear01 1;
+        N2 limitedLinear01 1;
+        H2O limitedLinear01 1;
+        CO2 limitedLinear01 1;
+        h limitedLinear 1;
+        e limitedLinear 1;
+    }
+    // div((muEff*dev2(T(grad(U)))))        Gauss linear;
+    div((viscousTerm*dev2(grad(U).T())))    Gauss linear;
+    div((thermo:mu*dev(grad(U).T())))       Gauss linear;
+    div((nuEff*dev2(T(grad(U)))))           Gauss linear;
+    div(phi,T)      Gauss limitedLinear 1;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected; //Gauss linear orthogonal;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(U)  linear;
+}
+
+snGradSchemes
+{
+    default         corrected; //orthogonal;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/system/fvSolution

@@ -0,0 +1,132 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-6;
+        relTol          0.01;
+        /*solver          GAMG;
+        tolerance       1e-06;
+        relTol          0.01;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration on;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1; */
+    }
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    "(rho|G)"
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-6;
+        relTol          0.01;
+        /*solver      smoothSolver;
+        smoother    symGaussSeidel;
+        tolerance   1e-6;
+        relTol      0.1;
+        nSweeps     3;*/
+    }
+
+    "(rho|G)Final"
+    {
+        $rho;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    "(U|h|e|R|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+        maxIter         100;
+    }
+
+    "(U|h|e|R|k|epsilon)Final"
+    {
+        $U;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    "(Yi|CO|CO2|H2|H2O|N2)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0.01;
+        maxIter         100;
+    }
+
+    "(Yi|CO|CO2|H2|H2O|N2)Final"
+    {
+        $Yi;
+        tolerance       1e-06;
+        relTol          0;
+        maxIter         100;
+    }
+
+    T
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor yes;
+    nOuterCorrectors 5;
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    rhoMin  rhoMin [ 1 -3 0 0 0 ] 0.4;
+    rhoMax  rhoMax [ 1 -3 0 0 0 ] 2.0;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        T   0.5;
+        "(Yi|CO|CO2|H2|H2O|N2)" 0.6;
+        p   0.6;
+    }
+    equations
+    {
+        ".*"            0.8;
+        "(h|e).*"       0.5;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/system/probesDict

@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \    /   O peration     | Version:  2.1.x                                 |
+|   \  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+	version		2.0;
+	format		ascii;
+	class		dictionary;
+	object		probesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+
+fields
+(
+	rho
+	p
+	T
+	Cpv
+	N2
+	CO2
+	CO
+	H2
+	H2O
+	ModSpeciesMassField_N2
+	ModSpeciesMassField_CO2
+	ModSpeciesMassField_CO
+	ModSpeciesMassField_H2
+	ModSpeciesMassField_H2O
+	X_CO
+	X_CO2
+	X_H2
+	X_N2
+	X_H2O
+	X_N2
+	
+);
+
+writeControl   timeStep;
+writeInterval  10000;
+
+
+// Locations to be probed.
+probeLocations
+(
+        (0.03 0.025 0.0)
+        (0.03 0.050 0.0)
+        (0.03 0.075 0.0)
+        (0.03 0.090 0.0)
+        (0.03 0.120 0.0)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/DEM/in.liggghts_init

@@ -0,0 +1,85 @@
+# Particle insertion into domain
+coarsegraining  7.0 model_check error
+atom_style      granular
+atom_modify     map array
+
+communicate     single vel yes
+boundary        f f f
+newton          off
+
+units           si
+processors      * 1 2
+
+region          reg block -0.034 0.034 0.000 0.15 -0.034 0.034
+create_box      2 reg
+
+neighbor        1e-4 bin
+neigh_modify    delay 0
+
+# Material properties required for granular pair styles
+fix         m1 all property/global youngsModulus peratomtype 5.e6 5.e6
+fix         m2 all property/global poissonsRatio peratomtype 0.45 0.45
+fix         m3 all property/global coefficientRestitution peratomtypepair 2 0.3 0.3 0.3 0.3
+fix         m4 all property/global coefficientFriction peratomtypepair 2 0.5 0.5 0.5 0.5
+
+# pair style
+pair_style  gran model hertz tangential history
+pair_coeff  * *
+
+# timestep, gravity
+timestep    1e-6
+fix         gravi all gravity 9.81 vector 0.0 -1.0 0.0
+
+# walls
+fix     zwalls1  all wall/gran model hertz tangential history primitive type 1 yplane 0.000
+fix     zwalls2  all wall/gran model hertz tangential history primitive type 1 yplane 0.15
+fix     cylwalls all wall/gran model hertz tangential history primitive type 1 ycylinder 0.034 0. 0.
+
+# particle distributions and insertion 
+region  insreg1 cylinder y 0.0 0.0 0.034 0.00 0.15 units box
+
+#insert Quartz particles
+fix     pts5 all particletemplate/sphere 61 atom_type 2 density constant 2630 radius constant 0.00025
+fix     pdd2 all particledistribution/discrete 1 1 pts5 1.0
+
+#insert iron-ore range
+fix     pts1 all particletemplate/sphere 21 atom_type 1 density constant 4300 radius constant 0.000375
+fix     pts2 all particletemplate/sphere 31 atom_type 1 density constant 4300 radius constant 0.000188
+fix     pts3 all particletemplate/sphere 41 atom_type 1 density constant 4300 radius constant 0.000125
+fix     pts4 all particletemplate/sphere 51 atom_type 1 density constant 4300 radius constant 6.25E-05
+
+fix     pdd1 all particledistribution/discrete 1 4 pts1 0.228 pts2 0.286 pts3 0.332 pts4 0.154
+
+fix     ts_check all check/timestep/gran 10 0.1 0.1
+
+# apply nve integration to all particles that are inserted as single particles
+fix     integr all nve/sphere
+
+fix     ins2 all insert/pack seed 5331 distributiontemplate pdd2 &
+          maxattempt 1000 insert_every once overlapcheck yes all_in yes vel constant 0.0 0.0 0.0 &
+          region insreg1 mass_in_region 0.2
+
+run     150000
+
+fix     ins1 all insert/pack seed 5330 distributiontemplate pdd1 &
+          maxattempt 1000 insert_every once overlapcheck yes all_in yes vel constant 0.0 0.0 0.0 &
+          region insreg1 mass_in_region 0.365
+
+# make groups
+group ore type 1
+group quartz type 2
+
+# screen output
+compute         rke all erotate/sphere
+thermo_style    custom step atoms ke c_rke vol
+thermo          1000
+thermo_modify   lost ignore norm no
+compute_modify  thermo_temp dynamic yes
+
+# insert the first particles so that dump is not empty
+run     1
+dump    dmp all custom 5000 post/dump.liggghts_init id type x y z vx vy vz fx fy fz radius mass
+
+run     300000 upto
+
+write_restart   ../DEM/post/restart/liggghts.restart

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/DEM/in.liggghts_run

@@ -0,0 +1,283 @@
+# chemistry in fluidized bed
+log             ../DEM/log.liggghts
+thermo_log      ../DEM/post/thermo.txt
+
+coarsegraining  7.0 model_check error
+
+atom_style      granular
+atom_modify     map array
+communicate     single vel yes
+
+boundary        f f f
+newton          off
+
+units           si
+processors      2 * 2
+
+# read the restart file
+read_restart    ../DEM/post/restart/liggghts.restart
+
+neighbor        1e-4 bin
+neigh_modify    delay 0 
+
+# Material properties required for granular pair styles
+fix         m1 all property/global youngsModulus peratomtype 5.e6 5.e6
+fix         m2 all property/global poissonsRatio peratomtype 0.45 0.45
+fix         m3 all property/global coefficientRestitution peratomtypepair 2 0.3 0.3 0.3 0.3
+fix         m4 all property/global coefficientFriction peratomtypepair 2 0.5 0.5 0.5 0.5
+
+# pair style
+pair_style  gran model hertz tangential history # Hertzian without cohesion
+pair_coeff  * *
+
+# timestep, gravity
+timestep    0.000005
+fix         gravi all gravity 9.81 vector 0.0 -1.0 0.0
+
+# walls
+fix     zwalls1  all wall/gran model hertz tangential history primitive type 1 yplane 0.000
+fix     zwalls2  all wall/gran model hertz tangential history primitive type 1 yplane 0.15
+fix     cylwalls all wall/gran model hertz tangential history primitive type 1 ycylinder 0.034 0. 0.
+
+# Thermal properties
+fix     ftco all property/global thermalConductivity peratomtype 1.4 1.4
+fix     ftca all property/global thermalCapacity peratomtype 3000 5000
+
+# re-make groups
+group   ore type 1
+group   quartz type 2
+
+###############################################
+fix ts_check all check/timestep/gran 10 0.1 0.1
+
+###############
+# cfd coupling
+fix     cfd all couple/cfd couple_every 50 mpi
+fix     cfd2 all couple/cfd/force
+
+# apply nve integration to all particles that are inserted as single particles
+fix     integr all nve/sphere 
+
+#compute kinetic energy of particles to see if they are moving changing
+compute  KinEn all ke
+variable ke_tot	equal c_KinEn
+
+# this one invokes heat transfer calculation, transfers per-particle temperature and adds convective heat flux to particles
+fix     tconv all couple/cfd/convection T0 1022
+
+###############
+# this should invoke chemistry
+fix     cfd3 ore couple/cfd/chemistry n_species 5 species_names H2 H2O CO CO2 N2 n_diff 2 diffusant_names CO H2
+
+# Activate for 3-layer unreacted core shrink model
+fix     cfd5 ore chem/shrink/core speciesA CO molMassA 0.02801 speciesC CO2 molMassC 0.04401 scale_reduction_rate 10.0 screen no
+fix     cfd6 ore chem/shrink/core speciesA H2 molMassA 0.00202 speciesC H2O molMassC 0.01801 scale_reduction_rate 10.0 screen no
+
+# Material properties for chemical reaction
+fix     k0_CO ore property/atom k0_cfd5 vector yes no no 17 25 2700
+fix     Ea_CO ore property/atom Ea_cfd5 vector yes no no 69488 95000 130940
+
+fix     k0_H2 ore property/atom k0_cfd6 vector yes no no 30 23 160
+fix     Ea_H2 ore property/atom Ea_cfd6 vector yes no no 63627 85000 105908
+
+# particle parameters
+fix     porosity        ore property/atom   porosity_       vector yes no no 0.61 0.34 0.19 0.17
+fix     tortuosity      ore property/global tortuosity_     scalar 3
+fix     pore_diameter   ore property/global pore_diameter_  scalar 1.e-7
+
+
+# define layer properties
+fix     LayerRelRadii ore property/atom relRadii vector yes no no  1.0 0.998 0.995 0.98
+
+### define fix for mass layer - initial testing
+fix     LayerMasses ore property/atom massLayer vector yes no no 0. 0. 0. 0.
+fix     LayerDensities ore property/global density_ore vector 7870. 5740. 5170. 5240.
+
+## define fix for rho_eff and fracRed - initialize as zero
+fix     fracRed ore property/atom fracRed_ore vector yes no no 0. 0. 0.
+fix     rhoeff  ore property/atom rhoeff      vector yes no no 0. 0. 0. 0.
+
+###############
+# Write data into files for post-processing
+group part type 1
+
+variable np equal count(part)
+variable WI equal 100000
+variable time equal time
+
+compute fr_1 all reduce sum f_fracRed[1] f_fracRed[2] f_fracRed[3]
+fix     fr_d1 all ave/time 1 1 1 c_fr_1[1] c_fr_1[2] c_fr_1[3]
+variable fr_d1_1 equal f_fr_d1[1]
+variable fr_d1_2 equal f_fr_d1[2]
+variable fr_d1_3 equal f_fr_d1[3]
+
+# Calculate total fractional reduction to be written into the dump file for postProcessing/graphical representation
+variable frOV_d1 equal 1/9*(f_fr_d1[3])+2/9*(f_fr_d1[2])+6/9*(f_fr_d1[1])
+
+fix printfr1 all print ${WI} "${time} ${fr_d1_1} ${fr_d1_2} ${fr_d1_3} ${frOV_d1}" file ../DEM/post/fr_d1.dat title "#time fr_d1_1 fr_d1_2 fr_d1_3 fr_overall_d1"
+
+variable fr_1_divided equal f_fr_d1[1]/${np}
+variable fr_2_divided equal f_fr_d1[2]/${np}
+variable fr_3_divided equal f_fr_d1[3]/${np}
+variable frOV_divided equal 1/9*v_fr_3_divided+2/9*v_fr_2_divided+6/9*v_fr_1_divided
+
+fix printfrdividied all print ${WI} "${time} ${fr_1_divided} ${fr_2_divided} ${fr_3_divided} ${frOV_divided}" file ../DEM/post/frDivided.dat title "#time fr_1 fr_2 fr_3 frOV"
+
+### Reactant gas mass change ###
+compute dma_CO all reduce sum f_dmA_cfd5[1] f_dmA_cfd5[2] f_dmA_cfd5[3]
+fix     dmA_f_CO all ave/time 1 1 1 c_dma_CO[1] c_dma_CO[2] c_dma_CO[3]
+variable dmA_CO_1 equal f_dmA_f_CO[1]
+variable dmA_CO_2 equal f_dmA_f_CO[2]
+variable dmA_CO_3 equal f_dmA_f_CO[3]
+
+compute dma_H2 all reduce sum f_dmA_cfd6[1] f_dmA_cfd6[2] f_dmA_cfd6[3]
+fix     dmA_f_H2 all ave/time 1 1 1 c_dma_H2[1] c_dma_H2[2] c_dma_H2[3]
+variable dmA_H2_1 equal f_dmA_f_H2[1]
+variable dmA_H2_2 equal f_dmA_f_H2[2]
+variable dmA_H2_3 equal f_dmA_f_H2[3]
+
+### Write Mass change files to file ###
+fix printdmA all print ${WI} "${time} ${dmA_CO_1} ${dmA_CO_2} ${dmA_CO_3} ${dmA_H2_1} ${dmA_H2_2} ${dmA_H2_3}" file ../DEM/post/dmA.dat title "#time dmA_CO_1 dmA_CO_2 dmA_CO_3 dmA_H2_1 dmA_H2_2 dmA_H2_3"
+
+###############
+
+compute Aterm_red_CO all reduce sum f_Aterm_cfd5[1] f_Aterm_cfd5[2] f_Aterm_cfd5[3]
+fix     Aterm1 all ave/time 10 1 10 c_Aterm_red_CO[1] c_Aterm_red_CO[2] c_Aterm_red_CO[3]
+variable a_CO_1 equal f_Aterm1[1]
+variable a_CO_2 equal f_Aterm1[2]
+variable a_CO_3 equal f_Aterm1[3]
+
+compute Aterm_red_H2 all reduce sum f_Aterm_cfd6[1] f_Aterm_cfd6[2] f_Aterm_cfd6[3]
+fix     Aterm2 all ave/time 10 1 10 c_Aterm_red_H2[1] c_Aterm_red_H2[2] c_Aterm_red_H2[3]
+variable a_H2_1 equal f_Aterm2[1]
+variable a_H2_2 equal f_Aterm2[2]
+variable a_H2_3 equal f_Aterm2[3]
+
+compute Bterm_red_CO all reduce sum f_Bterm_cfd5[1] f_Bterm_cfd5[2] f_Bterm_cfd5[3]
+fix     Bterm1 all ave/time 10 1 10 c_Bterm_red_CO[1] c_Bterm_red_CO[2] c_Bterm_red_CO[3]
+variable b_CO_1 equal f_Bterm1[1]
+variable b_CO_2 equal f_Bterm1[2]
+variable b_CO_3 equal f_Bterm1[3]
+
+compute Bterm_red_H2 all reduce sum f_Bterm_cfd6[1] f_Bterm_cfd6[2] f_Bterm_cfd6[3]
+fix     Bterm2 all ave/time 10 1 10 c_Bterm_red_H2[1] c_Bterm_red_H2[2] c_Bterm_red_H2[3]
+variable b_H2_1 equal f_Bterm2[1]
+variable b_H2_2 equal f_Bterm2[2]
+variable b_H2_3 equal f_Bterm2[3]
+
+compute Massterm_red_CO all reduce sum f_Massterm_cfd5
+fix     Massterm1 all ave/time 10 1 10 c_Massterm_red_CO
+variable mt_CO equal f_Massterm1
+
+compute Massterm_red_H2 all reduce sum f_Massterm_cfd6
+fix     Massterm2 all ave/time 10 1 10 c_Massterm_red_H2
+variable mt_H2 equal f_Massterm2
+
+fix printAterm all print ${WI} "${time} ${a_CO_1} ${a_CO_2} ${a_CO_3} ${a_H2_1} ${a_H2_2} ${a_H2_3}" file ../DEM/post/Aterm.dat title "#time a_CO_1 aCO_2 a_CO_3 a_H2_1 a_H2_2 a_H2_3"
+
+fix printBterm all print ${WI} "${time} ${b_CO_1} ${b_CO_2} ${b_CO_3} ${b_H2_1} ${b_H2_2} ${b_H2_3}" file ../DEM/post/Bterm.dat title "#time b_CO_1 b_CO_2 b_CO_3 b_H2_1 b_H2_2 b_H2_3"
+
+fix printMterm all print ${WI} "${time} ${mt_CO} ${mt_H2}" file ../DEM/post/Mterm.dat title "#time mt_CO mt_H2"
+
+###############
+compute mass_layer all reduce sum f_LayerMasses[1] f_LayerMasses[2] f_LayerMasses[3] f_LayerMasses[4]
+fix     massLayerPrintout all ave/time 1 1 1 c_mass_layer[1] c_mass_layer[2] c_mass_layer[3] c_mass_layer[4]
+variable mL1 equal f_massLayerPrintout[1]
+variable mL2 equal f_massLayerPrintout[2]
+variable mL3 equal f_massLayerPrintout[3]
+variable mL4 equal f_massLayerPrintout[4]
+
+compute dY_CO all reduce sum f_dY_cfd5[1] f_dY_cfd5[2] f_dY_cfd5[3]
+fix dY_CO_Output all ave/time 1 1 1 c_dY_CO[1] c_dY_CO[2] c_dY_CO[3]
+variable dY_CO_1 equal f_dY_CO_Output[1]
+variable dY_CO_2 equal f_dY_CO_Output[2]
+variable dY_CO_3 equal f_dY_CO_Output[3]
+
+compute dY_H2 all reduce sum f_dY_cfd6[1] f_dY_cfd6[2] f_dY_cfd6[3]
+fix dY_H2_Output all ave/time 1 1 1 c_dY_H2[1] c_dY_H2[2] c_dY_H2[3]
+variable dY_H2_1 equal f_dY_H2_Output[1]
+variable dY_H2_2 equal f_dY_H2_Output[2]
+variable dY_H2_3 equal f_dY_H2_Output[3]
+
+compute xA_CO all reduce sum f_X_CO
+fix     molarFractionA_CO all ave/time 1 1 1 c_xA_CO
+variable xA_CO_1 equal f_molarFractionA_CO
+
+compute xA_H2 all reduce sum f_X_H2
+fix     molarFractionA_H2 all ave/time 1 1 1 c_xA_H2
+variable xA_H2_1 equal f_molarFractionA_H2
+
+compute xC_CO2 all reduce sum f_X_CO2
+fix     molarFractionC_CO2 all ave/time 1 1 1 c_xC_CO2
+variable xC_CO2_1 equal f_molarFractionC_CO2
+
+compute xC_H2O all reduce sum f_X_H2O
+fix     molarFractionC_H2O all ave/time 1 1 1 c_xC_H2O
+variable xC_H2O_1 equal f_molarFractionC_H2O
+
+fix printMassLayer all print ${WI} "${time} ${mL1} ${mL2} ${mL3} ${mL4}" file ../DEM/post/MassLayers.dat title "#time mL_Fe mL_w mL_m mL_h"
+
+fix printdmYLayer all print ${WI} "${time} ${dY_CO_1} ${dY_CO_2} ${dY_CO_3} ${dY_H2_1} ${dY_H2_2} ${dY_H2_3}" file ../DEM/post/dmY.dat title "#time dY_CO_1 dY_CO_2 dY_CO_3 dY_H2_1 dY_H2_2 dY_H2_3"
+
+fix molarFractions all print ${WI} "${time} ${xA_CO_1} ${xA_H2_1} ${xC_CO2_1} ${xC_H2O_1}" file ../DEM/post/molarFractions.dat title "#title x_CO x_H2 x_CO2 x_H2O"
+
+###############
+
+compute effDiffBinary_red_CO all reduce sum f_effDiffBinary_cfd5[1] f_effDiffBinary_cfd5[2] f_effDiffBinary_cfd5[3]
+fix     effDiffBinary1 all ave/time 10 1 10 c_effDiffBinary_red_CO[1] c_effDiffBinary_red_CO[2] c_effDiffBinary_red_CO[3]
+variable dij_CO_1 equal f_effDiffBinary1[1]
+variable dij_CO_2 equal f_effDiffBinary1[2]
+variable dij_CO_3 equal f_effDiffBinary1[3]
+
+compute effDiffBinary_red_H2 all reduce sum f_effDiffBinary_cfd6[1] f_effDiffBinary_cfd6[2] f_effDiffBinary_cfd6[3]
+fix     effDiffBinary2 all ave/time 10 1 10 c_effDiffBinary_red_H2[1] c_effDiffBinary_red_H2[2] c_effDiffBinary_red_H2[3]
+variable dij_H2_1 equal f_effDiffBinary2[1]
+variable dij_H2_2 equal f_effDiffBinary2[2]
+variable dij_H2_3 equal f_effDiffBinary2[3]
+
+compute effDiffKnud_red_CO all reduce sum f_effDiffKnud_cfd5[1] f_effDiffKnud_cfd5[2] f_effDiffKnud_cfd5[3]
+fix     effDiffKnud1 all ave/time 10 1 10 c_effDiffKnud_red_CO[1] c_effDiffKnud_red_CO[2] c_effDiffKnud_red_CO[3]
+variable dik_CO_1 equal f_effDiffKnud1[1]
+variable dik_CO_2 equal f_effDiffKnud1[2]
+variable dik_CO_3 equal f_effDiffKnud1[3]
+
+compute effDiffKnud_red_H2 all reduce sum f_effDiffKnud_cfd6[1] f_effDiffKnud_cfd6[2] f_effDiffKnud_cfd6[3]
+fix     effDiffKnud2 all ave/time 10 1 10 c_effDiffKnud_red_H2[1] c_effDiffKnud_red_H2[2] c_effDiffKnud_red_H2[3]
+variable dik_H2_1 equal f_effDiffKnud2[1]
+variable dik_H2_2 equal f_effDiffKnud2[2]
+variable dik_H2_3 equal f_effDiffKnud2[3]
+
+fix printDiffTerms all print ${WI} "${time} ${dij_CO_1} ${dij_CO_2} ${dij_CO_3} ${dij_H2_1} ${dij_H2_2} ${dij_H2_3} ${dik_CO_1} ${dik_CO_2} ${dik_CO_3} ${dik_H2_1} ${dik_H2_2} ${dik_H2_3} " file ../DEM/post/DiffTerm.dat title "#time dij_CO_1 dij_CO_2 dij_CO_3 dij_H2_1 dij_H2_2 dij_H2_3 dik_CO_1 dik_CO_2 dik_CO_3 dik_H2_1 dik_H2_2 dik_H2_3"
+
+###############
+
+## Check Layer Radii and Write them to files ###
+compute layerRad1 all reduce sum f_LayerRelRadii[1]
+fix     redRad1 all ave/time 1 1 1 c_layerRad1
+variable rr1 equal f_redRad1
+
+compute layerRad2 all reduce sum f_LayerRelRadii[2]
+fix     redRad2 all ave/time 1 1 1 c_layerRad2
+variable rr2 equal f_redRad2
+
+compute layerRad3 all reduce sum f_LayerRelRadii[3]
+fix     redRad3 all ave/time 1 1 1 c_layerRad3
+variable rr3 equal f_redRad3
+
+compute layerRad4 all reduce sum f_LayerRelRadii[4]
+fix     redRad4 all ave/time 1 1 1 c_layerRad4
+variable rr4 equal f_redRad4
+
+fix printRelRadii all print ${WI} "${time} ${rr1} ${rr2} ${rr3} ${rr4}" file ../DEM/post/relRadii.dat title "#time relRad_1 relRad_2 relRad_3 relRad_4"
+
+###############
+# screen output
+compute         rke all erotate/sphere
+thermo_style    custom step atoms ke c_rke vol
+thermo          10
+thermo_modify   lost ignore norm no
+compute_modify  thermo_temp dynamic yes
+
+dump        dmp all custom ${WI} ../DEM/post/dump*.liggghts_run id type x y z vx vy vz fx fy fz radius mass f_fracRed[1] f_fracRed[2] f_fracRed[3]
+
+run         1

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/parCFDDEMrun.sh

@@ -0,0 +1,83 @@
+#!/bin/bash
+
+#===================================================================#
+# allrun script for testcase
+# M. Efe Kinaci - Sep 2018
+#===================================================================#
+
+#- source CFDEM env vars
+. ~/.bashrc
+
+#- include functions
+source $CFDEM_PROJECT_DIR/etc/functions.sh
+
+#--------------------------------------------------------------------------------#
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+logpath=$casePath
+headerText="R2_FluidBed"
+logfileName="log_$headerText"
+solverName="cfdemSolverRhoPimpleChem"   #"cfdemSolverPiso"  #
+nrProcs="4"
+machineFileName="none"   # yourMachinefileName | none
+debugMode="off"          # on | off| strict
+testHarnessPath="$CFDEM_TEST_HARNESS_PATH"
+runOctave="false"
+postproc="false"
+
+#--------------------------------------------------------------------------------#
+
+#- call function to run a parallel CFD-DEM case
+parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode
+
+
+if [ $runOctave == "true" ]
+    then
+        #------------------------------#
+        #  octave
+
+        #- change path
+        cd octave
+
+        #- rmove old graph
+        rm cfdemSolverPiso_ErgunTestMPI.eps
+
+        #- run octave
+        octave totalPressureDrop.m
+
+        #- show plot 
+        evince cfdemSolverPiso_ErgunTestMPI.eps
+
+        #- copy log file to test harness
+        cp ../../$logfileName $testHarnessPath
+        cp cfdemSolverPiso_ErgunTestMPI.eps $testHarnessPath
+fi
+    
+if [ $postproc == "true" ]
+  then
+
+    #- keep terminal open (if started in new terminal)
+    echo "simulation finished? ...press enter to proceed"
+    read
+
+    #- get VTK data from liggghts dump file
+    cd $casePath/DEM/post
+    python -i $CFDEM_LPP_DIR/lpp.py dump*.liggghts_run
+
+    #- get VTK data from CFD sim
+    cd $casePath/CFD
+    reconstructPar
+    foamToVTK                                                   #- serial run of foamToVTK
+    #source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh                       #- include functions
+    #pseudoParallelRun "foamToVTK" $nrPostProcProcessors          #- pseudo parallel run of foamToVTK
+
+    #- start paraview
+    paraview
+ 
+ #- keep terminal open (if started in new terminal)
+    echo "...press enter to clean up case"
+    echo "press Ctr+C to keep data"
+    read
+
+fi
+

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/parDEMrun.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+#===================================================================#
+# DEMrun script for testcase
+# M. Efe Kinaci - Sep 2018
+#===================================================================#
+
+#- source CFDEM env vars
+. ~/.bashrc
+
+#- include functions
+source $CFDEM_PROJECT_DIR/etc/functions.sh
+
+echo "starting DEM run in parallel..."
+#--------------------------------------------------------------------------------#
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+logpath="$casePath"
+headerText="run_liggghts_init_DEM"
+logfileName="log_$headerText"
+solverName="in.liggghts_init"
+nrProcs=4
+machineFileName="none"
+debugMode="off"
+#--------------------------------------------------------------------------------#
+
+#- call function to run DEM case
+parDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode
+

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/postrun.sh

@@ -0,0 +1,19 @@
+#!/bin/bash
+#- source CFDEM env vars
+. ~/.bashrc
+
+#- include functions
+source $CFDEM_PROJECT_DIR/etc/functions.sh
+
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+#- get VTK data from liggghts dump file
+cd $casePath/DEM/post
+python -i $CFDEM_LPP_DIR/lpp.py dump*.liggghts_run
+
+#- get VTK data from CFD sim
+cd $casePath/CFD
+reconstructPar
+foamToVTK
+
+

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/prerun.sh

@@ -0,0 +1,16 @@
+#!/bin/bash
+
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+# check if mesh was built
+if [ -f "$casePath/CFD/constant/polyMesh/boundary" ]; then
+    echo "mesh was built before - using old mesh"
+else
+    echo "mesh needs to be built"
+    cd $casePath/CFD
+    blockMesh
+fi
+
+# adapt settings for test harness run
+cp $casePath/CFD/system/controlDict_test $casePath/CFD/system/controlDict

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/run.config

@@ -0,0 +1,19 @@
+{
+  "type" : "CFDEMcoupling",
+  "runs" : [
+    {
+      "name" : "liggghts-init",
+      "input_script" : "DEM/in.liggghts_init",
+      "type" : "liggghts/mpi",
+      "nprocs" : 4
+    },
+    {
+      "name" : "cfdemrun",
+      "depends_on" : "liggghts-init",
+      "solver" : "cfdemSolverRhoPimpleChem",
+      "type" : "CFDEMcoupling/mpi",
+      "nprocs" : 4,
+      "pre_scripts" : ["prerun.sh"]
+    }
+  ]
+}