Skip to content

Commit 90c5a10

Browse files
TobiKattmannTobiKattmann
committed
Compute Species Variance now in a seperate function.
And not inflate the SetAnalyzeSurface funtion even more.
1 parent 3661da5 commit 90c5a10

3 files changed

Lines changed: 145 additions & 99 deletions

File tree

‎SU2_CFD/include/output/CFlowOutput.hpp‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,20 @@ class CFlowOutput : public CFVMOutput{
6767
*/
6868
void SetAnalyzeSurface(const CSolver* const* solver, const CGeometry *geometry, CConfig *config, bool output);
6969

70+
/*!
71+
* \brief Compute and Set flow species variance output field values
72+
* \param[in] solver - The container holding all solution data.
73+
* \param[in] geometry - Geometrical definition of the problem.
74+
* \param[in,out] config - Definition of the particular problem.
75+
* \param[in] Surface_Species_Total - Avg mass fraction of each species on all Marker_Analyze
76+
* \param[in] Surface_MassFlow_Abs_Total - Massflow on all Marker_Analyze
77+
* \param[in] Surface_Area_Total - Area of all Marker_Analyze
78+
*/
79+
void SetAnalyzeSurface_SpeciesVariance(const CSolver* const*solver, const CGeometry *geometry, CConfig *config,
80+
const su2activematrix Surface_Species_Total,
81+
const vector<su2double> Surface_MassFlow_Abs_Total,
82+
const vector<su2double> Surface_Area_Total);
83+
7084
/*!
7185
* \brief Add scalar (turbulence/species) history fields for the linear solver (FVMComp, FVMInc, FVMNEMO).
7286
*/

‎SU2_CFD/src/output/CFlowOutput.cpp‎

Lines changed: 130 additions & 98 deletions
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,7 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){
7575
AddHistoryOutput("SURFACE_SPECIES_" + std::to_string(iVar), "Avg_Species_" + std::to_string(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF", "Total average species " + std::to_string(iVar) + " on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT);
7676
}
7777
/// DESCRIPTION: Average Species
78-
AddHistoryOutput("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::FIXED, "SPECIES_COEFF", "Total species variance, measure for mixing quality. On all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT);
78+
AddHistoryOutput("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF", "Total species variance, measure for mixing quality. On all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT);
7979
/// END_GROUP
8080

8181
/// BEGIN_GROUP: AERO_COEFF_SURF, DESCRIPTION: Surface values on non-solid markers.
@@ -115,7 +115,7 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){
115115
AddHistoryOutputPerSurface("SURFACE_SPECIES_" + std::to_string(iVar), "Avg_Species_" + std::to_string(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT);
116116
}
117117
/// DESCRIPTION: Average Species
118-
AddHistoryOutputPerSurface("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::FIXED, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT);
118+
AddHistoryOutputPerSurface("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT);
119119
/// END_GROUP
120120
}
121121

@@ -544,28 +544,134 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry
544544
for (unsigned short iVar = 0; iVar < nSpecies; iVar++)
545545
SetHistoryOutputValue("SURFACE_SPECIES_" + std::to_string(iVar), Tot_Surface_Species[iVar]);
546546

547+
if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE)
548+
SetAnalyzeSurface_SpeciesVariance(solver, geometry, config, Surface_Species_Total, Surface_MassFlow_Abs_Total,
549+
Surface_Area_Total);
550+
551+
if ((rank == MASTER_NODE) && !config->GetDiscrete_Adjoint() && output) {
552+
553+
cout.precision(6);
554+
cout.setf(ios::scientific, ios::floatfield);
555+
cout << endl << "Computing surface mean values." << endl << endl;
556+
557+
for (iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
558+
cout << "Surface "<< config->GetMarker_Analyze_TagBound(iMarker_Analyze) << ":" << endl;
559+
560+
if (nDim == 3) { if (si_units) cout << setw(20) << "Area (m^2): "; else cout << setw(20) << "Area (ft^2): "; }
561+
else { if (si_units) cout << setw(20) << "Area (m): "; else cout << setw(20) << "Area (ft): "; }
562+
563+
if (si_units) cout << setw(15) << fabs(Surface_Area_Total[iMarker_Analyze]);
564+
else if (us_units) cout << setw(15) << fabs(Surface_Area_Total[iMarker_Analyze])*12.0*12.0;
565+
566+
cout << endl;
567+
568+
su2double MassFlow = config->GetSurface_MassFlow(iMarker_Analyze);
569+
if (si_units) cout << setw(20) << "Mf (kg/s): " << setw(15) << MassFlow;
570+
else if (us_units) cout << setw(20) << "Mf (lbs/s): " << setw(15) << MassFlow;
571+
572+
su2double NormalVelocity = config->GetSurface_NormalVelocity(iMarker_Analyze);
573+
if (si_units) cout << setw(20) << "Vn (m/s): " << setw(15) << NormalVelocity;
574+
else if (us_units) cout << setw(20) << "Vn (ft/s): " << setw(15) << NormalVelocity;
575+
576+
cout << endl;
577+
578+
su2double Uniformity = config->GetSurface_Uniformity(iMarker_Analyze);
579+
if (si_units) cout << setw(20) << "Uniformity (m/s): " << setw(15) << Uniformity;
580+
else if (us_units) cout << setw(20) << "Uniformity (ft/s): " << setw(15) << Uniformity;
581+
582+
su2double SecondaryStrength = config->GetSurface_SecondaryStrength(iMarker_Analyze);
583+
if (si_units) cout << setw(20) << "Secondary (m/s): " << setw(15) << SecondaryStrength;
584+
else if (us_units) cout << setw(20) << "Secondary (ft/s): " << setw(15) << SecondaryStrength;
585+
586+
cout << endl;
587+
588+
su2double MomentumDistortion = config->GetSurface_MomentumDistortion(iMarker_Analyze);
589+
cout << setw(20) << "Mom. Distortion: " << setw(15) << MomentumDistortion;
590+
591+
su2double SecondOverUniform = config->GetSurface_SecondOverUniform(iMarker_Analyze);
592+
cout << setw(20) << "Second/Uniform: " << setw(15) << SecondOverUniform;
593+
594+
cout << endl;
595+
596+
su2double Pressure = config->GetSurface_Pressure(iMarker_Analyze);
597+
if (si_units) cout << setw(20) << "P (Pa): " << setw(15) << Pressure;
598+
else if (us_units) cout << setw(20) << "P (psf): " << setw(15) << Pressure;
599+
600+
su2double TotalPressure = config->GetSurface_TotalPressure(iMarker_Analyze);
601+
if (si_units) cout << setw(20) << "PT (Pa): " << setw(15) <<TotalPressure;
602+
else if (us_units) cout << setw(20) << "PT (psf): " << setw(15) <<TotalPressure;
603+
604+
cout << endl;
605+
606+
su2double Mach = config->GetSurface_Mach(iMarker_Analyze);
607+
cout << setw(20) << "Mach: " << setw(15) << Mach;
608+
609+
su2double Density = config->GetSurface_Density(iMarker_Analyze);
610+
if (si_units) cout << setw(20) << "Rho (kg/m^3): " << setw(15) << Density;
611+
else if (us_units) cout << setw(20) << "Rho (lb/ft^3): " << setw(15) << Density*32.174;
612+
613+
cout << endl;
614+
615+
if (compressible || energy) {
616+
su2double Temperature = config->GetSurface_Temperature(iMarker_Analyze);
617+
if (si_units) cout << setw(20) << "T (K): " << setw(15) << Temperature;
618+
else if (us_units) cout << setw(20) << "T (R): " << setw(15) << Temperature;
619+
620+
su2double TotalTemperature = config->GetSurface_TotalTemperature(iMarker_Analyze);
621+
if (si_units) cout << setw(20) << "TT (K): " << setw(15) << TotalTemperature;
622+
else if (us_units) cout << setw(20) << "TT (R): " << setw(15) << TotalTemperature;
623+
624+
cout << endl;
625+
}
626+
627+
}
628+
cout.unsetf(ios_base::floatfield);
629+
630+
}
631+
632+
std::cout << std::resetiosflags(std::cout.flags());
633+
}
634+
635+
void CFlowOutput::SetAnalyzeSurface_SpeciesVariance(const CSolver* const*solver, const CGeometry *geometry, CConfig *config, su2activematrix Surface_Species_Total, vector<su2double> Surface_MassFlow_Abs_Total, vector<su2double> Surface_Area_Total) {
636+
637+
const unsigned short nMarker = config->GetnMarker_All();
638+
const unsigned short nDim = geometry->GetnDim();
639+
const unsigned short Kind_Average = config->GetKind_Average();
640+
641+
const bool species = config->GetKind_Species_Model() != SPECIES_MODEL::NONE;
642+
const auto nSpecies = config->GetnSpecies();
643+
644+
const bool axisymmetric = config->GetAxisymmetric();
645+
const unsigned short nMarker_Analyze = config->GetnMarker_Analyze();
646+
647+
const auto flow_nodes = solver[FLOW_SOL]->GetNodes();
648+
const CVariable* species_nodes = species ? solver[SPECIES_SOL]->GetNodes() : nullptr;
649+
547650
/*--- Compute Variance of species on the analyze markers. This is done after the rest as the average species value is
548651
* necessary. The variance is computed for all species together and not for each species alone. ---*/
549652
vector<su2double> Surface_SpeciesVariance(nMarker,0.0);
550653
su2double Tot_Surface_SpeciesVariance = 0.0;
551654

552655
/*--- sum += (Yj_i - mu_Yj)^2 * weight_i with i representing the node and j the species. ---*/
553-
for (iMarker = 0; iMarker < nMarker; iMarker++) {
656+
for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) {
554657

555658
if (config->GetMarker_All_Analyze(iMarker) == YES) {
556659

557660
/*--- Find iMarkerAnalyze to iMarker. As SpeciesAvg is accessed via iMarkerAnalyze. ---*/
558-
for (iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++)
661+
unsigned short iMarker_Analyze_Stored = std::numeric_limits<unsigned short>::max();
662+
for (unsigned short iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++)
559663
if (config->GetMarker_All_TagBound(iMarker) == config->GetMarker_Analyze_TagBound(iMarker_Analyze))
560-
break;
664+
iMarker_Analyze_Stored = iMarker_Analyze;
561665

562-
for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
563-
iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
666+
for (unsigned long iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
667+
const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
564668

565669
if (geometry->nodes->GetDomain(iPoint)) {
566670

671+
su2double Vector[3];
567672
geometry->vertex[iMarker][iVertex]->GetNormal(Vector);
568673

674+
su2double AxiFactor = 0.0;
569675
if (axisymmetric) {
570676
if (geometry->nodes->GetCoord(iPoint, 1) != 0.0)
571677
AxiFactor = 2.0*PI_NUMBER*geometry->nodes->GetCoord(iPoint, 1);
@@ -585,21 +691,25 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry
585691
AxiFactor = 1.0;
586692
}
587693

588-
Density = flow_nodes->GetDensity(iPoint);
589-
Area = 0.0; MassFlow = 0.0;
694+
const su2double Density = flow_nodes->GetDensity(iPoint);
695+
su2double Area = 0.0;
696+
su2double MassFlow = 0.0;
697+
su2double Velocity[3];
590698

591-
for (iDim = 0; iDim < nDim; iDim++) {
699+
for (unsigned short iDim = 0; iDim < nDim; iDim++) {
592700
Area += (Vector[iDim] * AxiFactor) * (Vector[iDim] * AxiFactor);
701+
Velocity[iDim] = flow_nodes->GetVelocity(iPoint,iDim);
593702
MassFlow += Vector[iDim] * AxiFactor * Density * Velocity[iDim];
594703
}
595704
Area= sqrt(Area);
596705

706+
su2double Weight;
597707
if (Kind_Average == AVERAGE_MASSFLUX) Weight = abs(MassFlow);
598708
else if (Kind_Average == AVERAGE_AREA) Weight = abs(Area);
599709
else Weight = 1.0;
600710

601711
for (unsigned short iVar = 0; iVar < nSpecies; iVar++)
602-
Surface_SpeciesVariance[iMarker] += pow(species_nodes->GetSolution(iPoint, iVar) - Surface_Species_Total(iMarker_Analyze, iVar), 2) * Weight;
712+
Surface_SpeciesVariance[iMarker] += pow(species_nodes->GetSolution(iPoint, iVar) - Surface_Species_Total(iMarker_Analyze_Stored, iVar), 2) * Weight;
603713
}
604714
}
605715
}
@@ -609,11 +719,11 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry
609719
vector<su2double> Surface_SpeciesVariance_Local(nMarker_Analyze,0.0);
610720
vector<su2double> Surface_SpeciesVariance_Total(nMarker_Analyze,0.0);
611721

612-
for (iMarker = 0; iMarker < nMarker; iMarker++) {
722+
for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) {
613723

614724
if (config->GetMarker_All_Analyze(iMarker) == YES) {
615725

616-
for (iMarker_Analyze= 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
726+
for (unsigned short iMarker_Analyze= 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
617727

618728
/*--- Add the Surface_MassFlow, and Surface_Area to the particular boundary ---*/
619729

@@ -623,11 +733,16 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry
623733
}
624734
}
625735
}
736+
737+
auto Allreduce = [](const vector<su2double>& src, vector<su2double>& dst) {
738+
SU2_MPI::Allreduce(src.data(), dst.data(), src.size(), MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
739+
};
626740
Allreduce(Surface_SpeciesVariance_Local, Surface_SpeciesVariance_Total);
627741

628742
/*--- Divide quantity by weight. ---*/
629-
for (iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
743+
for (unsigned short iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
630744

745+
su2double Weight;
631746
if (Kind_Average == AVERAGE_MASSFLUX) Weight = Surface_MassFlow_Abs_Total[iMarker_Analyze];
632747
else if (Kind_Average == AVERAGE_AREA) Weight = abs(Surface_Area_Total[iMarker_Analyze]);
633748
else Weight = 1.0;
@@ -641,96 +756,13 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry
641756
}
642757

643758
/*--- Set values on markers ---*/
644-
for (iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
759+
for (unsigned short iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
645760
su2double SpeciesVariance = Surface_SpeciesVariance_Total[iMarker_Analyze];
646761
SetHistoryOutputPerSurfaceValue("SURFACE_SPECIES_VARIANCE", SpeciesVariance, iMarker_Analyze);
647762
Tot_Surface_SpeciesVariance += SpeciesVariance;
648763
// Set value into config. Necessary to access as an OF.
649764
}
650765
SetHistoryOutputValue("SURFACE_SPECIES_VARIANCE", Tot_Surface_SpeciesVariance);
651-
652-
if ((rank == MASTER_NODE) && !config->GetDiscrete_Adjoint() && output) {
653-
654-
cout.precision(6);
655-
cout.setf(ios::scientific, ios::floatfield);
656-
cout << endl << "Computing surface mean values." << endl << endl;
657-
658-
for (iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
659-
cout << "Surface "<< config->GetMarker_Analyze_TagBound(iMarker_Analyze) << ":" << endl;
660-
661-
if (nDim == 3) { if (si_units) cout << setw(20) << "Area (m^2): "; else cout << setw(20) << "Area (ft^2): "; }
662-
else { if (si_units) cout << setw(20) << "Area (m): "; else cout << setw(20) << "Area (ft): "; }
663-
664-
if (si_units) cout << setw(15) << fabs(Surface_Area_Total[iMarker_Analyze]);
665-
else if (us_units) cout << setw(15) << fabs(Surface_Area_Total[iMarker_Analyze])*12.0*12.0;
666-
667-
cout << endl;
668-
669-
su2double MassFlow = config->GetSurface_MassFlow(iMarker_Analyze);
670-
if (si_units) cout << setw(20) << "Mf (kg/s): " << setw(15) << MassFlow;
671-
else if (us_units) cout << setw(20) << "Mf (lbs/s): " << setw(15) << MassFlow;
672-
673-
su2double NormalVelocity = config->GetSurface_NormalVelocity(iMarker_Analyze);
674-
if (si_units) cout << setw(20) << "Vn (m/s): " << setw(15) << NormalVelocity;
675-
else if (us_units) cout << setw(20) << "Vn (ft/s): " << setw(15) << NormalVelocity;
676-
677-
cout << endl;
678-
679-
su2double Uniformity = config->GetSurface_Uniformity(iMarker_Analyze);
680-
if (si_units) cout << setw(20) << "Uniformity (m/s): " << setw(15) << Uniformity;
681-
else if (us_units) cout << setw(20) << "Uniformity (ft/s): " << setw(15) << Uniformity;
682-
683-
su2double SecondaryStrength = config->GetSurface_SecondaryStrength(iMarker_Analyze);
684-
if (si_units) cout << setw(20) << "Secondary (m/s): " << setw(15) << SecondaryStrength;
685-
else if (us_units) cout << setw(20) << "Secondary (ft/s): " << setw(15) << SecondaryStrength;
686-
687-
cout << endl;
688-
689-
su2double MomentumDistortion = config->GetSurface_MomentumDistortion(iMarker_Analyze);
690-
cout << setw(20) << "Mom. Distortion: " << setw(15) << MomentumDistortion;
691-
692-
su2double SecondOverUniform = config->GetSurface_SecondOverUniform(iMarker_Analyze);
693-
cout << setw(20) << "Second/Uniform: " << setw(15) << SecondOverUniform;
694-
695-
cout << endl;
696-
697-
su2double Pressure = config->GetSurface_Pressure(iMarker_Analyze);
698-
if (si_units) cout << setw(20) << "P (Pa): " << setw(15) << Pressure;
699-
else if (us_units) cout << setw(20) << "P (psf): " << setw(15) << Pressure;
700-
701-
su2double TotalPressure = config->GetSurface_TotalPressure(iMarker_Analyze);
702-
if (si_units) cout << setw(20) << "PT (Pa): " << setw(15) <<TotalPressure;
703-
else if (us_units) cout << setw(20) << "PT (psf): " << setw(15) <<TotalPressure;
704-
705-
cout << endl;
706-
707-
su2double Mach = config->GetSurface_Mach(iMarker_Analyze);
708-
cout << setw(20) << "Mach: " << setw(15) << Mach;
709-
710-
su2double Density = config->GetSurface_Density(iMarker_Analyze);
711-
if (si_units) cout << setw(20) << "Rho (kg/m^3): " << setw(15) << Density;
712-
else if (us_units) cout << setw(20) << "Rho (lb/ft^3): " << setw(15) << Density*32.174;
713-
714-
cout << endl;
715-
716-
if (compressible || energy) {
717-
su2double Temperature = config->GetSurface_Temperature(iMarker_Analyze);
718-
if (si_units) cout << setw(20) << "T (K): " << setw(15) << Temperature;
719-
else if (us_units) cout << setw(20) << "T (R): " << setw(15) << Temperature;
720-
721-
su2double TotalTemperature = config->GetSurface_TotalTemperature(iMarker_Analyze);
722-
if (si_units) cout << setw(20) << "TT (K): " << setw(15) << TotalTemperature;
723-
else if (us_units) cout << setw(20) << "TT (R): " << setw(15) << TotalTemperature;
724-
725-
cout << endl;
726-
}
727-
728-
}
729-
cout.unsetf(ios_base::floatfield);
730-
731-
}
732-
733-
std::cout << std::resetiosflags(std::cout.flags());
734766
}
735767

736768
void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) {

‎TestCases/species_regression.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ def main():
4747
species2_primitiveVenturi.cfg_dir = "species_transport/venturi_primitive_3species"
4848
species2_primitiveVenturi.cfg_file = "species2_primitiveVenturi.cfg"
4949
species2_primitiveVenturi.test_iter = 50
50-
species2_primitiveVenturi.test_vals = [-5.957517, -5.187476, -5.037298, -5.851420, -1.511976, -6.046002, 5, -0.808614, 5 , -2.351161, 5, -0.247992, 0.000140, 0.000140, 0.000000, 0.000000]
50+
species2_primitiveVenturi.test_vals = [-5.957517, -5.187476, -5.037298, -5.851420, -1.511976, -6.046002, 5, -0.808614, 5 , -2.351161, 5, -0.247992, 0.000092, 0.000090, 0.000002, 0.000000]
5151
species2_primitiveVenturi.su2_exec = "mpirun -n 2 SU2_CFD"
5252
species2_primitiveVenturi.timeout = 1600
5353
species2_primitiveVenturi.new_output = True

0 commit comments

Comments
 (0)