Heat Source Implementation Issues in 3D-Viscous Cases

[Ihavenoidea5]

Hey everyone,

I've been able to implement an imposed heat source successfully in 2D inviscid & viscous cases, as well as the 3D inviscid case, but the heat source disappears or doesn't work at all when I try to implement the source in the 3D-RANS case using the SA Turb model. I should also mention that I am unable to impose the source if the mesh has irregular geometries in both the 2D-inviscid & viscous case (i.e. shock wedge/generator in rectangular mesh). Would anyone know why this might be happening?

I've modified the following files to get the heat source to work and I am not sure what else to do in order to get the desired results:

•  CDriver.cpp
  

• CConfig.hpp
• CConfig.cpp
• CEulerSolver.cpp
• flow_sources.hpp
• flow_sources.cpp

Please let me know if you have any ideas as to why this might be happening, thanks!

[WallyMaier]

Hi,
There shouldn't be any differences in the implementations. But it hard to understand without more details. Can you post any error messages you are seeing?

Can you post the implementation of source terms?

[pcarruscag]

Push your code to github so that people can compare the changes directly. You are asking for too much if you want someone to figure out what you changed

[Ihavenoidea5]

Yes, apologies for the hassle, I'm still pretty new to all of this.

The following bits of code are all additions to the source code in order implement the source term.

SU2/SU2_CFD/src/drivers/CDriver.cpp (Lines 1835-1837)

      }
      else if (config->GetImposedSource() == YES) {
        numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceImposed(nDim, nVar_Flow, config);
      }

SU2/Common/include/CConfig.hpp

(Lines 67-73)

  bool ImposedSourceTerm;   /*!< \brief Flag to know if there is an imposed source. */
  su2double Source_xorigin; /*!< \brief X-Position of Force.*/
  su2double Source_yorigin; /*!< \brief Y-Position of Force.*/
  su2double Source_zorigin; /*!< \brief Z-Position of Force.*/
  su2double Source_xscale;  /*!< \brief X-Direction Length scale.*/
  su2double Source_yscale;  /*!< \brief Y-Direction Length scale.*/
  su2double Source_zscale;  /*!< \brief Z-Direction Length scale.*/

(Lines 1421-1454)

  /*!
   * \brief Get imposed source term flag
   */
  bool GetImposedSource(void) const { return ImposedSourceTerm; }

  /*!
   * \brief Value of the x-position of the source term.
   */
  su2double GetSource_xorigin(void) const { return Source_xorigin; }

  /*!
  * \brief Value of the y-position of the source term.
  */
  su2double GetSource_yorigin(void) const { return Source_yorigin; }

  /*!
   * \brief Value of the z-position of the source term.
   */
  su2double GetSource_zorigin(void) const { return Source_zorigin; }

  /*!
   * \brief Value of the x-direction source term scale.
   */
  su2double GetSource_xscale(void) const { return Source_xscale; }

  /*!
   * \brief Value of the y-direction source term scale.
   */
  su2double GetSource_yscale(void) const { return Source_yscale; }

  /*!
   * \brief Value of the z-direction source term scale.
   */
  su2double GetSource_zscale(void) const { return Source_zscale; }

SU2/Common/src/CConfig.cpp (Lines 1073-1088)

/*--- Options related to the custom source term ---*/

  /*!\brief IMPOSED_SOURCE_TERM \n  DESCRIPTION: turn on imposed source term \ingroup Config*/
  addBoolOption("IMPOSED_SOURCE_TERM", ImposedSourceTerm, false);
   /* DESCRIPTION: X location at which source begins (meters) */
  addDoubleOption("SOURCE_X_ORIGIN", Source_xorigin, 0.0);
   /* DESCRIPTION: Y location at which source begins (meters) */
  addDoubleOption("SOURCE_Y_ORIGIN", Source_yorigin, 0.0);
   /* DESCRIPTION: Z location at which source begins (meters) */
  addDoubleOption("SOURCE_Z_ORIGIN", Source_zorigin, 0.0);
   /* DESCRIPTION: X-direction strength of source term */
  addDoubleOption("SOURCE_XSCALE", Source_xscale, 0.0);
   /* DESCRIPTION: Y-direction strength of source term */
  addDoubleOption("SOURCE_YSCALE", Source_yscale, 0.0);
   /* DESCRIPTION: Z-direction strength of source term */
  addDoubleOption("SOURCE_ZSCALE", Source_zscale, 0.0);

