Heat solver using scalar framework - Part 1

SU2_CFD/include/solvers/CHeatSolver.hpp

@@ -27,7 +27,7 @@
 
 #pragma once
 
-#include "CSolver.hpp"
+#include "CScalarSolver.hpp"
 #include "../variables/CHeatVariable.hpp"
 
 /*!
@@ -36,48 +36,38 @@
  * \author O. Burghardt
  * \version 7.4.0 "Blackbird"
  */
-class CHeatSolver final : public CSolver {
+class CHeatSolver final : public CScalarSolver<CHeatVariable> {
 protected:
   static constexpr size_t MAXNDIM = 3; /*!< \brief Max number of space dimensions, used in some static arrays. */
   static constexpr size_t MAXNVAR = 1; /*!< \brief Max number of variables, for static arrays. */
 
   const bool flow; /*!< \brief Use solver as a scalar transport equation of Temperature for the inc solver. */
   const bool heat_equation; /*!< \brief use solver for heat conduction in solids. */
 
-  unsigned short nVarFlow, nMarker;
+  su2double Global_Delta_Time = 0.0, Global_Delta_UnstTimeND = 0.0;
+
+  unsigned short nVarFlow;
   vector<vector<su2double> > HeatFlux;
   vector<su2double> HeatFlux_per_Marker;
   su2double Total_HeatFlux;
   su2double AllBound_HeatFlux;
   vector<su2double> AverageT_per_Marker;
   su2double Total_AverageT;
   su2double AllBound_AverageT;
-  vector<su2double>  Primitive_Flow_i;
+  vector<su2double> Primitive_Flow_i;
   vector<su2double> Primitive_Flow_j;
   vector<su2double> Surface_Areas;
   su2double Total_HeatFlux_Areas;
   su2double Total_HeatFlux_Areas_Monitor;
   vector<su2activematrix> ConjugateVar;
 
-  CHeatVariable* nodes = nullptr;  /*!< \brief The highest level in the variable hierarchy this solver can safely use. */
-
-  /*!
-   * \brief Return nodes to allow CSolver::base_nodes to be set.
-   */
-  inline CVariable* GetBaseClassPointerToNodes() override { return nodes; }
-
 public:
 
   /*!
    * \brief Constructor of the class.
    */
   CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CHeatSolver(void) override;
-
   /*!
    * \brief Restart residual and compute gradients.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -247,15 +237,6 @@ class CHeatSolver final : public CSolver {
                                   CConfig *config,
                                   unsigned short val_marker) override;
 
-  /*!
-   * \brief Impose a periodic boundary condition by summing contributions from the complete control volume.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] numerics - Description of the numerical method.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void BC_Periodic(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config) final;
-
   /*!
    * \brief Set the conjugate heat variables.
    * \param[in] val_marker        - marker index
@@ -330,26 +311,6 @@ class CHeatSolver final : public CSolver {
     }
   }
 
-  /*!
-   * \brief Update the solution using an implicit solver.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void ImplicitEuler_Iteration(CGeometry *geometry,
-                               CSolver **solver_container,
-                               CConfig *config) override;
-
-  /*!
-   * \brief Update the solution using an explicit solver.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void ExplicitEuler_Iteration(CGeometry *geometry,
-                               CSolver **solver_container,
-                               CConfig *config) override;
-
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CRadP1Solver.hpp

@@ -241,7 +241,7 @@ class CRadP1Solver final: public CRadSolver {
    * \param[in] Iteration - Index of the current iteration.
    */
   void SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config,
-                    unsigned short iMesh, unsigned long Iteration) override ;
+                    unsigned short iMesh, unsigned long Iteration) override;
 
   /*!
    * \brief Set the freestream temperature.

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -156,6 +156,8 @@ class CScalarSolver : public CSolver {
   void BC_Fluid_Interface_impl(const SolverSpecificNumericsFunc& SolverSpecificNumerics, CGeometry *geometry,
                                CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics,
                                CConfig *config) {
+    if (solver_container[FLOW_SOL] == nullptr) return;
+
     const auto nPrimVar = solver_container[FLOW_SOL]->GetnPrimVar();
     su2activevector PrimVar_j(nPrimVar);
     su2double solution_j[MAXNVAR] = {0.0};
@@ -436,6 +438,17 @@ class CScalarSolver : public CSolver {
     END_SU2_OMP_FOR
   }
 
+  /*!
+   * \brief This base implementation simply copies the time step of the flow solver.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   * \param[in] Iteration - Index of the current iteration.
+   */
+  void SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config,
+                    unsigned short iMesh, unsigned long Iteration) override;
+
   /*!
    * \brief Prepare an implicit iteration.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -466,7 +479,7 @@ class CScalarSolver : public CSolver {
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  void ImplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) override;
+  void ImplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) final;
 
   /*!
    * \brief Set the total residual adding the term that comes from the Dual Time-Stepping Strategy.
@@ -478,7 +491,7 @@ class CScalarSolver : public CSolver {
    * \param[in] RunTime_EqSystem - System of equations which is going to be solved.
    */
   void SetResidual_DualTime(CGeometry* geometry, CSolver** solver_container, CConfig* config, unsigned short iRKStep,
-                            unsigned short iMesh, unsigned short RunTime_EqSystem) final;
+                            unsigned short iMesh, unsigned short RunTime_EqSystem) override;
 
   /*!
    * \brief Load a solution from a restart file.

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -424,10 +424,21 @@ void CScalarSolver<VariableType>::BC_Far_Field(CGeometry* geometry, CSolver** so
 }
 
 template <class VariableType>
-void CScalarSolver<VariableType>::PrepareImplicitIteration(CGeometry* geometry, CSolver** solver_container,
-                                                           CConfig* config) {
+void CScalarSolver<VariableType>::SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config,
+                                               unsigned short iMesh, unsigned long Iteration) {
   const auto flowNodes = solver_container[FLOW_SOL]->GetNodes();
 
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    su2double dt = nodes->GetLocalCFL(iPoint) / flowNodes->GetLocalCFL(iPoint) * flowNodes->GetDelta_Time(iPoint);
+    nodes->SetDelta_Time(iPoint, dt);
+  }
+  END_SU2_OMP_FOR
+}
+
+template <class VariableType>
+void CScalarSolver<VariableType>::PrepareImplicitIteration(CGeometry* geometry, CSolver** solver_container,
+                                                           CConfig* config) {
   /*--- Set shared residual variables to 0 and declare
    *    local ones for current thread to work on. ---*/
 
