Explicit Euler in CScalarSolver

Common/src/CConfig.cpp

@@ -4332,6 +4332,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
   }
 
+  if(Kind_TimeIntScheme_Turb != EULER_IMPLICIT &&
+     Kind_TimeIntScheme_Turb != EULER_EXPLICIT){
+    SU2_MPI::Error("Only TIME_DISCRE_TURB = EULER_IMPLICIT, EULER_EXPLICIT have been implemented.", CURRENT_FUNCTION);
+  }
+
   if (nIntCoeffs == 0) {
     nIntCoeffs = 2;
     Int_Coeffs = new su2double[2]; Int_Coeffs[0] = 0.25; Int_Coeffs[1] = 0.5;

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -775,13 +775,8 @@ class CFVMFlowSolverBase : public CSolver {
     const su2double RK_FuncCoeff[] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
     const su2double RK_TimeCoeff[] = {0.5, 0.5, 1.0, 1.0};
 
-    /*--- Set shared residual variables to 0 and declare
-     *    local ones for current thread to work on. ---*/
-
-    SetResToZero();
-
+    /*--- Local residual variables for current thread ---*/
     su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0};
-    const su2double* coordMax[MAXNVAR] = {nullptr};
     unsigned long idxMax[MAXNVAR] = {0};
 
     /*--- Update the solution and residuals ---*/
@@ -832,23 +827,13 @@ class CFVMFlowSolverBase : public CSolver {
           }
 
           /*--- Update residual information for current thread. ---*/
-          resRMS[iVar] += Res*Res;
-          if (fabs(Res) > resMax[iVar]) {
-            resMax[iVar] = fabs(Res);
-            idxMax[iVar] = iPoint;
-            coordMax[iVar] = geometry->nodes->GetCoord(iPoint);
-          }
+          ResidualReductions_PerThread(iPoint, iVar, Res, resRMS, resMax, idxMax);
         }
       }
       END_SU2_OMP_FOR
       /*--- Reduce residual information over all threads in this rank. ---*/
-      SU2_OMP_CRITICAL
-      for (unsigned short iVar = 0; iVar < nVar; iVar++) {
-        Residual_RMS[iVar] += resRMS[iVar];
-        AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]);
-      }
-      END_SU2_OMP_CRITICAL
-      SU2_OMP_BARRIER
+      ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax);
+
     }
 
     /*--- MPI solution ---*/
@@ -857,9 +842,6 @@ class CFVMFlowSolverBase : public CSolver {
     CompleteComms(geometry, config, SOLUTION);
 
     if (!adjoint) {
-      /*--- Compute the root mean square residual ---*/
-      SetResidual_RMS(geometry, config);
-
       /*--- For verification cases, compute the global error metrics. ---*/
       ComputeVerificationError(geometry, config);
     }
@@ -893,12 +875,8 @@ class CFVMFlowSolverBase : public CSolver {
 
     const bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
 
-    /*--- Set shared residual variables to 0 and declare local ones for current thread to work on. ---*/
-
-    SetResToZero();
-
+    /*--- Local residual variables for current thread ---*/
     su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0};
-    const su2double* coordMax[MAXNVAR] = {nullptr};
     unsigned long idxMax[MAXNVAR] = {0};
 
     /*--- Add pseudotime term to Jacobian. ---*/
@@ -948,26 +926,15 @@ class CFVMFlowSolverBase : public CSolver {
         LinSysRes[total_index] = - (LinSysRes[total_index] + local_Res_TruncError[iVar]);
         LinSysSol[total_index] = 0.0;
 
-        su2double Res = fabs(LinSysRes[total_index]);
-        resRMS[iVar] += Res*Res;
-        if (Res > resMax[iVar]) {
-          resMax[iVar] = Res;
-          idxMax[iVar] = iPoint;
-          coordMax[iVar] = geometry->nodes->GetCoord(iPoint);
-        }
+        /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/
+        ResidualReductions_PerThread(iPoint, iVar, LinSysRes[total_index], resRMS, resMax, idxMax);
       }
     }
     END_SU2_OMP_FOR
-    SU2_OMP_CRITICAL
-    for (unsigned short iVar = 0; iVar < nVar; iVar++) {
-      Residual_RMS[iVar] += resRMS[iVar];
-      AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]);
-    }
-    END_SU2_OMP_CRITICAL
-    SU2_OMP_BARRIER
 
-    /*--- Compute the root mean square residual ---*/
-    SetResidual_RMS(geometry, config);
+    /*--- "Add" residuals from all threads to global residual variables. ---*/
+    ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax);
+
   }
 
   /*!

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -259,6 +259,14 @@ class CScalarSolver : public CSolver {
    */
   void CompleteImplicitIteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) final;
 
+  /*!
+   * \brief Update the solution using the explicit Euler scheme.
+   * \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) final;
+
   /*!
    * \brief Update the solution using an implicit solver.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -287,7 +295,8 @@ class CScalarSolver : public CSolver {
    * \param[in] val_iter - Current external iteration number.
    * \param[in] val_update_geo - Flag for updating coords and grid velocity.
    */
