diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 794622bb1a1..7becae80917 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -449,6 +449,11 @@ class CConfig { TURBOMACHINERY_TYPE *Kind_TurboMachinery; su2vector Kind_TurboInterface; + /* Turbomachinery objective functions */ + su2double *EntropyGeneration; + su2double *TotalPressureLoss; + su2double *KineticEnergyLoss; + /* Gradient smoothing options */ su2double SmoothingEps1; /*!< \brief Parameter for the identity part in gradient smoothing. */ su2double SmoothingEps2; /*!< \brief Parameter for the Laplace part in gradient smoothing. */ @@ -8147,6 +8152,27 @@ class CConfig { * \param[in] val_surface_species_variance - Value of the species variance. */ void SetSurface_Species_Variance(unsigned short val_marker, su2double val_surface_species_variance) { Surface_Species_Variance[val_marker] = val_surface_species_variance; } + + /*! + * \brief Set entropy generation for a turbomachinery zone + * \param[in] val_iZone - zone index + * \param[in] val_entropy_generation - value of entropy generation + */ + void SetEntropyGeneration(unsigned short val_iZone, su2double val_entropy_generation) { EntropyGeneration[val_iZone] = val_entropy_generation; } + + /*! + * \brief Get total pressure loss for a turbomachinery zone + * \param[in] val_iZone - zone index + * \param[in] val_total_pressure_loss - value of total pressure loss + */ + void SetTotalPressureLoss(unsigned short val_iZone, su2double val_total_pressure_loss) { TotalPressureLoss[val_iZone] = val_total_pressure_loss; } + + /*! + * \brief Get kinetic energy loss for a turbomachinery zone + * \param[in] val_iZone - zone index + * \param[in] val_kinetic_energy_loss - value of kinetic energy loss + */ + void SetKineticEnergyLoss(unsigned short val_iZone, su2double val_kinetic_energy_loss) { KineticEnergyLoss[val_iZone] = val_kinetic_energy_loss; } /*! * \brief Get the back pressure (static) at an outlet boundary. @@ -10009,6 +10035,13 @@ class CConfig { */ short FindInterfaceMarker(unsigned short iInterface) const; + /*! + * \brief Find the marker index (if any) that is part of a mixing plane interface pair. + * \param[in] nMarker - Number of the marker in a zone being tested, starting at 0. + * \return value > 1 if (on this mpi rank) the zone defined by config is part of the mixing plane. + */ + short FindMixingPlaneInterfaceMarker(unsigned short nMarker) const; + /*! * \brief Get whether or not to save solution data to libROM. * \return True if specified in config file. diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index f515e03e26e..7b428ae5c27 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -776,7 +776,7 @@ class CGeometry { inline virtual void GatherInOutAverageValues(CConfig* config, bool allocate) {} /*! - * \brief Store all the turboperformance in the solver in ZONE_0. + * \brief Store all the turboperformance in the solver in final zone. * \param[in] donor_geometry - Solution from the donor mesh. * \param[in] target_geometry - Solution from the target mesh. * \param[in] donorZone - counter of the donor solution diff --git a/Common/include/interface_interpolation/CInterpolator.hpp b/Common/include/interface_interpolation/CInterpolator.hpp index ae853554085..3a4a0c3a84c 100644 --- a/Common/include/interface_interpolation/CInterpolator.hpp +++ b/Common/include/interface_interpolation/CInterpolator.hpp @@ -101,6 +101,12 @@ class CInterpolator { }; vector > targetVertices; /*! \brief Donor information per marker per vertex of the target. */ + struct CSpanDonorInfo { + unsigned long donorSpan; // Refers to donor span + su2double coefficient; // Refers to coefficient + }; + vector> targetSpans; // > + /*! * \brief Constructor of the class. * \param[in] geometry_container - Geometrical definition of the problem. diff --git a/Common/include/interface_interpolation/CInterpolatorFactory.hpp b/Common/include/interface_interpolation/CInterpolatorFactory.hpp index 63c5eaf23e2..89172dc40e5 100644 --- a/Common/include/interface_interpolation/CInterpolatorFactory.hpp +++ b/Common/include/interface_interpolation/CInterpolatorFactory.hpp @@ -43,5 +43,5 @@ namespace CInterpolatorFactory { */ CInterpolator* CreateInterpolator(CGeometry**** geometry_container, const CConfig* const* config, const CInterpolator* transpInterpolator, unsigned iZone, unsigned jZone, - bool verbose = true); + bool mixing_plane, bool verbose = true); } // namespace CInterpolatorFactory diff --git a/Common/include/interface_interpolation/CMixingPlane.hpp b/Common/include/interface_interpolation/CMixingPlane.hpp new file mode 100644 index 00000000000..e396de0fd84 --- /dev/null +++ b/Common/include/interface_interpolation/CMixingPlane.hpp @@ -0,0 +1,43 @@ +/*! + * \file CMixingPlane.hpp + * \brief Header of mixing plane interpolation methods. + * \author J. Kelly + * \version 8.3.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + + #pragma once + #include "CInterpolator.hpp" + +/*! + * \brief Mixing plane interpolation. + * \note This contains several interpolation methods used in the mixing plane interpolation + * and enables the mixing state class structure for proper recording in AD mode + * \ingroup Interfaces + */ +class CMixingPlane final : public CInterpolator { + public: + CMixingPlane(CGeometry**** geometry_container, const CConfig* const* config, unsigned int iZone, + unsigned int jZone); + + void SetTransferCoeff(const CConfig* const* config) override; +}; \ No newline at end of file diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bf99257b138..cf0681846d3 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2067,6 +2067,9 @@ enum ENUM_OBJECTIVE { TOPOL_DISCRETENESS = 63, /*!< \brief Measure of the discreteness of the current topology. */ TOPOL_COMPLIANCE = 64, /*!< \brief Measure of the discreteness of the current topology. */ STRESS_PENALTY = 65, /*!< \brief Penalty function of VM stresses above a maximum value. */ + ENTROPY_GENERATION = 80, /*!< \brief Entropy generation turbomachinery objective function. */ + TOTAL_PRESSURE_LOSS = 81, /*!< \brief Total pressure loss turbomachinery objective function. */ + KINETIC_ENERGY_LOSS = 82 /*!< \breif Kinetic energy loss coefficient turbomachinery objective function. */ }; static const MapType Objective_Map = { MakePair("DRAG", DRAG_COEFFICIENT) @@ -2109,6 +2112,9 @@ static const MapType Objective_Map = { MakePair("TOPOL_DISCRETENESS", TOPOL_DISCRETENESS) MakePair("TOPOL_COMPLIANCE", TOPOL_COMPLIANCE) MakePair("STRESS_PENALTY", STRESS_PENALTY) + MakePair("ENTROPY_GENERATION", ENTROPY_GENERATION) + MakePair("TOTAL_PRESSURE_LOSS", TOTAL_PRESSURE_LOSS) + MakePair("KINETIC_ENERGY_LOSS", KINETIC_ENERGY_LOSS) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 0a9cca0f63f..956d01dfdac 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1049,6 +1049,11 @@ void CConfig::SetPointersNull() { nBlades = nullptr; FreeStreamTurboNormal = nullptr; + /*--- Turbomachinery Objective Functions ---*/ + EntropyGeneration = nullptr; + TotalPressureLoss = nullptr; + KineticEnergyLoss = nullptr; + top_optim_kernels = nullptr; top_optim_kernel_params = nullptr; top_optim_filter_radius = nullptr; @@ -3479,6 +3484,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250; } } + } else if (Multizone_Problem && DiscreteAdjoint) { + SU2_MPI::Error(string("OUTPUT_WRT_FREQ cannot be specified for this solver " + "(writing of restart and sensitivity files not possible for multizone discrete adjoint during runtime yet).\n" + "Please remove this option from the config file, output files will be written when solver finalizes.\n"), CURRENT_FUNCTION); } else if (nVolumeOutputFrequencies < nVolumeOutputFiles) { /*--- If there are fewer frequencies than files, repeat the last frequency. * This is useful to define 1 frequency for the restart file and 1 frequency for all the visualization files. ---*/ @@ -6172,6 +6181,11 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_CfgFile_ZoneInterface[iMarker_CfgFile] = YES; } + /*--- Allocate memory for turbomachinery objective functions ---*/ + EntropyGeneration = new su2double[nZone] (); + TotalPressureLoss = new su2double[nZone] (); + KineticEnergyLoss = new su2double[nZone] (); + /*--- Identification of Turbomachinery markers and flag them---*/ for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { @@ -6956,6 +6970,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case TOPOL_DISCRETENESS: cout << "Topology discreteness objective function." << endl; break; case TOPOL_COMPLIANCE: cout << "Topology compliance objective function." << endl; break; case STRESS_PENALTY: cout << "Stress penalty objective function." << endl; break; + case ENTROPY_GENERATION: cout << "Entropy generation objective function." << endl; break; } } else { @@ -8521,6 +8536,10 @@ CConfig::~CConfig() { delete [] nBlades; delete [] FreeStreamTurboNormal; + + delete [] EntropyGeneration; + delete [] TotalPressureLoss; + delete [] KineticEnergyLoss; } /*--- Input is the filename base, output is the completed filename. ---*/ @@ -8664,6 +8683,9 @@ string CConfig::GetObjFunc_Extension(string val_filename) const { case TOPOL_DISCRETENESS: AdjExt = "_topdisc"; break; case TOPOL_COMPLIANCE: AdjExt = "_topcomp"; break; case STRESS_PENALTY: AdjExt = "_stress"; break; + case ENTROPY_GENERATION: AdjExt = "_entg"; break; + case TOTAL_PRESSURE_LOSS: AdjExt = "_tot_press_loss"; break; + case KINETIC_ENERGY_LOSS: AdjExt = "_kin_en_loss"; break; } } else{ @@ -9942,6 +9964,23 @@ short CConfig::FindInterfaceMarker(unsigned short iInterface) const { return -1; } +short CConfig::FindMixingPlaneInterfaceMarker(unsigned short nMarker) const { + short mark; + for (auto iMarker = 0; iMarker < nMarker; iMarker++){ + /*--- If the tag GetMarker_All_MixingPlaneInterface equals the index we are looping at ---*/ + if (GetMarker_All_MixingPlaneInterface(iMarker)){ + /*--- We have identified the local index of the marker ---*/ + /*--- Store the identifier for the marker ---*/ + mark = iMarker; + /*--- Exit the for loop: we have found the local index for Mixing-Plane interface ---*/ + return mark; + } + /*--- If the tag hasn't matched any tag within the donor markers ---*/ + mark = -1; + } + return mark; +} + void CConfig::Tick(double *val_start_time) { #ifdef PROFILE diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 1ba13be36c7..581512daf4d 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -50,6 +50,8 @@ #include "../../include/geometry/primal_grid/CPrism.hpp" #include "../../include/geometry/primal_grid/CVertexMPI.hpp" +#include "../../../Common/include/tracy_structure.hpp" + #include #include #include @@ -6326,6 +6328,7 @@ void CPhysicalGeometry::GatherInOutAverageValues(CConfig* config, bool allocate) void CPhysicalGeometry::SetAvgTurboGeoValues(const CConfig* donor_config, CGeometry* donor_geometry, unsigned short donorZone) { + unsigned short iSpan; unsigned short nSpanMaxAllZones = donor_config->GetnSpanMaxAllZones(); diff --git a/Common/src/interface_interpolation/CInterpolatorFactory.cpp b/Common/src/interface_interpolation/CInterpolatorFactory.cpp index f1e1490b8d9..fc816c8efc9 100644 --- a/Common/src/interface_interpolation/CInterpolatorFactory.cpp +++ b/Common/src/interface_interpolation/CInterpolatorFactory.cpp @@ -32,11 +32,12 @@ #include "../../include/interface_interpolation/CNearestNeighbor.hpp" #include "../../include/interface_interpolation/CRadialBasisFunction.hpp" #include "../../include/interface_interpolation/CSlidingMesh.hpp" +#include "../../include/interface_interpolation/CMixingPlane.hpp" namespace CInterpolatorFactory { CInterpolator* CreateInterpolator(CGeometry**** geometry_container, const CConfig* const* config, const CInterpolator* transpInterpolator, unsigned iZone, unsigned jZone, - bool verbose) { + bool mixing_plane, bool verbose) { CInterpolator* interpolator = nullptr; /*--- Only print information on master node. ---*/ @@ -47,36 +48,41 @@ CInterpolator* CreateInterpolator(CGeometry**** geometry_container, const CConfi if (verbose) cout << " Setting coupling "; - /*--- Conservative interpolation is not applicable to the sliding - * mesh approach so that case is handled first. Then we either - * return a CMirror if the target requires conservative inter- - * polation, or the type of interpolator defined by "type". ---*/ + if (mixing_plane) { + if (verbose) cout << "using a mixing plane interpolation." << endl; + interpolator = new CMixingPlane(geometry_container, config, iZone, jZone); + } else { // Really awful thing to do + /*--- Conservative interpolation is not applicable to the sliding + * mesh approach so that case is handled first. Then we either + * return a CMirror if the target requires conservative inter- + * polation, or the type of interpolator defined by "type". ---*/ - if (type == INTERFACE_INTERPOLATOR::WEIGHTED_AVERAGE) { - if (verbose) cout << "using a sliding mesh approach." << endl; - interpolator = new CSlidingMesh(geometry_container, config, iZone, jZone); - } else if (config[jZone]->GetConservativeInterpolation()) { - if (verbose) cout << "using the mirror approach, \"transposing\" coefficients from opposite mesh." << endl; - interpolator = new CMirror(geometry_container, config, transpInterpolator, iZone, jZone); - } else { - switch (type) { - case INTERFACE_INTERPOLATOR::ISOPARAMETRIC: - if (verbose) cout << "using the isoparametric approach." << endl; - interpolator = new CIsoparametric(geometry_container, config, iZone, jZone); - break; + if (type == INTERFACE_INTERPOLATOR::WEIGHTED_AVERAGE) { + if (verbose) cout << "using a sliding mesh approach." << endl; + interpolator = new CSlidingMesh(geometry_container, config, iZone, jZone); + } else if (config[jZone]->GetConservativeInterpolation()) { + if (verbose) cout << "using the mirror approach, \"transposing\" coefficients from opposite mesh." << endl; + interpolator = new CMirror(geometry_container, config, transpInterpolator, iZone, jZone); + } else { + switch (type) { + case INTERFACE_INTERPOLATOR::ISOPARAMETRIC: + if (verbose) cout << "using the isoparametric approach." << endl; + interpolator = new CIsoparametric(geometry_container, config, iZone, jZone); + break; - case INTERFACE_INTERPOLATOR::NEAREST_NEIGHBOR: - if (verbose) cout << "using a nearest neighbor approach." << endl; - interpolator = new CNearestNeighbor(geometry_container, config, iZone, jZone); - break; + case INTERFACE_INTERPOLATOR::NEAREST_NEIGHBOR: + if (verbose) cout << "using a nearest neighbor approach." << endl; + interpolator = new CNearestNeighbor(geometry_container, config, iZone, jZone); + break; - case INTERFACE_INTERPOLATOR::RADIAL_BASIS_FUNCTION: - if (verbose) cout << "using a radial basis function approach." << endl; - interpolator = new CRadialBasisFunction(geometry_container, config, iZone, jZone); - break; + case INTERFACE_INTERPOLATOR::RADIAL_BASIS_FUNCTION: + if (verbose) cout << "using a radial basis function approach." << endl; + interpolator = new CRadialBasisFunction(geometry_container, config, iZone, jZone); + break; - default: - SU2_MPI::Error("Unknown type of interpolation.", CURRENT_FUNCTION); + default: + SU2_MPI::Error("Unknown type of interpolation.", CURRENT_FUNCTION); + } } } diff --git a/Common/src/interface_interpolation/CMixingPlane.cpp b/Common/src/interface_interpolation/CMixingPlane.cpp new file mode 100644 index 00000000000..cbae8e6987a --- /dev/null +++ b/Common/src/interface_interpolation/CMixingPlane.cpp @@ -0,0 +1,161 @@ +/*! + * \file CMixingPlane.cpp + * \brief Implementation of mixing plane interpolation methods. + * \author J. Kelly + * \version 8.3.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/interface_interpolation/CMixingPlane.hpp" +#include "../../include/CConfig.hpp" +#include "../../include/geometry/CGeometry.hpp" +#include "../../include/toolboxes/geometry_toolbox.hpp" + +CMixingPlane::CMixingPlane(CGeometry**** geometry_container, const CConfig* const* config, unsigned int iZone, + unsigned int jZone) + : CInterpolator(geometry_container, config, iZone, jZone) { + SetTransferCoeff(config); +} + +void CMixingPlane::SetTransferCoeff(const CConfig* const* config) { + const auto nMarkerInt = config[donorZone]->GetnMarker_MixingPlaneInterface() / 2; + const auto nDim = donor_geometry->GetnDim(); + + const auto donor_config = config[donorZone]; + const auto target_config = config[targetZone]; + + const auto nMarkerDonor = donor_geometry->GetnMarker(); // Number of markers on the interfacce + const auto nMarkerTarget = target_geometry->GetnMarker(); + //TODO turbo this approach only works if all the turboamchinery marker + // of all zones have the same amount of span wise sections. + //TODO turbo initialization needed for the MPI routine should be place somewhere else. + auto nSpanDonor = donor_config->GetnSpanWiseSections(); + auto nSpanTarget = target_config->GetnSpanWiseSections(); + + targetSpans.resize(config[donorZone]->GetnMarker_MixingPlaneInterface()); + + /*--- On the donor side ---*/ + for (auto iMarkerInt = 0u; iMarkerInt < nMarkerInt; iMarkerInt++){ + int markDonor = -1, markTarget = -1; + unsigned short donorFlag = 0, targetFlag = 0; + + markDonor = donor_config->FindMixingPlaneInterfaceMarker(donor_config->GetnMarker_All()); + donorFlag = (markDonor != -1) ? donor_config->GetMarker_All_MixingPlaneInterface(markDonor) : -1; + +#ifdef HAVE_MPI + auto buffMarkerDonor = new int[size]; + auto buffDonorFlag = new int[size]; + for (int iSize=0; iSize= 0.0) { + markDonor = buffMarkerDonor[iSize]; + donorFlag = buffDonorFlag[iSize]; + break; + } + } + delete [] buffMarkerDonor; + delete [] buffDonorFlag; +#endif + + markTarget = target_config->FindMixingPlaneInterfaceMarker(target_config->GetnMarker_All()); + targetFlag = (markTarget != -1) ? target_config->GetMarker_All_MixingPlaneInterface(markTarget) : -1; + + if (markTarget == -1 || markDonor == -1) continue; + + nSpanDonor = donor_config->GetnSpanWiseSections(); + nSpanTarget = target_config->GetnSpanWiseSections(); + + targetSpans[iMarkerInt].resize(nSpanTarget+1); + + const auto spanValuesDonor = donor_geometry->GetSpanWiseValue(donorFlag); + const auto spanValuesTarget = target_geometry->GetSpanWiseValue(targetFlag); + + /*--- Interpolation at hub, shroud & 1D values ---*/ + targetSpans[iMarkerInt][0].donorSpan = 0; + targetSpans[iMarkerInt][0].coefficient = 0.0; + if (nDim > 2) { + targetSpans[iMarkerInt][nSpanTarget-2].donorSpan = nSpanDonor-2; + targetSpans[iMarkerInt][nSpanTarget-2].coefficient = 0.0; + targetSpans[iMarkerInt][nSpanTarget-1].donorSpan = nSpanDonor-1; + targetSpans[iMarkerInt][nSpanTarget-1].coefficient = 0.0; + } + + for(auto iSpanTarget = 1; iSpanTarget < nSpanTarget - 2; iSpanTarget++){ + auto &targetSpan = targetSpans[iMarkerInt][iSpanTarget]; + + auto tSpan = 0; // Nearest donor span index + auto kSpan = 0; // Lower bound donor span for interpolation + su2double coeff = 0.0; // Interpolation coefficient + su2double minDist = 10E+06; + + switch(donor_config->GetKind_MixingPlaneInterface()){ + case MATCHING: + targetSpan.donorSpan = iSpanTarget; + targetSpan.coefficient = 0.0; + break; + + case NEAREST_SPAN: + // Find the nearest donor span + for (auto iSpanDonor = 0; iSpanDonor < nSpanDonor - 1; iSpanDonor++) { + const auto dist = abs(spanValuesTarget[iSpanTarget] - spanValuesDonor[iSpanDonor]); + if(dist < minDist){ + minDist = dist; + tSpan = iSpanDonor; + } + } + targetSpan.donorSpan = tSpan; + targetSpan.coefficient = 0.0; + break; + + case LINEAR_INTERPOLATION: + // Find the donor span interval that brackets the target span + for (auto iSpanDonor = 1; iSpanDonor < nSpanDonor - 2; iSpanDonor++) { + if(spanValuesTarget[iSpanTarget] >= spanValuesDonor[iSpanDonor] && + spanValuesTarget[iSpanTarget] <= spanValuesDonor[iSpanDonor + 1]){ + kSpan = iSpanDonor; + break; + } + } + // Calculate interpolation coefficient + coeff = (spanValuesTarget[iSpanTarget] - spanValuesDonor[kSpan]) / + (spanValuesDonor[kSpan + 1] - spanValuesDonor[kSpan]); + targetSpan.donorSpan = kSpan; + targetSpan.coefficient = coeff; + break; + + default: + SU2_MPI::Error("MixingPlane interface option not implemented yet", CURRENT_FUNCTION); + break; + } + } + } +} \ No newline at end of file diff --git a/Common/src/interface_interpolation/meson.build b/Common/src/interface_interpolation/meson.build index 8624b3ad4fa..6cd185ef873 100644 --- a/Common/src/interface_interpolation/meson.build +++ b/Common/src/interface_interpolation/meson.build @@ -4,4 +4,5 @@ common_src += files(['CInterpolatorFactory.cpp', 'CSlidingMesh.cpp', 'CIsoparametric.cpp', 'CNearestNeighbor.cpp', - 'CRadialBasisFunction.cpp']) + 'CRadialBasisFunction.cpp', + 'CMixingPlane.cpp']) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index fc7094dd29a..7b4eb794616 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -79,6 +79,8 @@ class CDriver : public CDriverBase { CInterface*** interface_container; /*!< \brief Definition of the interface of information and physics. */ bool dry_run; /*!< \brief Flag if SU2_CFD was started as dry-run via "SU2_CFD -d .cfg" */ + //std::shared_ptr TurbomachineryStagePerformance; /*!< \brief turbo stage performance calculator. */ + public: /*! * \brief Constructor of the class. @@ -198,13 +200,11 @@ class CDriver : public CDriverBase { * \param[in] config - Definition of the particular problem. * \param[in] solver - Container vector with all the solutions. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] interface_types - Type of coupling between the distinct (physical) zones. * \param[in] interface - Class defining the physical transfer of information. * \param[in] interpolation - Object defining the interpolation. */ void InitializeInterface(CConfig** config, CSolver***** solver, CGeometry**** geometry, - unsigned short** interface_types, CInterface*** interface, - vector>>& interpolation); + CInterface*** interface, vector>>& interpolation); /*! * \brief Definition and allocation of all solver classes. @@ -290,17 +290,14 @@ class CDriver : public CDriverBase { * \param[in] geometry - Geometrical definition of the problem. * \param[in] solver - Container vector with all the solutions. * \param[in] interface - Class defining the physical transfer of information. + * \param[in] iteration - Class defining the iteration strcuture. * \param[in] dummy - Definition of dummy driver */ void PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, CSolver***** solver, - CInterface*** interface, bool dummy); + CInterface*** interface, CIteration*** iteration, bool dummy); - /*! - * \brief Ramp some simulation settings for turbomachinery problems. - * \param[in] iter - Iteration for the ramp (can be outer or time depending on type of simulation). - * \note TODO This is not compatible with inner iterations because they are delegated to the iteration class. - */ - void RampTurbomachineryValues(unsigned long iter); + void PreprocessTurboVertex(CConfig** config, CGeometry**** geometry, CSolver***** solver, + CInterface*** interface, CIteration*** iteration, bool dummy); /*! * \brief A virtual member. diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp index 97772527cf3..1f3377a58c2 100644 --- a/SU2_CFD/include/drivers/CDriverBase.hpp +++ b/SU2_CFD/include/drivers/CDriverBase.hpp @@ -60,8 +60,7 @@ class CDriverBase { nZone, /*!< \brief Total number of zones in the problem. */ nDim, /*!< \brief Number of dimensions. */ iInst, /*!< \brief Iterator on instance levels. */ - *nInst, /*!< \brief Total number of instances in the problem (per zone). */ - **interface_types; /*!< \brief Type of coupling between the distinct (physical) zones. */ + *nInst; /*!< \brief Total number of instances in the problem (per zone). */ CConfig* driver_config = nullptr; /*!< \brief Definition of the driver configuration. */ COutput* driver_output = nullptr; /*!< \brief Definition of the driver output. */ diff --git a/SU2_CFD/include/drivers/CMultizoneDriver.hpp b/SU2_CFD/include/drivers/CMultizoneDriver.hpp index 414a13243af..64c17c6b632 100644 --- a/SU2_CFD/include/drivers/CMultizoneDriver.hpp +++ b/SU2_CFD/include/drivers/CMultizoneDriver.hpp @@ -83,15 +83,6 @@ class CMultizoneDriver : public CDriver { */ bool TransferData(unsigned short donorZone, unsigned short targetZone); - - /*! - * \brief Transfer the local turboperfomance quantities (for each blade row) from all the donorZones to the - * targetZone (ZONE_0). - * \note IMPORTANT: This approach of multi-zone performances rely upon the fact that turbomachinery markers follow - * the natural (stator-rotor) development of the real machine. - */ - void SetTurboPerformance(); - /*! * \brief Check the convergence at the outer level. */ diff --git a/SU2_CFD/include/interfaces/CInterface.hpp b/SU2_CFD/include/interfaces/CInterface.hpp index 65876d99eeb..8a741b2c7cb 100644 --- a/SU2_CFD/include/interfaces/CInterface.hpp +++ b/SU2_CFD/include/interfaces/CInterface.hpp @@ -65,6 +65,8 @@ class CInterface { su2double *Target_Variable = nullptr; bool valAggregated = false; + unsigned short InterfaceType; /*!< \brief The type of interface. */ + /*--- Mixing Plane interface variable ---*/ su2double *SpanValueCoeffTarget = nullptr; unsigned short *SpanLevelDonor = nullptr; @@ -146,6 +148,12 @@ class CInterface { for (auto iVar = 0u; iVar < nVar; iVar++) Target_Variable[iVar] += donorCoeff * bcastVariable[iVar]; } + inline virtual void RecoverTarget_Variable(const vector bcastVariable, unsigned long idx) { } + + inline virtual void RecoverTarget_Variable(const vector bcastVariable, unsigned long idx, su2double donorCoeff) { } + + + /*! * \brief A virtual member. * \param[in] target_solution - Solution from the target mesh. @@ -185,6 +193,21 @@ class CInterface { inline virtual void SetAverageValues(CSolver *donor_solution, CSolver *target_solution, unsigned short donorZone) { } + /*! + * \brief Interpolate data and broadcast it into all processors, for nonmatching meshes. + * \param[in] interpolator - Object defining the interpolation. + * \param[in] donor_solution - Solution from the donor mesh. + * \param[in] target_solution - Solution from the target mesh. + * \param[in] donor_geometry - Geometry of the donor mesh. + * \param[in] target_geometry - Geometry of the target mesh. + * \param[in] donor_config - Definition of the problem at the donor mesh. + * \param[in] target_config - Definition of the problem at the target mesh. + */ + inline virtual void BroadcastData_MixingPlane(const CInterpolator& interpolator, + CSolver *donor_solution, CSolver *target_solution, + CGeometry *donor_geometry, CGeometry *target_geometry, + const CConfig *donor_config, const CConfig *target_config) { }; + /*! * \brief Transfer pre-processing for the mixing plane inteface. * \param[in] donor_geometry - Geometry of the donor mesh. @@ -224,4 +247,15 @@ class CInterface { * \param[in] val_contact_resistance - Contact resistance value in m^2/W */ inline virtual void SetContactResistance(su2double val_contact_resistance) {}; + + /*! + * \brief Set the type of an interface + * \param[in] interface_type - The type of interface + */ + void SetInterfaceType(unsigned short interface_type) { InterfaceType = interface_type; } + + /*! + * \brief Get the type of an interface + */ + unsigned short GetInterfaceType(void) const { return InterfaceType; } }; diff --git a/SU2_CFD/include/interfaces/cfd/CMixingPlaneInterface.hpp b/SU2_CFD/include/interfaces/cfd/CMixingPlaneInterface.hpp index 2abd9be5d60..5ce053e3ba7 100644 --- a/SU2_CFD/include/interfaces/cfd/CMixingPlaneInterface.hpp +++ b/SU2_CFD/include/interfaces/cfd/CMixingPlaneInterface.hpp @@ -36,6 +36,7 @@ */ class CMixingPlaneInterface : public CInterface { public: + unsigned int nMixingVars; /*! * \overload * \param[in] val_nVar - Number of variables that need to be transferred. @@ -49,6 +50,21 @@ class CMixingPlaneInterface : public CInterface { */ void SetSpanWiseLevels(const CConfig *donor_config, const CConfig *target_config) override; + /*! + * \brief Interpolate data and broadcast it into all processors, for nonmatching meshes. + * \param[in] interpolator - Object defining the interpolation. + * \param[in] donor_solution - Solution from the donor mesh. + * \param[in] target_solution - Solution from the target mesh. + * \param[in] donor_geometry - Geometry of the donor mesh. + * \param[in] target_geometry - Geometry of the target mesh. + * \param[in] donor_config - Definition of the problem at the donor mesh. + * \param[in] target_config - Definition of the problem at the target mesh. + */ + void BroadcastData_MixingPlane(const CInterpolator& interpolator, + CSolver *donor_solution, CSolver *target_solution, + CGeometry *donor_geometry, CGeometry *target_geometry, + const CConfig *donor_config, const CConfig *target_config) override; + /*! * \brief Retrieve the variable that will be sent from donor mesh to target mesh. * \param[in] donor_solution - Solution from the donor mesh. @@ -61,6 +77,9 @@ class CMixingPlaneInterface : public CInterface { void GetDonor_Variable(CSolver *donor_solution, CGeometry *donor_geometry, const CConfig *donor_config, unsigned long Marker_Donor, unsigned long val_Span, unsigned long Point_Donor) override; + void InitializeTarget_Variable(CSolver *target_solution, unsigned long Marker_Target, + unsigned long Span_Target, unsigned short nDonorPoints) override; + /*! * \brief Set the variable that has been received from the target mesh into the target mesh. * \param[in] target_solution - Solution from the target mesh. @@ -73,6 +92,18 @@ class CMixingPlaneInterface : public CInterface { void SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config, unsigned long Marker_Target, unsigned long val_Span, unsigned long Point_Target) override; + inline void RecoverTarget_Variable(const vector bcastVariable, unsigned long iSpan) override{ + for (auto iVar = 0u; iVar < nMixingVars; iVar++) { + Target_Variable[iVar] = bcastVariable[iSpan * nMixingVars + iVar]; + } + } + + inline void RecoverTarget_Variable(const vector bcastVariable, unsigned long iSpan, su2double donorCoeff) override { + for (auto iVar = 0u; iVar < nMixingVars; iVar++) { + Target_Variable[iVar] = (1 - donorCoeff)*bcastVariable[iSpan * nMixingVars + iVar] + donorCoeff * bcastVariable[(iSpan + 1) * nMixingVars + iVar]; + } + } + /*! * \brief Store all the turboperformance in the solver in ZONE_0. * \param[in] donor_solution - Solution from the donor mesh. diff --git a/SU2_CFD/include/iteration/CFluidIteration.hpp b/SU2_CFD/include/iteration/CFluidIteration.hpp index 883540f65ac..ca472eda61e 100644 --- a/SU2_CFD/include/iteration/CFluidIteration.hpp +++ b/SU2_CFD/include/iteration/CFluidIteration.hpp @@ -115,12 +115,12 @@ class CFluidIteration : public CIteration { * \param[in] iZone - The current zone * \param[in] ramp_flag - Flag indicating type of ramp (grid or boundary) */ - void UpdateRamp(CGeometry**** geometry_container, CConfig** config_container, unsigned long iter, unsigned short iZone, RAMP_TYPE ramp_flag); + void UpdateRamps(CGeometry**** geometry_container, CConfig** config_container, unsigned long iter, unsigned short iZone, RAMP_TYPE ramp_flag); /*! * \brief Computes turboperformance. */ - void ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container, unsigned long ExtIter); + // void ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container); /*! * \brief Postprocesses the fluid system before heading to another physics system or the next iteration. diff --git a/SU2_CFD/include/iteration/CIteration.hpp b/SU2_CFD/include/iteration/CIteration.hpp index b9d68919ba6..d5af0faf5df 100644 --- a/SU2_CFD/include/iteration/CIteration.hpp +++ b/SU2_CFD/include/iteration/CIteration.hpp @@ -60,9 +60,7 @@ class CIteration { su2double StartTime{0.0}, /*!< \brief Tracking wall time. */ StopTime{0.0}, UsedTime{0.0}; - std::shared_ptr TurbomachineryPerformance; /*!< \brief turbo performance calculator. */ std::shared_ptr TurbomachineryStagePerformance; /*!< \brief turbo stage performance calculator. */ - public: /*! * \brief Constructor of the class. @@ -295,4 +293,24 @@ class CIteration { virtual void RegisterOutput(CSolver***** solver, CGeometry**** geometry, CConfig** config, unsigned short iZone, unsigned short iInst) {} + + /*! + * \brief Computes turboperformance. + */ + void ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container); + + /*! + * \brief Initialises turboperformance classes. + */ + void InitTurboPerformance(CGeometry *geometry, CConfig** config, CFluidModel *fluid, unsigned short val_iZone); + + inline vector> GetBladesPerformanceVector(CSolver***** solver, unsigned short nBladeRow){ + vector> bladePerformances; + bladePerformances.reserve(nBladeRow); + for (auto iBladeRow = 0u; iBladeRow < nBladeRow; iBladeRow++) { + auto const turboPerf = solver[iBladeRow][INST_0][MESH_0][FLOW_SOL]->GetTurboBladePerformance(); + bladePerformances.push_back(turboPerf); + } + return bladePerformances; + } }; diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp index 13ee01f7eda..d79108704b1 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp @@ -148,6 +148,8 @@ class CAvgGrad_Scalar : public CNumerics { FinishResidualCalc(config); AD::SetPreaccOut(Flux, nVar); + AD::SetPreaccOut(Jacobian_i, nVar, nVar); + AD::SetPreaccOut(Jacobian_j, nVar, nVar); AD::EndPreacc(); return ResidualType<>(Flux, Jacobian_i, Jacobian_j); diff --git a/SU2_CFD/include/output/CFlowCompOutput.hpp b/SU2_CFD/include/output/CFlowCompOutput.hpp index a0c2a368db5..f257a981048 100644 --- a/SU2_CFD/include/output/CFlowCompOutput.hpp +++ b/SU2_CFD/include/output/CFlowCompOutput.hpp @@ -53,6 +53,8 @@ class CFlowCompOutput final: public CFlowOutput { */ void LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) override; + void LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned short iZone); + /*! * \brief Set the available volume output fields * \param[in] config - Definition of the particular problem. @@ -74,6 +76,8 @@ class CFlowCompOutput final: public CFlowOutput { */ void SetHistoryOutputFields(CConfig *config) override; + void SetTurbomachineryObjectiveFunctions(CSolver *solver, CConfig *config); + /*! * \brief Check whether the base values for relative residuals should be initialized * \param[in] config - Definition of the particular problem. @@ -101,7 +105,7 @@ class CFlowCompOutput final: public CFlowOutput { * \param[in] OuterIter - Index of current outer iteration * \param[in] InnerIter - Index of current inner iteration */ - void SetTurboPerformance_Output(std::shared_ptr TurboPerf, CConfig *config, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter) override; + void SetTurboPerformance_Output(std::vector> TurboBladePerfs, CConfig *config, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter) override; /*! * \brief Sets the multizone turboperformacne screen output @@ -109,7 +113,7 @@ class CFlowCompOutput final: public CFlowOutput { * \param[in] TurboPerf - Turboperformance class * \param[in] config - Definition of the particular problem */ - void SetTurboMultiZonePerformance_Output(std::shared_ptr TurboStagePerf, std::shared_ptr TurboPerf, CConfig *config) override; + void SetTurboMultiZonePerformance_Output(std::shared_ptr TurboStagePerf, std::vector> TurboBladePerfs, CConfig *config) override; /*! * \brief Loads the turboperformacne history data @@ -117,7 +121,7 @@ class CFlowCompOutput final: public CFlowOutput { * \param[in] TurboPerf - Turboperformance class * \param[in] config - Definition of the particular problem */ - void LoadTurboHistoryData(std::shared_ptr TurboStagePerf, std::shared_ptr TurboPerf, CConfig *config) override; + void LoadTurboHistoryData(std::shared_ptr TurboStagePerf, std::vector> TurboBladePerfs, CConfig *config) override; /*! * \brief Write the kinematic and thermodynamic variables at each spanwise division @@ -126,6 +130,6 @@ class CFlowCompOutput final: public CFlowOutput { * \param[in] config - Descripiton of the particular problem * \param[in] val_iZone - Idientifier of current zone */ - void WriteTurboSpanwisePerformance(std::shared_ptr TurboPerf, CGeometry *geometry, CConfig **config, + void WriteTurboSpanwisePerformance(std::vector> TurboBladePerfs, CGeometry *geometry, CConfig **config, unsigned short val_iZone) override; }; diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index 3369b7e7a92..c8e5b3b7369 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -406,7 +406,7 @@ class COutput { */ void SetHistoryOutput(CGeometry ****geometry, CSolver *****solver_container, CConfig **config, std::shared_ptr TurboStagePerf, - std::shared_ptr TurboPerf, unsigned short val_iZone, + std::vector> TurboBladePerfs, unsigned short val_iZone, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter, unsigned short val_iInst); /*! @@ -977,7 +977,7 @@ class COutput { * \param[in] OuterIter - Index of current outer iteration * \param[in] InnerIter - Index of current inner iteration */ - inline virtual void SetTurboPerformance_Output(std::shared_ptr TurboPerf, CConfig *config, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter) {} + inline virtual void SetTurboPerformance_Output(std::vector> TurboBladePerfs, CConfig *config, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter) {} /*! * \brief Sets the multizone turboperformacne screen output @@ -985,7 +985,7 @@ class COutput { * \param[in] TurboPerf - Turboperformance class * \param[in] config - Definition of the particular problem */ - inline virtual void SetTurboMultiZonePerformance_Output(std::shared_ptr TurboStagePerf, std::shared_ptr TurboPerf, CConfig *config) {} + inline virtual void SetTurboMultiZonePerformance_Output(std::shared_ptr TurboStagePerf, std::vector> TurboPerf, CConfig *config) {} /*! * \brief Loads the turboperformacne history data @@ -993,7 +993,7 @@ class COutput { * \param[in] TurboPerf - Turboperformance class * \param[in] config - Definition of the particular problem */ - inline virtual void LoadTurboHistoryData(std::shared_ptr TurboStagePerf, std::shared_ptr TurboPerf, CConfig *config) {} + inline virtual void LoadTurboHistoryData(std::shared_ptr TurboStagePerf, std::vector> TurboPerf, CConfig *config) {} /*! * \brief Write the kinematic and thermodynamic variables at each spanwise division @@ -1002,7 +1002,7 @@ class COutput { * \param[in] config - Descripiton of the particular problem * \param[in] val_iZone - Idientifier of current zone */ - inline virtual void WriteTurboSpanwisePerformance(std::shared_ptr TurboPerf, CGeometry *geometry, CConfig **config, + inline virtual void WriteTurboSpanwisePerformance(std::vector> TurboBladePerfs, CGeometry *geometry, CConfig **config, unsigned short val_iZone) {}; /*! diff --git a/SU2_CFD/include/output/CTurboOutput.hpp b/SU2_CFD/include/output/CTurboOutput.hpp index 4fe06e70987..6b7e35bb83b 100644 --- a/SU2_CFD/include/output/CTurboOutput.hpp +++ b/SU2_CFD/include/output/CTurboOutput.hpp @@ -96,6 +96,11 @@ class CTurbomachineryState { CTurbomachineryState(unsigned short nDim, su2double area, su2double radius); + inline void SetZeroValues() { + Density = Pressure = Entropy = Enthalpy = Temperature = TotalTemperature = TotalPressure = TotalEnthalpy = 0.0; + AbsFlowAngle = FlowAngle = MassFlow = Rothalpy = TotalRelPressure = 0.0; + }; + void ComputeState(CFluidModel& fluidModel, const CTurbomachineryPrimitiveState& primitiveState); const su2double& GetDensity() const { return Density; } @@ -247,16 +252,16 @@ class CTurbomachineryStagePerformance { */ class CTurboOutput { private: - vector>> BladesPerformances; + vector> BladesPerformances; static void ComputePerBlade(vector> const bladePerformances, vector const bladePrimitives); static void ComputePerSpan(shared_ptr const spanPerformances, const CTurbomachineryCombinedPrimitiveStates& spanPrimitives); public: - CTurboOutput(CConfig** config, const CGeometry& geometry, CFluidModel& fluidModel); + CTurboOutput(CConfig** config, const CGeometry& geometry, CFluidModel& fluidModel, unsigned short iBladeRow); - const vector>>& GetBladesPerformances() const { return BladesPerformances; } + const vector>& GetBladesPerformances() const { return BladesPerformances; } - void ComputeTurbomachineryPerformance(vector> const primitives); + void ComputeTurbomachineryPerformance(vector const primitives, unsigned short iBladeRow); }; \ No newline at end of file diff --git a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp index 211613db2c7..cbfe7ec8dd1 100644 --- a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp @@ -106,6 +106,15 @@ class CDiscAdjSolver final : public CSolver { */ void RegisterOutput(CGeometry *geometry, CConfig *config) override; + /*! + * \brief performs the preprocessing of the adjoint AD-based solver. + * Registers the normals at turbomachinery boundaries of interest on the tape + * Called while tape is active. + * \param[in] geometry - the geometrical definition of the problem + * \param[in] config - the particular config + */ + void Register_VertexNormals(CGeometry *geometry, CConfig *config, bool input) override; + /*! * \brief Sets the adjoint values of the output of the flow (+turb.) iteration * before evaluation of the tape. diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index e3f029a3753..97e8a55268c 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -29,6 +29,7 @@ #include "CFVMFlowSolverBase.hpp" #include "../variables/CEulerVariable.hpp" +#include "../output/CTurboOutput.hpp" /*! * \class CEulerSolver @@ -123,20 +124,14 @@ class CEulerSolver : public CFVMFlowSolverBase AverageVelocity; vector AverageTurboVelocity; vector OldAverageTurboVelocity; - vector ExtAverageTurboVelocity; su2activematrix AveragePressure; su2activematrix OldAveragePressure; su2activematrix RadialEquilibriumPressure; - su2activematrix ExtAveragePressure; su2activematrix AverageDensity; su2activematrix OldAverageDensity; - su2activematrix ExtAverageDensity; su2activematrix AverageNu; su2activematrix AverageKine; su2activematrix AverageOmega; - su2activematrix ExtAverageNu; - su2activematrix ExtAverageKine; - su2activematrix ExtAverageOmega; su2activevector AverageMassFlowRate; su2activematrix DensityIn; @@ -154,8 +149,78 @@ class CEulerSolver : public CFVMFlowSolverBase > > CkInflow, CkOutflow1, CkOutflow2; + su2double EntropyGeneration; + su2double TotalPressureLoss; + su2double KineticEnergyLoss; + /*--- End of Turbomachinery Solver Variables ---*/ + vector> MixingState; // vector of vector of pointers... inner dim alloc'd elsewhere (welcome, to the night zone) + vector> MixingStateNodes; + + /*! + * \brief Get the outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_span - vertex index + * \param[in] val_state - requested state component + // * \param[in] donor_index- index of the donor node to get + */ + inline su2double GetMixingState(unsigned short val_marker, + unsigned long val_span, + unsigned short val_state + // unsigned long donor_index + ) const final { + return MixingState[val_marker][val_span][val_state]; + } + + /*! + * \brief Allocates the final pointer of MixingState depending on how many donor spans donate to it. That number is stored in MixingStateSpans[val_marker][val_Span]. + * \param[in] val_marker - marker index + * \param[in] val_span - vertex index + */ + inline void SetMixingStateStructure(unsigned short val_marker, unsigned long val_span) final { + assert(val_marker < MixingState.size()); + assert(val_span < MixingState[val_marker].size()); + if( MixingState[val_marker][val_span] != nullptr ) delete [] MixingState[val_marker][val_span]; + + MixingState[val_marker][val_span] = new su2double[8]; + } + + /*! + * \brief Set the outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_span - vertex index + * \param[in] val_state - requested state component + // * \param[in] donor_index - index of the donor node to set + * \param[in] component - set value + */ + inline void SetMixingState(unsigned short val_marker, + unsigned long val_span, + unsigned short val_state, + // unsigned long donor_index, + su2double component) final { + MixingState[val_marker][val_span][val_state] = component; + } + + /*! + * \brief Set the number of outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_span - vertex index + * \param[in] value - number of outer states + */ + inline void SetnMixingStates(unsigned short val_marker, + unsigned long val_span, + int value) final { MixingStateNodes[val_marker][val_span] = value; } + + /*! + * \brief Get the number of outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_vertex - vertex index + */ + inline int GetnMixingStates(unsigned short val_marker, unsigned long val_span) const final { + return MixingStateNodes[val_marker][val_span]; + } + /*! * \brief Preprocessing actions common to the Euler and NS solvers. * \param[in] geometry - Geometrical definition of the problem. @@ -1036,7 +1101,7 @@ class CEulerSolver : public CFVMFlowSolverBase GetTurboBladePerformance() const final { return TurbomachineryPerformance; } + /*! * \brief it take a velocity in the cartesian reference of framework and transform into the turbomachinery frame of reference. * \param[in] cartesianVelocity - cartesian components of velocity vector. @@ -1229,105 +1327,6 @@ class CEulerSolver : public CFVMFlowSolverBaseval_marker. - */ - inline su2double GetExtAverageNu(unsigned short valMarker, unsigned short valSpan) const final { - return ExtAverageNu[valMarker][valSpan]; - } - - /*! - * \brief Provide the average density at the boundary of interest. - * \param[in] val_marker - bound marker. - * \return Value of the Average turbulent Kine on the surface val_marker. - */ - inline su2double GetExtAverageKine(unsigned short valMarker, unsigned short valSpan) const final { - return ExtAverageKine[valMarker][valSpan]; - } - - /*! - * \brief Provide the average density at the boundary of interest. - * \param[in] val_marker - bound marker. - * \return Value of the Average turbulent Omega on the surface val_marker. - */ - inline su2double GetExtAverageOmega(unsigned short valMarker, unsigned short valSpan) const final { - return ExtAverageOmega[valMarker][valSpan]; - } - - /*! - * \brief Set the external average density at the boundary of interest. - * \param[in] val_marker - bound marker. - * \param[in] val_Span - value of the Span. - * \param[in] valDensity - value to set. - */ - inline void SetExtAverageDensity(unsigned short valMarker, - unsigned short valSpan, - su2double valDensity) final { - ExtAverageDensity[valMarker][valSpan] = valDensity; - } - - /*! - * \brief Set the external average density at the boundary of interest. - * \param[in] val_marker - bound marker. - * \param[in] val_Span - value of the Span. - * \param[in] valPressure - value to set. - */ - inline void SetExtAveragePressure(unsigned short valMarker, - unsigned short valSpan, - su2double valPressure) final { - ExtAveragePressure[valMarker][valSpan] = valPressure; - } - - /*! - * \brief Set the external the average turbo velocity average at the boundary of interest. - * \param[in] val_marker - bound marker. - * \return Value of the Average Total Pressure on the surface val_marker. - */ - inline void SetExtAverageTurboVelocity(unsigned short valMarker, - unsigned short valSpan, - unsigned short valIndex, - su2double valTurboVelocity) final { - ExtAverageTurboVelocity[valMarker][valSpan][valIndex] = valTurboVelocity; - } - - /*! - * \brief Set the external average turbulent Nu at the boundary of interest. - * \param[in] val_marker - bound marker. - * \param[in] val_Span - value of the Span. - * \param[in] valNu - value to set. - */ - inline void SetExtAverageNu(unsigned short valMarker, - unsigned short valSpan, - su2double valNu) final { - ExtAverageNu[valMarker][valSpan] = valNu; - } - - /*! - * \brief Set the external average turbulent Kine at the boundary of interest. - * \param[in] val_marker - bound marker. - * \param[in] val_Span - value of the Span. - * \param[in] valKine - value to set. - */ - inline void SetExtAverageKine(unsigned short valMarker, - unsigned short valSpan, - su2double valKine) final { - ExtAverageKine[valMarker][valSpan] = valKine; - } - - /*! - * \brief Set the external average turbulent Omega at the boundary of interest. - * \param[in] val_marker - bound marker. - * \param[in] val_Span - value of the Span. - * \param[in] valOmega - value to set. - */ - inline void SetExtAverageOmega(unsigned short valMarker, - unsigned short valSpan, - su2double valOmega) final { - ExtAverageOmega[valMarker][valSpan] = valOmega; - } - /*! * \brief Provide the inlet density to check convergence of conservative mixing-plane. * \param[in] inMarkerTP - bound marker. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index f0b3a4c493d..ae4344ca46d 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -59,6 +59,8 @@ #include "../limiters/CLimiterDetails.hpp" #include "../variables/CVariable.hpp" +#include "../output/CTurboOutput.hpp" + #ifdef HAVE_LIBROM #include "librom.h" #endif @@ -151,6 +153,8 @@ class CSolver { vector VertexTraction; /*- Temporary, this will be moved to a new postprocessing structure once in place -*/ vector VertexTractionAdjoint; /*- Also temporary -*/ + std::shared_ptr TurbomachineryPerformance; /*!< \brief turbo performance calculator. */ + string SolverName; /*!< \brief Store the name of the solver for output purposes. */ /*! @@ -1105,6 +1109,12 @@ class CSolver { CConfig *config, unsigned short val_marker) { } + inline virtual void SetTurboObjectiveFunction(short unsigned int ObjFunc, su2double val) { } + + inline virtual su2double GetTurboObjectiveFunction(short unsigned int ObjFunc, int bladeRow) const { return 0.0; } + + inline virtual std::shared_ptr GetTurboBladePerformance() const { return std::shared_ptr(nullptr); } + /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. @@ -1286,6 +1296,55 @@ class CSolver { */ virtual void Impose_Fixed_Values(const CGeometry *geometry, const CConfig *config) { } + /*! + * \brief Get the outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_span - vertex index + * \param[in] val_state - requested state component + * \param[in] donor_span - index of the donor span to get + */ + inline virtual su2double GetMixingState(unsigned short val_marker, + unsigned long val_span, + unsigned short val_state) const { return 0; } + + /*! + * \brief Allocates the final pointer of MixingState depending on how many donor vertex donate to it. That number is stored in MixingStateNodes[val_marker][val_vertex]. + * \param[in] val_marker - marker index + * \param[in] val_span - span index + */ + inline virtual void SetMixingStateStructure(unsigned short val_marker, unsigned long val_vertex) { } + + /*! + * \brief Set the outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_span - vertex index + * \param[in] val_state - requested state component + * \param[in] donor_index - index of the donor node to set + * \param[in] component - set value + */ + inline virtual void SetMixingState(unsigned short val_marker, + unsigned long val_span, + unsigned short val_state, + // unsigned long donor_span, // Do I care about where it comes from? + su2double component) { } + + /*! + * \brief Set the number of outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_span - span index + * \param[in] value - number of outer states + */ + inline virtual void SetnMixingStates(unsigned short val_marker, + unsigned long val_span, + int value) { } + + /*! + * \brief Get the number of outer state for mixing plane interface nodes. + * \param[in] val_marker - marker index + * \param[in] val_span - vertex index + */ + inline virtual int GetnMixingStates(unsigned short val_marker, unsigned long val_span) const { return 0; } + /*! * \brief Get the outer state for fluid interface nodes. * \param[in] val_marker - marker index @@ -3528,6 +3587,14 @@ class CSolver { */ inline virtual void RegisterOutput(CGeometry *geometry_container, CConfig *config) { } + /*! + * \brief A vritual member + * \param[in] geometry - the geometrical definition of the problem + * \param[in] config - the particular config + * \param[in] input - Boolean whether In- or Output should be registered + */ + inline virtual void Register_VertexNormals(CGeometry *geometry, CConfig *config, bool input) { }; + /*! * \brief A virtual member. * \param[in] geometry - The geometrical definition of the problem. @@ -3767,7 +3834,7 @@ class CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - inline virtual void InitTurboContainers(CGeometry *geometry, CConfig *config) { } + inline virtual void InitTurboContainers(CGeometry *geometry, CConfig **config, unsigned short iZone) { } /*! * \brief Get Primal variables for turbo performance computation @@ -3807,6 +3874,8 @@ class CSolver { */ inline virtual void GatherInOutAverageValues(CConfig *config, CGeometry *geometry) { } + inline virtual void ComputeTurboBladePerformance(CGeometry* geometry, CConfig* config, unsigned short iBlade) { }; + /*! * \brief A virtual member. * \param[in] val_marker - bound marker. @@ -3856,82 +3925,6 @@ class CSolver { */ inline virtual su2double GetAverageOmega(unsigned short valMarker, unsigned short valSpan) const { return 0.0; } - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Nu on the surface val_marker. - */ - inline virtual su2double GetExtAverageNu(unsigned short valMarker, unsigned short valSpan) const { return 0.0; } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Kine on the surface val_marker. - */ - inline virtual su2double GetExtAverageKine(unsigned short valMarker, unsigned short valSpan) const { return 0.0; } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Omega on the surface val_marker. - */ - inline virtual su2double GetExtAverageOmega(unsigned short valMarker, unsigned short valSpan) const { return 0.0; } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Density on the surface val_marker. - */ - inline virtual void SetExtAverageDensity(unsigned short valMarker, - unsigned short valSpan, - su2double valDensity) { } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Pressure on the surface val_marker. - */ - inline virtual void SetExtAveragePressure(unsigned short valMarker, - unsigned short valSpan, - su2double valPressure) { } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Total Pressure on the surface val_marker. - */ - inline virtual void SetExtAverageTurboVelocity(unsigned short valMarker, - unsigned short valSpan, - unsigned short valIndex, - su2double valTurboVelocity) { } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Nu on the surface val_marker. - */ - inline virtual void SetExtAverageNu(unsigned short valMarker, - unsigned short valSpan, - su2double valNu) { } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Kine on the surface val_marker. - */ - inline virtual void SetExtAverageKine(unsigned short valMarker, - unsigned short valSpan, - su2double valKine) { } - - /*! - * \brief A virtual member. - * \param[in] val_marker - bound marker. - * \return Value of the Average Omega on the surface val_marker. - */ - inline virtual void SetExtAverageOmega(unsigned short valMarker, - unsigned short valSpan, - su2double valOmega) { } - /*! * \brief A virtual member. * \param[in] inMarkerTP - bound marker. diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index d5f9acbd0a9..ce30c850889 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -838,6 +838,15 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { solvers[HEAT_SOL]->Heat_Fluxes(geometry, solvers, config); } + if(config->GetBoolTurbomachinery()){ + solvers[FLOW_SOL]->TurboAverageProcess(solvers, geometry, config, INFLOW); + solvers[FLOW_SOL]->TurboAverageProcess(solvers, geometry, config, OUTFLOW); + + /*--- Gather Inflow and Outflow quantities on the Master Node to compute performance ---*/ + solvers[FLOW_SOL]->GatherInOutAverageValues(config, geometry); + solvers[FLOW_SOL]->ComputeTurboBladePerformance(geometry, config, iZone); + } + direct_output[iZone]->SetHistoryOutput(geometry, solvers, config); ObjFunc += solvers[FLOW_SOL]->GetTotal_ComboObj(); break; diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index 59eb6f4929a..db02acc1f92 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -33,6 +33,7 @@ #include "../../include/iteration/CTurboIteration.hpp" #include "../../../Common/include/toolboxes/CQuasiNewtonInvLeastSquares.hpp" + CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator) : CSinglezoneDriver(confFile, diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 7e17ed9c601..85d509697af 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -226,6 +226,13 @@ CDriverBase(confFile, val_nZone, MPICommunicator), StopCalc(false), fsi(false), CGeometry::ComputeWallDistance(config_container, geometry_container); } + if (config_container[ZONE_0]->GetBoolTurbomachinery()){ + if (rank == MASTER_NODE) + cout << endl <<"---------------------- Turbo-Vertex Preprocessing ---------------------" << endl; + + PreprocessTurboVertex(config_container, geometry_container, solver_container, interface_container, iteration_container, dummy_geo); + } + /*--- Definition of the interface and transfer conditions between different zones. ---*/ if (nZone > 1) { @@ -233,7 +240,7 @@ CDriverBase(confFile, val_nZone, MPICommunicator), StopCalc(false), fsi(false), cout << endl <<"------------------- Multizone Interface Preprocessing -------------------" << endl; InitializeInterface(config_container, solver_container, geometry_container, - interface_types, interface_container, interpolator_container); + interface_container, interpolator_container); } if (fsi) { @@ -249,12 +256,9 @@ CDriverBase(confFile, val_nZone, MPICommunicator), StopCalc(false), fsi(false), if (rank == MASTER_NODE) cout << endl <<"---------------------- Turbomachinery Preprocessing ---------------------" << endl; - PreprocessTurbomachinery(config_container, geometry_container, solver_container, interface_container, dummy_geo); - } else { - mixingplane = false; + PreprocessTurbomachinery(config_container, geometry_container, solver_container, interface_container, iteration_container, dummy_geo); } - PreprocessPythonInterface(config_container, geometry_container, solver_container); @@ -314,7 +318,6 @@ void CDriver::InitializeContainers(){ grid_movement = nullptr; FFDBox = nullptr; interface_container = nullptr; - interface_types = nullptr; nInst = nullptr; /*--- Definition and of the containers for all possible zones. ---*/ @@ -330,16 +333,12 @@ void CDriver::InitializeContainers(){ FFDBox = new CFreeFormDefBox**[nZone] (); interpolator_container.resize(nZone); interface_container = new CInterface**[nZone] (); - interface_types = new unsigned short*[nZone] (); output_container = new COutput*[nZone] (); nInst = new unsigned short[nZone] (); driver_config = nullptr; driver_output = nullptr; - for (iZone = 0; iZone < nZone; iZone++) { - interface_types[iZone] = new unsigned short[nZone]; - nInst[iZone] = 1; - } + for (iZone = 0; iZone < nZone; iZone++) nInst[iZone] = 1; } @@ -417,12 +416,6 @@ void CDriver::Finalize() { if (rank == MASTER_NODE) cout << "Deleted CInterface container." << endl; } - if (interface_types != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) - delete [] interface_types[iZone]; - delete [] interface_types; - } - for (iZone = 0; iZone < nZone; iZone++) { if (geometry_container[iZone] != nullptr) { for (iInst = 0; iInst < nInst[iZone]; iInst++){ @@ -1286,14 +1279,14 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) { - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); } else { - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); } } else if (menter_sst) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); } } @@ -2408,8 +2401,7 @@ void CDriver::PreprocessDynamicMesh(CConfig *config, CGeometry **geometry, CSolv } -void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeometry**** geometry, - unsigned short** interface_types, CInterface ***interface, +void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeometry**** geometry, CInterface ***interface, vector > >& interpolation) { /*--- Setup interpolation and transfer for all possible donor/target pairs. ---*/ @@ -2418,19 +2410,9 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet for (auto donor = 0u; donor < nZone; donor++) { - /*--- Aliases to make code less verbose. ---*/ - auto& interface_type = interface_types[donor][target]; - - if (donor == target) { - interface_type = ZONES_ARE_EQUAL; - continue; - } - interface_type = NO_TRANSFER; - /*--- If there is a common interface setup the interpolation and transfer. ---*/ - - if (!CInterpolator::CheckZonesInterface(config[donor], config[target])) { - interface_type = NO_COMMON_INTERFACE; + if (!CInterpolator::CheckZonesInterface(config[donor], config[target]) || donor == target) { + continue; } else { /*--- Begin the creation of the communication pattern among zones. ---*/ @@ -2439,8 +2421,11 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet /*--- Setup the interpolation. ---*/ - interpolation[donor][target] = unique_ptr(CInterpolatorFactory::CreateInterpolator( - geometry, config, interpolation[target][donor].get(), donor, target)); + if (!config[donor]->GetBoolTurbomachinery()) { + interpolation[donor][target] = unique_ptr(CInterpolatorFactory::CreateInterpolator( + geometry, config, interpolation[target][donor].get(), donor, target, false)); + if (rank == MASTER_NODE) cout << " Transferring "; + } /*--- The type of variables transferred depends on the donor/target physics. ---*/ @@ -2454,10 +2439,8 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet /*--- Initialize the appropriate transfer strategy. ---*/ - if (rank == MASTER_NODE) cout << " Transferring "; - if (fluid_donor && structural_target) { - interface_type = FLOW_TRACTION; + auto nConst = 2; bool conservative = config[target]->GetConservativeInterpolation(); if(!config[ZONE_0]->GetDiscrete_Adjoint()) { @@ -2472,7 +2455,6 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet SU2_MPI::Error("Mesh deformation was not correctly specified for the fluid/heat zone.\n" "Use DEFORM_MESH=YES, and setup MARKER_DEFORM_MESH=(...)", CURRENT_FUNCTION); } - interface_type = BOUNDARY_DISPLACEMENTS; if (!config[donor]->GetTime_Domain()) interface[donor][target] = new CDisplacementsInterface(nDim, 0); else interface[donor][target] = new CDisplacementsInterface(2*nDim, 0); if (rank == MASTER_NODE) cout << "boundary displacements from the structural solver." << endl; @@ -2483,15 +2465,19 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet auto interfaceIndex = donor+target; // Here we assume that the interfaces at each side are the same kind switch (config[donor]->GetKind_TurboInterface(interfaceIndex)) { case TURBO_INTERFACE_KIND::MIXING_PLANE: { - interface_type = MIXING_PLANE; + interpolation[donor][target] = unique_ptr(CInterpolatorFactory::CreateInterpolator( + geometry, config, interpolation[target][donor].get(), donor, target, true)); + if (rank == MASTER_NODE) cout << " Transferring "; auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnVar(); interface[donor][target] = new CMixingPlaneInterface(nVar, 0); if (rank == MASTER_NODE) cout << "using a mixing-plane interface from donor zone " << donor << " to target zone " << target << "." << endl; break; } case TURBO_INTERFACE_KIND::FROZEN_ROTOR: { + interpolation[donor][target] = unique_ptr(CInterpolatorFactory::CreateInterpolator( + geometry, config, interpolation[target][donor].get(), donor, target, false)); + if (rank == MASTER_NODE) cout << " Transferring "; auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnPrimVar(); - interface_type = SLIDING_INTERFACE; interface[donor][target] = new CSlidingInterface(nVar, 0); if (rank == MASTER_NODE) cout << "using a fluid interface interface from donor zone " << donor << " to target zone " << target << "." << endl; } @@ -2499,30 +2485,31 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet } else{ auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnPrimVar(); - interface_type = SLIDING_INTERFACE; interface[donor][target] = new CSlidingInterface(nVar, 0); if (rank == MASTER_NODE) cout << "sliding interface." << endl; } } else if (heat_donor || heat_target) { + unsigned short interface_type; if (heat_donor && heat_target){ - interface_type = CONJUGATE_HEAT_SS; + interface_type = ENUM_TRANSFER::CONJUGATE_HEAT_SS; } else { const auto fluidZone = heat_target? donor : target; if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) || (config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET)) - interface_type = heat_target? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF; + interface_type = heat_target? ENUM_TRANSFER::CONJUGATE_HEAT_FS : ENUM_TRANSFER::CONJUGATE_HEAT_SF; else if (config[fluidZone]->GetWeakly_Coupled_Heat()) - interface_type = heat_target? CONJUGATE_HEAT_WEAKLY_FS : CONJUGATE_HEAT_WEAKLY_SF; + interface_type = heat_target? ENUM_TRANSFER::CONJUGATE_HEAT_WEAKLY_FS : ENUM_TRANSFER::CONJUGATE_HEAT_WEAKLY_SF; else - interface_type = NO_TRANSFER; + interface_type = ENUM_TRANSFER::NO_TRANSFER; } if (interface_type != NO_TRANSFER) { auto nVar = 4; interface[donor][target] = new CConjugateHeatInterface(nVar, 0); + interface[donor][target]->SetInterfaceType(interface_type); if (rank == MASTER_NODE) cout << "conjugate heat variables." << endl; } else { @@ -2534,7 +2521,6 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet SU2_MPI::Error("Could not determine the number of variables for transfer.", CURRENT_FUNCTION); auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnVar(); - interface_type = CONSERVATIVE_VARIABLES; interface[donor][target] = new CConservativeVarsInterface(nVar, 0); if (rank == MASTER_NODE) cout << "generic conservative variables." << endl; } @@ -2629,9 +2615,8 @@ void CDriver::PreprocessOutput(CConfig **config, CConfig *driver_config, COutput } - -void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, CSolver***** solver, - CInterface*** interface, bool dummy){ +void CDriver::PreprocessTurboVertex(CConfig** config, CGeometry**** geometry, CSolver***** solver, + CInterface*** interface, CIteration*** iteration, bool dummy){ unsigned short donorZone,targetZone, nMarkerInt, iMarkerInt; unsigned short nSpanMax = 0; @@ -2665,12 +2650,16 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, } } if (rank == MASTER_NODE) cout<<"Max number of span-wise sections among all zones: "<< nSpanMax<<"."<< endl; +} - - if (rank == MASTER_NODE) cout<<"Initialize solver containers for average quantities." << endl; - for (iZone = 0; iZone < nZone; iZone++) { - solver[iZone][INST_0][MESH_0][FLOW_SOL]->InitTurboContainers(geometry[iZone][INST_0][MESH_0],config[iZone]); - } +void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, CSolver***** solver, + CInterface*** interface, CIteration*** iteration, bool dummy){ + unsigned short donorZone,targetZone, nMarkerInt, iMarkerInt; + unsigned short nSpanMax = 0; + bool restart = (config[ZONE_0]->GetRestart() || config[ZONE_0]->GetRestart_Flow()); + mixingplane = config[ZONE_0]->GetBoolMixingPlaneInterface(); + bool discrete_adjoint = config[ZONE_0]->GetDiscrete_Adjoint(); + su2double areaIn, areaOut, nBlades, flowAngleIn, flowAngleOut; // TODO(turbo): make it general for turbo HB if (rank == MASTER_NODE) cout<<"Compute inflow and outflow average geometric quantities." << endl; @@ -2680,6 +2669,10 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, geometry[iZone][INST_0][MESH_0]->GatherInOutAverageValues(config[iZone], true); } + if (rank == MASTER_NODE) cout<<"Initialize solver containers for average quantities." << endl; + for (iZone = 0; iZone < nZone; iZone++) { + solver[iZone][INST_0][MESH_0][FLOW_SOL]->InitTurboContainers(geometry[iZone][INST_0][MESH_0],config, iZone); + } if(mixingplane){ if (rank == MASTER_NODE) cout << "Set span-wise sections between zones on Mixing-Plane interface." << endl; @@ -2692,10 +2685,6 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, } } - for (iZone = 0; iZone < nZone-1; iZone++) { - geometry[nZone-1][INST_0][MESH_0]->SetAvgTurboGeoValues(config[iZone],geometry[iZone][INST_0][MESH_0], iZone); - } - /*--- Transfer number of blade to ZONE_0 to correctly compute turbo performance---*/ for (iZone = 1; iZone < nZone; iZone++) { nBlades = config[iZone]->GetnBlades(iZone); @@ -2713,23 +2702,6 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, } } - - if(mixingplane){ - if (rank == MASTER_NODE) cout<<"Preprocessing of the Mixing-Plane Interface." << endl; - for (donorZone = 0; donorZone < nZone; donorZone++) { - nMarkerInt = config_container[donorZone]->GetnMarker_MixingPlaneInterface()/2; - for (iMarkerInt = 1; iMarkerInt <= nMarkerInt; iMarkerInt++){ - for (targetZone = 0; targetZone < nZone; targetZone++) { - if (interface_types[donorZone][targetZone]==MIXING_PLANE){ - interface[donorZone][targetZone]->PreprocessAverage(geometry[donorZone][INST_0][MESH_0], geometry[targetZone][INST_0][MESH_0], - config[donorZone], config[targetZone], - iMarkerInt); - } - } - } - } - } - if(!restart && !discrete_adjoint){ if (rank == MASTER_NODE) cout<<"Initialize turbomachinery solution quantities." << endl; for(iZone = 0; iZone < nZone; iZone++) { @@ -2756,7 +2728,6 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, } } - } CDriver::~CDriver() = default; diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index a570755d29f..cbfc4576df0 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -547,14 +547,20 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar config_container[targetZone]); }; - switch (interface_types[donorZone][targetZone]) { + // Zones are equal or unconnected + if(donorZone == targetZone || interface_container[donorZone][targetZone] == nullptr) return UpdateMesh; + + switch (interface_container[donorZone][targetZone]->GetInterfaceType()) { case SLIDING_INTERFACE: BroadcastData(FLOW_SOL, FLOW_SOL); /*--- Additional transfer for turbulence variables. ---*/ if (config_container[targetZone]->GetKind_Solver() == MAIN_SOLVER::RANS || - config_container[targetZone]->GetKind_Solver() == MAIN_SOLVER::INC_RANS) { + config_container[targetZone]->GetKind_Solver() == MAIN_SOLVER::INC_RANS || + config_container[targetZone]->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_RANS || + config_container[targetZone]->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_INC_RANS + ) { BroadcastData(TURB_SOL, TURB_SOL); } @@ -587,28 +593,16 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar break; case MIXING_PLANE: { - const auto nMarkerInt = config_container[donorZone]->GetnMarker_MixingPlaneInterface() / 2; - - /*--- Transfer the average value from the donorZone to the targetZone ---*/ - /*--- Loops over the mixing planes defined in the config file to find the correct mixing plane for the donor-target combination ---*/ - for (auto iMarkerInt = 1; iMarkerInt <= nMarkerInt; iMarkerInt++) { - interface_container[donorZone][targetZone]->AllgatherAverage(solver_container[donorZone][INST_0][MESH_0][FLOW_SOL],solver_container[targetZone][INST_0][MESH_0][FLOW_SOL], - geometry_container[donorZone][INST_0][MESH_0],geometry_container[targetZone][INST_0][MESH_0], - config_container[donorZone], config_container[targetZone], iMarkerInt ); - } - - /*--- Set average value donorZone->targetZone ---*/ - interface_container[donorZone][targetZone]->SetAverageValues(solver_container[donorZone][INST_0][MESH_0][FLOW_SOL],solver_container[targetZone][INST_0][MESH_0][FLOW_SOL], donorZone); - - /*--- Set average geometrical properties FROM donorZone IN targetZone ---*/ - geometry_container[targetZone][INST_0][MESH_0]->SetAvgTurboGeoValues(config_container[iZone],geometry_container[iZone][INST_0][MESH_0], iZone); - + interface_container[donorZone][targetZone]->BroadcastData_MixingPlane( + *interpolator_container[donorZone][targetZone].get(), + solver_container[donorZone][INST_0][MESH_0][FLOW_SOL], + solver_container[targetZone][INST_0][MESH_0][FLOW_SOL], + geometry_container[donorZone][INST_0][MESH_0], + geometry_container[targetZone][INST_0][MESH_0], + config_container[donorZone], + config_container[targetZone]); break; } - case NO_TRANSFER: - case ZONES_ARE_EQUAL: - case NO_COMMON_INTERFACE: - break; default: if(rank == MASTER_NODE) cout << "WARNING: One of the intended interface transfer routines is not " @@ -619,16 +613,6 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar return UpdateMesh; } - - -void CMultizoneDriver::SetTurboPerformance() { - for (auto donorZone = 1u; donorZone < nZone; donorZone++) { - interface_container[donorZone][ZONE_0]->SetAverageValues(solver_container[donorZone][INST_0][MESH_0][FLOW_SOL], - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL], - donorZone); - } -} - bool CMultizoneDriver::Monitor(unsigned long TimeIter) { /*--- Check whether the inner solver has converged --- */ diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index 039fad933a4..ed33ce184bc 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -85,19 +85,9 @@ void CIntegration::Space_Integration(CGeometry *geometry, if (config->GetBoolGiles() && config->GetSpatialFourier()){ solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, conv_bound_numerics, INFLOW); - solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, conv_bound_numerics, OUTFLOW); } - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - if (config->GetBoolTurbomachinery()){ - /*--- Average quantities at the inflow and outflow boundaries ---*/ - solver_container[MainSolver]->TurboAverageProcess(solver_container, geometry,config,INFLOW); - solver_container[MainSolver]->TurboAverageProcess(solver_container, geometry, config, OUTFLOW); - } - } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - /*--- Weak boundary conditions ---*/ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index f2a39186909..1ad6c6f8b1c 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -668,20 +668,6 @@ void CMultiGridIntegration::NonDimensional_Parameters(CGeometry **geometry, CSol solver_container[FinestMesh][FLOW_SOL]->Momentum_Forces(geometry[FinestMesh], config); solver_container[FinestMesh][FLOW_SOL]->Friction_Forces(geometry[FinestMesh], config); - /*--- Calculate the turbo performance ---*/ - if (config->GetBoolTurbomachinery()){ - - /*--- Average quantities at the inflow and outflow boundaries ---*/ - - solver_container[FinestMesh][FLOW_SOL]->TurboAverageProcess(solver_container[FinestMesh], geometry[FinestMesh],config,INFLOW); - solver_container[FinestMesh][FLOW_SOL]->TurboAverageProcess(solver_container[FinestMesh], geometry[FinestMesh], config, OUTFLOW); - - /*--- Gather Inflow and Outflow quantities on the Master Node to compute performance ---*/ - - solver_container[FinestMesh][FLOW_SOL]->GatherInOutAverageValues(config, geometry[FinestMesh]); - - } - break; case RUNTIME_ADJFLOW_SYS: diff --git a/SU2_CFD/src/interfaces/CInterface.cpp b/SU2_CFD/src/interfaces/CInterface.cpp index 95b95b4a18a..1fa0d1d8755 100644 --- a/SU2_CFD/src/interfaces/CInterface.cpp +++ b/SU2_CFD/src/interfaces/CInterface.cpp @@ -199,443 +199,3 @@ void CInterface::BroadcastData(const CInterpolator& interpolator, } } } - -void CInterface::PreprocessAverage(CGeometry *donor_geometry, CGeometry *target_geometry, - const CConfig *donor_config, const CConfig *target_config, - unsigned short iMarkerInt){ - - unsigned short nMarkerDonor, nMarkerTarget; // Number of markers on the interface, donor and target side - unsigned short iMarkerDonor, iMarkerTarget; // Variables for iteration over markers - unsigned short iSpan,jSpan, tSpan = 0, kSpan = 0, nSpanDonor, nSpanTarget, Donor_Flag = 0, Target_Flag = 0; - int Marker_Donor = -1, Marker_Target = -1; - - const su2double *SpanValuesDonor, *SpanValuesTarget; - su2double dist, test, dist2, test2; - - nMarkerDonor = donor_geometry->GetnMarker(); - nMarkerTarget = target_geometry->GetnMarker(); - //TODO turbo this approach only works if all the turboamchinery marker - // of all zones have the same amount of span wise sections. - //TODO turbo initialization needed for the MPI routine should be place somewhere else. - nSpanDonor = donor_config->GetnSpanWiseSections(); - nSpanTarget = target_config->GetnSpanWiseSections(); - - /*--- On the donor side ---*/ - for (iMarkerDonor = 0; iMarkerDonor < nMarkerDonor; iMarkerDonor++){ - /*--- If the tag GetMarker_All_MixingPlaneInterface equals the index we are looping at ---*/ - if ( donor_config->GetMarker_All_MixingPlaneInterface(iMarkerDonor) == iMarkerInt ){ - /*--- We have identified the local index of the Donor marker ---*/ - /*--- Now we are going to store the average values that belong to Marker_Donor on each processor ---*/ - /*--- Store the identifier for the structural marker ---*/ - Marker_Donor = iMarkerDonor; - Donor_Flag = donor_config->GetMarker_All_TurbomachineryFlag(iMarkerDonor); - /*--- Exit the for loop: we have found the local index for Mixing-Plane interface ---*/ - break; - } - /*--- If the tag hasn't matched any tag within the donor markers ---*/ - Marker_Donor = -1; - Donor_Flag = -1; - - } - -#ifdef HAVE_MPI - auto BuffMarkerDonor = new int[size]; - auto BuffDonorFlag = new int[size]; - for (int iSize=0; iSize= 0.0){ - Marker_Donor = BuffMarkerDonor[iSize]; - Donor_Flag = BuffDonorFlag[iSize]; - break; - } - } - delete [] BuffMarkerDonor; - delete [] BuffDonorFlag; -#endif - - /*--- On the target side we have to identify the marker as well ---*/ - - for (iMarkerTarget = 0; iMarkerTarget < nMarkerTarget; iMarkerTarget++){ - /*--- If the tag GetMarker_All_MixingPlaneInterface(iMarkerTarget) equals the index we are looping at ---*/ - if ( target_config->GetMarker_All_MixingPlaneInterface(iMarkerTarget) == iMarkerInt ){ - /*--- Store the identifier for the fluid marker ---*/ - - // here i should then store it in the target zone - - Marker_Target = iMarkerTarget; - Target_Flag = target_config->GetMarker_All_TurbomachineryFlag(iMarkerTarget); - /*--- Exit the for loop: we have found the local index for iMarkerFSI on the FEA side ---*/ - break; - } - /*--- If the tag hasn't matched any tag within the Flow markers ---*/ - Marker_Target = -1; - Target_Flag = -1; - } - - if (Marker_Target != -1 && Marker_Donor != -1){ - - SpanValuesDonor = donor_geometry->GetSpanWiseValue(Donor_Flag); - SpanValuesTarget = target_geometry->GetSpanWiseValue(Target_Flag); - - - for(iSpan = 1; iSpan SpanValuesDonor[jSpan]){ - dist = test; - kSpan = jSpan; - } - if(test2 < dist2){ - dist2 = test2; - tSpan = jSpan; - } - - } - switch(donor_config->GetKind_MixingPlaneInterface()){ - case MATCHING: - SpanLevelDonor[iSpan] = iSpan; - SpanValueCoeffTarget[iSpan] = 0.0; - break; - case NEAREST_SPAN: - SpanLevelDonor[iSpan] = tSpan; - SpanValueCoeffTarget[iSpan] = 0.0; - break; - case LINEAR_INTERPOLATION: - SpanLevelDonor[iSpan] = kSpan; - SpanValueCoeffTarget[iSpan] = (SpanValuesTarget[iSpan] - SpanValuesDonor[kSpan]) - /(SpanValuesDonor[kSpan + 1] - SpanValuesDonor[kSpan]); - break; - default: - SU2_MPI::Error("MixingPlane interface option not implemented yet", CURRENT_FUNCTION); - break; - - } - } - } - -} - - -void CInterface::AllgatherAverage(CSolver *donor_solution, CSolver *target_solution, - CGeometry *donor_geometry, CGeometry *target_geometry, - const CConfig *donor_config, const CConfig *target_config, unsigned short iMarkerInt){ - - unsigned short nMarkerDonor, nMarkerTarget; // Number of markers on the interface, donor and target side - unsigned short iMarkerDonor, iMarkerTarget; // Variables for iteration over markers - unsigned short iSpan, nSpanDonor, nSpanTarget; - int Marker_Donor = -1, Marker_Target = -1; - su2double *avgPressureDonor = nullptr, *avgDensityDonor = nullptr, *avgNormalVelDonor = nullptr, - *avgTangVelDonor = nullptr, *avg3DVelDonor = nullptr, *avgNuDonor = nullptr, - *avgOmegaDonor = nullptr, *avgKineDonor = nullptr; - su2double *avgPressureTarget = nullptr, *avgDensityTarget = nullptr, *avgNormalVelTarget = nullptr, - *avg3DVelTarget = nullptr, *avgTangVelTarget = nullptr, *avgNuTarget = nullptr, - *avgOmegaTarget = nullptr, *avgKineTarget = nullptr; - -#ifdef HAVE_MPI - int iSize; - su2double *BuffAvgPressureDonor = nullptr, *BuffAvgDensityDonor = nullptr, *BuffAvgNormalVelDonor = nullptr, - *BuffAvg3DVelDonor = nullptr, *BuffAvgTangVelDonor = nullptr, *BuffAvgNuDonor = nullptr, - *BuffAvgKineDonor = nullptr, *BuffAvgOmegaDonor = nullptr; - int nSpanSize, *BuffMarkerDonor; -#endif - - - nMarkerTarget = target_geometry->GetnMarker(); - nMarkerDonor = donor_geometry->GetnMarker(); - nSpanDonor = donor_config->GetnSpanWiseSections() +1; - nSpanTarget = target_config->GetnSpanWiseSections() +1; - - - avgDensityDonor = new su2double[nSpanDonor]; - avgPressureDonor = new su2double[nSpanDonor]; - avgNormalVelDonor = new su2double[nSpanDonor]; - avgTangVelDonor = new su2double[nSpanDonor]; - avg3DVelDonor = new su2double[nSpanDonor]; - avgNuDonor = new su2double[nSpanDonor]; - avgKineDonor = new su2double[nSpanDonor]; - avgOmegaDonor = new su2double[nSpanDonor]; - - for (iSpan = 0; iSpan < nSpanDonor; iSpan++){ - avgDensityDonor[iSpan] = -1.0; - avgPressureDonor[iSpan] = -1.0; - avgNormalVelDonor[iSpan] = -1.0; - avgTangVelDonor[iSpan] = -1.0; - avg3DVelDonor[iSpan] = -1.0; - avgNuDonor[iSpan] = -1.0; - avgKineDonor[iSpan] = -1.0; - avgOmegaDonor[iSpan] = -1.0; - } - - avgDensityTarget = new su2double[nSpanTarget]; - avgPressureTarget = new su2double[nSpanTarget]; - avgNormalVelTarget = new su2double[nSpanTarget]; - avgTangVelTarget = new su2double[nSpanTarget]; - avg3DVelTarget = new su2double[nSpanTarget]; - avgNuTarget = new su2double[nSpanTarget]; - avgKineTarget = new su2double[nSpanTarget]; - avgOmegaTarget = new su2double[nSpanTarget]; - - - for (iSpan = 0; iSpan < nSpanTarget; iSpan++){ - avgDensityTarget[iSpan] = -1.0; - avgPressureTarget[iSpan] = -1.0; - avgNormalVelTarget[iSpan] = -1.0; - avgTangVelTarget[iSpan] = -1.0; - avg3DVelTarget[iSpan] = -1.0; - avgNuTarget[iSpan] = -1.0; - avgKineTarget[iSpan] = -1.0; - avgOmegaTarget[iSpan] = -1.0; - } - - /*--- Outer loop over the markers on the Mixing-Plane interface: compute one by one ---*/ - /*--- The tags are always an integer greater than 1: loop from 1 to nMarkerMixingPlane ---*/ - Marker_Donor = -1; - Marker_Target = -1; - - /*--- The donor and target markers are tagged with the same index. - *--- This is independent of the MPI domain decomposition. - *--- We need to loop over all markers on both sides ---*/ - - /*--- On the donor side ---*/ - - for (iMarkerDonor = 0; iMarkerDonor < nMarkerDonor; iMarkerDonor++){ - /*--- If the tag GetMarker_All_MixingPlaneInterface equals the index we are looping at ---*/ - if ( donor_config->GetMarker_All_MixingPlaneInterface(iMarkerDonor) == iMarkerInt ){ - /*--- We have identified the local index of the Donor marker ---*/ - /*--- Now we are going to store the average values that belong to Marker_Donor on each processor ---*/ - /*--- Store the identifier for the structural marker ---*/ - Marker_Donor = iMarkerDonor; - /*--- Exit the for loop: we have found the local index for Mixing-Plane interface ---*/ - break; - } - /*--- If the tag hasn't matched any tag within the donor markers ---*/ - Marker_Donor = -1; - - } - /*--- Here we want to make available the quantities for all the processors and collect them in a buffer - * for each span of the donor the span-wise height vector also so - * that then we can interpolate on the target side ---*/ - if (Marker_Donor != -1){ - for(iSpan = 0; iSpan < nSpanDonor; iSpan++){ - GetDonor_Variable(donor_solution, donor_geometry, donor_config, Marker_Donor, iSpan, rank); - avgDensityDonor[iSpan] = Donor_Variable[0]; - avgPressureDonor[iSpan] = Donor_Variable[1]; - avgNormalVelDonor[iSpan] = Donor_Variable[2]; - avgTangVelDonor[iSpan] = Donor_Variable[3]; - avg3DVelDonor[iSpan] = Donor_Variable[4]; - avgNuDonor[iSpan] = Donor_Variable[5]; - avgKineDonor[iSpan] = Donor_Variable[6]; - avgOmegaDonor[iSpan] = Donor_Variable[7]; - } - } - -#ifdef HAVE_MPI - nSpanSize = size*nSpanDonor; - BuffAvgDensityDonor = new su2double[nSpanSize]; - BuffAvgPressureDonor = new su2double[nSpanSize]; - BuffAvgNormalVelDonor = new su2double[nSpanSize]; - BuffAvgTangVelDonor = new su2double[nSpanSize]; - BuffAvg3DVelDonor = new su2double[nSpanSize]; - BuffAvgNuDonor = new su2double[nSpanSize]; - BuffAvgKineDonor = new su2double[nSpanSize]; - BuffAvgOmegaDonor = new su2double[nSpanSize]; - BuffMarkerDonor = new int[size]; - - for (iSpan=0;iSpan 0.0){ - for (iSpan = 0; iSpan < nSpanDonor; iSpan++){ - avgDensityDonor[iSpan] = BuffAvgDensityDonor[nSpanDonor*iSize + iSpan]; - avgPressureDonor[iSpan] = BuffAvgPressureDonor[nSpanDonor*iSize + iSpan]; - avgNormalVelDonor[iSpan] = BuffAvgNormalVelDonor[nSpanDonor*iSize + iSpan]; - avgTangVelDonor[iSpan] = BuffAvgTangVelDonor[nSpanDonor*iSize + iSpan]; - avg3DVelDonor[iSpan] = BuffAvg3DVelDonor[nSpanDonor*iSize + iSpan]; - avgNuDonor[iSpan] = BuffAvgNuDonor[nSpanDonor*iSize + iSpan]; - avgKineDonor[iSpan] = BuffAvgKineDonor[nSpanDonor*iSize + iSpan]; - avgOmegaDonor[iSpan] = BuffAvgOmegaDonor[nSpanDonor*iSize + iSpan]; - } - Marker_Donor = BuffMarkerDonor[iSize]; - break; - } - } - delete [] BuffAvgDensityDonor; - delete [] BuffAvgPressureDonor; - delete [] BuffAvgNormalVelDonor; - delete [] BuffAvgTangVelDonor; - delete [] BuffAvg3DVelDonor; - delete [] BuffAvgNuDonor; - delete [] BuffAvgKineDonor; - delete [] BuffAvgOmegaDonor; - delete [] BuffMarkerDonor; -#endif - - /*--- On the target side we have to identify the marker as well ---*/ - for (iMarkerTarget = 0; iMarkerTarget < nMarkerTarget; iMarkerTarget++){ - /*--- If the tag GetMarker_All_MixingPlaneInterface(iMarkerTarget) equals the index we are looping at ---*/ - if ( target_config->GetMarker_All_MixingPlaneInterface(iMarkerTarget) == iMarkerInt ){ - /*--- Store the identifier for the fluid marker ---*/ - Marker_Target = iMarkerTarget; - /*--- Exit the for loop: we have found the local index for iMarkerFSI on the FEA side ---*/ - break; - } - /*--- If the tag hasn't matched any tag within the Flow markers ---*/ - Marker_Target = -1; - - } - - - if (Marker_Target != -1 && Marker_Donor != -1){ - - /*--- linear interpolation of the average value of for the internal span-wise levels ---*/ - for(iSpan = 1; iSpan < nSpanTarget -2 ; iSpan++){ - avgDensityTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avgDensityDonor[SpanLevelDonor[iSpan] + 1] - - avgDensityDonor[SpanLevelDonor[iSpan]]); - avgDensityTarget[iSpan] += avgDensityDonor[SpanLevelDonor[iSpan]]; - avgPressureTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avgPressureDonor[SpanLevelDonor[iSpan] + 1] - - avgPressureDonor[SpanLevelDonor[iSpan]]); - avgPressureTarget[iSpan] += avgPressureDonor[SpanLevelDonor[iSpan]]; - avgNormalVelTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avgNormalVelDonor[SpanLevelDonor[iSpan] + 1] - - avgNormalVelDonor[SpanLevelDonor[iSpan]]); - avgNormalVelTarget[iSpan] += avgNormalVelDonor[SpanLevelDonor[iSpan]]; - avgTangVelTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avgTangVelDonor[SpanLevelDonor[iSpan] + 1] - - avgTangVelDonor[SpanLevelDonor[iSpan]]); - avgTangVelTarget[iSpan] += avgTangVelDonor[SpanLevelDonor[iSpan]]; - avg3DVelTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avg3DVelDonor[SpanLevelDonor[iSpan] + 1] - - avg3DVelDonor[SpanLevelDonor[iSpan]]); - avg3DVelTarget[iSpan] += avg3DVelDonor[SpanLevelDonor[iSpan]]; - avgNuTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avgNuDonor[SpanLevelDonor[iSpan] + 1] - - avgNuDonor[SpanLevelDonor[iSpan]]); - avgNuTarget[iSpan] += avgNuDonor[SpanLevelDonor[iSpan]]; - avgKineTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avgKineDonor[SpanLevelDonor[iSpan] + 1] - - avgKineDonor[SpanLevelDonor[iSpan]]); - avgKineTarget[iSpan] += avgKineDonor[SpanLevelDonor[iSpan]]; - avgOmegaTarget[iSpan] = SpanValueCoeffTarget[iSpan]*(avgOmegaDonor[SpanLevelDonor[iSpan] + 1] - - avgOmegaDonor[SpanLevelDonor[iSpan] ]); - avgOmegaTarget[iSpan] += avgOmegaDonor[SpanLevelDonor[iSpan]]; - } - - - /*--- transfer values at the hub ---*/ - avgDensityTarget[0] = avgDensityDonor[0]; - avgPressureTarget[0] = avgPressureDonor[0]; - avgNormalVelTarget[0] = avgNormalVelDonor[0]; - avgTangVelTarget[0] = avgTangVelDonor[0]; - avg3DVelTarget[0] = avg3DVelDonor[0]; - avgNuTarget[0] = avgNuDonor[0]; - avgKineTarget[0] = avgKineDonor[0]; - avgOmegaTarget[0] = avgOmegaDonor[0]; - - /*--- transfer values at the shroud ---*/ - avgDensityTarget[nSpanTarget - 2] = avgDensityDonor[nSpanDonor - 2]; - avgPressureTarget[nSpanTarget - 2] = avgPressureDonor[nSpanDonor - 2]; - avgNormalVelTarget[nSpanTarget - 2] = avgNormalVelDonor[nSpanDonor - 2]; - avgTangVelTarget[nSpanTarget - 2] = avgTangVelDonor[nSpanDonor - 2]; - avg3DVelTarget[nSpanTarget - 2] = avg3DVelDonor[nSpanDonor - 2]; - avgNuTarget[nSpanTarget - 2] = avgNuDonor[nSpanDonor - 2]; - avgKineTarget[nSpanTarget - 2] = avgKineDonor[nSpanDonor - 2]; - avgOmegaTarget[nSpanTarget - 2] = avgOmegaDonor[nSpanDonor - 2]; - - /*--- transfer 1D values ---*/ - avgDensityTarget[nSpanTarget - 1] = avgDensityDonor[nSpanDonor - 1]; - avgPressureTarget[nSpanTarget - 1] = avgPressureDonor[nSpanDonor - 1]; - avgNormalVelTarget[nSpanTarget - 1] = avgNormalVelDonor[nSpanDonor - 1]; - avgTangVelTarget[nSpanTarget - 1] = avgTangVelDonor[nSpanDonor - 1]; - avg3DVelTarget[nSpanTarget - 1] = avg3DVelDonor[nSpanDonor - 1]; - avgNuTarget[nSpanTarget - 1] = avgNuDonor[nSpanDonor - 1]; - avgKineTarget[nSpanTarget - 1] = avgKineDonor[nSpanDonor - 1]; - avgOmegaTarget[nSpanTarget - 1] = avgOmegaDonor[nSpanDonor - 1]; - - - /*---finally, the interpolated value is sent to the target zone ---*/ - for(iSpan = 0; iSpan < nSpanTarget ; iSpan++){ - Target_Variable[0] = avgDensityTarget[iSpan]; - Target_Variable[1] = avgPressureTarget[iSpan]; - Target_Variable[2] = avgNormalVelTarget[iSpan]; - Target_Variable[3] = avgTangVelTarget[iSpan]; - Target_Variable[4] = avg3DVelTarget[iSpan]; - Target_Variable[5] = avgNuTarget[iSpan]; - Target_Variable[6] = avgKineTarget[iSpan]; - Target_Variable[7] = avgOmegaTarget[iSpan]; - - - SetTarget_Variable(target_solution, target_geometry, target_config, Marker_Target, iSpan, rank); - } - } - - delete [] avgDensityDonor; - delete [] avgPressureDonor; - delete [] avgNormalVelDonor; - delete [] avgTangVelDonor; - delete [] avg3DVelDonor; - delete [] avgNuDonor; - delete [] avgKineDonor; - delete [] avgOmegaDonor; - - - delete [] avgDensityTarget; - delete [] avgPressureTarget; - delete [] avgNormalVelTarget; - delete [] avgTangVelTarget; - delete [] avg3DVelTarget; - delete [] avgNuTarget; - delete [] avgKineTarget; - delete [] avgOmegaTarget; -} diff --git a/SU2_CFD/src/interfaces/cfd/CConservativeVarsInterface.cpp b/SU2_CFD/src/interfaces/cfd/CConservativeVarsInterface.cpp index e5b0c8e2b21..df308fb2f21 100644 --- a/SU2_CFD/src/interfaces/cfd/CConservativeVarsInterface.cpp +++ b/SU2_CFD/src/interfaces/cfd/CConservativeVarsInterface.cpp @@ -33,6 +33,7 @@ CConservativeVarsInterface::CConservativeVarsInterface(unsigned short val_nVar, unsigned short val_nConst) : CInterface(val_nVar, val_nConst) { + InterfaceType = ENUM_TRANSFER::CONSERVATIVE_VARIABLES; } void CConservativeVarsInterface::GetDonor_Variable(CSolver *donor_solution, CGeometry *donor_geometry, diff --git a/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp b/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp index 84cf6a2b803..309a01ba70f 100644 --- a/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp +++ b/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp @@ -27,14 +27,17 @@ */ #include "../../../include/interfaces/cfd/CMixingPlaneInterface.hpp" +#include "../../../Common/include/interface_interpolation/CInterpolator.hpp" #include "../../../../Common/include/CConfig.hpp" #include "../../../../Common/include/geometry/CGeometry.hpp" #include "../../../include/solvers/CSolver.hpp" CMixingPlaneInterface::CMixingPlaneInterface(unsigned short val_nVar, unsigned short val_nConst){ - nVar = val_nVar; - Donor_Variable = new su2double[nVar + 5](); - Target_Variable = new su2double[nVar + 5](); + nVar = val_nVar; // Solver vars + nMixingVars = 8; // Always 8 vars in turbo MP + Donor_Variable = new su2double[nMixingVars](); + Target_Variable = new su2double[nMixingVars](); + InterfaceType = ENUM_TRANSFER::MIXING_PLANE; } void CMixingPlaneInterface::SetSpanWiseLevels(const CConfig *donor_config, const CConfig *target_config){ @@ -53,68 +56,157 @@ void CMixingPlaneInterface::SetSpanWiseLevels(const CConfig *donor_config, const } } -void CMixingPlaneInterface::GetDonor_Variable(CSolver *donor_solution, CGeometry *donor_geometry, - const CConfig *donor_config, unsigned long Marker_Donor, - unsigned long iSpan, unsigned long rank) { +void CMixingPlaneInterface::BroadcastData_MixingPlane(const CInterpolator& interpolator, + CSolver *donor_solution, CSolver *target_solution, + CGeometry *donor_geometry, CGeometry *target_geometry, + const CConfig *donor_config, const CConfig *target_config) { + static_assert(su2activematrix::Storage == StorageType::RowMajor,""); + + /*--- Loop over interface markers. ---*/ + + for (auto iMarkerInt = 0u; iMarkerInt < donor_config->GetnMarker_MixingPlaneInterface()/2; iMarkerInt++) { + + /*--- Find the markers containing the interface ---*/ + short markDonor = donor_config->FindMixingPlaneInterfaceMarker(donor_geometry->GetnMarker()); + short markTarget= target_config->FindMixingPlaneInterfaceMarker(target_geometry->GetnMarker()); + + /*--- Check if this interface connects the two zones, if not continue. ---*/ + if(!CInterpolator::CheckInterfaceBoundary(markDonor, markTarget)) continue; + + // The number of spans is available on every rank + const auto nSpanDonor = donor_config->GetnSpanWiseSections(); + + /*--- Fill send buffers. ---*/ + + vector sendDonorMarker(nSpanDonor); + vector sendDonorVar(nSpanDonor * nMixingVars); + + if (markDonor != -1) { + // cout << "markDonor Identified!" << endl; + for (auto iSpan = 0ul; iSpan < nSpanDonor; iSpan++) { + GetDonor_Variable(donor_solution, donor_geometry, donor_config, markDonor, iSpan, 0); + for (auto iVar = 0u; iVar < nMixingVars; iVar++) sendDonorVar[iSpan * nMixingVars + iVar] = Donor_Variable[iVar]; + sendDonorMarker[iSpan] = markDonor; + } + } +#ifdef HAVE_MPI + /*--- Gather data. ---*/ + const int nTotalDonors = nSpanDonor * size; + const int nSpanDonorVars = nSpanDonor * nMixingVars; + vector buffDonorMarker(nTotalDonors); + vector buffDonorVar(nTotalDonors * nMixingVars); + + SU2_MPI::Allgather(sendDonorMarker.data(), nSpanDonor, MPI_UNSIGNED_SHORT, + buffDonorMarker.data(), nSpanDonor, MPI_UNSIGNED_SHORT, + SU2_MPI::GetComm()); + + SU2_MPI::Allgather(sendDonorVar.data(), nSpanDonorVars, MPI_DOUBLE, + buffDonorVar.data(), nSpanDonorVars, MPI_DOUBLE, + SU2_MPI::GetComm()); + + for (auto iSize = 0; iSize < size; iSize++){ + if (buffDonorVar[nSpanDonorVars * iSize] > 0.0) { + for (auto iSpan = 0ul; iSpan < nSpanDonor; iSpan++){ + for (auto iVar = 0u; iVar < nMixingVars; iVar++) sendDonorVar[iSpan * nMixingVars + iVar] = buffDonorVar[iSize * nSpanDonorVars + iVar]; // This could be wrong in 3D + } + markDonor = buffDonorMarker[iSize * nSpanDonor]; + } + } +#endif + + /*--- This rank does not need to do more work. ---*/ + if (!(markTarget != -1 && markDonor != -1)) continue; + + /*--- Loop over target spans. ---*/ + auto nTargetSpan = target_config->GetnSpanWiseSections() + 1; + + for (auto iTargetSpan = 0ul; iTargetSpan < nTargetSpan; iTargetSpan++) { - unsigned short nDim = nVar - 2; - bool turbulent = (donor_config->GetKind_Turb_Model() != TURB_MODEL::NONE); + auto& targetSpan = interpolator.targetSpans[iMarkerInt][markTarget]; + InitializeTarget_Variable(target_solution, markTarget, iTargetSpan, nSpanDonor); + if ((iTargetSpan == 0) || (iTargetSpan > nTargetSpan - 3)) { + /*--- Transfer values at hub, shroud and 1D values ---*/ + unsigned long donorSpan; + if (iTargetSpan == 0) donorSpan = 0; + else if (iTargetSpan == nTargetSpan - 2) donorSpan = nSpanDonor - 2; + else if (iTargetSpan == nTargetSpan - 1) donorSpan = nSpanDonor - 1; - Donor_Variable[0] = donor_solution->GetAverageDensity(Marker_Donor, iSpan); - Donor_Variable[1] = donor_solution->GetAveragePressure(Marker_Donor, iSpan); - Donor_Variable[2] = donor_solution->GetAverageTurboVelocity(Marker_Donor, iSpan)[0]; - Donor_Variable[3] = donor_solution->GetAverageTurboVelocity(Marker_Donor, iSpan)[1]; + RecoverTarget_Variable(sendDonorVar, donorSpan); - if(nDim == 3){ - Donor_Variable[4] = donor_solution->GetAverageTurboVelocity(Marker_Donor, iSpan)[2]; + SetTarget_Variable(target_solution, target_geometry, target_config, markTarget, iTargetSpan, 0); + } + else { + /*--- Get the global index of the donor and the interpolation coefficient. ---*/ + + const auto donorSpan = targetSpan.donorSpan; + const auto donorCoeff = targetSpan.coefficient; + + /*--- Recover the Target_Variable from the buffer of variables. ---*/ + RecoverTarget_Variable(sendDonorVar, donorSpan, donorCoeff); + + SetTarget_Variable(target_solution, target_geometry, target_config, markTarget, iTargetSpan, 0); + } + } + } +} + +void CMixingPlaneInterface::GetDonor_Variable(CSolver *donor_solution, CGeometry *donor_geometry, + const CConfig *donor_config, unsigned long Marker_Donor, + unsigned long Span_Donor, unsigned long Point_Donor) { + + Donor_Variable[0] = donor_solution->GetAverageDensity(Marker_Donor, Span_Donor); + Donor_Variable[1] = donor_solution->GetAveragePressure(Marker_Donor, Span_Donor); + Donor_Variable[2] = donor_solution->GetAverageTurboVelocity(Marker_Donor, Span_Donor)[0]; + Donor_Variable[3] = donor_solution->GetAverageTurboVelocity(Marker_Donor, Span_Donor)[1]; + + if(donor_geometry->GetnDim() == 3){ + Donor_Variable[4] = donor_solution->GetAverageTurboVelocity(Marker_Donor, Span_Donor)[2]; } else{ Donor_Variable[4] = -1.0; } - if(turbulent){ - Donor_Variable[5] = donor_solution->GetAverageNu(Marker_Donor, iSpan); - Donor_Variable[6] = donor_solution->GetAverageKine(Marker_Donor, iSpan); - Donor_Variable[7] = donor_solution->GetAverageOmega(Marker_Donor, iSpan); + if(donor_config->GetKind_Turb_Model() != TURB_MODEL::NONE){ + Donor_Variable[5] = donor_solution->GetAverageNu(Marker_Donor, Span_Donor); + Donor_Variable[6] = donor_solution->GetAverageKine(Marker_Donor, Span_Donor); + Donor_Variable[7] = donor_solution->GetAverageOmega(Marker_Donor, Span_Donor); } else{ Donor_Variable[5] = -1.0; Donor_Variable[6] = -1.0; Donor_Variable[7] = -1.0; } - } +void CMixingPlaneInterface::InitializeTarget_Variable(CSolver *target_solution, unsigned long Marker_Target, + unsigned long Span_Target, unsigned short nSpanDonor) { -void CMixingPlaneInterface::SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, - const CConfig *target_config, unsigned long Marker_Target, - unsigned long iSpan, unsigned long rank) { + target_solution->SetnMixingStates(Marker_Target, Span_Target, nSpanDonor); // This is to allocate + target_solution->SetMixingStateStructure(Marker_Target, Span_Target); + target_solution->SetnMixingStates(Marker_Target, Span_Target, 0); // Reset counter to 0 - unsigned short nDim = nVar - 2; - bool turbulent = (target_config->GetKind_Turb_Model() != TURB_MODEL::NONE); +} +void CMixingPlaneInterface::SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, + const CConfig *target_config, unsigned long Marker_Target, + unsigned long Span_Target, unsigned long Point_Target) { - target_solution->SetExtAverageDensity(Marker_Target, iSpan, Target_Variable[0]); - target_solution->SetExtAveragePressure(Marker_Target, iSpan, Target_Variable[1]); - target_solution->SetExtAverageTurboVelocity(Marker_Target, iSpan, 0, Target_Variable[2]); - target_solution->SetExtAverageTurboVelocity(Marker_Target, iSpan, 1, Target_Variable[3]); + unsigned short iVar, iDonorSpan; + /*--- Set the mixing plane solution with the value of the Target Variable ---*/ - if(nDim == 3){ - target_solution->SetExtAverageTurboVelocity(Marker_Target, iSpan, 2, Target_Variable[4]); - } + iDonorSpan = target_solution->GetnMixingStates(Marker_Target, Span_Target); - if(turbulent){ - target_solution->SetExtAverageNu(Marker_Target, iSpan, Target_Variable[5]); - target_solution->SetExtAverageKine(Marker_Target, iSpan, Target_Variable[6]); - target_solution->SetExtAverageOmega(Marker_Target, iSpan, Target_Variable[7]); - } + for (iVar = 0; iVar < nMixingVars; iVar++) + target_solution->SetMixingState(Marker_Target, Span_Target, iVar, Target_Variable[iVar]); + target_solution->SetnMixingStates( Marker_Target, Span_Target, iDonorSpan + 1 ); } void CMixingPlaneInterface::SetAverageValues(CSolver *donor_solution, CSolver *target_solution, unsigned short donorZone){ + unsigned short iSpan; for(iSpan = 0; iSpanGetFrozen_Visc_Disc()) { solvers0[ADJTURB_SOL]->RegisterSolution(geometry0, config[iZone]); } + if (config[iZone]->GetBoolTurbomachinery()) { + geometry0->RegisterCoordinates(); + solvers0[ADJFLOW_SOL]->Register_VertexNormals(geometry0, config[iZone], true); + } if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { solvers0[ADJSPECIES_SOL]->RegisterSolution(geometry0, config[iZone]); } @@ -458,12 +463,15 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, CConfig** config, unsigned short iZone, unsigned short iInst, RECORDING kind_recording) { + auto solvers0 = solver[iZone][iInst][MESH_0]; auto geometry0 = geometry[iZone][iInst][MESH_0]; if ((kind_recording == RECORDING::MESH_COORDS) || (kind_recording == RECORDING::CLEAR_INDICES) || - (kind_recording == RECORDING::SOLUTION_AND_MESH)) { + (kind_recording == RECORDING::SOLUTION_AND_MESH) || + (kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH) || + (kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH)) { /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ SU2_OMP_PARALLEL @@ -480,16 +488,23 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** solvers0[FLOW_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); - if (config[iZone]->GetBoolTurbomachinery()) { - solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], INFLOW); - solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], OUTFLOW); - } if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) { solvers0[TURB_SOL]->Postprocessing(geometry0, solvers0, config[iZone], MESH_0); solvers0[TURB_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); solvers0[TURB_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); } + if (config[iZone]->GetBoolTurbomachinery()) { + solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], INFLOW); + solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], OUTFLOW); + solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], INFLOW); + solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], OUTFLOW); + if (config[iZone]->GetBoolGiles() && config[iZone]->GetSpatialFourier()){ + auto conv_bound_numerics = numerics[iZone][iInst][MESH_0][FLOW_SOL][CONV_BOUND_TERM + omp_get_thread_num()*MAX_TERMS]; + solvers0[FLOW_SOL]->PreprocessBC_Giles(geometry0, config[iZone], conv_bound_numerics, INFLOW); + solvers0[FLOW_SOL]->PreprocessBC_Giles(geometry0, config[iZone], conv_bound_numerics, OUTFLOW); + } + } if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { solvers0[SPECIES_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); solvers0[SPECIES_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); @@ -528,6 +543,9 @@ void CDiscAdjFluidIteration::RegisterOutput(CSolver***** solver, CGeometry**** g if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) { solvers0[ADJTURB_SOL]->RegisterOutput(geometry0, config[iZone]); } + if (config[iZone]->GetBoolTurbomachinery()) { + solvers0[ADJFLOW_SOL]->Register_VertexNormals(geometry0, config[iZone], false); + } if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { solvers0[ADJSPECIES_SOL]->RegisterOutput(geometry0, config[iZone]); } diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index 529b3e12453..92dc60c408a 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -226,20 +226,21 @@ bool CFluidIteration::Monitor(COutput* output, CIntegration**** integration, CGe /*--- Turbomachinery Specific Montior ---*/ if (config[ZONE_0]->GetBoolTurbomachinery()){ if (val_iZone == config[ZONE_0]->GetnZone()-1) { - ComputeTurboPerformance(solver, geometry, config, config[val_iZone]->GetnInner_Iter()); + ComputeTurboPerformance(solver, geometry, config); + auto TurbomachineryBladePerformances = GetBladesPerformanceVector(solver, config[val_iZone]->GetnZone()); - output->SetHistoryOutput(geometry, solver, - config, TurbomachineryStagePerformance, TurbomachineryPerformance, val_iZone, config[val_iZone]->GetTimeIter(), config[val_iZone]->GetOuterIter(), - config[val_iZone]->GetInnerIter(), val_iInst); + output->SetHistoryOutput(geometry, solver, config, TurbomachineryStagePerformance, TurbomachineryBladePerformances, + val_iZone, config[val_iZone]->GetTimeIter(), config[val_iZone]->GetOuterIter(), + config[val_iZone]->GetInnerIter(), val_iInst); } /*--- Update ramps, grid first then outlet boundary ---*/ if (config[val_iZone]->GetRampMotionFrame()) - UpdateRamp(geometry, config, config[val_iZone]->GetInnerIter(), val_iZone, RAMP_TYPE::GRID); + UpdateRamps(geometry, config, config[val_iZone]->GetInnerIter(), val_iZone, RAMP_TYPE::GRID); } // Outside turbo scope as Riemann boundaries can be ramped (pressure only) if (config[val_iZone]->GetRampOutflow()) - UpdateRamp(geometry, config, config[val_iZone]->GetInnerIter(), val_iZone, RAMP_TYPE::BOUNDARY); + UpdateRamps(geometry, config, config[val_iZone]->GetInnerIter(), val_iZone, RAMP_TYPE::BOUNDARY); output->SetHistoryOutput(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], config[val_iZone], config[val_iZone]->GetTimeIter(), config[val_iZone]->GetOuterIter(), @@ -257,7 +258,7 @@ bool CFluidIteration::Monitor(COutput* output, CIntegration**** integration, CGe return StopCalc; } -void CFluidIteration::UpdateRamp(CGeometry**** geometry_container, CConfig** config_container, unsigned long iter, unsigned short iZone, RAMP_TYPE ramp_flag) { +void CFluidIteration::UpdateRamps(CGeometry**** geometry_container, CConfig** config_container, unsigned long iter, unsigned short iZone, RAMP_TYPE ramp_flag) { /*--- Generic function for handling ramps ---*/ // Grid updates (i.e. rotation/translation) handled seperately to boundary (i.e. pressure/mass flow) updates auto* config = config_container[iZone]; @@ -332,39 +333,6 @@ void CFluidIteration::UpdateRamp(CGeometry**** geometry_container, CConfig** con } } -void CFluidIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container, unsigned long ExtIter) { - unsigned short nDim = geometry_container[ZONE_0][INST_0][MESH_0]->GetnDim(); - unsigned short nBladesRow = config_container[ZONE_0]->GetnMarker_Turbomachinery(); - unsigned short iBlade=0, iSpan; - vector TurboPrimitiveIn, TurboPrimitiveOut; - std::vector> bladesPrimitives; - - if (rank == MASTER_NODE) { - for (iBlade = 0; iBlade < nBladesRow; iBlade++){ - /* Blade Primitive initialized per blade */ - std::vector bladePrimitives; - auto nSpan = config_container[iBlade]->GetnSpanWiseSections(); - for (iSpan = 0; iSpan < nSpan + 1; iSpan++) { - TurboPrimitiveIn= solver[iBlade][INST_0][MESH_0][FLOW_SOL]->GetTurboPrimitive(iBlade, iSpan, true); - TurboPrimitiveOut= solver[iBlade][INST_0][MESH_0][FLOW_SOL]->GetTurboPrimitive(iBlade, iSpan, false); - auto spanInletPrimitive = CTurbomachineryPrimitiveState(TurboPrimitiveIn, nDim, geometry_container[iBlade][INST_0][MESH_0]->GetTangGridVelIn(iBlade, iSpan)); - auto spanOutletPrimitive = CTurbomachineryPrimitiveState(TurboPrimitiveOut, nDim, geometry_container[iBlade][INST_0][MESH_0]->GetTangGridVelOut(iBlade, iSpan)); - auto spanCombinedPrimitive = CTurbomachineryCombinedPrimitiveStates(spanInletPrimitive, spanOutletPrimitive); - bladePrimitives.push_back(spanCombinedPrimitive); - } - bladesPrimitives.push_back(bladePrimitives); - } - TurbomachineryPerformance->ComputeTurbomachineryPerformance(bladesPrimitives); - - auto nSpan = config_container[ZONE_0]->GetnSpanWiseSections(); - auto InState = TurbomachineryPerformance->GetBladesPerformances().at(ZONE_0).at(nSpan)->GetInletState(); - nSpan = config_container[nZone-1]->GetnSpanWiseSections(); - auto OutState = TurbomachineryPerformance->GetBladesPerformances().at(nZone-1).at(nSpan)->GetOutletState(); - - TurbomachineryStagePerformance->ComputePerformanceStage(InState, OutState, config_container[nZone-1]); - } -} - void CFluidIteration::Postprocess(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, @@ -408,6 +376,9 @@ void CFluidIteration::Solve(COutput* output, CIntegration**** integration, CGeom Iterate(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, INST_0); + /*--- Postprocessing Step ---*/ + Postprocess(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, val_iInst); + /*--- Monitor the pseudo-time ---*/ StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, INST_0); diff --git a/SU2_CFD/src/iteration/CIteration.cpp b/SU2_CFD/src/iteration/CIteration.cpp index a9df4be5397..08e27fb190b 100644 --- a/SU2_CFD/src/iteration/CIteration.cpp +++ b/SU2_CFD/src/iteration/CIteration.cpp @@ -209,3 +209,23 @@ void CIteration::Output(COutput* output, CGeometry**** geometry, CSolver***** so output->SetResultFiles(geometry[val_iZone][INST_0][MESH_0], config[val_iZone], solver[val_iZone][INST_0][MESH_0], InnerIter); } + +void CIteration::InitTurboPerformance(CGeometry* geometry, CConfig** config, CFluidModel* fluid, unsigned short val_iZone) { + TurbomachineryStagePerformance = std::make_shared(*fluid); +} + +void CIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container) { + // Computes the turboperformance per blade in zone iBlade + const auto nDim = geometry_container[ZONE_0][INST_0][MESH_0]->GetnDim(); + const auto nBladesRow = config_container[ZONE_0]->GetnMarker_Turbomachinery(); + const auto nZone = config_container[ZONE_0]->GetnZone(); + + auto TurbomachineryBladePerformances = GetBladesPerformanceVector(solver, nZone); + + auto nSpan = config_container[ZONE_0]->GetnSpanWiseSections(); + auto InState = TurbomachineryBladePerformances.at(ZONE_0)->GetBladesPerformances().at(nSpan)->GetInletState(); + nSpan = config_container[nZone-1]->GetnSpanWiseSections(); + auto OutState = TurbomachineryBladePerformances.at(nZone-1)->GetBladesPerformances().at(nSpan)->GetOutletState(); + + TurbomachineryStagePerformance->ComputePerformanceStage(InState, OutState, config_container[nZone-1]); +} diff --git a/SU2_CFD/src/iteration/CTurboIteration.cpp b/SU2_CFD/src/iteration/CTurboIteration.cpp index 2116576d9bf..1ddd23fdb0e 100644 --- a/SU2_CFD/src/iteration/CTurboIteration.cpp +++ b/SU2_CFD/src/iteration/CTurboIteration.cpp @@ -33,34 +33,37 @@ void CTurboIteration::Preprocess(COutput* output, CIntegration**** integration, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { - /*--- Average quantities at the inflow and outflow boundaries ---*/ - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( - solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], INFLOW); - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( - solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], OUTFLOW); - if (config[val_iZone]->GetBoolTurbomachinery()) { - InitTurboPerformance(geometry[val_iZone][INST_0][MESH_0], config, - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->GetFluidModel()); - } + /*--- Average quantities at the inflow and outflow boundaries ---*/ + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( + solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], INFLOW); + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( + solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], OUTFLOW); + + InitTurboPerformance(geometry[val_iZone][val_iInst][MESH_0], config, solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->GetFluidModel()); + } void CTurboIteration::Postprocess(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { - /*--- Average quantities at the inflow and outflow boundaries ---*/ - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( - solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], INFLOW); - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( - solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], OUTFLOW); + /*--- Average quantities at the inflow and outflow boundaries ---*/ + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( + solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], INFLOW); + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->TurboAverageProcess( + solver[val_iZone][val_iInst][MESH_0], geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], OUTFLOW); + + /*--- Gather Inflow and Outflow quantities on the Master Node to compute performance ---*/ + + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->GatherInOutAverageValues(config[val_iZone], geometry[val_iZone][val_iInst][MESH_0]); + + /*--- Compute the turboperformance ---*/ + + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->ComputeTurboBladePerformance(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], val_iZone); - /*--- Gather Inflow and Outflow quantities on the Master Node to compute performance ---*/ - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->GatherInOutAverageValues(config[val_iZone], - geometry[val_iZone][val_iInst][MESH_0]); } void CTurboIteration::InitTurboPerformance(CGeometry* geometry, CConfig** config, CFluidModel* fluid) { - TurbomachineryPerformance = std::make_shared(config, *geometry, *fluid); - TurbomachineryStagePerformance = std::make_shared(*fluid); -} \ No newline at end of file + TurbomachineryStagePerformance = std::make_shared(*fluid); +} diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 324b13c8775..6a6266d1c8d 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -244,6 +244,8 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi } AD::SetPreaccOut(Flux, nVar); + AD::SetPreaccOut(Jacobian_i, nVar, nVar); + AD::SetPreaccOut(Jacobian_j, nVar, nVar); AD::EndPreacc(); return ResidualType<>(Flux, Jacobian_i, Jacobian_j); diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index dff4abf037f..287bc223407 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -479,6 +479,8 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) } AD::SetPreaccOut(Proj_Flux_Tensor, nVar); + AD::SetPreaccOut(Jacobian_i, nVar, nVar); + AD::SetPreaccOut(Jacobian_j, nVar, nVar); AD::EndPreacc(); return ResidualType<>(Proj_Flux_Tensor, Jacobian_i, Jacobian_j); diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 2476498e84d..aba523129cb 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -503,7 +503,7 @@ bool CFlowCompOutput::WriteHistoryFileOutput(const CConfig *config) { return !config->GetFinite_Difference_Mode() && COutput::WriteHistoryFileOutput(config); } -void CFlowCompOutput::SetTurboPerformance_Output(std::shared_ptr TurboPerf, +void CFlowCompOutput::SetTurboPerformance_Output(std::vector> TurboBladePerfs, CConfig *config, unsigned long TimeIter, unsigned long OuterIter, @@ -515,8 +515,6 @@ void CFlowCompOutput::SetTurboPerformance_Output(std::shared_ptr T curInnerIter = InnerIter; stringstream TurboInOutTable, TurboPerfTable; - auto BladePerformance = TurboPerf->GetBladesPerformances(); - /*-- Table for Turbomachinery Performance Values --*/ PrintingToolbox::CTablePrinter TurboInOut(&TurboInOutTable); @@ -529,7 +527,7 @@ void CFlowCompOutput::SetTurboPerformance_Output(std::shared_ptr T for (unsigned short iZone = 0; iZone <= config->GetnZone()-1; iZone++) { auto nSpan = config->GetnSpan_iZones(iZone); - const auto& BladePerf = BladePerformance.at(iZone).at(nSpan); + const auto& BladePerf = TurboBladePerfs.at(iZone)->GetBladesPerformances().at(nSpan); TurboInOut<<" BLADE ROW INDEX "< T cout< TurboStagePerf, std::shared_ptr TurboPerf, CConfig *config) { +void CFlowCompOutput::SetTurboMultiZonePerformance_Output(std::shared_ptr TurboStagePerf, std::vector> TurboPerf, CConfig *config) { stringstream TurboMZPerf; @@ -579,11 +577,10 @@ void CFlowCompOutput::SetTurboMultiZonePerformance_Output(std::shared_ptr TurboStagePerf, std::shared_ptr TurboPerf, CConfig *config) { - auto BladePerformance = TurboPerf->GetBladesPerformances(); +void CFlowCompOutput::LoadTurboHistoryData(std::shared_ptr TurboStagePerf, std::vector> TurboBladePerfs, CConfig *config) { for (unsigned short iZone = 0; iZone <= config->GetnZone()-1; iZone++) { auto nSpan = config->GetnSpan_iZones(iZone); - const auto& BladePerf = BladePerformance.at(iZone).at(nSpan); + const auto& BladePerf = TurboBladePerfs.at(iZone)->GetBladesPerformances().at(nSpan); stringstream tag; tag << iZone + 1; @@ -625,7 +622,7 @@ void CFlowCompOutput::LoadTurboHistoryData(std::shared_ptrGetTotalPressureLoss()); } -void CFlowCompOutput::WriteTurboSpanwisePerformance(std::shared_ptr TurboPerf, CGeometry *geometry, CConfig **config, unsigned short val_iZone) { +void CFlowCompOutput::WriteTurboSpanwisePerformance(std::vector> TurboBladePerfs, CGeometry *geometry, CConfig **config, unsigned short val_iZone) { string inMarker_Tag, outMarker_Tag, inMarkerTag_Mix; unsigned short nZone = config[val_iZone]->GetnZone(); @@ -637,8 +634,6 @@ void CFlowCompOutput::WriteTurboSpanwisePerformance(std::shared_ptrGetBladesPerformances(); - /*--- Start of write file turboperformance spanwise ---*/ SpanWiseValuesIn = geometry->GetSpanWiseValue(INFLOW); SpanWiseValuesOut = geometry->GetSpanWiseValue(OUTFLOW); @@ -670,7 +665,7 @@ void CFlowCompOutput::WriteTurboSpanwisePerformance(std::shared_ptrGetnSpanWiseSections(); iSpan++){ - const auto& BladePerf = BladePerformance.at(val_iZone).at(iSpan); + const auto& BladePerf = TurboBladePerfs.at(val_iZone)->GetBladesPerformances().at(iSpan); file.width(30); file << SpanWiseValuesIn[iSpan]; file.width(15); file << iSpan; @@ -714,7 +709,7 @@ void CFlowCompOutput::WriteTurboSpanwisePerformance(std::shared_ptrGetnSpanWiseSections(); iSpan++){ - const auto& BladePerf = BladePerformance.at(val_iZone).at(iSpan); + const auto& BladePerf = TurboBladePerfs.at(val_iZone)->GetBladesPerformances().at(iSpan); file.width(30); file << SpanWiseValuesOut[iSpan]; file.width(15); file << iSpan; @@ -764,7 +759,7 @@ void CFlowCompOutput::WriteTurboSpanwisePerformance(std::shared_ptrGetnSpanWiseSections(); iSpan++){ - const auto& BladePerf = BladePerformance.at(val_iZone).at(iSpan); + const auto& BladePerf = TurboBladePerfs.at(val_iZone)->GetBladesPerformances().at(iSpan); file.width(30); file << SpanWiseValuesIn[iSpan]; file.width(15); file << iSpan; @@ -828,7 +823,7 @@ void CFlowCompOutput::WriteTurboSpanwisePerformance(std::shared_ptrGetnSpanWiseSections(); iSpan++){ - const auto& BladePerf = BladePerformance.at(val_iZone).at(iSpan); + const auto& BladePerf = TurboBladePerfs.at(val_iZone)->GetBladesPerformances().at(iSpan); file.width(30); file << SpanWiseValuesOut[iSpan]; file.width(15); file << iSpan; diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp index a3a2289191b..e081b0b7f9d 100644 --- a/SU2_CFD/src/output/COutput.cpp +++ b/SU2_CFD/src/output/COutput.cpp @@ -228,7 +228,7 @@ void COutput::SetHistoryOutput(CGeometry *geometry, } -void COutput::SetHistoryOutput(CGeometry ****geometry, CSolver *****solver, CConfig **config, std::shared_ptr(TurboStagePerf), std::shared_ptr TurboPerf, unsigned short val_iZone, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter, unsigned short val_iInst){ +void COutput::SetHistoryOutput(CGeometry ****geometry, CSolver *****solver, CConfig **config, std::shared_ptr(TurboStagePerf), std::vector> TurboBladePerfs, unsigned short val_iZone, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter, unsigned short val_iInst){ unsigned long Iter= InnerIter; @@ -237,19 +237,19 @@ void COutput::SetHistoryOutput(CGeometry ****geometry, CSolver *****solver, CCon /*--- Turbomachinery Performance Screen summary output---*/ if (Iter%100 == 0 && rank == MASTER_NODE) { - SetTurboPerformance_Output(TurboPerf, config[val_iZone], TimeIter, OuterIter, InnerIter); - SetTurboMultiZonePerformance_Output(TurboStagePerf, TurboPerf, config[val_iZone]); + SetTurboPerformance_Output(TurboBladePerfs, config[val_iZone], TimeIter, OuterIter, InnerIter); //Blade-row index scree + SetTurboMultiZonePerformance_Output(TurboStagePerf, TurboBladePerfs, config[val_iZone]); //Stage performance screen } for (int iZone = 0; iZone < config[ZONE_0]->GetnZone(); iZone ++){ if (rank == MASTER_NODE) { - WriteTurboSpanwisePerformance(TurboPerf, geometry[iZone][val_iInst][MESH_0], config, iZone); + WriteTurboSpanwisePerformance(TurboBladePerfs, geometry[iZone][val_iInst][MESH_0], config, iZone); //Spanwise files } } /*--- Update turboperformance history file*/ if (rank == MASTER_NODE){ - LoadTurboHistoryData(TurboStagePerf, TurboPerf, config[val_iZone]); + LoadTurboHistoryData(TurboStagePerf, TurboBladePerfs, config[val_iZone]); //History files } } diff --git a/SU2_CFD/src/output/CTurboOutput.cpp b/SU2_CFD/src/output/CTurboOutput.cpp index 5ecb0ac542f..a64d680c133 100644 --- a/SU2_CFD/src/output/CTurboOutput.cpp +++ b/SU2_CFD/src/output/CTurboOutput.cpp @@ -39,8 +39,7 @@ CTurbomachineryCombinedPrimitiveStates::CTurbomachineryCombinedPrimitiveStates( : InletPrimitiveState(inletPrimitiveState), OutletPrimitiveState(outletPrimitiveState) {} CTurbomachineryState::CTurbomachineryState() { - Density = Pressure = Entropy = Enthalpy = Temperature = TotalTemperature = TotalPressure = TotalEnthalpy = 0.0; - AbsFlowAngle = FlowAngle = MassFlow = Rothalpy = TotalRelPressure = 0.0; + SetZeroValues(); Area = Radius = 0.0; } @@ -160,51 +159,46 @@ void CPropellorBladePerformance::ComputePerformance(const CTurbomachineryCombine // TODO: to be implemented } -CTurboOutput::CTurboOutput(CConfig** config, const CGeometry& geometry, CFluidModel& fluidModel) { - unsigned short nBladesRow = config[ZONE_0]->GetnMarker_Turbomachinery(); +CTurboOutput::CTurboOutput(CConfig** config, const CGeometry& geometry, CFluidModel& fluidModel, unsigned short iBladeRow) { unsigned short nDim = geometry.GetnDim(); - for (unsigned short iBladeRow = 0; iBladeRow < nBladesRow; iBladeRow++) { - vector> bladeSpanPerformances; - unsigned short nSpan = config[iBladeRow]->GetnSpanWiseSections(); - for (unsigned short iSpan = 0; iSpan < nSpan + 1; iSpan++) { - su2double areaIn = geometry.GetSpanAreaIn(iBladeRow, iSpan); - su2double areaOut = geometry.GetSpanAreaOut(iBladeRow, iSpan); - su2double radiusIn = geometry.GetTurboRadiusIn(iBladeRow, iSpan); - su2double radiusOut = geometry.GetTurboRadiusOut(iBladeRow, iSpan); - - /* Switch between the Turbomachinery Performance Kind */ - switch (config[ZONE_0]->GetKind_TurboPerf(iBladeRow)) { - case TURBO_PERF_KIND::TURBINE: - bladeSpanPerformances.push_back( - make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); - break; - - case TURBO_PERF_KIND::COMPRESSOR: - bladeSpanPerformances.push_back( - make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); - break; - - case TURBO_PERF_KIND::PROPELLOR: - bladeSpanPerformances.push_back( - make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); - break; - - default: - bladeSpanPerformances.push_back( - make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); - break; - } + vector> bladeSpanPerformances; + unsigned short nSpan = config[iBladeRow]->GetnSpanWiseSections(); + for (unsigned short iSpan = 0; iSpan < nSpan + 1; iSpan++) { + auto areaIn = geometry.GetSpanAreaIn(iBladeRow, iSpan); + auto areaOut = geometry.GetSpanAreaOut(iBladeRow, iSpan); + auto radiusIn = geometry.GetTurboRadiusIn(iBladeRow, iSpan); + auto radiusOut = geometry.GetTurboRadiusOut(iBladeRow, iSpan); + + /* Switch between the Turbomachinery Performance Kind */ + switch (config[ZONE_0]->GetKind_TurboPerf(iBladeRow)) { + case TURBO_PERF_KIND::TURBINE: + bladeSpanPerformances.push_back( + make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); + break; + + case TURBO_PERF_KIND::COMPRESSOR: + bladeSpanPerformances.push_back( + make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); + break; + + case TURBO_PERF_KIND::PROPELLOR: + bladeSpanPerformances.push_back( + make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); + break; + + default: + bladeSpanPerformances.push_back( + make_shared(fluidModel, nDim, areaIn, radiusIn, areaOut, radiusOut)); + break; } - BladesPerformances.push_back(bladeSpanPerformances); } + BladesPerformances = bladeSpanPerformances; } void CTurboOutput::ComputeTurbomachineryPerformance( - vector> const bladesPrimitives) { - for (unsigned i = 0; i < BladesPerformances.size(); ++i) { - ComputePerBlade(BladesPerformances[i], bladesPrimitives[i]); - } + vector const bladePrimitives, unsigned short iBladeRow) { + ComputePerBlade(BladesPerformances, bladePrimitives); } void CTurboOutput::ComputePerBlade(vector> const bladePerformances, @@ -265,6 +259,8 @@ void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CT fluidModel.SetTDState_Ps(OutState.GetPressure(), InState.GetEntropy()); su2double enthalpyOutIs = fluidModel.GetStaticEnergy() + OutState.GetPressure() / fluidModel.GetDensity(); su2double totEnthalpyOutIs = enthalpyOutIs + 0.5 * OutState.GetVelocityValue() * OutState.GetVelocityValue(); + su2double tangVel = OutState.GetTangVelocity(); + su2double relVelOutIs2 = 2 * (OutState.GetRothalpy() - enthalpyOutIs) + tangVel * tangVel; /*--- Compute compressor stage performance ---*/ NormEntropyGen = (OutState.GetEntropy() - InState.GetEntropy()) / InState.GetEntropy(); @@ -273,4 +269,7 @@ void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CT TotalTotalEfficiency = (totEnthalpyOutIs - InState.GetTotalEnthalpy()) / EulerianWork; TotalStaticPressureRatio = OutState.GetPressure() / InState.GetTotalPressure(); TotalTotalPressureRatio = OutState.GetTotalPressure() / InState.GetTotalPressure(); + TotalPressureLoss = (InState.GetTotalRelPressure() - OutState.GetTotalRelPressure()) / + (InState.GetTotalRelPressure() - InState.GetPressure()); + KineticEnergyLoss = 2 * (OutState.GetEnthalpy() - enthalpyOutIs) / relVelOutIs2; } \ No newline at end of file diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 8bae5e6ca36..605d648cbdd 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -221,14 +221,14 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && config->GetBoolTurbomachinery()){ - BPressure = config->GetPressureOut_BC(); - Temperature = config->GetTotalTemperatureIn_BC(); - if (!reset){ AD::RegisterInput(BPressure); AD::RegisterInput(Temperature); } + BPressure = config->GetPressureOut_BC(); + Temperature = config->GetTotalTemperatureIn_BC(); + config->SetPressureOut_BC(BPressure); config->SetTotalTemperatureIn_BC(Temperature); } @@ -304,6 +304,24 @@ void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) { direct_solver->RegisterSolutionExtra(false, config); } +void CDiscAdjSolver::Register_VertexNormals(CGeometry *geometry, CConfig *config, bool input) { + // Only need vertex normals for boundaries declared as turbomachinery markers? + + for (auto iMarker=0ul; iMarker < nMarker; iMarker++) { + for (auto iMarkerTP = 1; iMarkerTP < config->GetnMarker_Turbomachinery() + 1; iMarkerTP++) { + if (config->GetMarker_All_Turbomachinery(iMarker) == iMarkerTP) { + for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + const auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + for (auto iDim = 0u; iDim < nDim; iDim++) { + if (input) AD::RegisterInput(Normal[iDim]); + else AD::RegisterOutput(Normal[iDim]); + } + } + } + } + } +} + void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) { const bool time_n1_needed = config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND; diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index dcca14edd19..c6984afac55 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -35,9 +35,10 @@ #include "../../include/fluid/CDataDrivenFluid.hpp" #include "../../include/fluid/CCoolProp.hpp" #include "../../include/numerics_simd/CNumericsSIMD.hpp" +#include "../../include/limiters/CLimiterDetails.hpp" +#include "../../include/output/COutput.hpp" #include "../../include/output/CTurboOutput.hpp" - CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, const bool navier_stokes) : CFVMFlowSolverBase(*geometry, *config) { @@ -358,7 +359,9 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, } CEulerSolver::~CEulerSolver() { - + for (auto& vec : MixingState) { + for (auto ptr : vec) delete [] ptr; + } for(auto& model : FluidModel) delete model; } @@ -383,37 +386,31 @@ void CEulerSolver::InstantiateEdgeNumerics(const CSolver* const* solver_containe END_SU2_OMP_SAFE_GLOBAL_ACCESS } -void CEulerSolver::InitTurboContainers(CGeometry *geometry, CConfig *config){ +void CEulerSolver::InitTurboContainers(CGeometry *geometry, CConfig **config_container, unsigned short iZone){ - /*--- Initialize quantities for the average process for internal flow ---*/ + auto config = config_container[iZone]; + /*--- Initialize quantities for the average process for internal flow ---*/ const auto nSpanWiseSections = config->GetnSpanWiseSections(); AverageVelocity.resize(nMarker); AverageTurboVelocity.resize(nMarker); OldAverageTurboVelocity.resize(nMarker); - ExtAverageTurboVelocity.resize(nMarker); AverageFlux.resize(nMarker); SpanTotalFlux.resize(nMarker); AveragePressure.resize(nMarker,nSpanWiseSections+1) = su2double(0.0); OldAveragePressure = AveragePressure; RadialEquilibriumPressure = AveragePressure; - ExtAveragePressure = AveragePressure; AverageDensity = AveragePressure; OldAverageDensity = AveragePressure; - ExtAverageDensity = AveragePressure; AverageNu = AveragePressure; AverageKine = AveragePressure; AverageOmega = AveragePressure; - ExtAverageNu = AveragePressure; - ExtAverageKine = AveragePressure; - ExtAverageOmega = AveragePressure; for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) { AverageVelocity[iMarker].resize(nSpanWiseSections+1,nDim) = su2double(0.0); AverageTurboVelocity[iMarker].resize(nSpanWiseSections+1,nDim) = su2double(0.0); OldAverageTurboVelocity[iMarker].resize(nSpanWiseSections+1,nDim) = su2double(0.0); - ExtAverageTurboVelocity[iMarker].resize(nSpanWiseSections+1,nDim) = su2double(0.0); AverageFlux[iMarker].resize(nSpanWiseSections+1,nVar) = su2double(0.0); SpanTotalFlux[iMarker].resize(nSpanWiseSections+1,nVar) = su2double(0.0); } @@ -455,6 +452,20 @@ void CEulerSolver::InitTurboContainers(CGeometry *geometry, CConfig *config){ CkOutflow2[iMarker] = CkOutflow1[iMarker]; } } + + auto nMarkerInterface = config->GetnMarker_ZoneInterface(); + + MixingState.resize(nMarkerInterface); + MixingStateNodes.resize(nMarkerInterface); + + for (unsigned long iMarker = 0; iMarker < nMarkerInterface; iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == GILES_BOUNDARY) { + MixingState[iMarker].resize(nSpanWiseSections+1); + MixingStateNodes[iMarker].resize(nSpanWiseSections+1); + } + } + + TurbomachineryPerformance = std::make_shared(config_container, *geometry, *GetFluidModel(), iZone); } void CEulerSolver::Set_MPI_ActDisk(CSolver **solver_container, CGeometry *geometry, CConfig *config) { @@ -4767,6 +4778,9 @@ void CEulerSolver::Evaluate_ObjFunc(const CConfig *config, CSolver**) { case SURFACE_MACH: Total_ComboObj+=Weight_ObjFunc*config->GetSurface_Mach(0); break; + case ENTROPY_GENERATION: + Total_ComboObj+=Weight_ObjFunc*EntropyGeneration; + break; default: break; } @@ -5567,6 +5581,18 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain for (iDim = 0; iDim < nDim; iDim++) ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; + auto nDonorSpan = GetnMixingStates(val_marker, iSpan); + su2double donorAverages[8] = {0.0}; + for (auto mixVar = 0u; mixVar < 8; mixVar++) { + donorAverages[mixVar] = GetMixingState(val_marker, iSpan, mixVar); + } + const auto ExtAverageDensity = donorAverages[0]; + const auto ExtAveragePressure = donorAverages[1]; + su2double ExtAverageTurboVelocity[3] = {0.0}; + ExtAverageTurboVelocity[0] = donorAverages[2]; + ExtAverageTurboVelocity[1] = donorAverages[3]; + ExtAverageTurboVelocity[2] = donorAverages[4]; + /*--- Build the external state u_e from boundary data and internal node ---*/ switch(config->GetKind_Data_Riemann(Marker_Tag)) @@ -5607,15 +5633,15 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain case MIXING_IN: /* --- compute total averaged quantities ---*/ - GetFluidModel()->SetTDState_Prho(ExtAveragePressure[val_marker][iSpan], ExtAverageDensity[val_marker][iSpan]); - AverageEnthalpy = GetFluidModel()->GetStaticEnergy() + ExtAveragePressure[val_marker][iSpan]/ExtAverageDensity[val_marker][iSpan]; + GetFluidModel()->SetTDState_Prho(ExtAveragePressure, ExtAverageDensity); + AverageEnthalpy = GetFluidModel()->GetStaticEnergy() + ExtAveragePressure/ExtAverageDensity; AverageEntropy = GetFluidModel()->GetEntropy(); FlowDirMixMag = 0; for (iDim = 0; iDim < nDim; iDim++) - FlowDirMixMag += ExtAverageTurboVelocity[val_marker][iSpan][iDim]*ExtAverageTurboVelocity[val_marker][iSpan][iDim]; + FlowDirMixMag += ExtAverageTurboVelocity[iDim]*ExtAverageTurboVelocity[iDim]; for (iDim = 0; iDim < nDim; iDim++){ - FlowDirMix[iDim] = ExtAverageTurboVelocity[val_marker][iSpan][iDim]/sqrt(FlowDirMixMag); + FlowDirMix[iDim] = ExtAverageTurboVelocity[iDim]/sqrt(FlowDirMixMag); } @@ -5643,7 +5669,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain case MIXING_OUT: /*--- Retrieve the static pressure for this boundary. ---*/ - Pressure_e = ExtAveragePressure[val_marker][iSpan]; + Pressure_e = ExtAveragePressure; Density_e = Density_i; /* --- Compute the boundary state u_e --- */ @@ -6214,7 +6240,6 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu deltaSpan = SpanWiseValues[nSpanWiseSections-1]*spanPercent; coeffrelfacAvg = (relfacAvgCfg - extrarelfacAvg)/deltaSpan; } - for (iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ /*--- Compute under relaxation for the Hub and Shroud Avg and Fourier Coefficient---*/ if(nDim == 3){ @@ -6253,6 +6278,27 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu kend_max = geometry->GetnFreqSpanMax(config->GetMarker_All_TurbomachineryFlag(val_marker)); conv_numerics->GetRMatrix(AverageSoundSpeed, AverageDensity[val_marker][iSpan], R_Matrix); + su2double ExtAverageDensity, ExtAveragePressure; + su2double ExtAverageTurboVelocity[3] = {0.0}; + su2double donorAverages[8] = {0.0}; + switch (config->GetKind_Data_Giles(Marker_Tag)){ + case MIXING_IN: case MIXING_IN_1D: case MIXING_OUT: case MIXING_OUT_1D: + // for (auto iMarkerGeo = 0u; iMarkeGeo < geometry->GetnMarker(); iMarkerGeo++){ + // if (config->Get) + // } + for (auto mixVar = 0u; mixVar < 8; mixVar++) { + donorAverages[mixVar] = GetMixingState(val_marker, iSpan, mixVar); + } + ExtAverageDensity = donorAverages[0]; + ExtAveragePressure = donorAverages[1]; + ExtAverageTurboVelocity[0] = donorAverages[2]; + ExtAverageTurboVelocity[1] = donorAverages[3]; + ExtAverageTurboVelocity[2] = donorAverages[4]; + break; + default: + break; + } + switch(config->GetKind_Data_Giles(Marker_Tag)){ case TOTAL_CONDITIONS_PT: @@ -6370,16 +6416,16 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu case MIXING_IN: case MIXING_OUT: /* --- Compute average jump of primitive at the mixing-plane interface--- */ - deltaprim[0] = ExtAverageDensity[val_marker][iSpan] - AverageDensity[val_marker][iSpan]; - deltaprim[1] = ExtAverageTurboVelocity[val_marker][iSpan][0] - AverageTurboVelocity[val_marker][iSpan][0]; - deltaprim[2] = ExtAverageTurboVelocity[val_marker][iSpan][1] - AverageTurboVelocity[val_marker][iSpan][1]; + deltaprim[0] = ExtAverageDensity - AverageDensity[val_marker][iSpan]; + deltaprim[1] = ExtAverageTurboVelocity[0] - AverageTurboVelocity[val_marker][iSpan][0]; + deltaprim[2] = ExtAverageTurboVelocity[1] - AverageTurboVelocity[val_marker][iSpan][1]; if (nDim == 2){ - deltaprim[3] = ExtAveragePressure[val_marker][iSpan] - AveragePressure[val_marker][iSpan]; + deltaprim[3] = ExtAveragePressure - AveragePressure[val_marker][iSpan]; } else { - deltaprim[3] = ExtAverageTurboVelocity[val_marker][iSpan][2] - AverageTurboVelocity[val_marker][iSpan][2]; - deltaprim[4] = ExtAveragePressure[val_marker][iSpan] - AveragePressure[val_marker][iSpan]; + deltaprim[3] = ExtAverageTurboVelocity[2] - AverageTurboVelocity[val_marker][iSpan][2]; + deltaprim[4] = ExtAveragePressure - AveragePressure[val_marker][iSpan]; } @@ -6392,16 +6438,16 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu case MIXING_IN_1D: case MIXING_OUT_1D: /* --- Compute average jump of primitive at the mixing-plane interface--- */ - deltaprim[0] = ExtAverageDensity[val_marker][nSpanWiseSections] - AverageDensity[val_marker][nSpanWiseSections]; - deltaprim[1] = ExtAverageTurboVelocity[val_marker][nSpanWiseSections][0] - AverageTurboVelocity[val_marker][nSpanWiseSections][0]; - deltaprim[2] = ExtAverageTurboVelocity[val_marker][nSpanWiseSections][1] - AverageTurboVelocity[val_marker][nSpanWiseSections][1]; + deltaprim[0] = ExtAverageDensity - AverageDensity[val_marker][nSpanWiseSections]; + deltaprim[1] = ExtAverageTurboVelocity[0] - AverageTurboVelocity[val_marker][nSpanWiseSections][0]; + deltaprim[2] = ExtAverageTurboVelocity[1] - AverageTurboVelocity[val_marker][nSpanWiseSections][1]; if (nDim == 2){ - deltaprim[3] = ExtAveragePressure[val_marker][nSpanWiseSections] - AveragePressure[val_marker][nSpanWiseSections]; + deltaprim[3] = ExtAveragePressure - AveragePressure[val_marker][nSpanWiseSections]; } else { - deltaprim[3] = ExtAverageTurboVelocity[val_marker][nSpanWiseSections][2] - AverageTurboVelocity[val_marker][nSpanWiseSections][2]; - deltaprim[4] = ExtAveragePressure[val_marker][nSpanWiseSections] - AveragePressure[val_marker][nSpanWiseSections]; + deltaprim[3] = ExtAverageTurboVelocity[2] - AverageTurboVelocity[val_marker][nSpanWiseSections][2]; + deltaprim[4] = ExtAveragePressure - AveragePressure[val_marker][nSpanWiseSections]; } /* --- Compute average jump of charachteristic variable at the mixing-plane interface--- */ @@ -8801,11 +8847,16 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon const auto nSpanWiseSections = config->GetnSpanWiseSections(); const auto iZone = config->GetiZone(); + const bool spalart_allmaras = (config->GetKind_Turb_Model() == TURB_MODEL::SA); + const bool menter_sst = (config->GetKind_Turb_Model() == TURB_MODEL::SST); for (auto iSpan= 0u; iSpan < nSpanWiseSections; iSpan++){ su2double TotalAreaVelocity[MAXNDIM]={0.0}, TotalAreaPressure{0}, - TotalAreaDensity{0}; + TotalAreaDensity{0}, + TotalAreaNu{0}, + TotalAreaKine{0}, + TotalAreaOmega{0}; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++){ for (auto iMarkerTP=1; iMarkerTP < config->GetnMarker_Turbomachinery()+1; iMarkerTP++){ if (config->GetMarker_All_Turbomachinery(iMarker) == iMarkerTP){ @@ -8822,6 +8873,16 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon auto Pressure = nodes->GetPressure(iPoint); auto Density = nodes->GetDensity(iPoint); + /*--- This is in Euler, however we also need to average the turbulent variables so we do it here too ---*/ + su2double Kine{0}, Omega{0}, Nu{0}; + if(menter_sst){ + Kine = solver[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + Omega = solver[TURB_SOL]->GetNodes()->GetSolution(iPoint,1); + } + if(spalart_allmaras){ + Nu = solver[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + } + su2double UnitNormal[MAXNDIM]={0}, TurboNormal[MAXNDIM]={0}, TurboVelocity[MAXNDIM], @@ -8842,6 +8903,9 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon TotalAreaDensity += Area*Density; for (auto iDim = 0u; iDim < nDim; iDim++) TotalAreaVelocity[iDim] += Area*Velocity[iDim]; + TotalAreaNu += Area*Nu; + TotalAreaKine += Area*Kine; + TotalAreaOmega += Area*Omega; } } } @@ -8855,9 +8919,15 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon su2double MyTotalAreaDensity = TotalAreaDensity; su2double MyTotalAreaPressure = TotalAreaPressure; + su2double MyTotalAreaNu = TotalAreaNu; + su2double MyTotalAreaKine = TotalAreaKine; + su2double MyTotalAreaOmega = TotalAreaOmega; SU2_MPI::Allreduce(&MyTotalAreaDensity, &TotalAreaDensity, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); SU2_MPI::Allreduce(&MyTotalAreaPressure, &TotalAreaPressure, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&MyTotalAreaNu, &TotalAreaNu, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&MyTotalAreaKine, &TotalAreaKine, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&MyTotalAreaOmega, &TotalAreaOmega, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); auto* MyTotalAreaVelocity = new su2double[nDim]; @@ -8889,6 +8959,10 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon for (auto iDim = 0u; iDim < nDim; iDim++) AverageVelocity[iMarker][iSpan][iDim] = TotalAreaVelocity[iDim] / TotalArea; + AverageNu[iMarker][iSpan] = TotalAreaNu / TotalArea; + AverageKine[iMarker][iSpan] = TotalAreaKine / TotalArea; + AverageOmega[iMarker][iSpan] = TotalAreaOmega / TotalArea; + /* --- compute static averaged quantities ---*/ ComputeTurboVelocity(AverageVelocity[iMarker][iSpan], AverageTurboNormal , AverageTurboVelocity[iMarker][iSpan], marker_flag, config->GetKind_TurboMachinery(iZone)); @@ -8919,6 +8993,10 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon for (auto iDim = 0u; iDim < nDim; iDim++) AverageVelocity[iMarker][nSpanWiseSections][iDim] = AverageVelocity[iMarker][nSpanWiseSections/2][iDim]; + AverageNu[iMarker][nSpanWiseSections] = AverageNu[iMarker][nSpanWiseSections/2]; + AverageKine[iMarker][nSpanWiseSections] = AverageKine[iMarker][nSpanWiseSections/2]; + AverageOmega[iMarker][nSpanWiseSections] = AverageOmega[iMarker][nSpanWiseSections/2]; + /* --- compute static averaged quantities ---*/ ComputeTurboVelocity(AverageVelocity[iMarker][nSpanWiseSections], AverageTurboNormal , AverageTurboVelocity[iMarker][nSpanWiseSections], marker_flag, config->GetKind_TurboMachinery(iZone)); @@ -9539,3 +9617,33 @@ void CEulerSolver::GatherInOutAverageValues(CConfig *config, CGeometry *geometry } } } + +void CEulerSolver::ComputeTurboBladePerformance(CGeometry* geometry, CConfig* config, unsigned short iBlade) { + // Computes the turboperformance per blade in zone iBlade + const auto nDim = geometry->GetnDim(); + vector TurboPrimitiveIn, TurboPrimitiveOut; + if (rank == MASTER_NODE) { + /* Blade Primitive initialized per blade */ + std::vector bladePrimitives; + auto nSpan = config->GetnSpanWiseSections(); + for (auto iSpan = 0; iSpan < nSpan + 1; iSpan++) { + TurboPrimitiveIn= GetTurboPrimitive(iBlade, iSpan, true); + TurboPrimitiveOut= GetTurboPrimitive(iBlade, iSpan, false); + auto spanInletPrimitive = CTurbomachineryPrimitiveState(TurboPrimitiveIn, nDim, geometry->GetTangGridVelIn(iBlade, iSpan)); + auto spanOutletPrimitive = CTurbomachineryPrimitiveState(TurboPrimitiveOut, nDim, geometry->GetTangGridVelOut(iBlade, iSpan)); + auto spanCombinedPrimitive = CTurbomachineryCombinedPrimitiveStates(spanInletPrimitive, spanOutletPrimitive); + bladePrimitives.push_back(spanCombinedPrimitive); + } + TurbomachineryPerformance->ComputeTurbomachineryPerformance(bladePrimitives, iBlade); + + auto BladePerf = TurbomachineryPerformance->GetBladesPerformances().at(nSpan); + + EntropyGeneration = BladePerf->GetEntropyGen(); + TotalPressureLoss = BladePerf->GetTotalPressureLoss(); + KineticEnergyLoss = BladePerf->GetKineticEnergyLoss(); + } +} + +// void CEulerSolver::RegisterSolutionExtra(bool input, CConfig* config) { +// nodes->RegisterGridVelocity(); +// } \ No newline at end of file diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 7649165d0be..a1e17b0d276 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -1075,7 +1075,10 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c /*--- Loop over all the vertices on this boundary marker ---*/ for (auto iSpan = 0u; iSpan < nSpanWiseSections ; iSpan++){ - su2double extAverageNu = solver_container[FLOW_SOL]->GetExtAverageNu(val_marker, iSpan); + auto nDonorSpan = solver_container[FLOW_SOL]->GetnMixingStates(val_marker, iSpan); + const auto extAverageNu = solver_container[FLOW_SOL]->GetMixingState(val_marker, iSpan, 5); + + // su2double extAverageNu = solver_container[FLOW_SOL]->GetExtAverageNu(val_marker, iSpan); /*--- Loop over all the vertices on this boundary marker ---*/ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index b36f5c2b876..3b18c17b49f 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -803,8 +803,9 @@ void CTurbSSTSolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_ for (auto iSpan = 0u; iSpan < nSpanWiseSections ; iSpan++){ - su2double extAverageKine = solver_container[FLOW_SOL]->GetExtAverageKine(val_marker, iSpan); - su2double extAverageOmega = solver_container[FLOW_SOL]->GetExtAverageOmega(val_marker, iSpan); + auto nDonorSpan = solver_container[FLOW_SOL]->GetnMixingStates(val_marker, iSpan); + const auto extAverageKine = solver_container[FLOW_SOL]->GetMixingState(val_marker, iSpan, 6); + const auto extAverageOmega = solver_container[FLOW_SOL]->GetMixingState(val_marker, iSpan, 7); su2double solution_j[] = {extAverageKine, extAverageOmega}; /*--- Loop over all the vertices on this boundary marker ---*/ diff --git a/TestCases/disc_adj_turbomachinery/axial_stage_2D/Axial_stage2D.cfg b/TestCases/disc_adj_turbomachinery/axial_stage_2D/Axial_stage2D.cfg new file mode 100755 index 00000000000..d87f092e6a7 --- /dev/null +++ b/TestCases/disc_adj_turbomachinery/axial_stage_2D/Axial_stage2D.cfg @@ -0,0 +1,211 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D Discrete Adjoint Axial stage % +% Author: J. Kelly % +% Institution: University of Liverpool. % +% Date: Oct 11th, 2025 % +% File Version 8.3.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +MULTIZONE= YES +CONFIG_LIST= (zone_1.cfg, zone_2.cfg) +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= RANS +KIND_TURB_MODEL= SST +READ_BINARY_RESTART=NO +RESTART_SOL= NO +% +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 0.05 +AOA= 0.0 +FREESTREAM_PRESSURE= 140000.0 +FREESTREAM_TEMPERATURE= 300.0 +FREESTREAM_DENSITY= 1.7418 +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_TURBULENCEINTENSITY = 0.03 +FREESTREAM_TURB2LAMVISCRATIO = 100.0 +INIT_OPTION= TD_CONDITIONS +% +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.00 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 +REF_DIMENSIONALIZATION= DIMENSIONAL +% +% ------------------------------ EQUATION OF STATE ----------------------------% +% +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.058 +% +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_REF= 1.716E-5 +MU_T_REF= 273.15 +SUTHERLAND_CONSTANT= 110.4 +% +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall1, 0.0, wall2, 0.0) +% +MARKER_PERIODIC= ( periodic1, periodic2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04463756775, 0.0, periodic3, periodic4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04463756775, 0.0) +% +%-------- INFLOW/OUTFLOW BOUNDARY CONDITION SPECIFIC FOR TURBOMACHINERY --------% +% +% Inflow and Outflow markers must be specified, for each blade (zone), following the natural groth of the machine (i.e, from the first blade to the last) +MARKER_TURBOMACHINERY= (inflow, outmix, inmix, outflow) +% +MARKER_ZONE_INTERFACE= (outmix, inmix) +% Mixing-plane interface markers must be specified to activate the transfer of information between zones +MARKER_MIXINGPLANE_INTERFACE= (outmix, inmix) +% +MARKER_GILES= (inflow, TOTAL_CONDITIONS_PT, 169623.33, 305.76, 1.0, 0.0, 0.0,1.0,1.0, outmix, MIXING_OUT, 0.0, 0.0, 0.0, 0.0, 0.0,1.0,1.0, inmix, MIXING_IN, 0.0, 0.0, 0.0, 0.0, 0.0,1.0, 1.0 outflow, STATIC_PRESSURE, 99741.00, 0.0, 0.0, 0.0, 0.0,1.0,1.0) +% +SPATIAL_FOURIER= YES +% +% +%---------------------------- TURBOMACHINERY SIMULATION -----------------------------% +% +TURBOMACHINERY_KIND= AXIAL AXIAL +TURBO_PERF_KIND = TURBINE TURBINE +TURBULENT_MIXINGPLANE= YES +RAMP_OUTLET_PRESSURE= NO +RAMP_OUTLET_PRESSURE_COEFF= (140000.0, 10.0, 2000) +AVERAGE_PROCESS_KIND= MIXEDOUT +PERFORMANCE_AVERAGE_PROCESS_KIND= MIXEDOUT +MIXEDOUT_COEFF= (1.0, 1.0E-05, 15) +AVERAGE_MACH_LIMIT= 0.05 +% +% ------------------------ SURFACES IDENTIFICATION ----------------------------% +% +MARKER_PLOTTING= wall1, wall2 +MARKER_MONITORING= wall1, wall2 +MARKER_DESIGNING= wall1, wall2 +MARKER_ANALYZE= wall1, wall2 +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 10.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.3, 1.2, 1.0, 10.0) +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 1E-4 +LINEAR_SOLVER_ITER= 20 +% +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +% +VENKAT_LIMITER_COEFF= 0.05 +LIMITER_ITER= 999999 +% +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= VAN_ALBADA_EDGE +ENTROPY_FIX_COEFF= 0.1 +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= NO +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +TIME_DISCRE_TURB= EULER_IMPLICIT +CFL_REDUCTION_TURB= 1.0 +% +% ----------------------- DESIGN VARIABLE PARAMETERS --------------------------% +% +DV_KIND= FFD_CONTROL_POINT_2D, FFD_CONTROL_POINT_2D +DV_MARKER= ( wall1, wall2 ) +DV_PARAM=( STATOR, 0, 0, 1.0, 0.0 ); ( ROTOR, 2, 2, 1.0, 0.0 ) +DV_VALUE=0.0,0.0 +% +% FFD_CONTROL_POINT_2D +DEFINITION_DV= ( 19, 1.0 | wall1 | STATOR, 0, 0, 1.0, 0.0 ); ( 19, 1.0 | wall2 | ROTOR, 2, 2, 1.0, 0.0 ); +% +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, JACOBI) +DEFORM_LINEAR_SOLVER_PREC= ILU +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_SOLVER_ITER= 1000 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Minimum residual criteria for the linear solver convergence of grid deformation +DEFORM_LINEAR_SOLVER_ERROR= 1E-14 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% -------------------- FREE-FORM DEFORMATION PARAMETERS -----------------------% +% +FFD_TOLERANCE= 1E-10 +FFD_ITERATIONS= 500 +FFD_DEFINITION= (STATOR, 0.139, -0.0038, 0.0, 0.186, -0.0567, 0.0, 0.193, -0.031, 0.0, 0.156, 0.010, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); (ROTOR, 0.200, -0.040, 0.0, 0.259, -0.001, 0.0, 0.259, 0.032, 0.0, 0.200, -0.007, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) +FFD_DEGREE= ( 2 , 2 , 0 ); (2, 2, 0) +FFD_CONTINUITY= 2ND_DERIVATIVE +% +% --------------------- OBJECTIVE FUNCTION DEFINITION -----------------------% +% +OBJECTIVE_FUNCTION= ENTROPY_GENERATION +% +% ------------------------- CONVERGENCE PARAMETERS --------------------------% +% +OUTER_ITER= 101 +CONV_RESIDUAL_MINVAL= -16 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], RMS_DENSITY[0], RMS_ENERGY[0], RMS_DENSITY[1], RMS_ENERGY[1], COMBO +HISTORY_OUTPUT= COMBO +MESH_FILENAME= axial_stage_2D.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out.su2 +SOLUTION_FILENAME= solution_flow +SOLUTION_ADJ_FILENAME= solution_adj +TABULAR_FORMAT= CSV +OUTPUT_FILES= RESTART_ASCII, TECPLOT, SURFACE_TECPLOT +CONV_FILENAME= history +RESTART_FILENAME= restart_flow +RESTART_ADJ_FILENAME= restart_adj +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad.dat +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +SURFACE_SENS_FILENAME= surface_sens +WRT_ZONE_CONV= NO +WRT_ZONE_HIST= YES +OUTPUT_PRECISION= 16 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index c9ffc63fef6..0e41a02c3ce 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -234,6 +234,15 @@ def main(): discadj_trans_stator.test_vals_aarch64 = [79.000000, 0.696755, 0.485950, 0.569475, -0.990065] test_list.append(discadj_trans_stator) + # Axial stage 2D + discadj_axial_stage = TestCase('axial_stage_2D') + discadj_axial_stage.cfg_dir = "disc_adj_turbomachinery/axial_stage_2D" + discadj_axial_stage.cfg_file = "Axial_stage2D.cfg" + discadj_axial_stage.test_iter = 79 + discadj_axial_stage.test_vals = [79, 0.668058, 0.483608, 0.518789, -1.013227] + discadj_axial_stage.test_vals_aarch64 = [79, 0.668058, 0.483608, 0.518789, -1.013227] + test_list.append(discadj_axial_stage) + ################################### ### Structural Adjoint ### ###################################