Add option to evolve temperature eqn

Docs/sphinx_documentation/source/FluidEquations.rst

@@ -12,12 +12,16 @@ Fluid Variables
    +-----------------------+--------------------------------------------------+
    | :math:`\mu_s`         | scalar diffusivity                               |
    +-----------------------+--------------------------------------------------+
+   | :math:`\mu_T`         | thermal conductivity                             |
+   +-----------------------+--------------------------------------------------+
    | :math:`{\bf g}`       | Gravitational acceleration                       |
    +-----------------------+--------------------------------------------------+
    | :math:`{\bf H}_U`     | :math:`= (H_x , H_y , H_z )`, External Forces    |
    +-----------------------+--------------------------------------------------+
    | :math:`H_s`           | External sources                                 |
    +-----------------------+--------------------------------------------------+
+   | :math:`H_T`           | External heat sources                            |
+   +-----------------------+--------------------------------------------------+
 
 Fluid Equations
 ===============
@@ -43,5 +47,9 @@ or, for non-conservative,
 
 .. math:: \frac{\partial s}{\partial t} + U \cdot \nabla s  =  \nabla \cdot \mu_s \nabla s + H_s
 
-By default, :math:`H_s = 0` and :math:`{\bf H}_U = {\bf 0}`.
+Optional temperature:
+
+.. math:: \rho c_p \left( \frac{\partial T}{\partial t} + U \cdot \nabla T \right)  =  \nabla \cdot \mu_T \nabla T + H_T
+
+By default, :math:`H_s = 0`,  :math:`H_T = 0`, and :math:`{\bf H}_U = {\bf 0}`.
 If gravity is set during runtime, then :math:`{\bf H}_U` defaults to :math:`\rho {\bf g}`.

Docs/sphinx_documentation/source/InputsAlgorithm.rst

@@ -18,6 +18,8 @@ The following inputs must be preceded by "incflo."
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  ntrac               |  number of tracers                                                    |  int        |  1           |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
+|  use_temperature     |  include temperature equation?                                        |  bool       |  false       |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  constant_density    |  Only evolve the continuity equation if false                         |  bool       |  true        |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  rho_0               |  density (if constant)                                                |  Real       |  1.0         |
@@ -28,6 +30,8 @@ The following inputs must be preceded by "incflo."
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  mu_s                |  scalar diffusivity                                                   |  Real(s)    |  0.0         |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
+|  mu_T                |  thermal conductivity                                                 |  Real(s)    |  0.0         |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  use_cc_proj         |  Use cell-centered rather than nodal pressure; this changes the       |  bool       | false        |
 |                      |  which approximate projection we use                                  |             |              |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+

Docs/sphinx_documentation/source/InputsProblemDefinition.rst

@@ -50,15 +50,18 @@ Setting basic boundary conditions can be specified by inputs preceded by "xlo",
 |                    | * 'po'  or 'pressure_outflow'                                             |             |           |
 |                    | * 'mi'  or 'mass_inflow'                                                  |             |           |
 |                    | * 'nsw' or 'no_slip_wall'                                                 |             |           |
+|                    | * 'sw' or 'slip_wall'                                                     |             |           |
 |                    | * 'mixed'                                                                 |             |           |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
 | pressure           | Sets boundary pressure for pressure inflows, outflows and mass inflows    |    Real     |  None     |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
-| velocity           | Sets boundary velocity for mass inflows                                   |    Reals    | (0, 0, 0) |
+| velocity           | Sets Dirichlet boundary condition for velocity                            |    Reals    | (0, 0, 0) |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
-| density            | Sets boundary density for mass inflows                                    |    Real     |  1.0      |
+| density            | Sets Dirichlet boundary condition for density                             |    Real     |  1.0      |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
-| tracer             | Sets boundary tracer for mass inflows                                     |    Real     |  0.0      |
+| tracer             | Sets Dirichlet boundary condition for tracer                              |    Real     |  0.0      |
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
+| temperature        | Sets Dirichlet boundary condition for temperature                         |    Real     |  1.0      |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
 
    The 'mixed' boundary type allows for inflow and outflow on the same domain face.

src/CMakeLists.txt

@@ -16,6 +16,8 @@ target_sources(incflo
    incflo_update_density.cpp
    incflo_update_tracer.cpp
    incflo_update_velocity.cpp
+   incflo_update_temperature.cpp
+   incflo_specific_heat.cpp
    incflo_utils.cpp
    main.cpp
    )

src/Make.package

@@ -10,6 +10,8 @@ CEXE_sources += incflo_tagging.cpp
 CEXE_sources += incflo_update_density.cpp
 CEXE_sources += incflo_update_tracer.cpp
 CEXE_sources += incflo_update_velocity.cpp
+CEXE_sources += incflo_update_temperature.cpp
+CEXE_sources += incflo_specific_heat.cpp
 CEXE_sources += incflo_utils.cpp
 CEXE_sources += main.cpp
 

