diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index b8a9dcd20082..47d9b6f673ae 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2254,6 +2254,8 @@ struct StreamwisePeriodicValues { su2double Streamwise_Periodic_MassFlow; /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */ su2double Streamwise_Periodic_IntegratedHeatFlow; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */ su2double Streamwise_Periodic_InletTemperature; /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */ + su2double Streamwise_Periodic_BoundaryArea; /*!< \brief Global Surface area of the streamwise periodic interface. */ + su2double Streamwise_Periodic_AvgDensity; /*!< \brief Area avg density on the periodic interface. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 054ecbd42ec1..b7f9cff7a24b 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4803,8 +4803,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION); if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0) SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION); - if (DiscreteAdjoint && Kind_Streamwise_Periodic == ENUM_STREAMWISE_PERIODIC::MASSFLOW) - SU2_MPI::Error("Discrete Adjoint currently not validated for prescribed MASSFLOW.", CURRENT_FUNCTION); if (Ref_Inc_NonDim != DIMENSIONAL) SU2_MPI::Error("Streamwise Periodicity only works with \"INC_NONDIM= DIMENSIONAL\", the nondimensionalization with source terms doesn;t work in general.", CURRENT_FUNCTION); if (Axisymmetric) diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index 6a5c333b9b5c..980d74e7e817 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -39,7 +39,7 @@ class CIncEulerSolver : public CFVMFlowSolverBase { protected: vector FluidModel; /*!< \brief fluid model used in the solver. */ - StreamwisePeriodicValues SPvals; + StreamwisePeriodicValues SPvals, SPvalsUpdated; /*! * \brief Preprocessing actions common to the Euler and NS solvers. @@ -395,8 +395,30 @@ class CIncEulerSolver : public CFVMFlowSolverBase AD_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + VectorType SolutionExtra; /*!< \brief Stores adjoint solution for extra solution variables. + Currently only streamwise periodic pressure-drop for massflow prescribed flows. */ + VectorType ExternalExtra; /*!< \brief External storage for the adjoint value (i.e. for the OF mainly */ + + VectorType SolutionExtra_BGS_k; /*!< \brief Intermediate storage, enables cross term extraction as that is also pushed to Solution. */ + + protected: unsigned long nPoint = 0; /*!< \brief Number of points in the domain. */ unsigned long nDim = 0; /*!< \brief Number of dimension of the problem. */ unsigned long nVar = 0; /*!< \brief Number of variables of the problem. */ @@ -406,6 +413,31 @@ class CVariable { for(unsigned long iVar = 0; iVar < nVar; iVar++) External(iPoint,iVar) += val_sol[iVar]; } + /*! + * \brief Store the adjoint solution of the extra adjoint into the external container. + */ + void Set_ExternalExtra_To_SolutionExtra() { + assert(SolutionExtra.size() == ExternalExtra.size()); + for (auto iEntry = 0ul; iEntry < SolutionExtra.size(); iEntry++) + ExternalExtra[iEntry] = SolutionExtra[iEntry]; + } + + /*! + * \brief Add the external contribution to the solution for the extra adjoint solutions. + */ + void Add_ExternalExtra_To_SolutionExtra() { + assert(SolutionExtra.size() == ExternalExtra.size()); + for (auto iEntry = 0ul; iEntry < SolutionExtra.size(); iEntry++) + SolutionExtra[iEntry] += ExternalExtra[iEntry]; + } + + /*! + * \brief Return the extra adjoint solution. + * \return Reference to extra adjoint solution. + */ + inline VectorType& GetSolutionExtra() { return SolutionExtra; } + inline const VectorType& GetSolutionExtra() const { return SolutionExtra; } + /*! * \brief Update the variables using a conservative format. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index aa96d12f8f5e..02f25b1ed1f1 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -430,7 +430,7 @@ void CDiscAdjMultizoneDriver::Run() { /*--- Compute residual from Solution and Solution_BGS_k and update the latter. ---*/ SetResidual_BGS(iZone); - + Set_BGSSolution_k_To_Solution(iZone); } /*--- Set the multizone output. ---*/ diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 07eb87561337..beac491a0390 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -403,8 +403,10 @@ bool CMultizoneDriver::OuterConvergence(unsigned long OuterIter) { auto solvers = solver_container[iZone][INST_0][MESH_0]; for (unsigned short iSol = 0; iSol < MAX_SOLS; iSol++){ - if (solvers[iSol] != nullptr) + if (solvers[iSol] != nullptr) { solvers[iSol]->ComputeResidual_Multizone(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + solvers[iSol]->GetNodes()->Set_BGSSolution_k(); + } } /*--- Make sure that everything is loaded into the output container. ---*/ diff --git a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp index de2a560baced..d9a398b3f2b6 100644 --- a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp @@ -113,6 +113,10 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: Maximum residual of the temperature. AddHistoryOutput("RMS_ADJ_TEMPERATURE", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { + AddHistoryOutput("ADJOINT_SOLEXTRA", "Adjoint_SolExtra", ScreenOutputFormat::FIXED, "ADJOINT_SOLEXTRA", "Adjoint value of the first extra Solution.", HistoryFieldType::COEFFICIENT); + } + AddHistoryOutputFields_AdjScalarRMS_RES(config); if (config->AddRadiation()){ @@ -205,6 +209,11 @@ void CAdjFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS if (config->AddRadiation()){ SetHistoryOutputValue("RMS_ADJ_RAD_ENERGY", log10(adjrad_solver->GetRes_RMS(0))); } + + if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { + SetHistoryOutputValue("ADJOINT_SOLEXTRA", adjflow_solver->GetNodes()->GetSolutionExtra()[0]); + } + SetHistoryOutputValue("MAX_ADJ_PRESSURE", log10(adjflow_solver->GetRes_Max(0))); SetHistoryOutputValue("MAX_ADJ_VELOCITY-X", log10(adjflow_solver->GetRes_Max(1))); SetHistoryOutputValue("MAX_ADJ_VELOCITY-Y", log10(adjflow_solver->GetRes_Max(2))); diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 57962aa82278..ecac01b70cb7 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -29,20 +29,19 @@ #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" -CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *direct_solver, +CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *direct_sol, unsigned short Kind_Solver, unsigned short iMesh) : CSolver() { + /*-- Store some information about direct solver ---*/ + + KindDirect_Solver = Kind_Solver; + direct_solver = direct_sol; + adjoint = true; nVar = direct_solver->GetnVar(); nDim = geometry->GetnDim(); - /*--- Initialize arrays to NULL ---*/ - - /*-- Store some information about direct solver ---*/ - this->KindDirect_Solver = Kind_Solver; - this->direct_solver = direct_solver; - nMarker = config->GetnMarker_All(); nPoint = geometry->GetnPoint(); nPointDomain = geometry->GetnPointDomain(); @@ -86,6 +85,11 @@ CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *di nodes = new CDiscAdjVariable(Solution.data(), nPoint, nDim, nVar, config); SetBaseClassPointerToNodes(); + /*--- Allocate extra solution variables, if any are in use. ---*/ + + const auto nVarExtra = direct_solver->RegisterSolutionExtra(true, config); + nodes->AllocateAdjointSolutionExtra(nVarExtra); + switch(KindDirect_Solver){ case RUNTIME_FLOW_SYS: SolverName = "ADJ.FLOW"; @@ -163,6 +167,8 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) { /*--- Boolean true indicates that an input is registered ---*/ direct_solver->GetNodes()->RegisterSolution(true); + direct_solver->RegisterSolutionExtra(true, config); + if (time_n_needed) direct_solver->GetNodes()->RegisterSolution_time_n(); @@ -298,6 +304,8 @@ void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) { /*--- Register variables as output of the solver iteration. Boolean false indicates that an output is registered ---*/ direct_solver->GetNodes()->RegisterSolution(false); + + direct_solver->RegisterSolutionExtra(false, config); } void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) { @@ -337,6 +345,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi } END_SU2_OMP_FOR + direct_solver->ExtractAdjoint_SolutionExtra(nodes->GetSolutionExtra(), config); + /*--- Residuals and time_n terms are not needed when evaluating multizone cross terms. ---*/ if (CrossTerm) return; @@ -468,6 +478,8 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) { direct_solver->GetNodes()->SetAdjointSolution(iPoint,Solution); } END_SU2_OMP_FOR + + direct_solver->SetAdjoint_SolutionExtra(nodes->GetSolutionExtra(), config); } void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*) { diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 13002535ebd4..908cb8bf461b 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -100,6 +100,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned string filename_ = "flow"; filename_ = config->GetFilename(filename_, ".meta", Unst_RestartIter); Read_SU2_Restart_Metadata(geometry, config, adjoint, filename_); + if (rank==MASTER_NODE) cout << "Setting streamwise periodic pressure drop from restart metadata file." << endl; } /*--- Set the gamma value ---*/ @@ -1614,8 +1615,32 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont } // if turbulent - /*--- Set delta_p, m_dot, inlet_T, integrated_heat ---*/ - numerics->SetStreamwisePeriodicValues(SPvals); + if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { + /*---------------------------------------------------------------------------------------------*/ + /*--- Update the Pressure Drop [Pa] for the Momentum source term if Massflow is prescribed. ---*/ + /*--- The Pressure drop is iteratively adapted to result in the prescribed Target-Massflow. ---*/ + /*---------------------------------------------------------------------------------------------*/ + + /*--- Compute update to Delta p based on massflow-difference ---*/ + const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity; + const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea; + const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref()); + const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow; + const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2)); + + /*--- Store updated pressure difference ---*/ + const su2double damping_factor = config->GetInc_Outlet_Damping(); + SPvalsUpdated = SPvals; + SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP; + if (!config->GetDiscrete_Adjoint()) + SPvals = SPvalsUpdated; + + /*--- Set delta_p, m_dot, inlet_T, integrated_heat ---*/ + numerics->SetStreamwisePeriodicValues(SPvalsUpdated); + } + else { + numerics->SetStreamwisePeriodicValues(SPvals); + } AD::StartNoSharedReading(); @@ -3069,14 +3094,10 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf } void CIncEulerSolver::SetFreeStream_Solution(const CConfig *config){ - - unsigned long iPoint; - unsigned short iDim; - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPoint; iPoint++){ + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ nodes->SetSolution(iPoint,0, Pressure_Inf); - for (iDim = 0; iDim < nDim; iDim++){ + for (unsigned short iDim = 0; iDim < nDim; iDim++){ nodes->SetSolution(iPoint,iDim+1, Velocity_Inf[iDim]); } nodes->SetSolution(iPoint,nDim+1, Temperature_Inf); @@ -3084,3 +3105,24 @@ void CIncEulerSolver::SetFreeStream_Solution(const CConfig *config){ END_SU2_OMP_FOR } + +unsigned long CIncEulerSolver::RegisterSolutionExtra(bool input, const CConfig* config) { + if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { + if (input) AD::RegisterInput(SPvals.Streamwise_Periodic_PressureDrop); + else AD::RegisterOutput(SPvalsUpdated.Streamwise_Periodic_PressureDrop); + return 1; + } + return 0; +} + +void CIncEulerSolver::SetAdjoint_SolutionExtra(const su2activevector& adj_sol, const CConfig* config) { + if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { + SU2_TYPE::SetDerivative(SPvalsUpdated.Streamwise_Periodic_PressureDrop, SU2_TYPE::GetValue(adj_sol[0])); + } +} + +void CIncEulerSolver::ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, const CConfig* config) { + if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { + adj_sol[0] = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop); + } +} diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 57ca1d5544e0..2e3e2e81de62 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -52,6 +52,12 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short default: break; } + + /*--- Set the initial Streamwise periodic pressure drop value. ---*/ + + if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) + // Note during restarts, the flow.meta is read first. But that sets the cfg-value so we are good here. + SPvals.Streamwise_Periodic_PressureDrop = config->GetStreamwise_Periodic_PressureDrop(); } void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, @@ -132,10 +138,6 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry, /*--- Pressure-Drop update in case of a prescribed massflow. ---*/ /*-------------------------------------------------------------------------------------------------*/ - const auto nZone = geometry->GetnZone(); - const auto InnerIter = config->GetInnerIter(); - const auto OuterIter = config->GetOuterIter(); - su2double Area_Local = 0.0, MassFlow_Local = 0.0, Average_Density_Local = 0.0, @@ -189,44 +191,8 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry, /*--- Set solver variables ---*/ SPvals.Streamwise_Periodic_MassFlow = MassFlow_Global; SPvals.Streamwise_Periodic_InletTemperature = Temperature_Global; - - /*--- As deltaP changes with prescribed massflow the const config value should only be used once. ---*/ - if((nZone==1 && InnerIter==0) || - (nZone>1 && OuterIter==0 && InnerIter==0)) { - SPvals.Streamwise_Periodic_PressureDrop = config->GetStreamwise_Periodic_PressureDrop() / config->GetPressure_Ref(); - } - - if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { - /*------------------------------------------------------------------------------------------------*/ - /*--- 2. Update the Pressure Drop [Pa] for the Momentum source term if Massflow is prescribed. ---*/ - /*--- The Pressure drop is iteratively adapted to result in the prescribed Target-Massflow. ---*/ - /*------------------------------------------------------------------------------------------------*/ - - /*--- Load/define all necessary variables ---*/ - const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref()); - const su2double damping_factor = config->GetInc_Outlet_Damping(); - su2double Pressure_Drop_new, ddP; - - /*--- Compute update to Delta p based on massflow-difference ---*/ - ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2)); - - /*--- Store updated pressure difference ---*/ - Pressure_Drop_new = SPvals.Streamwise_Periodic_PressureDrop + damping_factor*ddP; - /*--- During restarts, this routine GetStreamwise_Periodic_Properties can get called multiple times - (e.g. 4x for INC_RANS restart). Each time, the pressure drop gets updated. For INC_RANS restarts - it gets called 2x before the restart files are read such that the current massflow is - Area*inital-velocity which can be way off! - With this there is still a slight inconsistency wrt to a non-restarted simulation: The restarted "zero-th" - iteration does not get a pressure-update but the continuing simulation would have an update here. This can be - fully neglected if the pressure drop is converged. And for all other cases it should be minor difference at - best. ---*/ - if((nZone==1 && InnerIter>0) || - (nZone>1 && OuterIter>0)) { - SPvals.Streamwise_Periodic_PressureDrop = Pressure_Drop_new; - } - - } // if massflow - + SPvals.Streamwise_Periodic_BoundaryArea = Area_Global; + SPvals.Streamwise_Periodic_AvgDensity = Average_Density_Global; if (config->GetEnergy_Equation()) { /*---------------------------------------------------------------------------------------------*/ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 5d9b81c3dc5b..bed7f10ffc72 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -2198,12 +2198,16 @@ void CSolver::Add_External_To_Solution() { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { base_nodes->AddSolution(iPoint, base_nodes->Get_External(iPoint)); } + + base_nodes->Add_ExternalExtra_To_SolutionExtra(); } void CSolver::Add_Solution_To_External() { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { base_nodes->Add_External(iPoint, base_nodes->GetSolution(iPoint)); } + + base_nodes->Set_ExternalExtra_To_SolutionExtra(); } void CSolver::Update_Cross_Term(CConfig *config, su2passivematrix &cross_term) { @@ -4166,8 +4170,6 @@ void CSolver::ComputeResidual_Multizone(const CGeometry *geometry, const CConfig for (unsigned short iVar = 0; iVar < nVar; iVar++) { const su2double Res = (base_nodes->Get_BGSSolution(iPoint,iVar) - base_nodes->Get_BGSSolution_k(iPoint,iVar))*domain; - base_nodes->Set_BGSSolution_k(iPoint,iVar, base_nodes->Get_BGSSolution(iPoint,iVar)); - /*--- Update residual information for current thread. ---*/ resRMS[iVar] += Res*Res; if (fabs(Res) > resMax[iVar]) { diff --git a/SU2_CFD/src/variables/CDiscAdjVariable.cpp b/SU2_CFD/src/variables/CDiscAdjVariable.cpp index 69bc25dbd848..affa11b1f7ae 100644 --- a/SU2_CFD/src/variables/CDiscAdjVariable.cpp +++ b/SU2_CFD/src/variables/CDiscAdjVariable.cpp @@ -49,3 +49,11 @@ void CDiscAdjVariable::Set_External_To_DualTimeDer() { assert(External.size() == DualTime_Derivative.size()); parallelCopy(External.size(), DualTime_Derivative.data(), External.data()); } + +void CDiscAdjVariable::AllocateAdjointSolutionExtra(unsigned long nVarExtra) { + if (nVarExtra == 0) return; + SolutionExtra.resize(nVarExtra) = su2double(1e-16); + /*--- These are only for multizone, but since nVarExtra is small we allocate by default. ---*/ + SolutionExtra_BGS_k.resize(nVarExtra) = su2double(1e-16); + ExternalExtra.resize(nVarExtra) = su2double(0.0); +} diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 80d574ceecc6..21d72d120e85 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -101,11 +101,17 @@ void CVariable::Set_Solution_time_n1() { void CVariable::Set_BGSSolution_k() { assert(Solution_BGS_k.size() == Solution.size()); parallelCopy(Solution.size(), Solution.data(), Solution_BGS_k.data()); + + assert(SolutionExtra_BGS_k.size() == SolutionExtra.size()); + parallelCopy(SolutionExtra.size(), SolutionExtra.data(), SolutionExtra_BGS_k.data()); } void CVariable::Restore_BGSSolution_k() { assert(Solution.size() == Solution_BGS_k.size()); parallelCopy(Solution_BGS_k.size(), Solution_BGS_k.data(), Solution.data()); + + assert(SolutionExtra.size() == SolutionExtra_BGS_k.size()); + parallelCopy(SolutionExtra_BGS_k.size(), SolutionExtra_BGS_k.data(), SolutionExtra.data()); } void CVariable::SetExternalZero() { parallelSet(External.size(), 0.0, External.data()); } diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configFluid.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configFluid.cfg new file mode 100644 index 000000000000..e4085f06cdba --- /dev/null +++ b/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configFluid.cfg @@ -0,0 +1,141 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Unit Cell flow around pin array (fluid) % +% Author: T. Kattmann % +% Date: 2022.02.10 % +% File Version 7.3.0 "Blackbird" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_RANS +% +KIND_TURB_MODEL= SST +% +%__OF_AVGT__OBJECTIVE_FUNCTION= AVG_TEMPERATURE +%__OF_AVGT__OBJECTIVE_WEIGHT= 0.0 +OBJECTIVE_FUNCTION= SURFACE_PRESSURE_DROP +OBJECTIVE_WEIGHT= 1.0 +%__OF_MDOT__OBJECTIVE_FUNCTION= SURFACE_MASSFLOW +%__OF_MDOT__OBJECTIVE_WEIGHT= 1.0 +%__OF_DRAG__OBJECTIVE_FUNCTION= DRAG +%__OF_DRAG__OBJECTIVE_WEIGHT= 1.0 +% +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_INIT= 1045.0 +INC_VELOCITY_INIT= ( 0.1, 0.0, 0.0 ) +% +INC_ENERGY_EQUATION = YES +INC_TEMPERATURE_INIT= 338.0 +INC_NONDIM= DIMENSIONAL +SPECIFIC_HEAT_CP= 3540.0 +% +FREESTREAM_TURBULENCEINTENSITY= 0.05 +FREESTREAM_TURB2LAMVISCRATIO= 10.0 +% +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001385 +% +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +% Pr_lam = mu_lam [Pa*s] * c_p [J/(kg*K)] / lambda[W/(m*K)] +% = 1.385e-3 * 3540 / 0.42 +% = 11.7 +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +PRANDTL_LAM= 11.7 +% +TURBULENT_CONDUCTIVITY_MODEL= CONSTANT_PRANDTL_TURB +PRANDTL_TURB= 0.90 +% +% --------------------- STREAMWISE PERIODICITY DEFINITION ---------------------% +% +KIND_STREAMWISE_PERIODIC= MASSFLOW +STREAMWISE_PERIODIC_PRESSURE_DROP= 208.023676 +STREAMWISE_PERIODIC_MASSFLOW= 0.85 +INC_OUTLET_DAMPING= 0.1 +% +STREAMWISE_PERIODIC_TEMPERATURE= NO +% +% inner pin length 0.00376991 m = (0.00322-0.00262)*2*pi +% with 5e5 W/m that is Q = -1884.95559215 +STREAMWISE_PERIODIC_OUTLET_HEAT= -1884.95559215 +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_SYM= ( fluid_symmetry_lower, fluid_symmetry_upper ) +% +MARKER_PERIODIC= ( fluid_inlet, fluid_outlet, 0.0,0.0,0.0, 0.0,0.0,0.0, 0.0111544,0.0,0.0 ) +% +% Alternative to periodic simulation with velocity inlet and pressure outlet +%MARKER_HEATFLUX= ( fluid_pin1_interface, 5e5, \ +% fluid_pin2_interface, 5e5, \ +% fluid_pin3_interface, 5e5 ) +% +% Alternative options for non-periodic flow +%INC_INLET_TYPE= VELOCITY_INLET +%MARKER_INLET= ( fluid_inlet, 338.0, 0.75, 1.0, 0.0, 0.0 ) +% Test vals to hinder outlet backflow +%MARKER_INLET= ( fluid_inlet, 338.0, 0.3, 1.0, 0.0, 0.0 ) +% +%INC_OUTLET_TYPE= PRESSURE_OUTLET +%MARKER_OUTLET= ( fluid_outlet, 0.0 ) +% +% ------------------------ SURFACES IDENTIFICATION ----------------------------% +% +MARKER_MONITORING= fluid_pin2_interface +% +%__DIRECT__MARKER_ANALYZE= fluid_outlet, fluid_inlet, fluid_pin2_interface +%__OF_AVGT__MARKER_ANALYZE= fluid_pin2_interface +%__OF_MDOT__MARKER_ANALYZE= fluid_inlet +MARKER_ANALYZE= fluid_outlet, fluid_inlet +MARKER_ANALYZE_AVERAGE= AREA +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +%ITER= 3500 +NUM_METHOD_GRAD= GREEN_GAUSS +%__DIRECT__CFL_NUMBER= 1e3 +% Optimal CFL numbers differ from OF to OF +%__OF_AVGT__CFL_NUMBER= 1e3 +%__OF_DRAG__CFL_NUMBER= 1e2 +CFL_NUMBER= 1e2 +%__OF_MDOT__CFL_NUMBER= 1e2 +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1e-15 +%__DIRECT__LINEAR_SOLVER_ITER= 10 +%__OF_AVGT__LINEAR_SOLVER_ITER= 10 +%__OF_DRAG__LINEAR_SOLVER_ITER= 15 +LINEAR_SOLVER_ITER= 15 +%__OF_MDOT__LINEAR_SOLVER_ITER= 15 +% +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= NO +SLOPE_LIMITER_TURB= NONE +TIME_DISCRE_TURB= EULER_IMPLICIT +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Note: SURFACE_MASSFLOW does only converge to ~1e-7 here +CONV_FIELD= SURFACE_STATIC_TEMPERATURE, DRAG, SURFACE_PRESSURE_DROP +CONV_CAUCHY_EPS= 1e-15 +CONV_CAUCHY_ELEMS= 100 +CONV_STARTITER= 10 diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configMaster.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configMaster.cfg new file mode 100644 index 000000000000..13084c9b96b4 --- /dev/null +++ b/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configMaster.cfg @@ -0,0 +1,114 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D cylinder array with CHT coupling % +% Author: T. Kattmann % +% Institution: None % +% Date: 2022.02.10 % +% File Version 7.3.0 "Blackbird" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= MULTIPHYSICS +% +CONFIG_LIST= configFluid.cfg, configSolid.cfg +% +MARKER_ZONE_INTERFACE= fluid_pin1_interface, solid_pin1_interface, fluid_pin2_interface, solid_pin2_interface, fluid_pin3_interface, solid_pin3_interface +MARKER_CHT_INTERFACE= fluid_pin1_interface, solid_pin1_interface, fluid_pin2_interface, solid_pin2_interface, fluid_pin3_interface, solid_pin3_interface +% +% Number of total iterations +%__DIRECT__OUTER_ITER= 4000 +OUTER_ITER= 22000 +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= AVG_TEMPERATURE[1], SURFACE_PRESSURE_DROP[0], DRAG[0] +CONV_CAUCHY_EPS= 1e-15 +CONV_CAUCHY_ELEMS= 100 +CONV_STARTITER= 10 +% +% -------------------------------------- OUTPUT -------------------------------% +% +%__DIRECT__SCREEN_OUTPUT= OUTER_ITER, RMS_PRESSURE[0], RMS_VELOCITY-X[0], RMS_VELOCITY-Y[0], RMS_TEMPERATURE[0], RMS_TKE[0], RMS_DISSIPATION[0], RMS_TEMPERATURE[1], STREAMWISE_MASSFLOW[0], STREAMWISE_DP[0], AVG_TEMPERATURE[1] +SCREEN_OUTPUT= OUTER_ITER, RMS_ADJ_PRESSURE[0], RMS_ADJ_VELOCITY-X[0], RMS_ADJ_VELOCITY-Y[0], RMS_ADJ_TEMPERATURE[0], RMS_ADJ_TKE[0], RMS_ADJ_DISSIPATION[0], RMS_ADJ_TEMPERATURE[1], ADJOINT_SOLEXTRA[0] +SCREEN_WRT_FREQ_OUTER= 100 +% +%__DIRECT__HISTORY_OUTPUT= ITER, BGS_RES[0], BGS_RES[1], RMS_RES[0], RMS_RES[1], STREAMWISE_PERIODIC[0], FLOW_COEFF[0], FLOW_COEFF_SURF[0], AERO_COEFF[0], HEAT[0], HEAT[1], LINSOL[0], LINSOL[1] +HISTORY_OUTPUT= ITER, BGS_RES[0], BGS_RES[1], RMS_RES[0], RMS_RES[1], ADJOINT_SOLEXTRA[0] +OUTPUT_PRECISION= 16 +% +OUTPUT_FILES= RESTART, PARAVIEW_MULTIBLOCK +OUTPUT_WRT_FREQ= 10000 +% +SOLUTION_FILENAME= restart +SOLUTION_ADJ_FILENAME= restart_adj +GRAD_OBJFUNC_FILENAME= of_grad.csv +% +MESH_FILENAME= 2D-PinArray.su2 +% +% -------------------- FREE-FORM DEFORMATION PARAMETERS -----------------------% +% +FFD_TOLERANCE= 1E-10 +FFD_ITERATIONS= 500 +% +% FFD box definition: 2D case (FFD_BoxTag, X1, Y1, 0.0, X2, Y2, 0.0, X3, Y3, 0.0, X4, Y4, 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_DEFINITION= (BOX, 0.0029772, 0.0, 0.0, 0.0081772, 0.0, 0.0, 0.0081772, 0.00322, 0.0, 0.0029772, 0.00322, 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 ); \ + (FRONTBOX, -0.0026, 0.00322, 0.0, 0.0026, 0.00322, 0.0, 0.0026, 0.0, 0.0, -0.0026, 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,0.0 ); \ + (BACKBOX, 0.0085544, 0.00322, 0.0, 0.0137544, 0.00322, 0.0, 0.0137544, 0.0, 0.0, 0.0085544, 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,0.0 ) +% +% FFD box degree: 2D case (x_degree, y_degree, 0) +FFD_DEGREE= (8, 8, 0); (8, 8, 0); (8, 8, 0) +% +% Surface grid continuity at the intersection with the faces of the FFD boxes. +% To keep a particular level of surface continuity, SU2 automatically freezes the right +% number of control point planes (NO_DERIVATIVE, 1ST_DERIVATIVE, 2ND_DERIVATIVE, USER_INPUT) +FFD_CONTINUITY= USER_INPUT +% +% BEZIER, BSPLINE_UNIFORM +%FFD_BLENDING= BEZIER +FFD_BLENDING= BSPLINE_UNIFORM +FFD_BSPLINE_ORDER= 4,3,2 +% +% ----------------------- DESIGN VARIABLE PARAMETERS --------------------------% +% +% Config options for writing the FFD-box into the mesh. +% Comment these options if they appear elsewhere in the .cfg file. +%DV_KIND= FFD_SETTING +%DV_PARAM= ( 1.0 ) +%DV_VALUE= 1.0 + +DV_KIND= FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D +% +% Marker of the surface in which we are going apply the shape deformation +DV_MARKER= fluid_pin2_interface, solid_pin2_interface, fluid_pin1_interface, solid_pin1_interface, fluid_pin3_interface, solid_pin3_interface, fluid_inlet, solid_pin1_inlet, fluid_outlet, solid_pin3_outlet +% +% Parameters of the shape deformation +% - FFD_SETTING ( 1.0 ) +% - FFD_CONTROL_POINT_2D ( FFD_BoxTag, i_Ind, j_Ind, x_Disp, y_Disp ) +DV_PARAM= \ +( BOX, 1, 5, 0.0, 1.0);\ +( BOX, 2, 5, 0.0, 1.0);\ +( BOX, 3, 5, 0.0, 1.0);\ +( BOX, 4, 5, 0.0, 1.0);\ +( BOX, 5, 5, 0.0, 1.0);\ +( BOX, 6, 5, 0.0, 1.0);\ +( BOX, 7, 5, 0.0, 1.0) +% +% Value of the shape deformation +DV_VALUE=0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 +% +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +DEFORM_LINEAR_SOLVER= FGMRES +DEFORM_LINEAR_SOLVER_PREC= ILU +DEFORM_LINEAR_SOLVER_ERROR= 1E-14 +% +DEFORM_NONLINEAR_ITER= 1 +DEFORM_LINEAR_SOLVER_ITER= 1000 +% +DEFORM_CONSOLE_OUTPUT= YES +DEFORM_STIFFNESS_TYPE= WALL_DISTANCE +% +% Deformation coefficient (linear elasticity limits from -1.0 to 0.5, a larger value is also possible) +DEFORM_COEFF = 0.1 diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configSolid.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configSolid.cfg new file mode 100644 index 000000000000..3692dae7f268 --- /dev/null +++ b/TestCases/incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d/configSolid.cfg @@ -0,0 +1,69 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Unit Cell flow around pin array (solid) % +% Author: T. Kattmann % +% Institution: None % +% Date: 2022.02.10 % +% File Version 7.3.0 "Blackbird" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= HEAT_EQUATION +% +%__OF_AVGT__OBJECTIVE_FUNCTION= AVG_TEMPERATURE +%__OF_AVGT__OBJECTIVE_WEIGHT= 1.0 +OBJECTIVE_FUNCTION= AVG_TEMPERATURE +OBJECTIVE_WEIGHT= 0.0 +%__OF_MDOT__OBJECTIVE_FUNCTION= AVG_TEMPERATURE +%__OF_MDOT__OBJECTIVE_WEIGHT= 0.0 +%__OF_DRAG__OBJECTIVE_FUNCTION= AVG_TEMPERATURE +%__OF_DRAG__OBJECTIVE_WEIGHT= 0.0 +% +% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------% +% +INC_NONDIM= DIMENSIONAL +FREESTREAM_TEMPERATURE= 345.0 +MATERIAL_DENSITY= 2719 +SPECIFIC_HEAT_CP = 871.0 +THERMAL_CONDUCTIVITY_CONSTANT= 200 +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_SYM= (solid_pin1_walls, solid_pin2_walls, solid_pin3_walls) +% +MARKER_HEATFLUX= (solid_pin1_inner, 5e5, \ + solid_pin2_inner, 5e5, \ + solid_pin3_inner, 5e5, \ + solid_pin1_inlet, 0.0, \ + solid_pin3_outlet, 0.0 ) +% +% ------------------------ SURFACES IDENTIFICATION ----------------------------% +% +MARKER_MONITORING = ( solid_pin2_inner ) +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +%__DIRECT__CFL_NUMBER= 1e4 +CFL_NUMBER= 1e3 +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-15 +LINEAR_SOLVER_ITER= 20 +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= AVG_TEMPERATURE +CONV_CAUCHY_EPS= 1e-15 +CONV_CAUCHY_ELEMS= 100 +CONV_STARTITER= 10 +% +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +TIME_DISCRE_HEAT= EULER_IMPLICIT diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 41a5286cd5fc..c66fa9308946 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -346,6 +346,18 @@ def main(): da_sp_pinArray_cht_2d_dp_hf.multizone = True test_list.append(da_sp_pinArray_cht_2d_dp_hf) + # 2D DA cht streamwise periodic case, 2 zones, PressureDrop objective, additional pressure drop adjoint equation + da_sp_pinArray_cht_2d_mf = TestCase('da_sp_pinArray_cht_2d_mf') + da_sp_pinArray_cht_2d_mf.cfg_dir = "incomp_navierstokes/streamwise_periodic/dp-adjoint_chtPinArray_2d" + da_sp_pinArray_cht_2d_mf.cfg_file = "configMaster.cfg" + da_sp_pinArray_cht_2d_mf.test_iter = 100 + da_sp_pinArray_cht_2d_mf.test_vals = [-4.609357, -1.273838, -1.502734, -18.503852, -0.834358, -5.813324, -19.074376, -48.287655] + da_sp_pinArray_cht_2d_mf.su2_exec = "mpirun -n 2 SU2_CFD_AD" + da_sp_pinArray_cht_2d_mf.timeout = 1600 + da_sp_pinArray_cht_2d_mf.tol = 0.00001 + da_sp_pinArray_cht_2d_mf.multizone = True + test_list.append(da_sp_pinArray_cht_2d_mf) + # 2D unsteady CHT vortex shedding at RE=200. TAVG_Temperature OF da_unsteadyCHT_cylinder = TestCase('da_unsteadyCHT_cylinder') da_unsteadyCHT_cylinder.cfg_dir = "coupled_cht/disc_adj_unsteadyCHT_cylinder" diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index fa40656d3510..6f9574dbc709 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -76,7 +76,7 @@ def main(): sp_pinArray_2d_mf_hf.cfg_dir = "../Tutorials/incompressible_flow/Inc_Streamwise_Periodic" sp_pinArray_2d_mf_hf.cfg_file = "sp_pinArray_2d_mf_hf.cfg" sp_pinArray_2d_mf_hf.test_iter = 25 - sp_pinArray_2d_mf_hf.test_vals = [-4.600340, 1.470386, -0.778623, 266.569743] #last 4 lines + sp_pinArray_2d_mf_hf.test_vals = [-4.626384, 1.444465, -0.750978, 241.757337] #last 4 lines sp_pinArray_2d_mf_hf.su2_exec = "mpirun -n 2 SU2_CFD" sp_pinArray_2d_mf_hf.timeout = 1600 sp_pinArray_2d_mf_hf.tol = 0.00001