-  virtual void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, bool val_update_geo) override = 0;
+  void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter,
+                           bool val_update_geo) override = 0;
 
   /*!
    * \brief Scalar solvers support OpenMP+MPI.

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -79,7 +79,8 @@ CScalarSolver<VariableType>::~CScalarSolver() {
 
 template <class VariableType>
 void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver** solver_container,
-                                               CNumerics** numerics_container, CConfig* config, unsigned short iMesh) {
+                                                  CNumerics** numerics_container, CConfig* config,
+                                                  unsigned short iMesh) {
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool muscl = config->GetMUSCL_Turb();
   const bool limiter = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER);
@@ -256,7 +257,7 @@ void CScalarSolver<VariableType>::SumEdgeFluxes(CGeometry* geometry) {
 
 template <class VariableType>
 void CScalarSolver<VariableType>::BC_Periodic(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics,
-                                           CConfig* config) {
+                                              CConfig* config) {
   /*--- Complete residuals for periodic boundary conditions. We loop over
    the periodic BCs in matching pairs so that, in the event that there are
    adjacent periodic markers, the repeated points will have their residuals
@@ -271,7 +272,7 @@ void CScalarSolver<VariableType>::BC_Periodic(CGeometry* geometry, CSolver** sol
 
 template <class VariableType>
 void CScalarSolver<VariableType>::PrepareImplicitIteration(CGeometry* geometry, CSolver** solver_container,
-                                                        CConfig* config) {
+                                                           CConfig* config) {
   const auto flowNodes = solver_container[FLOW_SOL]->GetNodes();
 
   /*--- Set shared residual variables to 0 and declare
@@ -280,7 +281,6 @@ void CScalarSolver<VariableType>::PrepareImplicitIteration(CGeometry* geometry,
   SetResToZero();
 
   su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0};
-  const su2double* coordMax[MAXNVAR] = {nullptr};
   unsigned long idxMax[MAXNVAR] = {0};
 
   /*--- Build implicit system ---*/
@@ -308,31 +308,19 @@ void CScalarSolver<VariableType>::PrepareImplicitIteration(CGeometry* geometry,
       LinSysRes[total_index] = -LinSysRes[total_index];
       LinSysSol[total_index] = 0.0;
 
-      su2double Res = fabs(LinSysRes[total_index]);
-      resRMS[iVar] += Res * Res;
-      if (Res > resMax[iVar]) {
-        resMax[iVar] = Res;
-        idxMax[iVar] = iPoint;
-        coordMax[iVar] = geometry->nodes->GetCoord(iPoint);
-      }
+      /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/
+      ResidualReductions_PerThread(iPoint, iVar, LinSysRes[total_index], resRMS, resMax, idxMax);
     }
   }
   END_SU2_OMP_FOR
-  SU2_OMP_CRITICAL
-  for (unsigned short iVar = 0; iVar < nVar; iVar++) {
-    Residual_RMS[iVar] += resRMS[iVar];
-    AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]);
-  }
-  END_SU2_OMP_CRITICAL
-  SU2_OMP_BARRIER
 
-  /*--- Compute the root mean square residual ---*/
-  SetResidual_RMS(geometry, config);
+  /*--- "Add" residuals from all threads to global residual variables. ---*/
+  ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax);
 }
 
 template <class VariableType>
 void CScalarSolver<VariableType>::CompleteImplicitIteration(CGeometry* geometry, CSolver** solver_container,
-                                                         CConfig* config) {
+                                                            CConfig* config) {
   const bool compressible = (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE);
 
   const auto flowNodes = solver_container[FLOW_SOL]->GetNodes();
@@ -380,7 +368,7 @@ void CScalarSolver<VariableType>::CompleteImplicitIteration(CGeometry* geometry,
 
 template <class VariableType>
 void CScalarSolver<VariableType>::ImplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container,
-                                                       CConfig* config) {
+                                                          CConfig* config) {
   PrepareImplicitIteration(geometry, solver_container, config);
 
   /*--- Solve or smooth the linear system. ---*/
@@ -404,10 +392,41 @@ void CScalarSolver<VariableType>::ImplicitEuler_Iteration(CGeometry* geometry, C
   CompleteImplicitIteration(geometry, solver_container, config);
 }
 
+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 Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
+
+    for (auto iVar = 0u; iVar < nVar; iVar++) {
+      /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/
+      ResidualReductions_PerThread(iPoint, iVar, LinSysRes(iPoint, iVar), resRMS, resMax, idxMax);
+      /*--- Explicit Euler step: ---*/
+      LinSysSol(iPoint, iVar) = -dt / Vol * LinSysRes(iPoint, iVar);
+    }
+  }
+  END_SU2_OMP_FOR
+
+  /*--- "Add" residuals from all threads to global residual variables. ---*/
+  ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax);
+
+  /*--- Use LinSysSol for solution update. ---*/
+  CompleteImplicitIteration(geometry, solver_container, config);
+}
+
 template <class VariableType>
 void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSolver** solver_container, CConfig* config,
-                                                    unsigned short iRKStep, unsigned short iMesh,
-                                                    unsigned short RunTime_EqSystem) {
+                                                       unsigned short iRKStep, unsigned short iMesh,
+                                                       unsigned short RunTime_EqSystem) {
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST);
   const bool second_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);