src/boundary_conditions/README.org

@@ -1,17 +1,19 @@
 * fillpatch
 
-|        | pi       | po       | mi       | nsw                  | sw                                    | dir_dep              |
-|--------+----------+----------+----------+----------------------+---------------------------------------|-----------------------
-| v_n    | foextrap | foextrap | ext_dir  | ext_dir (0)          | ext_dir (0)                           | ext_dir if inflowing |
-|        |          |          |          |                      |                                       | foextrap otherwise   |
-| v_t    | foextrap | foextrap | ext_dir  | ext_dir (0)          | hoextrap                              | ext_dir if inflowing |
-|        |          |          |          |                      |                                       | foextrap otherwise   |
-| rho    | foextrap | foextrap | ext_dir  | foextrap             | if (advection_type == "BDS") foextrap | ext_dir if inflowing |
-|        |          |          |          |                      | else hoextrap                         | foextrap otherwise   |
-| tracer | foextrap | foextrap | ext_dir  | ext_dir if specified | ext_dir if specified, otherwise       | ext_dir if inflowing |
-|        |          |          |          | foextrap otherwise   | if (advection_type == "BDS") foextrap | foextrap otherwise   |
-|        |          |          |          |                      | else hoextrap                         |                      |
-| force  | foextrap | foextrap | foextrap | foextrap             | foextrap                              | foextrap             |
+|             | pi       | po       | mi       | nsw                    | sw                                    | dir_dep              |
+|-------------+----------+----------+----------+------------------------+---------------------------------------|-----------------------
+| v_n         | foextrap | foextrap | ext_dir  | ext_dir (0)            | ext_dir (0)                           | ext_dir if inflowing |
+|             |          |          |          |                        |                                       | foextrap otherwise   |
+| v_t         | foextrap | foextrap | ext_dir  | ext_dir (0)            | hoextrap                              | ext_dir if inflowing |
+|             |          |          |          |                        |                                       | foextrap otherwise   |
+| rho         | foextrap | foextrap | ext_dir  | foextrap               | if (advection_type == "BDS") foextrap | ext_dir if inflowing |
+|             |          |          |          |                        | else hoextrap                         | foextrap otherwise   |
+| tracer      | foextrap | foextrap | ext_dir  | ext_dir if specified   | ext_dir if specified, otherwise       | ext_dir if inflowing |
+|             |          |          |          | foextrap otherwise     | if (advection_type == "BDS") foextrap | foextrap otherwise   |
+|             |          |          |          |                        | else hoextrap                         |                      |
+| temperature | foextrap | foextrap | ext_dir  | ext_dir if specified   | ext_dir if specified, otherwise       | ext_dir if inflowing |
+|             |          |          |          | reflect_even otherwise | reflect_even                          | foextrap otherwise   |
+| force       | foextrap | foextrap | foextrap | foextrap               | foextrap                              | foextrap             |
 
 * projection
 
@@ -32,5 +34,5 @@
 |---------+---------+---------+-----------+------------------------+------------------------|------------
 | v_n     | Neumann | Neumann | Dirichlet | Dirichlet              | Dirichlet              | Dirichlet |
 | v_t     | Neumann | Neumann | Dirichlet | Dirichlet              | Neumann                | Dirichlet |
-| tracer  | Neumann | Neumann | Dirichlet | Dirchelet if specified | Dirichlet if specified | Dirichlet |
+| scalars | Neumann | Neumann | Dirichlet | Dirchelet if specified | Dirichlet if specified | Dirichlet |
 |         |         |         |           | Neumann otherwise      | Neumann otherwise      |           |

src/boundary_conditions/boundary_conditions.cpp

@@ -13,6 +13,7 @@ void incflo::init_bcs ()
     m_bcrec_velocity.resize(AMREX_SPACEDIM);
     m_bcrec_density.resize(1);
     if (m_ntrac > 0) { m_bcrec_tracer.resize(m_ntrac); }
+    m_bcrec_temperature.resize(1);
 
     auto f = [this] (std::string const& bcid, Orientation ori)
     {
@@ -21,6 +22,7 @@ void incflo::init_bcs ()
                      m_bc_velocity[ori][1] = 0.0;,
                      m_bc_velocity[ori][2] = 0.0;);
         m_bc_tracer[ori].resize(m_ntrac,0.0);
+        m_bc_temperature[ori] = 1.0;
 
         ParmParse pp(bcid);
         std::string bc_type_in = "null";
@@ -41,6 +43,7 @@ void incflo::init_bcs ()
                          m_bcrec_velocity[2].set(ori, BCType::foextrap););
             m_bcrec_density[0].set(ori, BCType::foextrap);
             for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