@@ -440,11 +451,8 @@ void CScalarSolver<VariableType>::PrepareImplicitIteration(CGeometry* geometry,
 
   SU2_OMP_FOR_(schedule(static, omp_chunk_size) SU2_NOWAIT)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
-    /// TODO: This could be the SetTime_Step of this solver.
-    su2double dt = nodes->GetLocalCFL(iPoint) / flowNodes->GetLocalCFL(iPoint) * flowNodes->GetDelta_Time(iPoint);
-    nodes->SetDelta_Time(iPoint, dt);
-
     /*--- Modify matrix diagonal to improve diagonal dominance. ---*/
+    const su2double dt = nodes->GetDelta_Time(iPoint);
 
     if (dt != 0.0) {
       su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
@@ -476,8 +484,6 @@ void CScalarSolver<VariableType>::CompleteImplicitIteration(CGeometry* geometry,
                                                             CConfig* config) {
   const bool compressible = (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE);
 
-  const auto flowNodes = solver_container[FLOW_SOL]->GetNodes();
-
   ComputeUnderRelaxationFactor(config);
 
   /*--- Update solution (system written in terms of increments) ---*/
@@ -486,6 +492,8 @@ void CScalarSolver<VariableType>::CompleteImplicitIteration(CGeometry* geometry,
     /*--- Update the scalar solution. For transport equations, where Solution is not equivalent with the transported
      * quantity, multiply the respective factor.  ---*/
     if (Conservative) {
+      const auto flowNodes = solver_container[FLOW_SOL]->GetNodes();
+
       SU2_OMP_FOR_STAT(omp_chunk_size)
       for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
         /*--- Multiply the Solution var with density to get the conservative transported quantity, if necessary. ---*/
@@ -547,16 +555,13 @@ void CScalarSolver<VariableType>::ImplicitEuler_Iteration(CGeometry* geometry, C
 template <class VariableType>
 void CScalarSolver<VariableType>::ExplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container,
                                                           CConfig* config) {
-  const auto flowNodes = solver_container[FLOW_SOL]->GetNodes();
-
   /*--- Local residual variables for current thread ---*/
   su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0};
   unsigned long idxMax[MAXNVAR] = {0};
 
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
-    const su2double dt = nodes->GetLocalCFL(iPoint) / flowNodes->GetLocalCFL(iPoint) * flowNodes->GetDelta_Time(iPoint);
-    nodes->SetDelta_Time(iPoint, dt);
+    const su2double dt = nodes->GetDelta_Time(iPoint);
     const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
 
     for (auto iVar = 0u; iVar < nVar; iVar++) {

SU2_CFD/include/variables/CHeatVariable.hpp

@@ -27,20 +27,18 @@
 
 #pragma once
 
-#include "CVariable.hpp"
+#include "CScalarVariable.hpp"
 
 /*!
  * \class CHeatVariable
  * \brief Class for defining the variables of the finite-volume heat equation solver.
  * \author O. Burghardt
  * \version 7.4.0 "Blackbird"
  */
-class CHeatVariable final : public CVariable {
-protected:
-  CVectorOfMatrix& Gradient_Reconstruction;  /*!< \brief Reference to the gradient of the primitive variables for MUSCL reconstruction for the convective term */
-  CVectorOfMatrix Gradient_Aux;              /*!< \brief Auxiliary structure to store a second gradient for reconstruction, if required. */
-
+class CHeatVariable final : public CScalarVariable {
 public:
+  static constexpr size_t MAXNVAR = 1; /*!< \brief Max number of variables, for static arrays. */
+
   /*!
    * \brief Constructor of the class.
    * \param[in] heat - Values of the Heat solution (initialization value).
@@ -51,29 +49,10 @@ class CHeatVariable final : public CVariable {
    */
   CHeatVariable(su2double heat, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CHeatVariable() override = default;
-
-  /*!
-   * \brief Get the array of the reconstruction variables gradient at a node.
-   * \param[in] iPoint - Index of the current node.
-   * \return Array of the reconstruction variables gradient at a node.
-   */
-  inline CMatrixView<su2double> GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; }
-
-  /*!
-   * \brief Get the reconstruction gradient for primitive variable at all points.
-   * \return Reference to variable reconstruction gradient.
-   */
-  inline CVectorOfMatrix& GetGradient_Reconstruction() final { return Gradient_Reconstruction; }
-  inline const CVectorOfMatrix& GetGradient_Reconstruction() const final { return Gradient_Reconstruction; }
-
   /*!
    * \brief Get the temperature of the point.
    * \return Value of the temperature of the point.
    */
-  inline su2double GetTemperature(unsigned long iPoint) const final { return Solution(iPoint,0); }
+  inline su2double GetTemperature(unsigned long iPoint) const final { return Solution(iPoint, 0); }
 
 };