SU2/SU2_CFD/src/solvers/CEulerSolver.cpp
(Line 1996)
const bool imposed_source = config->GetImposedSource();

(Lines 2007-2030)

if (imposed_source) {

    /*--- loop over points ---*/
    SU2_OMP_FOR_DYN(omp_chunk_size)
    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {

      /*--- Set solution  ---*/
      numerics->SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint));

      /*--- Set control volume ---*/
      numerics->SetVolume(geometry->nodes->GetVolume(iPoint));

      /*--- Set coordinates ---*/
      numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));

      /*--- Compute Source term Residual ---*/
      auto residual = numerics->ComputeResidual(config);

      /*--- Add Residual ---*/
      LinSysRes.AddBlock(iPoint, residual);

    }

  }

SU2/SU2_CFD/include/numerics/flow/flow_sources.hpp (Lines 422-448)

/*!
 * \class CSourceImposed
 * \brief Class for the source term integration of the imposed source term.
 * \ingroup SourceDiscr
 * \author A.D
 */
class CSourceImposed final : public CSourceBase_Flow {
private:
  bool compressible;

public:
  /*!
   * \param[in] val_nDim - Number of dimensions of the problem.
   * \param[in] val_nVar - Number of variables of the problem.
   * \param[in] config - Definition of the particular problem.
   */
  CSourceImposed(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);

  /*!
   * \brief Source term integration for the poissonal potential.
   * \param[in] config - Definition of the particular problem.
   * \return Lightweight const-view of residual and Jacobian.
   */
  ResidualType<> ComputeResidual(const CConfig* config) override;

};

SU2/SU2_CFD/src/numerics/flow/flow_sources.cpp (Lines 821-870)

CSourceImposed::CSourceImposed(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config) :
		CSourceBase_Flow(val_nDim, val_nVar, config) {

  bool compressible = (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE);

}

CNumerics::ResidualType<> CSourceImposed::ComputeResidual(const CConfig *config) {
  su2double xorigin = config->GetSource_xorigin();
  su2double yorigin = config->GetSource_yorigin();
  su2double zorigin = config->GetSource_zorigin();

  su2double xscale = config->GetSource_xscale();
  su2double yscale = config->GetSource_yscale();
  su2double zscale = config->GetSource_zscale();

  const su2double heating_magnitude = 10000000.0; // Total power; should be an input.

  if (compressible && nDim == 3) {
    // shifted and scaled coordinates
    su2double x = (Coord_i[0] - xorigin)/xscale;
    su2double y = (Coord_i[1] - yorigin)/yscale;
    su2double z = (Coord_i[2] - zorigin)/zscale;

    su2double spatial_distribution = exp(-(x*x + y*y + z*z)) / (xscale*yscale*zscale*pow(M_PI,1.5));
    su2double power_volume = heating_magnitude * spatial_distribution;

    residual[0] = 0.0;
    residual[1] = 0.0;
    residual[2] = 0.0;
    residual[3] = 0.0;
    residual[4] = -Volume * power_volume;
  }

  if (compressible && nDim == 2) {
    // shifted and scaled coordinates
    su2double x = (Coord_i[0] - xorigin)/xscale;
    su2double y = (Coord_i[1] - yorigin)/yscale;

    su2double spatial_distribution = exp(-(x*x + y*y)) / (xscale*yscale*M_PI);
    su2double power_volume = heating_magnitude * spatial_distribution;

    residual[0] = 0.0;
    residual[1] = 0.0;
    residual[2] = 0.0;
    residual[3] = -Volume * power_volume;
  }

  return ResidualType<>(residual, jacobian, nullptr);
}

Thank you for your patience and assistance, please let me know if you would like the code presented in a different manner.

[WallyMaier]

Generally, you can either push you code via a branch on Github or make a pull request if you want your contribution to in the main code.

But from what I see here, it looks likes everything should be fine. I would check that the value of your residuals are what you expect