+            m_bcrec_temperature[0].set(ori, BCType::foextrap);
         }
         else if (bc_type == "pressure_outflow" || bc_type == "po")
         {
@@ -56,6 +59,7 @@ void incflo::init_bcs ()
                          m_bcrec_velocity[2].set(ori, BCType::foextrap););
             m_bcrec_density[0].set(ori, BCType::foextrap);
             for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
+            m_bcrec_temperature[0].set(ori, BCType::foextrap);
         }
         else if (bc_type == "mass_inflow" || bc_type == "mi")
         {
@@ -72,13 +76,15 @@ void incflo::init_bcs ()
 
             pp.query("density", m_bc_density[ori]);
             pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
+            pp.query("temperature", m_bc_temperature[ori]);
 
             // Set mathematical BCs
             AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::ext_dir);,
                          m_bcrec_velocity[1].set(ori, BCType::ext_dir);,
                          m_bcrec_velocity[2].set(ori, BCType::ext_dir););
             m_bcrec_density[0].set(ori, BCType::ext_dir);
             for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::ext_dir); }
+            m_bcrec_temperature[0].set(ori, BCType::ext_dir);
         }
         else if (bc_type == "direction_dependent" || bc_type == "dd" )
         {
@@ -97,12 +103,14 @@ void incflo::init_bcs ()
 
             pp.query("density", m_bc_density[ori]);
             pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
+            pp.query("temperature", m_bc_temperature[ori]);
 
             AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::direction_dependent);,
                          m_bcrec_velocity[1].set(ori, BCType::direction_dependent);,
                          m_bcrec_velocity[2].set(ori, BCType::direction_dependent););
             m_bcrec_density[0].set(ori, BCType::direction_dependent);
             for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::direction_dependent); }
+            m_bcrec_temperature[0].set(ori, BCType::direction_dependent);
         }
         else if (bc_type == "no_slip_wall" || bc_type == "nsw")
         {
@@ -126,6 +134,7 @@ void incflo::init_bcs ()
             // We potentially read in values at no-slip walls in the event that the
             // tracer has Dirichlet bcs
             pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
+            pp.query("temperature", m_bc_temperature[ori]);
 
             // Set mathematical BCs
             AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::ext_dir);,
@@ -137,6 +146,11 @@ void incflo::init_bcs ()
             } else {
                 for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
             }
+            if ( pp.contains("temperature") ) {
+                for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::ext_dir); }
+            } else {
+                for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::reflect_even); }
+            }
         }
         else if (bc_type == "slip_wall" || bc_type == "sw")
         {
@@ -151,6 +165,7 @@ void incflo::init_bcs ()
             // We potentially read in values at slip walls in the event that the
             // tracer has Dirichlet bcs
             pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
+            pp.query("temperature", m_bc_temperature[ori]);
 
             // Tangential directions have hoextrap
             AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::hoextrap);,
@@ -175,6 +190,11 @@ void incflo::init_bcs ()
                     }
                 }
             }
+            if ( pp.contains("temperature") ) {
+                for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::ext_dir); }
+            } else {
+                for (auto& b : m_bcrec_temperature) { b.set(ori, BCType::reflect_even); }
+            }
         }
         else if (bc_type == "mixed" )
         {
@@ -208,13 +228,15 @@ void incflo::init_bcs ()
 
             pp.query("density", m_bc_density[ori]);
             pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
+            pp.query("temperature", m_bc_temperature[ori]);
 
             // Set mathematical BCs. BC_mask will handle Dirichlet part.
             AMREX_D_TERM(m_bcrec_velocity[0].set(ori, BCType::foextrap);,
                          m_bcrec_velocity[1].set(ori, BCType::foextrap);,
                          m_bcrec_velocity[2].set(ori, BCType::foextrap););
             m_bcrec_density[0].set(ori, BCType::foextrap);
             for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::foextrap); }
+            m_bcrec_temperature[0].set(ori, BCType::foextrap);
         }
         else
         {
@@ -231,6 +253,7 @@ void incflo::init_bcs ()
                              m_bcrec_velocity[2].set(ori, BCType::int_dir););
                 m_bcrec_density[0].set(ori, BCType::int_dir);
                 for (auto& b : m_bcrec_tracer) { b.set(ori, BCType::int_dir); }
+                m_bcrec_temperature[0].set(ori, BCType::int_dir);
             } else {
                 amrex::Abort("Wrong BC type for periodic boundary");
             }
@@ -298,6 +321,16 @@ void incflo::init_bcs ()
             (m_bcrec_tracer_d.data(), m_bcrec_tracer.data(), sizeof(BCRec)*m_ntrac);
     }
 
+    if (m_use_temperature) {
+        m_bcrec_temperature_d.resize(1);
+#ifdef AMREX_USE_GPU
+        Gpu::htod_memcpy
+#else
+        std::memcpy
+#endif
+            (m_bcrec_temperature_d.data(), m_bcrec_temperature.data(), sizeof(BCRec));
+    }
+
     // force
     {
         const int ncomp = std::max(m_ntrac, AMREX_SPACEDIM);