From 02b466c8ada5ec31e8e1456cea0a1119dce97303 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 4 Apr 2023 16:19:07 +0100 Subject: [PATCH 01/20] use the normal unsteady options for FEA problems instead of DYN_** --- Common/include/CConfig.hpp | 32 ------------------- Common/include/option_structure.hpp | 12 ------- Common/src/CConfig.cpp | 21 +++++------- .../interfaces/fsi/CFlowTractionInterface.cpp | 2 +- .../src/output/output_structure_legacy.cpp | 2 +- SU2_CFD/src/solvers/CFEASolver.cpp | 15 +++++---- .../fea_fsi/DynBeam_2d/configBeam_2d.cfg | 19 +++++------ 7 files changed, 28 insertions(+), 75 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index ba5190f4b24..d790a723604 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -150,7 +150,6 @@ class CConfig { su2double CL_Target; /*!< \brief Fixed Cl mode Target Cl. */ su2double Confinement_Param; /*!< \brief Confinement paramenter for Vorticity Confinement method. */ TIME_MARCHING TimeMarching; /*!< \brief Steady or unsteady (time stepping or dual time stepping) computation. */ - unsigned short Dynamic_Analysis; /*!< \brief Static or dynamic structural analysis. */ su2double FixAzimuthalLine; /*!< \brief Fix an azimuthal line due to misalignments of the nearfield. */ su2double **DV_Value; /*!< \brief Previous value of the design variable. */ su2double Venkat_LimiterCoeff; /*!< \brief Limiter coefficient */ @@ -168,9 +167,6 @@ class CConfig { su2double HarmonicBalance_Period; /*!< \brief Period of oscillation to be used with harmonic balance computations. */ su2double Delta_UnstTime, /*!< \brief Time step for unsteady computations. */ Delta_UnstTimeND; /*!< \brief Time step for unsteady computations (non dimensional). */ - su2double Delta_DynTime, /*!< \brief Time step for dynamic structural computations. */ - Total_DynTime, /*!< \brief Total time for dynamic structural computations. */ - Current_DynTime; /*!< \brief Global time of the dynamic structural computations. */ su2double Total_UnstTime, /*!< \brief Total time for unsteady computations. */ Total_UnstTimeND; /*!< \brief Total time for unsteady computations (non dimensional). */ su2double Current_UnstTime, /*!< \brief Global time of the unsteady simulation. */ @@ -8723,34 +8719,6 @@ class CConfig { */ bool GetSteadyRestart(void) const { return SteadyRestart; } - /*! - * \brief Provides information about the time integration of the structural analysis, and change the write in the output - * files information about the iteration. - * \return The kind of time integration: Static or dynamic analysis - */ - unsigned short GetDynamic_Analysis(void) const { return Dynamic_Analysis; } - - /*! - * \brief If we are prforming an unsteady simulation, there is only - * one value of the time step for the complete simulation. - * \return Value of the time step in an unsteady simulation (non dimensional). - */ - su2double GetDelta_DynTime(void) const { return Delta_DynTime; } - - /*! - * \brief If we are prforming an unsteady simulation, there is only - * one value of the time step for the complete simulation. - * \return Value of the time step in an unsteady simulation (non dimensional). - */ - su2double GetTotal_DynTime(void) const { return Total_DynTime; } - - /*! - * \brief If we are prforming an unsteady simulation, there is only - * one value of the time step for the complete simulation. - * \return Value of the time step in an unsteady simulation (non dimensional). - */ - su2double GetCurrent_DynTime(void) const { return Current_DynTime; } - /*! * \brief Get the current instance. * \return Current instance identifier. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 4c934f0612b..57f89986b66 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2368,18 +2368,6 @@ enum class RECORDING { SOLUTION_AND_MESH, }; -/*! - * \brief Types of schemes for dynamic structural computations - */ -enum ENUM_DYNAMIC { - STATIC = 0, /*!< \brief A static structural computation. */ - DYNAMIC = 1 /*!< \brief Use a time stepping strategy for dynamic computations. */ -}; -static const MapType Dynamic_Map = { - MakePair("NO", STATIC) - MakePair("YES", DYNAMIC) -}; - /*! * \brief Types of input file formats */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index ea45d2d0159..5b3eb0c10a0 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2434,12 +2434,6 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Temporary: pseudo static analysis (no density in dynamic analysis) * Options: NO, YES \ingroup Config */ addBoolOption("PSEUDO_STATIC", PseudoStatic, false); - /* DESCRIPTION: Dynamic or static structural analysis */ - addEnumOption("DYNAMIC_ANALYSIS", Dynamic_Analysis, Dynamic_Map, STATIC); - /* DESCRIPTION: Time Step for dynamic analysis (s) */ - addDoubleOption("DYN_TIMESTEP", Delta_DynTime, 0.0); - /* DESCRIPTION: Total Physical Time for dual time stepping simulations (s) */ - addDoubleOption("DYN_TIME", Total_DynTime, 1.0); /* DESCRIPTION: Parameter alpha for Newmark scheme (s) */ addDoubleOption("NEWMARK_BETA", Newmark_beta, 0.25); /* DESCRIPTION: Parameter delta for Newmark scheme (s) */ @@ -3032,6 +3026,12 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){ newString.append("SOLID_TEMPERATURE_INIT is deprecated. Use FREESTREAM_TEMPERATURE instead.\n\n"); else if (!option_name.compare("SA_QCR")) newString.append("SA_QCR is deprecated. Use SA_OPTIONS=QCR2000 instead.\n\n"); + else if (!option_name.compare("DYN_TIMESTEP")) + newString.append("DYN_TIMESTEP is deprecated. Use TIME_STEP instead.\n\n"); + else if (!option_name.compare("DYN_TIME")) + newString.append("DYN_TIME is deprecated. Use MAX_TIME instead.\n\n"); + else if (!option_name.compare("DYNAMIC_ANALYSIS")) + newString.append("DYNAMIC_ANALYSIS is deprecated. Use TIME_DOMAIN instead.\n\n"); else { /*--- Find the most likely candidate for the unrecognized option, based on the length of start and end character sequences shared by candidates and the option. ---*/ @@ -3707,7 +3707,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Time_Domain){ Delta_UnstTime = Time_Step; - Delta_DynTime = Time_Step; if (TimeMarching == TIME_MARCHING::TIME_STEPPING){ InnerIter = 1; } @@ -6840,7 +6839,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { else { if (Time_Domain) { cout << "Dynamic structural analysis."<< endl; - cout << "Time step provided by the user for the dynamic analysis(s): "<< Delta_DynTime << "." << endl; + cout << "Time step provided by the user for the dynamic analysis(s): "<< Time_Step << "." << endl; } else { cout << "Static structural analysis." << endl; } @@ -6915,7 +6914,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "Generalized-alpha method." << endl; break; case STRUCT_TIME_INT::NEWMARK_IMPLICIT: - if (Dynamic_Analysis) cout << "Newmark implicit method for the structural time integration." << endl; + if (Time_Domain) cout << "Newmark implicit method for the structural time integration." << endl; switch (Kind_Linear_Solver) { case BCGSTAB: cout << "BCGSTAB is used for solving the linear system." << endl; @@ -8461,9 +8460,6 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, case MAIN_SOLVER::FEM_ELASTICITY: case MAIN_SOLVER::DISC_ADJ_FEM: - - Current_DynTime = static_cast(TimeIter)*Delta_DynTime; - if (val_system == RUNTIME_FEA_SYS) { SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE); SetKind_TimeIntScheme(NONE); @@ -10040,7 +10036,6 @@ void CConfig::SetMultizone(const CConfig *driver_config, const CConfig* const* c /*--- Fix the Time Step for all subdomains, for the case of time-dependent problems ---*/ if (driver_config->GetTime_Domain()){ Delta_UnstTime = driver_config->GetTime_Step(); - Delta_DynTime = driver_config->GetTime_Step(); Time_Domain = true; } diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index d95bb5a9841..1aaa0d422fe 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -132,7 +132,7 @@ void CFlowTractionInterface::GetPhysical_Constants(CSolver *flow_solution, CSolv /*--- Apply a ramp to the transfer of the fluid loads ---*/ su2double ModAmpl = 0.0; - su2double CurrentTime = struct_config->GetCurrent_DynTime(); + su2double CurrentTime = struct_config->GetCurrent_UnstTime(); bool Ramp_Load = struct_config->GetRamp_Load(); su2double Ramp_Time = struct_config->GetRamp_Time(); diff --git a/SU2_CFD/src/output/output_structure_legacy.cpp b/SU2_CFD/src/output/output_structure_legacy.cpp index 40e495d6b66..896923965aa 100644 --- a/SU2_CFD/src/output/output_structure_legacy.cpp +++ b/SU2_CFD/src/output/output_structure_legacy.cpp @@ -1670,7 +1670,7 @@ void COutputLegacy::SetConvHistoryBody(ofstream *ConvHist_file, } else if (fem && !fsi) { if (dynamic) { - cout << endl << "Simulation time: " << config[val_iZone]->GetCurrent_DynTime() << ". Time step: " << config[val_iZone]->GetDelta_DynTime() << "."; + //cout << endl << "Simulation time: " << config[val_iZone]->GetCurrent_DynTime() << ". Time step: " << config[val_iZone]->GetDelta_DynTime() << "."; } } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 4f30cb91450..96d2c6e0a72 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -50,6 +50,7 @@ CFEASolver::CFEASolver(LINEAR_SOLVER_MODE mesh_deform_mode) : CFEASolverBase(mes CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(geometry, config) { bool dynamic = (config->GetTime_Domain()); + config->SetDelta_UnstTimeND(config->GetDelta_UnstTime()); /*--- Test whether we consider dielectric elastomers ---*/ bool de_effects = config->GetDE_Effects(); @@ -1282,7 +1283,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, if (outputReactions) { - bool dynamic = (config->GetDynamic_Analysis() == DYNAMIC); + const bool dynamic = config->GetTime_Domain(); ofstream myfile; myfile.open ("Reactions.txt"); @@ -1481,7 +1482,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, con void CFEASolver::Compute_IntegrationConstants(const CConfig *config) { - su2double Delta_t= config->GetDelta_DynTime(); + su2double Delta_t= config->GetDelta_UnstTime(); su2double gamma = config->GetNewmark_gamma(), beta = config->GetNewmark_beta(); @@ -1652,7 +1653,7 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, const CConfig *config, unsigned const su2double *DispDirLocal = config->GetDisp_Dir(TagBound); su2double DispDirMod = Norm(nDim, DispDirLocal); - su2double CurrentTime = config->GetCurrent_DynTime(); + su2double CurrentTime = config->GetCurrent_UnstTime(); su2double RampTime = config->GetRamp_Time(); su2double ModAmpl = Compute_LoadCoefficient(CurrentTime, RampTime, config); @@ -1818,7 +1819,7 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, const CConfig *config, unsi /*--- Retrieve the normal pressure and the application conditions for the considered boundary. ---*/ - su2double CurrentTime = config->GetCurrent_DynTime(); + su2double CurrentTime = config->GetCurrent_UnstTime(); su2double Ramp_Time = config->GetRamp_Time(); su2double ModAmpl = Compute_LoadCoefficient(CurrentTime, Ramp_Time, config); @@ -1914,7 +1915,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, const CConfig *config, unsigne /*--- Compute the norm of the vector that was passed in the config file. ---*/ su2double LoadNorm = Norm(nDim, Load_Dir_Local); - su2double CurrentTime=config->GetCurrent_DynTime(); + su2double CurrentTime=config->GetCurrent_UnstTime(); su2double Ramp_Time = config->GetRamp_Time(); su2double ModAmpl = Compute_LoadCoefficient(CurrentTime, Ramp_Time, config); @@ -2036,7 +2037,7 @@ su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double R /*--- This offset introduces the ramp load in dynamic cases starting from the restart point. ---*/ bool offset = (restart && fsi && (!stat_fsi)); - su2double DeltaT = config->GetDelta_DynTime(); + su2double DeltaT = config->GetDelta_UnstTime(); su2double OffsetTime = offset? DeltaT * (config->GetRestart_Iter()-1) : su2double(0.0); /*--- Polynomial functions from https://en.wikipedia.org/wiki/Smoothstep ---*/ @@ -2507,7 +2508,7 @@ void CFEASolver::Solve_System(CGeometry *geometry, CConfig *config) { void CFEASolver::PredictStruct_Displacement(CGeometry *geometry, const CConfig *config) { const unsigned short predOrder = config->GetPredictorOrder(); - const su2double Delta_t = config->GetDelta_DynTime(); + const su2double Delta_t = config->GetDelta_UnstTime(); const bool dynamic = config->GetTime_Domain(); if(predOrder > 2 && rank == MASTER_NODE) diff --git a/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg b/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg index c4181d6339d..1df1677212b 100644 --- a/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg +++ b/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg @@ -18,7 +18,7 @@ RESTART_ITER= 1 ELASTICITY_MODULUS=3E7 POISSON_RATIO=0.3 MATERIAL_DENSITY=7854 -FORMULATION_ELASTICITY_2D = PLANE_STRESS +FORMULATION_ELASTICITY_2D= PLANE_STRESS TIME_DOMAIN=YES TIME_STEP=0.01 MAX_TIME= 0.1 @@ -27,17 +27,18 @@ TIME_ITER=7 TIME_DISCRE_FEA= NEWMARK_IMPLICIT NEWMARK_BETA=0.2601 NEWMARK_GAMMA=0.52 -MARKER_CLAMPED = ( left , right ) -MARKER_PRESSURE= ( lower, 0) -MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0) -LINEAR_SOLVER= FGMRES -LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-8 -LINEAR_SOLVER_ITER= 1000 +MARKER_CLAMPED= ( left , right ) +MARKER_PRESSURE= ( lower, 0 ) +MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0 ) +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 100 MESH_FORMAT= SU2 TABULAR_FORMAT= CSV VOLUME_FILENAME= beam RESTART_FILENAME= restart_beam.dat SOLUTION_FILENAME= solution_beam.dat OUTPUT_WRT_FREQ= 1 -OUTPUT_FILES= (RESTART_ASCII) +OUTPUT_FILES= RESTART_ASCII, PARAVIEW +SCREEN_OUTPUT= CUR_TIME, INNER_ITER, RMS_RES, LINSOL From 3430b3a39b7a39ce1aa0b1936f1667d058851027 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 8 Apr 2023 10:20:09 -0700 Subject: [PATCH 02/20] fix relaxation issue and add testcase --- SU2_CFD/src/iteration/CFEAIteration.cpp | 17 ++--- .../py_wrapper/custom_load_fea/config.cfg | 32 +++++++++ TestCases/py_wrapper/custom_load_fea/run.py | 70 +++++++++++++++++++ 3 files changed, 108 insertions(+), 11 deletions(-) create mode 100644 TestCases/py_wrapper/custom_load_fea/config.cfg create mode 100755 TestCases/py_wrapper/custom_load_fea/run.py diff --git a/SU2_CFD/src/iteration/CFEAIteration.cpp b/SU2_CFD/src/iteration/CFEAIteration.cpp index c827f28467e..8390fdab76a 100644 --- a/SU2_CFD/src/iteration/CFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CFEAIteration.cpp @@ -184,21 +184,16 @@ void CFEAIteration::Update(COutput* output, CIntegration**** integration, CGeome CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { - const bool dynamic = (config[val_iZone]->GetTime_Domain()); - const bool fsi = config[val_iZone]->GetFSI_Simulation(); - CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; /*----------------- Update structural solver ----------------------*/ - if (dynamic) { - integration[val_iZone][val_iInst][FEA_SOL]->SetDualTime_Solver( - geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0][FEA_SOL], config[val_iZone], - MESH_0); - - } else if (fsi) { - /*--- For FSI problems, output the relaxed result, which is the one transferred into the fluid domain (for restart - * purposes) ---*/ + if (config[val_iZone]->GetTime_Domain()) { + integration[val_iZone][val_iInst][FEA_SOL]->SetDualTime_Solver(geometry[val_iZone][val_iInst][MESH_0], feaSolver, + config[val_iZone], MESH_0); + } else if (config[val_iZone]->GetFSI_Simulation() && config[val_iZone]->GetRelaxation()) { + /*--- For FSI problems with relaxation, output the relaxed result, which is the one transferred into the fluid + * domain (for consistent restart purposes). ---*/ if (config[val_iZone]->GetKind_TimeIntScheme_FEA() == STRUCT_TIME_INT::NEWMARK_IMPLICIT) { feaSolver->ImplicitNewmark_Relaxation(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } diff --git a/TestCases/py_wrapper/custom_load_fea/config.cfg b/TestCases/py_wrapper/custom_load_fea/config.cfg new file mode 100644 index 00000000000..c0352190d77 --- /dev/null +++ b/TestCases/py_wrapper/custom_load_fea/config.cfg @@ -0,0 +1,32 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: 2D Beam with custom load via Python wrapper % +% File Version 7.5.1 "Blackbird" % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= ELASTICITY +GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS +FORMULATION_ELASTICITY_2D= PLANE_STRESS + +MATERIAL_MODEL= NEO_HOOKEAN +ELASTICITY_MODULUS= 10000 +POISSON_RATIO= 0.3 + +MARKER_CLAMPED= ( x_minus ) +MARKER_FLUID_LOAD= ( x_plus, y_minus, y_plus ) + +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 100 + +MESH_FORMAT= RECTANGLE +MESH_BOX_SIZE= ( 17, 5, 0 ) +MESH_BOX_LENGTH= ( 0.5, 0.05, 0 ) + +OUTPUT_FILES= RESTART_ASCII, PARAVIEW +SCREEN_OUTPUT= INNER_ITER, RMS_RES, LINSOL, VMS + +INNER_ITER= 20 +CONV_FIELD= REL_RMS_RTOL +CONV_RESIDUAL_MINVAL= -6 diff --git a/TestCases/py_wrapper/custom_load_fea/run.py b/TestCases/py_wrapper/custom_load_fea/run.py new file mode 100755 index 00000000000..c1ac289b4de --- /dev/null +++ b/TestCases/py_wrapper/custom_load_fea/run.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python + +## \file run.py +# \brief FEA case with custom load. +# \version 7.5.1 "Blackbird" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2023, 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 . + +import pysu2 +from mpi4py import MPI +import numpy as np + +def main(): + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + + # Initialize the corresponding driver of SU2, this includes solver preprocessing. + try: + SU2Driver = pysu2.CSinglezoneDriver('config.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2.CDriver : ', exception) + raise + + # Get the ID of the marker we want to deform. + AllMarkerIDs = SU2Driver.GetMarkerIndices() + MarkerName = 'y_minus' + MarkerID = AllMarkerIDs[MarkerName] if MarkerName in AllMarkerIDs else -1 + + # Number of vertices on the specified marker (per rank). + nVertex = SU2Driver.GetNumberMarkerNodes(MarkerID) if MarkerID >= 0 else 0 + + # Apply a load based on the coordinates. + MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) + L = 0.5 + dx = L / 16 # known from mesh settings in this case. + for iVertex in range(nVertex): + x = MarkerCoords(iVertex, 0) + nodalForce = (2 * x / L) * dx + # Half load due to half dx on first and last node. + if abs(x) < 1e-6 or abs(x - L) < 1e-6: + nodalForce = nodalForce / 2 + SU2Driver.SetMarkerCustomFEALoad(MarkerID, iVertex, (0, nodalForce)) + + # Solve. + SU2Driver.StartSolver() + + # Finalize the solver and exit cleanly. + SU2Driver.Finalize() + + +if __name__ == '__main__': + main() From 6f2424d419c60dac0ac850ea7b6adfb9a65f9b0f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 8 Apr 2023 10:50:48 -0700 Subject: [PATCH 03/20] parallel regression --- TestCases/parallel_regression.py | 9 ++++++++ TestCases/py_wrapper/custom_load_fea/run.py | 23 +++++++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 1b359311076..6973bdfb075 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1300,6 +1300,15 @@ def main(): pywrapper_aeroelastic.unsteady = True test_list.append(pywrapper_aeroelastic) + # Custom FEA load + pywrapper_custom_fea_load = TestCase('pywrapper_custom_fea_load') + pywrapper_custom_fea_load.cfg_dir = "py_wrapper/custom_load_fea" + pywrapper_custom_fea_load.cfg_file = "config.cfg" + pywrapper_custom_fea_load.test_iter = 13 + pywrapper_custom_fea_load.test_vals = [-7.263559, -4.946814, -14.165142, 34, -6.385218, 3.2058e+02] + pywrapper_custom_fea_load.command = TestCase.Command("mpirun -np 2", "python", "run.py") + test_list.append(pywrapper_custom_fea_load) + # FSI, 2d pywrapper_fsi2d = TestCase('pywrapper_fsi2d') pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" diff --git a/TestCases/py_wrapper/custom_load_fea/run.py b/TestCases/py_wrapper/custom_load_fea/run.py index c1ac289b4de..58733c9c504 100755 --- a/TestCases/py_wrapper/custom_load_fea/run.py +++ b/TestCases/py_wrapper/custom_load_fea/run.py @@ -62,6 +62,29 @@ def main(): # Solve. SU2Driver.StartSolver() + # Find the tip displacement. + MarkerName = 'x_plus' + MarkerID = AllMarkerIDs[MarkerName] if MarkerName in AllMarkerIDs else -1 + nVertex = SU2Driver.GetNumberMarkerNodes(MarkerID) if MarkerID >= 0 else 0 + MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) + + SolverID = SU2Driver.GetSolverIndices()["FEA"] + Solution = SU2Driver.MarkerSolution(SolverID, MarkerID) + DispID = SU2Driver.GetFEASolutionIndices()["DISPLACEMENT_Y"] + + Disp = 0 + NodeFound = False + for iVertex in range(nVertex): + y = MarkerCoords(iVertex, 1) + if abs(y - 0.025) < 1e-6: + Disp = Solution(iVertex, DispID) + NodeFound = True + + if NodeFound: + print(f"Vertical displacement of tip: {Disp}") + # Test the value against expected. + assert abs(Disp / 0.095439 - 1) < 1e-5, "Test FAILED" + # Finalize the solver and exit cleanly. SU2Driver.Finalize() From 30f6c1b2833b8302e9b15a97252963c0dbd79710 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 8 Apr 2023 13:26:39 -0700 Subject: [PATCH 04/20] cleanup configs --- TestCases/fea_fsi/MixElemsKnowles/config.cfg | 1 - TestCases/fea_fsi/SquareCyl_Beam/config.cfg | 2 +- TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg | 3 +-- TestCases/fea_fsi/stat_fsi/configFEA.cfg | 2 +- TestCases/fea_topology/config.cfg | 1 - TestCases/fea_topology/quick_start/settings.cfg | 1 - TestCases/fea_topology/quick_start/settings_compliance.cfg | 1 - TestCases/fea_topology/quick_start/settings_volfrac.cfg | 1 - 8 files changed, 3 insertions(+), 9 deletions(-) diff --git a/TestCases/fea_fsi/MixElemsKnowles/config.cfg b/TestCases/fea_fsi/MixElemsKnowles/config.cfg index 94e5328e4b9..ffa4573286a 100644 --- a/TestCases/fea_fsi/MixElemsKnowles/config.cfg +++ b/TestCases/fea_fsi/MixElemsKnowles/config.cfg @@ -11,7 +11,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DIRECT -DYNAMIC_ANALYSIS= NO % % Optimization OBJECTIVE_FUNCTION= REFERENCE_NODE diff --git a/TestCases/fea_fsi/SquareCyl_Beam/config.cfg b/TestCases/fea_fsi/SquareCyl_Beam/config.cfg index ab740fe4255..d4033e5f4b1 100644 --- a/TestCases/fea_fsi/SquareCyl_Beam/config.cfg +++ b/TestCases/fea_fsi/SquareCyl_Beam/config.cfg @@ -67,7 +67,7 @@ NONLINEAR_FEM_INT_ITER = 20 % -------------------------- DYNAMIC SIMULATION -------------------------------% % -DYNAMIC_ANALYSIS= YES +TIME_DOMAIN= YES TIME_DISCRE_FEA= NEWMARK_IMPLICIT NEWMARK_BETA=0.2601 NEWMARK_GAMMA=0.52 diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg index 1aee4275ed5..e3b82a9e601 100644 --- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg +++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg @@ -6,7 +6,7 @@ % Date: 2016.02.01 % % File Version 7.5.1 "Blackbird" % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + SOLVER= ELASTICITY MATH_PROBLEM= DIRECT GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS @@ -15,7 +15,6 @@ MESH_FILENAME= meshBeam_3d.su2 ELASTICITY_MODULUS=3E7 POISSON_RATIO=0.3 MATERIAL_DENSITY=7854 -DYNAMIC_ANALYSIS= NO MARKER_CLAMPED = ( left , right ) MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0) MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0) diff --git a/TestCases/fea_fsi/stat_fsi/configFEA.cfg b/TestCases/fea_fsi/stat_fsi/configFEA.cfg index 3ce2b13322c..0f4e66ad1f7 100755 --- a/TestCases/fea_fsi/stat_fsi/configFEA.cfg +++ b/TestCases/fea_fsi/stat_fsi/configFEA.cfg @@ -59,7 +59,7 @@ FORMULATION_ELASTICITY_2D = PLANE_STRAIN POISSON_RATIO=0.4 % -------------------------- DYNAMIC SIMULATION -------------------------------% -DYNAMIC_ANALYSIS= NO +TIME_DOMAIN= NO TIME_DISCRE_FEA= NEWMARK_IMPLICIT % -------------------------- STRUCTURAL SOLVER --------------------------------% diff --git a/TestCases/fea_topology/config.cfg b/TestCases/fea_topology/config.cfg index cf12a17353a..021e878895d 100644 --- a/TestCases/fea_topology/config.cfg +++ b/TestCases/fea_topology/config.cfg @@ -66,7 +66,6 @@ ITER=1 % Physics SOLVER= ELASTICITY MATH_PROBLEM= DISCRETE_ADJOINT -DYNAMIC_ANALYSIS= NO RESTART_SOL= NO % SOLUTION_FILENAME=solution_structure.dat diff --git a/TestCases/fea_topology/quick_start/settings.cfg b/TestCases/fea_topology/quick_start/settings.cfg index 44d8b3eafaf..8a0edb1041a 100644 --- a/TestCases/fea_topology/quick_start/settings.cfg +++ b/TestCases/fea_topology/quick_start/settings.cfg @@ -1,7 +1,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DIRECT -DYNAMIC_ANALYSIS= NO % ITER= 1 % diff --git a/TestCases/fea_topology/quick_start/settings_compliance.cfg b/TestCases/fea_topology/quick_start/settings_compliance.cfg index d9cf275df75..57241e0615f 100644 --- a/TestCases/fea_topology/quick_start/settings_compliance.cfg +++ b/TestCases/fea_topology/quick_start/settings_compliance.cfg @@ -1,7 +1,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DISCRETE_ADJOINT -DYNAMIC_ANALYSIS= NO % ITER= 1 % diff --git a/TestCases/fea_topology/quick_start/settings_volfrac.cfg b/TestCases/fea_topology/quick_start/settings_volfrac.cfg index a3ec2ebc3cb..c7da2fe14f5 100644 --- a/TestCases/fea_topology/quick_start/settings_volfrac.cfg +++ b/TestCases/fea_topology/quick_start/settings_volfrac.cfg @@ -1,7 +1,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DISCRETE_ADJOINT -DYNAMIC_ANALYSIS= NO % ITER= 1 % From 610fa17e7b2415dad7deb7078ac28beb272908db Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 8 Apr 2023 19:41:20 -0700 Subject: [PATCH 05/20] fix testcases --- .../fea_fsi/DynBeam_2d/configBeam_2d.cfg | 3 +- TestCases/parallel_regression.py | 2 +- TestCases/py_wrapper/custom_load_fea/run.py | 46 +++++++++---------- 3 files changed, 25 insertions(+), 26 deletions(-) diff --git a/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg b/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg index 1df1677212b..1303f5457be 100644 --- a/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg +++ b/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg @@ -6,7 +6,7 @@ % Date: 2016.02.01 % % File Version 7.5.1 "Blackbird" % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + SOLVER= ELASTICITY MATH_PROBLEM= DIRECT GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS @@ -41,4 +41,3 @@ RESTART_FILENAME= restart_beam.dat SOLUTION_FILENAME= solution_beam.dat OUTPUT_WRT_FREQ= 1 OUTPUT_FILES= RESTART_ASCII, PARAVIEW -SCREEN_OUTPUT= CUR_TIME, INNER_ITER, RMS_RES, LINSOL diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 6973bdfb075..4154ebc1a99 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1305,7 +1305,7 @@ def main(): pywrapper_custom_fea_load.cfg_dir = "py_wrapper/custom_load_fea" pywrapper_custom_fea_load.cfg_file = "config.cfg" pywrapper_custom_fea_load.test_iter = 13 - pywrapper_custom_fea_load.test_vals = [-7.263559, -4.946814, -14.165142, 34, -6.385218, 3.2058e+02] + pywrapper_custom_fea_load.test_vals = [-7.263559, -4.946814, -14.165142, 34, -6.380144, 320.58] pywrapper_custom_fea_load.command = TestCase.Command("mpirun -np 2", "python", "run.py") test_list.append(pywrapper_custom_fea_load) diff --git a/TestCases/py_wrapper/custom_load_fea/run.py b/TestCases/py_wrapper/custom_load_fea/run.py index 58733c9c504..24881ef85df 100755 --- a/TestCases/py_wrapper/custom_load_fea/run.py +++ b/TestCases/py_wrapper/custom_load_fea/run.py @@ -26,11 +26,9 @@ import pysu2 from mpi4py import MPI -import numpy as np def main(): comm = MPI.COMM_WORLD - rank = comm.Get_rank() # Initialize the corresponding driver of SU2, this includes solver preprocessing. try: @@ -48,16 +46,17 @@ def main(): nVertex = SU2Driver.GetNumberMarkerNodes(MarkerID) if MarkerID >= 0 else 0 # Apply a load based on the coordinates. - MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) - L = 0.5 - dx = L / 16 # known from mesh settings in this case. - for iVertex in range(nVertex): - x = MarkerCoords(iVertex, 0) - nodalForce = (2 * x / L) * dx - # Half load due to half dx on first and last node. - if abs(x) < 1e-6 or abs(x - L) < 1e-6: - nodalForce = nodalForce / 2 - SU2Driver.SetMarkerCustomFEALoad(MarkerID, iVertex, (0, nodalForce)) + if nVertex > 0: + MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) + L = 0.5 + dx = L / 16 # known from mesh settings in this case. + for iVertex in range(nVertex): + x = MarkerCoords(iVertex, 0) + nodalForce = (2 * x / L) * dx + # Half load due to half dx on first and last node. + if abs(x) < 1e-6 or abs(x - L) < 1e-6: + nodalForce = nodalForce / 2 + SU2Driver.SetMarkerCustomFEALoad(MarkerID, iVertex, (0, nodalForce)) # Solve. SU2Driver.StartSolver() @@ -66,19 +65,20 @@ def main(): MarkerName = 'x_plus' MarkerID = AllMarkerIDs[MarkerName] if MarkerName in AllMarkerIDs else -1 nVertex = SU2Driver.GetNumberMarkerNodes(MarkerID) if MarkerID >= 0 else 0 - MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) - - SolverID = SU2Driver.GetSolverIndices()["FEA"] - Solution = SU2Driver.MarkerSolution(SolverID, MarkerID) - DispID = SU2Driver.GetFEASolutionIndices()["DISPLACEMENT_Y"] - Disp = 0 NodeFound = False - for iVertex in range(nVertex): - y = MarkerCoords(iVertex, 1) - if abs(y - 0.025) < 1e-6: - Disp = Solution(iVertex, DispID) - NodeFound = True + + if nVertex > 0: + MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) + SolverID = SU2Driver.GetSolverIndices()["FEA"] + Solution = SU2Driver.MarkerSolution(SolverID, MarkerID) + DispID = SU2Driver.GetFEASolutionIndices()["DISPLACEMENT_Y"] + + for iVertex in range(nVertex): + y = MarkerCoords(iVertex, 1) + if abs(y - 0.025) < 1e-6: + Disp = Solution(iVertex, DispID) + NodeFound = True if NodeFound: print(f"Vertical displacement of tip: {Disp}") From 6b6fccb26f8ddf23af96a20ab5b9f28025d13011 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 9 Apr 2023 18:54:32 -0700 Subject: [PATCH 06/20] do not call Newmark relaxation also for transient when relaxation is off --- SU2_CFD/include/solvers/CFEASolver.hpp | 23 ++----------------- SU2_CFD/include/solvers/CSolver.hpp | 21 +---------------- .../integration/CStructuralIntegration.cpp | 13 +++-------- SU2_CFD/src/solvers/CFEASolver.cpp | 14 +++++------ 4 files changed, 13 insertions(+), 58 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index ed0aee60236..86eb0813e6a 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -682,28 +682,9 @@ class CFEASolver : public CFEASolverBase { inline su2double GetFSI_ConvValue(unsigned short val_index) const final { return FSI_Conv[val_index]; } /*! - * \brief Retrieve the value of the dynamic Aitken relaxation factor. - * \return Value of the dynamic Aitken relaxation factor. + * \brief Store the value of the last Aitken relaxation factor in the current time step. */ - inline su2double GetWAitken_Dyn(void) const final { return WAitken_Dyn; } - - /*! - * \brief Retrieve the value of the last Aitken relaxation factor in the previous time step. - * \return Value of the last Aitken relaxation factor in the previous time step. - */ - inline su2double GetWAitken_Dyn_tn1(void) const final { return WAitken_Dyn_tn1; } - - /*! - * \brief Set the value of the dynamic Aitken relaxation factor - * \param[in] Value of the dynamic Aitken relaxation factor - */ - inline void SetWAitken_Dyn(su2double waitk) final { WAitken_Dyn = waitk; } - - /*! - * \brief Set the value of the last Aitken relaxation factor in the current time step. - * \param[in] Value of the last Aitken relaxation factor in the current time step. - */ - inline void SetWAitken_Dyn_tn1(su2double waitk_tn1) final { WAitken_Dyn_tn1 = waitk_tn1; } + inline void SetWAitken_Dyn_tn1() final { WAitken_Dyn_tn1 = WAitken_Dyn; } /*! * \brief Set the value of the load increment for nonlinear structural analysis diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index cdfa1946151..892552859b0 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3698,27 +3698,8 @@ class CSolver { /*! * \brief A virtual member. - * \return Value of the dynamic Aitken relaxation factor */ - inline virtual su2double GetWAitken_Dyn(void) const { return 0; } - - /*! - * \brief A virtual member. - * \return Value of the last Aitken relaxation factor in the previous time step. - */ - inline virtual su2double GetWAitken_Dyn_tn1(void) const { return 0; } - - /*! - * \brief A virtual member. - * \param[in] Value of the dynamic Aitken relaxation factor - */ - inline virtual void SetWAitken_Dyn(su2double waitk) { } - - /*! - * \brief A virtual member. - * \param[in] Value of the last Aitken relaxation factor in the previous time step. - */ - inline virtual void SetWAitken_Dyn_tn1(su2double waitk_tn1) { } + inline virtual void SetWAitken_Dyn_tn1() { } /*! * \brief A virtual member. diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp index 1dd5573923c..04b207f6fa0 100644 --- a/SU2_CFD/src/integration/CStructuralIntegration.cpp +++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp @@ -193,7 +193,7 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * void CStructuralIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver, const CConfig *config, unsigned short iMesh) { - bool fsi = config->GetFSI_Simulation(); + const bool fsi = config->GetFSI_Simulation(); /*--- Update the solution according to the integration scheme used ---*/ @@ -201,7 +201,7 @@ void CStructuralIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolv case (STRUCT_TIME_INT::CD_EXPLICIT): break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): - if (fsi) solver->ImplicitNewmark_Relaxation(geometry, config); + if (fsi && config->GetRelaxation()) solver->ImplicitNewmark_Relaxation(geometry, config); break; case (STRUCT_TIME_INT::GENERALIZED_ALPHA): solver->GeneralizedAlpha_UpdateSolution(geometry, config); @@ -215,12 +215,5 @@ void CStructuralIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolv /*--- If FSI problem, save the last Aitken relaxation parameter of the previous time step ---*/ - if (fsi) { - - su2double WAitk=0.0; - - WAitk = solver->GetWAitken_Dyn(); - solver->SetWAitken_Dyn_tn1(WAitk); - - } + if (fsi) solver->SetWAitken_Dyn_tn1(); } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index c7a1b4b08e4..74870622ea0 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2576,26 +2576,26 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c if (RelaxMethod_FSI == BGS_RELAXATION::NONE) { - SetWAitken_Dyn(1.0); + WAitken_Dyn = 1.0; } else if (RelaxMethod_FSI == BGS_RELAXATION::FIXED) { - SetWAitken_Dyn(config->GetAitkenStatRelax()); + WAitken_Dyn = config->GetAitkenStatRelax(); } else if (RelaxMethod_FSI == BGS_RELAXATION::AITKEN) { if (iOuterIter == 0) { - WAitkDyn_tn1 = GetWAitken_Dyn_tn1(); + WAitkDyn_tn1 = WAitken_Dyn_tn1; WAitkDyn_Max = config->GetAitkenDynMaxInit(); WAitkDyn_Min = config->GetAitkenDynMinInit(); WAitkDyn = min(WAitkDyn_tn1, WAitkDyn_Max); WAitkDyn = max(WAitkDyn, WAitkDyn_Min); - SetWAitken_Dyn(WAitkDyn); + WAitken_Dyn = WAitkDyn; } else { @@ -2626,7 +2626,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c SU2_MPI::Allreduce(&sbuf_numAitk, &rbuf_numAitk, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); SU2_MPI::Allreduce(&sbuf_denAitk, &rbuf_denAitk, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - WAitkDyn = GetWAitken_Dyn(); + WAitkDyn = WAitken_Dyn; if (rbuf_denAitk > EPS) { WAitkDyn = - 1.0 * WAitkDyn * rbuf_numAitk / rbuf_denAitk ; @@ -2635,7 +2635,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c WAitkDyn = max(WAitkDyn, 0.1); WAitkDyn = min(WAitkDyn, 1.0); - SetWAitken_Dyn(WAitkDyn); + WAitken_Dyn = WAitkDyn; } @@ -2648,7 +2648,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c void CFEASolver::SetAitken_Relaxation(CGeometry *geometry, const CConfig *config) { - const su2double WAitken = GetWAitken_Dyn(); + const su2double WAitken = WAitken_Dyn; const bool dynamic = config->GetTime_Domain(); /*--- To nPoint to avoid communication. ---*/ From 0acec215cdb9eb3248126020799c288daee8ab24 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 10 Apr 2023 09:15:02 -0700 Subject: [PATCH 07/20] cleanup limiter docs --- config_template.cfg | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/config_template.cfg b/config_template.cfg index de994a0ce60..a99ec1a8922 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -733,7 +733,7 @@ CONV_NUM_METHOD_SPECIES= SCALAR_UPWIND % Required for 2nd order upwind schemes (NO, YES) MUSCL_SPECIES= NO % -% Slope limiter for species equations (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE) +% Slope limiter for species equations (same as SLOPE_LIMITER_TURB) SLOPE_LIMITER_SPECIES = NONE % % Time discretization for species equations (EULER_IMPLICIT, EULER_EXPLICIT) @@ -1095,36 +1095,37 @@ CUSTOM_OBJFUNC= 'DRAG + 10 * pow(fmax(0.4-LIFT, 0), 2)' % Required for 2nd order upwind schemes (NO, YES) MUSCL_FLOW= YES % -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE) -% +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE, +% NISHIKAWA_R3, NISHIKAWA_R4, NISHIKAWA_R5) SLOPE_LIMITER_FLOW= VENKATAKRISHNAN % -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. -% Required for 2nd order upwind schemes (NO, YES) +% Same as MUSCL_FLOW but for turbulence. +% MUSCL_TURB= NO % -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN) +% Slope limiter (same as SLOPE_LIMITER_FLOW except VAN_ALBADA_EDGE) % SLOPE_LIMITER_TURB= VENKATAKRISHNAN % -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the adjoint flow equations. -% Required for 2nd order upwind schemes (NO, YES) +% Same as MUSCL_FLOW but for the continuous adjoint equations. +% MUSCL_ADJFLOW= YES % -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE, -% SHARP_EDGES, WALL_DISTANCE, NISHIKAWA_R3, NISHIKAWA_R4, NISHIKAWA_R5) +% Slope limiter (same as SLOPE_LIMITER_FLOW plus SHARP_EDGES, WALL_DISTANCE) +% SLOPE_LIMITER_ADJFLOW= VENKATAKRISHNAN % -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence adjoint equations. -% Required for 2nd order upwind schemes (NO, YES) +% Same as MUSCL_FLOW but for continuous adjoint turbulence equations. +% MUSCL_ADJTURB= NO % % Slope limiter (see SLOPE_LIMITER_ADJFLOW) +% SLOPE_LIMITER_ADJTURB= VENKATAKRISHNAN % -% Coefficient for the Venkat's limiter (upwind scheme). A larger values decrease -% the extent of limiting, values approaching zero cause -% lower-order approximation to the solution (0.05 by default) +% Coefficient for Venkatakrishnan-type limiters (upwind scheme). +% A larger value decreases the extent of limiting, values approaching zero +% cause lower-order approximation to the solution (0.05 by default) VENKAT_LIMITER_COEFF= 0.05 % % Reference coefficient for detecting sharp edges (3.0 by default). From 9782af92b3e71921e3c6f9330d9a3c9c1839340c Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Apr 2023 15:02:45 -0700 Subject: [PATCH 08/20] cleanup unsused code --- Common/include/CConfig.hpp | 6 ------ SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp | 12 +++--------- 2 files changed, 3 insertions(+), 15 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index b01d4d1c5ca..f7b57a2445b 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -9251,12 +9251,6 @@ class CConfig { */ void SetnTime_Iter(unsigned long val_iter) { nTimeIter = val_iter; } - /*! - * \brief Get the number of pseudo-time iterations - * \return Number of pseudo-time steps run for the single-zone problem - */ - unsigned long GetnIter(void) const { return nIter; } - /*! * \brief Get the restart iteration * \return Iteration for the restart of multizone problems diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index 4a6e61dfc2c..75ad91681f3 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -34,7 +34,7 @@ CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(c fem_iteration = new CFEAIteration(config); // TEMPORARY output only for standalone structural problems - if ((!config->GetFSI_Simulation()) && (rank == MASTER_NODE)) { + if (config->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) { bool de_effects = config->GetDE_Effects(); unsigned short iVar; @@ -56,9 +56,7 @@ CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(c for (iVar = 0; iVar < config->GetnElectric_Field(); iVar++) myfile_res << "Sens_EField_" << iVar << "\t"; } - myfile_res << endl; - - myfile_res.close(); + myfile_res << '\n'; } } @@ -365,25 +363,21 @@ void CDiscAdjFEAIteration::Postprocess(COutput* output, CIntegration**** integra if (config[iZone]->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) { /*--- Header of the temporary output file ---*/ ofstream myfile_res; - bool outputDVFEA = false; + bool outputDVFEA = true; switch (config[iZone]->GetDV_FEA()) { case YOUNG_MODULUS: myfile_res.open("grad_young.opt"); - outputDVFEA = true; break; case POISSON_RATIO: myfile_res.open("grad_poisson.opt"); - outputDVFEA = true; break; case DENSITY_VAL: case DEAD_WEIGHT: myfile_res.open("grad_density.opt"); - outputDVFEA = true; break; case ELECTRIC_FIELD: myfile_res.open("grad_efield.opt"); - outputDVFEA = true; break; default: outputDVFEA = false; From 45243a7cc48ca4f185f5bccadbcab855c7f29918 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Apr 2023 18:47:11 -0700 Subject: [PATCH 09/20] fix structural unsteady adjoints --- .../src/drivers/CDiscAdjMultizoneDriver.cpp | 14 ++-- .../src/drivers/CDiscAdjSinglezoneDriver.cpp | 18 +++-- .../integration/CStructuralIntegration.cpp | 6 +- .../src/iteration/CDiscAdjFEAIteration.cpp | 65 ++++++------------- SU2_CFD/src/iteration/CFEAIteration.cpp | 36 ++++------ SU2_CFD/src/solvers/CFEASolver.cpp | 18 +++-- 6 files changed, 63 insertions(+), 94 deletions(-) diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 586c9cabfdc..0a2a77ae31f 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -763,15 +763,15 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { } void CDiscAdjMultizoneDriver::SetAdjObjFunction() { - - const auto IterAvg_Obj = config_container[ZONE_0]->GetIter_Avg_Objective(); su2double seeding = 1.0; - if (config_container[ZONE_0]->GetTime_Marching() != TIME_MARCHING::STEADY){ - if (TimeIter < IterAvg_Obj){ - // Default behavior (in case no specific window is chosen) is to use Square-Windowing, i.e. the numerator equals 1.0 + if (config_container[ZONE_0]->GetTime_Domain()) { + const auto IterAvg_Obj = config_container[ZONE_0]->GetIter_Avg_Objective(); + if (TimeIter < IterAvg_Obj) { + /*--- Default behavior when no window is chosen is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ auto windowEvaluator = CWindowingTools(); - su2double weight = windowEvaluator.GetWndWeight(config_container[ZONE_0]->GetKindWindow(), TimeIter, IterAvg_Obj-1); + const su2double weight = + windowEvaluator.GetWndWeight(config_container[ZONE_0]->GetKindWindow(), TimeIter, IterAvg_Obj - 1); seeding = weight / IterAvg_Obj; } else { @@ -780,6 +780,8 @@ void CDiscAdjMultizoneDriver::SetAdjObjFunction() { } if (rank == MASTER_NODE) { AD::SetDerivative(ObjFunc_Index, SU2_TYPE::GetValue(seeding)); + } else { + AD::SetDerivative(ObjFunc_Index, 0.0); } } diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index 5a42be040d7..f9cfb269a8e 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -319,23 +319,21 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ } void CDiscAdjSinglezoneDriver::SetAdjObjFunction(){ - - const auto IterAvg_Obj = config->GetIter_Avg_Objective(); su2double seeding = 1.0; - CWindowingTools windowEvaluator = CWindowingTools(); - - if (config->GetTime_Marching() != TIME_MARCHING::STEADY){ - if (TimeIter < IterAvg_Obj){ - /*--- Default behavior (in case no specific window is chosen) is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ - seeding = windowEvaluator.GetWndWeight(config->GetKindWindow(),TimeIter, IterAvg_Obj-1)/ (static_cast(IterAvg_Obj)); + if (config->GetTime_Domain()) { + const auto IterAvg_Obj = config->GetIter_Avg_Objective(); + if (TimeIter < IterAvg_Obj) { + /*--- Default behavior when no window is chosen is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ + auto windowEvaluator = CWindowingTools(); + const su2double weight = windowEvaluator.GetWndWeight(config->GetKindWindow(), TimeIter, IterAvg_Obj - 1); + seeding = weight / IterAvg_Obj; } else { seeding = 0.0; } } - - if (rank == MASTER_NODE){ + if (rank == MASTER_NODE) { SU2_TYPE::SetDerivative(ObjFunc, SU2_TYPE::GetValue(seeding)); } else { SU2_TYPE::SetDerivative(ObjFunc, 0.0); diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp index 04b207f6fa0..276cf4a01bf 100644 --- a/SU2_CFD/src/integration/CStructuralIntegration.cpp +++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp @@ -135,8 +135,6 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * switch (config->GetKind_TimeIntScheme_FEA()) { case (STRUCT_TIME_INT::CD_EXPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, numerics, config); - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, numerics, config); break; @@ -161,7 +159,7 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * } } - /*--- Solver linear system ---*/ + /*--- Solve linear system ---*/ solver_container[MainSolver]->Solve_System(geometry, config); @@ -169,8 +167,6 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * switch (config->GetKind_TimeIntScheme_FEA()) { case (STRUCT_TIME_INT::CD_EXPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Update(geometry, config); - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): solver_container[MainSolver]->ImplicitNewmark_Update(geometry, config); break; diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index 75ad91681f3..7decbd05b57 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -56,7 +56,7 @@ CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(c for (iVar = 0; iVar < config->GetnElectric_Field(); iVar++) myfile_res << "Sens_EField_" << iVar << "\t"; } - myfile_res << '\n'; + myfile_res << "\n"; } } @@ -66,7 +66,6 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short iZone, unsigned short iInst) { - unsigned long iPoint; auto solvers0 = solver[iZone][iInst][MESH_0]; auto geometry0 = geometry[iZone][iInst][MESH_0]; auto dirNodes = solvers0[FEA_SOL]->GetNodes(); @@ -91,19 +90,12 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat /*--- Load solution timestep n ---*/ LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter); + } - /*--- Store FEA solution also in the adjoint solver in order to be able to reset it later ---*/ - - for (iPoint = 0; iPoint < geometry0->GetnPoint(); iPoint++) { - adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); - } + /*--- Store FEA solution also in the adjoint solver in order to be able to reset it later. ---*/ - } else { - /*--- Store FEA solution also in the adjoint solver in order to be able to reset it later ---*/ - - for (iPoint = 0; iPoint < geometry0->GetnPoint(); iPoint++) { - adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); - } + for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) { + adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); } solvers0[ADJFEA_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, 0, RUNTIME_ADJFEA_SYS, false); @@ -113,27 +105,19 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat void CDiscAdjFEAIteration::LoadDynamic_Solution(CGeometry**** geometry, CSolver***** solver, CConfig** config, unsigned short iZone, unsigned short iInst, int val_DirectIter) { - unsigned short iVar; - unsigned long iPoint; - bool update_geo = false; // TODO: check + /*--- Set to false to prevent updating Solution_time_n when loading primal solutions of unsteady cases. ---*/ + const bool update_geo = false; + auto*** solvers = solver[iZone][iInst]; if (val_DirectIter >= 0) { if (rank == MASTER_NODE && iZone == ZONE_0) - cout << " Loading FEA solution from direct iteration " << val_DirectIter << "." << endl; - solver[iZone][iInst][MESH_0][FEA_SOL]->LoadRestart( - geometry[iZone][iInst], solver[iZone][iInst], config[iZone], val_DirectIter, update_geo); + cout << " Loading FEA solution from direct iteration " << val_DirectIter << ".\n"; + solvers[MESH_0][FEA_SOL]->LoadRestart(geometry[iZone][iInst], solvers, config[iZone], val_DirectIter, update_geo); } else { - /*--- If there is no solution file we set the freestream condition ---*/ + /*--- If there is no solution file we set the initial conditions. ---*/ if (rank == MASTER_NODE && iZone == ZONE_0) - cout << " Setting static conditions at direct iteration " << val_DirectIter << "." << endl; - /*--- Push solution back to correct array ---*/ - for (iPoint = 0; iPoint < geometry[iZone][iInst][MESH_0]->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < solver[iZone][iInst][MESH_0][FEA_SOL]->GetnVar(); iVar++) { - solver[iZone][iInst][MESH_0][FEA_SOL]->GetNodes()->SetSolution(iPoint, iVar, 0.0); - solver[iZone][iInst][MESH_0][FEA_SOL]->GetNodes()->SetSolution_Accel(iPoint, iVar, 0.0); - solver[iZone][iInst][MESH_0][FEA_SOL]->GetNodes()->SetSolution_Vel(iPoint, iVar, 0.0); - } - } + cout << " Setting static conditions at direct iteration " << val_DirectIter << ".\n"; + solvers[MESH_0][FEA_SOL]->SetInitialCondition(geometry[iZone][iInst], solvers, config[iZone], val_DirectIter); } } @@ -354,9 +338,7 @@ void CDiscAdjFEAIteration::Postprocess(COutput* output, CIntegration**** integra myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iVar) << "\t"; } - myfile_res << endl; - - myfile_res.close(); + myfile_res << "\n"; } // TEST: for implementation of python framework in standalone structural problems @@ -385,23 +367,14 @@ void CDiscAdjFEAIteration::Postprocess(COutput* output, CIntegration**** integra } if (outputDVFEA) { - unsigned short iDV; - unsigned short nDV = solvers0[ADJFEA_SOL]->GetnDVFEA(); - - myfile_res << "INDEX" - << "\t" - << "GRAD" << endl; - + const auto nDV = solvers0[ADJFEA_SOL]->GetnDVFEA(); + myfile_res << "INDEX\tGRAD\n"; myfile_res.precision(15); - for (iDV = 0; iDV < nDV; iDV++) { - myfile_res << iDV; - myfile_res << "\t"; - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iDV); - myfile_res << endl; + for (auto iDV = 0u; iDV < nDV; iDV++) { + myfile_res << iDV << "\t"; + myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iDV) << "\n"; } - - myfile_res.close(); } } } diff --git a/SU2_CFD/src/iteration/CFEAIteration.cpp b/SU2_CFD/src/iteration/CFEAIteration.cpp index 973f3600c6d..91b9e715681 100644 --- a/SU2_CFD/src/iteration/CFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CFEAIteration.cpp @@ -69,8 +69,8 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom } else if (nonlinear && !incremental_load) { /*--- THIS IS THE DIRECT APPROACH (NO INCREMENTAL LOAD APPLIED) ---*/ - /*--- Keep the current inner iter, we need to restore it in discrete adjoint cases as file output depends on it - * ---*/ + /*--- Keep the current inner iter, we need to restore it in discrete adjoint cases + * because file output depends on it. ---*/ const auto CurIter = config[val_iZone]->GetInnerIter(); /*--- Newton-Raphson subiterations ---*/ @@ -84,11 +84,11 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom if (disc_adj_fem) { config[val_iZone]->SetInnerIter(CurIter); break; - } StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, INST_0); + } + StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, + FFDBox, val_iZone, INST_0); - if (StopCalc && (IntIter > 0)) break; - + if (StopCalc && (IntIter > 0)) break; } } else { @@ -101,13 +101,16 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom /*--- Run two nonlinear iterations to check if incremental loading can be skipped ---*/ - for (IntIter = 0; IntIter < 2; ++IntIter) { + auto Iterate = [&](unsigned long IntIter) { config[val_iZone]->SetInnerIter(IntIter); - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst); StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, INST_0); + }; + + for (IntIter = 0; IntIter < 2; ++IntIter) { + Iterate(IntIter); } /*--- Early return if we already meet the convergence criteria. ---*/ @@ -125,13 +128,7 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom /*--- Newton-Raphson subiterations ---*/ for (IntIter = 2; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) { - config[val_iZone]->SetInnerIter(IntIter); - - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst); - - StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, INST_0); - + Iterate(IntIter); if (StopCalc) break; } @@ -162,14 +159,7 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom /*--- Newton-Raphson subiterations ---*/ for (IntIter = 0; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) { - config[val_iZone]->SetInnerIter(IntIter); - - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, - val_iInst); - - StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, INST_0); - + Iterate(IntIter); if (StopCalc && (IntIter > 0)) break; } } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 74870622ea0..07018b95233 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -615,8 +615,8 @@ void CFEASolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_con SU2_OMP_PARALLEL { + su2double zeros[MAXNVAR] = {0.0}; if (!config->GetPrestretch()) { - su2double zeros[MAXNVAR] = {0.0}; SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) nodes->SetSolution(iPoint, zeros); @@ -628,6 +628,14 @@ void CFEASolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_con nodes->SetSolution(iPoint, nodes->GetPrestretch(iPoint)); END_SU2_OMP_FOR } + if (config->GetTime_Domain()) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + nodes->SetSolution_Vel(iPoint, zeros); + nodes->SetSolution_Accel(iPoint, zeros); + } + END_SU2_OMP_FOR + } } END_SU2_OMP_PARALLEL } @@ -3108,8 +3116,8 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c /*--- Detect a wrong solution file. ---*/ if (counter != nPointDomain) { - SU2_MPI::Error(string("The solution file ") + filename + string(" doesn't match with the mesh file!\n") + - string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); + SU2_MPI::Error("The solution file " + filename + " doesn't match with the mesh file!\n" + "It could be empty lines at the end of the file.", CURRENT_FUNCTION); } /*--- MPI. If dynamic, we also need to communicate the old solution. ---*/ @@ -3117,7 +3125,9 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c InitiateComms(geometry[MESH_0], config, SOLUTION_FEA); CompleteComms(geometry[MESH_0], config, SOLUTION_FEA); - if (dynamic) nodes->Set_Solution_time_n(); + /*--- It's important to not push back the solution when this function is used to load solutions for + * unsteady discrete adjoints, otherwise we overwrite one of the two solutions needed. ---*/ + if (dynamic && val_update_geo) nodes->Set_Solution_time_n(); if (fluid_structure) { for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { From 8d2b09697a3b39ae33a5ae04ae566be5fe33d275 Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Sat, 15 Apr 2023 21:00:27 -0700 Subject: [PATCH 10/20] Update SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp --- SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 0a2a77ae31f..6277f409b8d 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -780,8 +780,6 @@ void CDiscAdjMultizoneDriver::SetAdjObjFunction() { } if (rank == MASTER_NODE) { AD::SetDerivative(ObjFunc_Index, SU2_TYPE::GetValue(seeding)); - } else { - AD::SetDerivative(ObjFunc_Index, 0.0); } } From 355d5e00d3a60e22e6efde8026469c49487d6506 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 16 Apr 2023 18:14:31 -0700 Subject: [PATCH 11/20] add unsteady FEA ad testcase --- .../py_wrapper/custom_load_fea/run_ad.py | 221 ++++++++++++++++++ TestCases/serial_regression_AD.py | 13 ++ 2 files changed, 234 insertions(+) create mode 100644 TestCases/py_wrapper/custom_load_fea/run_ad.py diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py new file mode 100644 index 00000000000..c39e05234d8 --- /dev/null +++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py @@ -0,0 +1,221 @@ +#!/usr/bin/env python + +## \file run.py +# \brief Unsteady adjoint FEA case with custom load. +# \version 7.5.1 "Blackbird" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2023, 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 . + +import pysu2ad +from mpi4py import MPI + +common_settings = """ +SOLVER= ELASTICITY +GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS +FORMULATION_ELASTICITY_2D= PLANE_STRESS + +TIME_DOMAIN= YES +TIME_STEP=0.01 +TIME_DISCRE_FEA= NEWMARK_IMPLICIT +NEWMARK_GAMMA= 0.5 +NEWMARK_BETA= 0.25 + +MATERIAL_MODEL= NEO_HOOKEAN +ELASTICITY_MODULUS= 10000 +POISSON_RATIO= 0.3 +MATERIAL_DENSITY= __DENSITY__ + +MARKER_CLAMPED= ( x_minus ) +MARKER_FLUID_LOAD= ( x_plus, y_minus, y_plus ) + +LINEAR_SOLVER= CONJUGATE_GRADIENT +DISCADJ_LIN_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +DISCADJ_LIN_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-5 +LINEAR_SOLVER_ITER= 100 + +MESH_FORMAT= RECTANGLE +MESH_BOX_SIZE= ( 17, 5, 0 ) +MESH_BOX_LENGTH= ( 0.5, 0.05, 0 ) + +OUTPUT_FILES= RESTART, PARAVIEW +OUTPUT_WRT_FREQ= 1 +OBJECTIVE_FUNCTION= STRESS_PENALTY +STRESS_PENALTY_PARAM= ( 500, 20 ) + +INNER_ITER= 20 +CONV_RESIDUAL_MINVAL= -4 +CONV_STARTITER= 5 +MAX_TIME= 0.2 +TIME_ITER= 21 +""" + +primal_settings = """ +MATH_PROBLEM= DIRECT +SCREEN_OUTPUT= TIME_ITER, CUR_TIME, INNER_ITER, RMS_RES, LINSOL, VMS, STRESS_PENALTY, TAVG_STRESS_PENALTY +HISTORY_OUTPUT= ITER, RMS_RES, STRUCT_COEFF, TAVG_STRUCT_COEFF +CONV_FIELD= REL_RMS_RTOL +""" + +adjoint_settings = """ +MATH_PROBLEM= DISCRETE_ADJOINT +SCREEN_OUTPUT= TIME_ITER, CUR_TIME, INNER_ITER, ADJOINT_DISP_X, ADJOINT_DISP_Y, LINSOL, SENSITIVITY +CONV_FIELD= ADJOINT_DISP_X, ADJOINT_DISP_Y +FEA_ADVANCED_MODE= YES +UNST_ADJOINT_ITER= 21 +ITER_AVERAGE_OBJ= 0 +SOLUTION_FILENAME= restart.dat +""" + + +def ApplyLoad(driver, marker_id, peak_load): + """ + Apply a load based on the coordinates and return the derivatives + of the nodal forces with respect to the peak load. + """ + derivatives = [] + if marker_id < 0: return + + marker_coords = driver.MarkerCoordinates(marker_id) + l = 0.5 + dx = l / 16 # known from mesh settings in this case. + for i_vertex in range(driver.GetNumberMarkerNodes(marker_id)): + x = marker_coords(i_vertex, 0) + nodal_force = (peak_load * x / l) * dx + # Half load due to half dx on first and last node. + if abs(x) < 1e-6 or abs(x - l) < 1e-6: + nodal_force = nodal_force / 2 + driver.SetMarkerCustomFEALoad(marker_id, i_vertex, (0, nodal_force)) + derivatives.append(nodal_force / peak_load) + + return derivatives + + +def RunPrimal(density, peak_load): + """Runs the primal solver for a given density and peak load and returns the time average objective function.""" + comm = MPI.COMM_WORLD + + with open('config_unsteady.cfg', 'w') as f: + f.write(common_settings.replace('__DENSITY__', str(density)) + primal_settings) + + # Initialize the primal driver of SU2, this includes solver preprocessing. + try: + driver = pysu2ad.CSinglezoneDriver('config_unsteady.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2ad.CSinglezoneDriver : ', exception) + raise + + # Get the ID of the marker where the load is applied. + all_marker_ids = driver.GetMarkerIndices() + marker_name = 'y_minus' + marker_id = all_marker_ids[marker_name] if marker_name in all_marker_ids else -1 + + # Apply a load based on the coordinates. + ApplyLoad(driver, marker_id, peak_load) + + # Solve. + driver.StartSolver() + + # Get the time average + tavg_stress_penalty = driver.GetOutputValue('TAVG_STRESS_PENALTY') + + # Finalize the solver and exit cleanly. + driver.Finalize() + + return tavg_stress_penalty + + +def RunAdjoint(density, peak_load): + """Runs the adjoint solver and returns the sensitivity of the objective function to the peak load.""" + comm = MPI.COMM_WORLD + + with open('config_unsteady_ad.cfg', 'w') as f: + f.write(common_settings.replace('__DENSITY__', str(density)) + adjoint_settings) + + # Initialize the adjoint driver of SU2, this includes solver preprocessing. + try: + driver = pysu2ad.CDiscAdjSinglezoneDriver('config_unsteady_ad.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2ad.CDiscAdjSinglezoneDriver : ', exception) + raise + + # Get the ID of the marker where the load is applied. + all_marker_ids = driver.GetMarkerIndices() + marker_name = 'y_minus' + marker_id = all_marker_ids[marker_name] if marker_name in all_marker_ids else -1 + + # Apply the same load that was used in the primal problem. + derivatives = ApplyLoad(driver, marker_id, peak_load) + + n_vertex = driver.GetNumberMarkerNodes(marker_id) if marker_id >= 0 else 0 + load_sens = 0.0 + + # Run the time loop in python to extract sensitivities at each step. + for time_iter in range(driver.GetNumberTimeIter()): + # Preprocess adjoint iteration (AD recording). + driver.Preprocess(time_iter) + + # Run one time iteration. + driver.Run() + driver.Postprocess() + driver.Update() + + # Accumulate sensitivies. + for i_vertex in range(n_vertex): + load_sens += derivatives[i_vertex] * driver.GetMarkerFEALoadSensitivity(marker_id, i_vertex)[1] + + # Monitor the solver and output solution to file if required. + driver.Monitor(time_iter) + driver.Output(time_iter) + + # Finalize the solver and exit cleanly. + driver.Finalize() + + return load_sens + + +def main(): + # Run the primal with 2 loads to compute the sensitivity via finite differences. + obj_pert_load = RunPrimal(1, 2.002) + obj_pert_rho = RunPrimal(1.0001, 2) + obj = RunPrimal(1, 2) + sens_load_fd = (obj_pert_load - obj) / 0.002 + sens_rho_fd = (obj_pert_rho - obj) / 0.0001 + + sens_load = RunAdjoint(1, 2) + # Read the density sensitivity from the special output file. + with open('Results_Reverse_Adjoint.txt') as f: + sens_rho = float(f.readlines()[-1].strip().split('\t')[-1]) + + print(" Finite Differences\tDiscrete Adjoint") + print(f"Load {sens_load_fd}\t{sens_load}") + print(f"Rho {sens_rho_fd}\t{sens_rho}") + + assert abs(sens_load / sens_load_fd - 1) < 1e-4, "Error in load derivative." + assert abs(sens_rho / sens_rho_fd - 1) < 1e-3, "Error in density derivative." + + # Print results for the regression script to check. + print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho) + + +if __name__ == '__main__': + main() diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index d9885b1dcd1..dce2b58b25a 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -343,6 +343,19 @@ def main(): test_list.append(pywrapper_FEA_AD_FlowLoad) pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test()) + # FEA unsteady AD Load Sensitivity + pywrapper_Unst_FEA_AD = TestCase('pywrapper_Unst_FEA_AD') + pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" + pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" + pywrapper_Unst_FEA_AD.test_iter = 100 + pywrapper_Unst_FEA_AD.test_vals = [0.2569256889872473, 0.256926360722213, 0.03166391031539373, 0.0316910814265723] + pywrapper_Unst_FEA_AD.command = TestCase.Command(exec = "python", param = "run_ad.py") + pywrapper_Unst_FEA_AD.timeout = 1600 + pywrapper_Unst_FEA_AD.tol = 0.00001 + pywrapper_Unst_FEA_AD.new_output = False + test_list.append(pywrapper_Unst_FEA_AD) + pass_list.append(pywrapper_Unst_FEA_AD.run_test()) + # Flow AD Mesh Displacement Sensitivity pywrapper_CFD_AD_MeshDisp = TestCase('pywrapper_CFD_AD_MeshDisp') pywrapper_CFD_AD_MeshDisp.cfg_dir = "py_wrapper/disc_adj_flow/mesh_disp_sens" From ee9ef29328909c55d589e18123c203b01fcbc274 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 16 Apr 2023 18:21:41 -0700 Subject: [PATCH 12/20] fix warning --- TestCases/py_wrapper/custom_load_fea/run_ad.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py index c39e05234d8..9e1f8e4a7b6 100644 --- a/TestCases/py_wrapper/custom_load_fea/run_ad.py +++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py @@ -93,7 +93,7 @@ def ApplyLoad(driver, marker_id, peak_load): of the nodal forces with respect to the peak load. """ derivatives = [] - if marker_id < 0: return + if marker_id < 0: return derivatives marker_coords = driver.MarkerCoordinates(marker_id) l = 0.5 From b06ebdb2cc369b57af6ba363fb1bb9439af830fd Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 16 Apr 2023 18:30:17 -0700 Subject: [PATCH 13/20] remove CD explicit as it is not implemented --- Common/include/option_structure.hpp | 2 -- Common/src/CConfig.cpp | 3 --- SU2_CFD/src/integration/CStructuralIntegration.cpp | 4 ---- SU2_CFD/src/solvers/CFEASolver.cpp | 6 ------ 4 files changed, 15 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index f08925fd40b..e661b451d1b 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1434,12 +1434,10 @@ static const MapType Heat_TimeStep_Map = { * \brief Type of time integration schemes */ enum class STRUCT_TIME_INT { - CD_EXPLICIT, /*!< \brief Support for implementing an explicit method. */ NEWMARK_IMPLICIT, /*!< \brief Implicit Newmark integration definition. */ GENERALIZED_ALPHA, /*!< \brief Support for implementing another implicit method. */ }; static const MapType Time_Int_Map_FEA = { - MakePair("CD_EXPLICIT", STRUCT_TIME_INT::CD_EXPLICIT) MakePair("NEWMARK_IMPLICIT", STRUCT_TIME_INT::NEWMARK_IMPLICIT) MakePair("GENERALIZED_ALPHA", STRUCT_TIME_INT::GENERALIZED_ALPHA) }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index e2d1e37cee0..4458e846616 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6870,9 +6870,6 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { if (fea) { switch (Kind_TimeIntScheme_FEA) { - case STRUCT_TIME_INT::CD_EXPLICIT: - cout << "Explicit time integration (NOT IMPLEMENTED YET)." << endl; - break; case STRUCT_TIME_INT::GENERALIZED_ALPHA: cout << "Generalized-alpha method." << endl; break; diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp index 276cf4a01bf..5fc18babb3f 100644 --- a/SU2_CFD/src/integration/CStructuralIntegration.cpp +++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp @@ -134,7 +134,6 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * /*--- Set the Jacobian according to the different time integration methods ---*/ switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, numerics, config); break; @@ -166,7 +165,6 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * /*--- Update solution ---*/ switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): solver_container[MainSolver]->ImplicitNewmark_Update(geometry, config); break; @@ -194,8 +192,6 @@ void CStructuralIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolv /*--- Update the solution according to the integration scheme used ---*/ switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): if (fsi && config->GetRelaxation()) solver->ImplicitNewmark_Relaxation(geometry, config); break; diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 07018b95233..e4877179706 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1340,9 +1340,6 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, else if (dynamic) { switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - cout << "NOT IMPLEMENTED YET" << endl; - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): /*--- Loop over all points, and set aux vector TimeRes_Aux = a0*U+a2*U'+a3*U'' ---*/ @@ -1495,9 +1492,6 @@ void CFEASolver::Compute_IntegrationConstants(const CConfig *config) { su2double gamma = config->GetNewmark_gamma(), beta = config->GetNewmark_beta(); switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - cout << "NOT IMPLEMENTED YET" << endl; - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): /*--- Integration constants for Newmark scheme ---*/ From bd95febb7e369379834efd7a87e175e121432b09 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 17 Apr 2023 11:58:07 -0700 Subject: [PATCH 14/20] parallel test instead of serial --- TestCases/parallel_regression_AD.py | 13 +++++++++++++ TestCases/py_wrapper/custom_load_fea/run_ad.py | 16 +++++++++++----- TestCases/serial_regression_AD.py | 13 ------------- 3 files changed, 24 insertions(+), 18 deletions(-) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index a67de6139db..1ea499b38b6 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -408,6 +408,19 @@ def main(): test_list.append(pywrapper_FEA_AD_FlowLoad) pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test()) + # FEA unsteady AD Load Sensitivity + pywrapper_Unst_FEA_AD = TestCase('pywrapper_Unst_FEA_AD') + pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" + pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" + pywrapper_Unst_FEA_AD.test_iter = 100 + pywrapper_Unst_FEA_AD.test_vals = [0.2569257, 0.2569264, 0.03166391, 0.03169108] + pywrapper_Unst_FEA_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py") + pywrapper_Unst_FEA_AD.timeout = 1600 + pywrapper_Unst_FEA_AD.tol = 0.00001 + pywrapper_Unst_FEA_AD.new_output = False + test_list.append(pywrapper_Unst_FEA_AD) + pass_list.append(pywrapper_Unst_FEA_AD.run_test()) + # Flow AD Mesh Displacement Sensitivity pywrapper_CFD_AD_MeshDisp = TestCase('pywrapper_CFD_AD_MeshDisp') pywrapper_CFD_AD_MeshDisp.cfg_dir = "py_wrapper/disc_adj_flow/mesh_disp_sens" diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py index 9e1f8e4a7b6..92e454083b3 100644 --- a/TestCases/py_wrapper/custom_load_fea/run_ad.py +++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py @@ -190,10 +190,13 @@ def RunAdjoint(density, peak_load): # Finalize the solver and exit cleanly. driver.Finalize() - return load_sens + return comm.allreduce(load_sens) def main(): + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + # Run the primal with 2 loads to compute the sensitivity via finite differences. obj_pert_load = RunPrimal(1, 2.002) obj_pert_rho = RunPrimal(1.0001, 2) @@ -206,15 +209,18 @@ def main(): with open('Results_Reverse_Adjoint.txt') as f: sens_rho = float(f.readlines()[-1].strip().split('\t')[-1]) - print(" Finite Differences\tDiscrete Adjoint") - print(f"Load {sens_load_fd}\t{sens_load}") - print(f"Rho {sens_rho_fd}\t{sens_rho}") + if rank == 0: + print(" Finite Differences\tDiscrete Adjoint") + print(f"Load {sens_load_fd}\t{sens_load}") + print(f"Rho {sens_rho_fd}\t{sens_rho}") assert abs(sens_load / sens_load_fd - 1) < 1e-4, "Error in load derivative." assert abs(sens_rho / sens_rho_fd - 1) < 1e-3, "Error in density derivative." # Print results for the regression script to check. - print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho) + if rank == 0: + print("\n") + print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho) if __name__ == '__main__': diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index dce2b58b25a..d9885b1dcd1 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -343,19 +343,6 @@ def main(): test_list.append(pywrapper_FEA_AD_FlowLoad) pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test()) - # FEA unsteady AD Load Sensitivity - pywrapper_Unst_FEA_AD = TestCase('pywrapper_Unst_FEA_AD') - pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" - pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" - pywrapper_Unst_FEA_AD.test_iter = 100 - pywrapper_Unst_FEA_AD.test_vals = [0.2569256889872473, 0.256926360722213, 0.03166391031539373, 0.0316910814265723] - pywrapper_Unst_FEA_AD.command = TestCase.Command(exec = "python", param = "run_ad.py") - pywrapper_Unst_FEA_AD.timeout = 1600 - pywrapper_Unst_FEA_AD.tol = 0.00001 - pywrapper_Unst_FEA_AD.new_output = False - test_list.append(pywrapper_Unst_FEA_AD) - pass_list.append(pywrapper_Unst_FEA_AD.run_test()) - # Flow AD Mesh Displacement Sensitivity pywrapper_CFD_AD_MeshDisp = TestCase('pywrapper_CFD_AD_MeshDisp') pywrapper_CFD_AD_MeshDisp.cfg_dir = "py_wrapper/disc_adj_flow/mesh_disp_sens" From 9cd53506214563d07731839081c53bad55f9165c Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 17 Apr 2023 12:37:35 -0700 Subject: [PATCH 15/20] update values --- TestCases/parallel_regression_AD.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 1ea499b38b6..772b7012db7 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -413,7 +413,7 @@ def main(): pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" pywrapper_Unst_FEA_AD.test_iter = 100 - pywrapper_Unst_FEA_AD.test_vals = [0.2569257, 0.2569264, 0.03166391, 0.03169108] + pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.031988, 0.032015] pywrapper_Unst_FEA_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py") pywrapper_Unst_FEA_AD.timeout = 1600 pywrapper_Unst_FEA_AD.tol = 0.00001 From 9b28e74216891cf914d5a5bb6a6a9f2efcdf8e9d Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 19 Apr 2023 18:55:30 -0700 Subject: [PATCH 16/20] cleanup the time integration of fea material sensitivities --- Common/include/CConfig.hpp | 30 +---- Common/src/CConfig.cpp | 5 + SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp | 52 +++------ SU2_CFD/include/solvers/CSolver.hpp | 36 ------ .../src/iteration/CDiscAdjFEAIteration.cpp | 109 +++++------------- .../numerics/elasticity/CFEAElasticity.cpp | 8 +- SU2_CFD/src/output/CAdjElasticityOutput.cpp | 75 +++++++----- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 11 +- SU2_CFD/src/solvers/CFEASolver.cpp | 4 +- SU2_PY/SU2/io/historyMap.py | 8 +- TestCases/disc_adj_fsi/Airfoil_2d/config.cfg | 2 +- TestCases/disc_adj_fsi/config.cfg | 2 +- .../py_wrapper/custom_load_fea/run_ad.py | 14 ++- 13 files changed, 125 insertions(+), 231 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 14acb60ec69..626b8c30f9a 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -2108,17 +2108,11 @@ class CConfig { su2double GetElasticyMod(unsigned short id_val) const { return ElasticityMod[id_val]; } /*! - * \brief Decide whether to apply DE effects to the model. - * \return TRUE if the DE effects are to be applied, FALSE otherwise. - */ + * \brief Decide whether to apply DE effects to the model. + * \return TRUE if the DE effects are to be applied, FALSE otherwise. + */ bool GetDE_Effects(void) const { return DE_Effects; } - /*! - * \brief Decide whether to predict the DE effects for the next time step. - * \return TRUE if the DE effects are to be applied, FALSE otherwise. - */ - bool GetDE_Predicted(void); - /*! * \brief Get the number of different electric constants. * \return Value of the DE modulus. @@ -8749,22 +8743,10 @@ class CConfig { unsigned short GetnIntCoeffs(void) const { return nIntCoeffs; } /*! - * \brief Get the number of different values for the elasticity modulus. - * \return Number of different values for the elasticity modulus. - */ - unsigned short GetnElasticityMod(void) const { return nElasticityMod; } - - /*! - * \brief Get the number of different values for the Poisson ratio. - * \return Number of different values for the Poisson ratio. - */ - unsigned short GetnPoissonRatio(void) const { return nPoissonRatio; } - - /*! - * \brief Get the number of different values for the Material density. - * \return Number of different values for the Material density. + * \brief Get the number of different materials for the elasticity solver. + * \return Number of different materials. */ - unsigned short GetnMaterialDensity(void) const { return nMaterialDensity; } + unsigned short GetnElasticityMat(void) const { return nElasticityMod; } /*! * \brief Get the integration coefficients for the Generalized Alpha - Newmark integration integration scheme. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 4458e846616..e705b7f38d4 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4648,6 +4648,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i MaterialDensity = new su2double[1]; MaterialDensity[0] = 7854; } + if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity) { + SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, and MATERIAL_DENSITY need to have the same number " + "of entries (the number of materials).", CURRENT_FUNCTION); + } + if (nElectric_Constant == 0) { nElectric_Constant = 1; Electric_Constant = new su2double[1]; Electric_Constant[0] = 0.0; diff --git a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp index abd8fcddfb0..446306efdc3 100644 --- a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp @@ -50,7 +50,8 @@ class CDiscAdjFEASolver final : public CSolver { su2double* val = nullptr; /*!< \brief Value of the variable. */ su2double* LocalSens = nullptr; /*!< \brief Local sensitivity (domain). */ su2double* GlobalSens = nullptr; /*!< \brief Global sensitivity (mpi). */ - su2double* TotalSens = nullptr; /*!< \brief Total sensitivity (time domain). */ + su2double* OldSens = nullptr; /*!< \brief Previous global sensitivity, used to update the total. */ + su2double* TotalSens = nullptr; /*!< \brief Total sensitivity (integrated over time). */ su2double& operator[] (unsigned short i) { return val[i]; } const su2double& operator[] (unsigned short i) const { return val[i]; } @@ -61,6 +62,7 @@ class CDiscAdjFEASolver final : public CSolver { val = new su2double[n](); LocalSens = new su2double[n](); GlobalSens = new su2double[n](); + OldSens = new su2double[n](); TotalSens = new su2double[n](); } @@ -69,6 +71,7 @@ class CDiscAdjFEASolver final : public CSolver { delete [] val; delete [] LocalSens; delete [] GlobalSens; + delete [] OldSens; delete [] TotalSens; } @@ -80,10 +83,19 @@ class CDiscAdjFEASolver final : public CSolver { for (auto i = 0u; i < size; ++i) LocalSens[i] = SU2_TYPE::GetDerivative(val[i]); SU2_MPI::Allreduce(LocalSens, GlobalSens, size, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + + for (auto i = 0u; i < size; ++i) { + /*--- Update the total by subtracting the old and adding the new value. + * Then update the old value for the next call to this function. ---*/ + TotalSens[i] += GlobalSens[i] - OldSens[i]; + OldSens[i] = GlobalSens[i]; + } } void UpdateTotal() { - for (auto i = 0u; i < size; ++i) TotalSens[i] += GlobalSens[i]; + /*--- Clears the old values such that on the next time step the total is + * incremented instead of updated. ---*/ + for (auto i = 0u; i < size; ++i) OldSens[i] = 0.0; } ~SensData() { clear(); } @@ -213,42 +225,6 @@ class CDiscAdjFEASolver final : public CSolver { */ inline su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const override { return DV.TotalSens[iDVFEA]; } - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Young Modulus E - */ - inline su2double GetGlobal_Sens_E(unsigned short iVal) const override { return E.GlobalSens[iVal]; } - - /*! - * \brief A virtual member. - * \return Value of the Mach sensitivity for the Poisson's ratio Nu - */ - inline su2double GetGlobal_Sens_Nu(unsigned short iVal) const override { return Nu.GlobalSens[iVal]; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Electric Field in the region iEField - */ - inline su2double GetGlobal_Sens_EField(unsigned short iEField) const override { return EField.GlobalSens[iEField]; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the FEA DV in the region iDVFEA - */ - inline su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const override { return DV.GlobalSens[iDVFEA]; } - - /*! - * \brief Get the total sensitivity for the structural density - * \return Value of the structural density sensitivity - */ - inline su2double GetGlobal_Sens_Rho(unsigned short iVal) const override { return Rho.GlobalSens[iVal]; } - - /*! - * \brief Get the total sensitivity for the structural weight - * \return Value of the structural weight sensitivity - */ - inline su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const override { return Rho_DL.GlobalSens[iVal]; } - /*! * \brief Get the value of the Young modulus from the adjoint solver * \return Value of the Young modulus from the adjoint solver diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 892552859b0..ba28bf8b193 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3201,42 +3201,6 @@ class CSolver { */ inline virtual su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const { return 0.0; } - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Young Modulus E - */ - inline virtual su2double GetGlobal_Sens_E(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Poisson's ratio Nu - */ - inline virtual su2double GetGlobal_Sens_Nu(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the structural density sensitivity - */ - inline virtual su2double GetGlobal_Sens_Rho(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the structural weight sensitivity - */ - inline virtual su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Electric Field in the region iEField - */ - inline virtual su2double GetGlobal_Sens_EField(unsigned short iEField) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the FEA DV in the region iDVFEA - */ - inline virtual su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const { return 0.0; } - /*! * \brief A virtual member. * \return Value of the Young modulus from the adjoint solver diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index 7decbd05b57..a37971f0aeb 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -32,32 +32,6 @@ CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(config), CurrentRecording(NONE) { fem_iteration = new CFEAIteration(config); - - // TEMPORARY output only for standalone structural problems - if (config->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) { - bool de_effects = config->GetDE_Effects(); - unsigned short iVar; - - /*--- Header of the temporary output file ---*/ - ofstream myfile_res; - myfile_res.open("Results_Reverse_Adjoint.txt"); - - myfile_res << "Obj_Func" - << " "; - for (iVar = 0; iVar < config->GetnElasticityMod(); iVar++) myfile_res << "Sens_E_" << iVar << "\t"; - - for (iVar = 0; iVar < config->GetnPoissonRatio(); iVar++) myfile_res << "Sens_Nu_" << iVar << "\t"; - - if (config->GetTime_Domain()) { - for (iVar = 0; iVar < config->GetnMaterialDensity(); iVar++) myfile_res << "Sens_Rho_" << iVar << "\t"; - } - - if (de_effects) { - for (iVar = 0; iVar < config->GetnElectric_Field(); iVar++) myfile_res << "Sens_EField_" << iVar << "\t"; - } - - myfile_res << "\n"; - } } CDiscAdjFEAIteration::~CDiscAdjFEAIteration() = default; @@ -71,31 +45,49 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat auto dirNodes = solvers0[FEA_SOL]->GetNodes(); auto adjNodes = solvers0[ADJFEA_SOL]->GetNodes(); - /*--- For the dynamic adjoint, load direct solutions from restart files. ---*/ + auto StoreDirectSolution = [&]() { + for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) { + adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); + } + }; + + /*--- For the dynamic adjoint, load direct solutions from restart files. + * For steady, store the direct solution to be able to reset it later. ---*/ if (config[iZone]->GetTime_Domain()) { const int TimeIter = config[iZone]->GetTimeIter(); const int Direct_Iter = SU2_TYPE::Int(config[iZone]->GetUnst_AdjointIter()) - TimeIter - 1; - /*--- We want to load the already converged solution at timesteps n and n-1 ---*/ - - /*--- Load solution at timestep n-1 ---*/ + /*--- We need the already converged solution at timesteps n and n-1. On the first + * adjoint time step we load both, then n-1 becomes "n" and we only load n-2. ---*/ - LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1); + if (TimeIter > 0) { + /*--- Save n-1 to become n. ---*/ + for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) { + adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution_time_n(iPoint)); + } + } - /*--- Push solution back to correct array ---*/ + /*--- Load solution at timestep n-1 and push back to correct array. ---*/ + LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1); dirNodes->Set_Solution_time_n(); - /*--- Load solution timestep n ---*/ + /*--- Load or set solution at timestep n. ---*/ - LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter); - } - - /*--- Store FEA solution also in the adjoint solver in order to be able to reset it later. ---*/ + if (TimeIter == 0) { + LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter); + StoreDirectSolution(); + } else { + /*--- Set n-1 as n. ---*/ + for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) + for (auto iVar = 0u; iVar < solvers0[ADJFEA_SOL]->GetnVar(); iVar++) + dirNodes->SetSolution(iPoint, iVar, adjNodes->GetSolution_Direct(iPoint)[iVar]); + } - for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) { - adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); + } else { + /*--- Steady. ---*/ + StoreDirectSolution(); } solvers0[ADJFEA_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, 0, RUNTIME_ADJFEA_SYS, false); @@ -178,7 +170,7 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge const int mat_knowles = MAT_KNOWLES+offset; const int de_term = DE_TERM+offset; - for (unsigned short iProp = 0; iProp < config[iZone]->GetnElasticityMod(); iProp++) { + for (unsigned short iProp = 0; iProp < config[iZone]->GetnElasticityMat(); iProp++) { su2double E = adj_solver->GetVal_Young(iProp); su2double nu = adj_solver->GetVal_Poisson(iProp); su2double rho = adj_solver->GetVal_Rho(iProp); @@ -302,45 +294,8 @@ void CDiscAdjFEAIteration::Postprocess(COutput* output, CIntegration**** integra CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short iZone, unsigned short iInst) { - const bool dynamic = (config[iZone]->GetTime_Domain()); auto solvers0 = solver[iZone][iInst][MESH_0]; - // TEMPORARY output only for standalone structural problems - if (config[iZone]->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) { - unsigned short iVar; - - const bool de_effects = config[iZone]->GetDE_Effects(); - - /*--- Header of the temporary output file ---*/ - ofstream myfile_res; - myfile_res.open("Results_Reverse_Adjoint.txt", ios::app); - - myfile_res.precision(15); - - myfile_res << config[iZone]->GetTimeIter() << "\t"; - - solvers0[FEA_SOL]->Evaluate_ObjFunc(config[iZone], solvers0); - myfile_res << scientific << solvers0[FEA_SOL]->GetTotal_ComboObj() << "\t"; - - for (iVar = 0; iVar < config[iZone]->GetnElasticityMod(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_E(iVar) << "\t"; - for (iVar = 0; iVar < config[iZone]->GetnPoissonRatio(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_Nu(iVar) << "\t"; - if (dynamic) { - for (iVar = 0; iVar < config[iZone]->GetnMaterialDensity(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_Rho(iVar) << "\t"; - } - if (de_effects) { - for (iVar = 0; iVar < config[iZone]->GetnElectric_Field(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_EField(iVar) << "\t"; - } - for (iVar = 0; iVar < solvers0[ADJFEA_SOL]->GetnDVFEA(); iVar++) { - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iVar) << "\t"; - } - - myfile_res << "\n"; - } - // TEST: for implementation of python framework in standalone structural problems if (config[iZone]->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) { /*--- Header of the temporary output file ---*/ diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index 6fd3a35042c..fc87cbc19a5 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -38,23 +38,19 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, bool body_forces = config->GetDeadLoad(); // Body forces (dead loads). bool pseudo_static = config->GetPseudoStatic(); - unsigned short iVar, nProp; + unsigned short iVar; /*--- Initialize vector structures for multiple material definition ---*/ - nProp = config->GetnElasticityMod(); + const auto nProp = config->GetnElasticityMat(); E_i = new su2double[nProp]; for (iVar = 0; iVar < nProp; iVar++) E_i[iVar] = config->GetElasticyMod(iVar); - nProp = config->GetnPoissonRatio(); - Nu_i = new su2double[nProp]; for (iVar = 0; iVar < nProp; iVar++) Nu_i[iVar] = config->GetPoissonRatio(iVar); - nProp = config->GetnMaterialDensity(); - Rho_s_i = new su2double[nProp]; // For inertial effects Rho_s_DL_i = new su2double[nProp]; // For dead loads diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp index 8ed7017aba2..f55e441c2d6 100644 --- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp +++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp @@ -27,6 +27,7 @@ #include "../../include/output/CAdjElasticityOutput.hpp" +#include #include "../../../Common/include/geometry/CGeometry.hpp" #include "../../include/solvers/CSolver.hpp" @@ -50,8 +51,8 @@ CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim) requestedScreenFields.emplace_back("INNER_ITER"); requestedScreenFields.emplace_back("ADJOINT_DISP_X"); requestedScreenFields.emplace_back("ADJOINT_DISP_Y"); - requestedScreenFields.emplace_back("SENS_E"); - requestedScreenFields.emplace_back("SENS_NU"); + requestedScreenFields.emplace_back("SENS_E_0"); + requestedScreenFields.emplace_back("SENS_NU_0"); nRequestedScreenFields = requestedScreenFields.size(); } @@ -62,9 +63,9 @@ CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim) nRequestedVolumeFields = requestedVolumeFields.size(); } - if (find(requestedVolumeFields.begin(), requestedVolumeFields.end(), string("SENSITIVITY")) == requestedVolumeFields.end()) { + if (find(requestedVolumeFields.begin(), requestedVolumeFields.end(), "SENSITIVITY") == requestedVolumeFields.end()) { requestedVolumeFields.emplace_back("SENSITIVITY"); - nRequestedVolumeFields ++; + nRequestedVolumeFields++; } stringstream ss; @@ -97,47 +98,69 @@ CAdjElasticityOutput::~CAdjElasticityOutput() = default; void CAdjElasticityOutput::SetHistoryOutputFields(CConfig *config){ - // Residuals + /*--- Residuals ---*/ AddHistoryOutput("ADJOINT_DISP_X", "rms[Ux_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the X displacements.", HistoryFieldType::RESIDUAL); AddHistoryOutput("ADJOINT_DISP_Y", "rms[Uy_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the Y displacements.", HistoryFieldType::RESIDUAL); - AddHistoryOutput("ADJOINT_DISP_Z", "rms[Uz_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the Z displacements.", HistoryFieldType::RESIDUAL); + if (nVar_FEM == 3) { + AddHistoryOutput("ADJOINT_DISP_Z", "rms[Uz_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the Z displacements.", HistoryFieldType::RESIDUAL); + } - //Sensitivities - AddHistoryOutput("SENS_E", "Sens[E]", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Elasticity modulus"); - AddHistoryOutput("SENS_NU", "Sens[Nu]", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Poisson ratio"); + /*--- Sensitivities ---*/ + for (auto iVar = 0u; iVar < config->GetnElasticityMat(); iVar++) { + const auto iVarS = std::to_string(iVar); + AddHistoryOutput("SENS_E_" + iVarS, "Sens[E" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Elasticity modulus"); + AddHistoryOutput("SENS_NU_" + iVarS, "Sens[Nu" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Poisson ratio"); + if (config->GetTime_Domain() && !config->GetPseudoStatic()) { + AddHistoryOutput("SENS_RHO_" + iVarS, "Sens[Rho" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Material density"); + } + if (config->GetDeadLoad()) { + AddHistoryOutput("SENS_RHO_DL_" + iVarS, "Sens[RhoDL" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Dead load density"); + } + } + if (config->GetDE_Effects()) { + for (auto iVar = 0u; iVar < config->GetnElectric_Field(); iVar++) { + const auto iVarS = std::to_string(iVar); + AddHistoryOutput("SENS_EFIELD_" + iVarS, "Sens[EField" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Electric field"); + } + } AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); - AddHistoryOutput("BGS_ADJ_DISP_X", "bgs[A_Ux]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint X displacement.", HistoryFieldType::RESIDUAL); - AddHistoryOutput("BGS_ADJ_DISP_Y", "bgs[A_Uy]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Y displacement.", HistoryFieldType::RESIDUAL); - AddHistoryOutput("BGS_ADJ_DISP_Z", "bgs[A_Uz]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Z displacement.", HistoryFieldType::RESIDUAL); - + if (multiZone) { + AddHistoryOutput("BGS_ADJ_DISP_X", "bgs[A_Ux]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint X displacement.", HistoryFieldType::RESIDUAL); + AddHistoryOutput("BGS_ADJ_DISP_Y", "bgs[A_Uy]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Y displacement.", HistoryFieldType::RESIDUAL); + if (nVar_FEM == 3) { + AddHistoryOutput("BGS_ADJ_DISP_Z", "bgs[A_Uz]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Z displacement.", HistoryFieldType::RESIDUAL); + } + } } inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { SetHistoryOutputValue("ADJOINT_DISP_X", log10(solver[ADJFEA_SOL]->GetRes_RMS(0))); SetHistoryOutputValue("ADJOINT_DISP_Y", log10(solver[ADJFEA_SOL]->GetRes_RMS(1))); - if (nVar_FEM == 3){ + if (nVar_FEM == 3) { SetHistoryOutputValue("ADJOINT_DISP_Z", log10(solver[ADJFEA_SOL]->GetRes_RMS(2))); } - su2double Total_SensE = 0.0; su2double Total_SensNu = 0.0; - if (config->GetnElasticityMod() == 1) { - Total_SensE = solver[ADJFEA_SOL]->GetGlobal_Sens_E(0); - Total_SensNu = solver[ADJFEA_SOL]->GetGlobal_Sens_Nu(0); + for (unsigned short iVar = 0; iVar < config->GetnElasticityMat(); iVar++) { + const auto iVarS = std::to_string(iVar); + SetHistoryOutputValue("SENS_E_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_E(iVar)); + SetHistoryOutputValue("SENS_NU_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Nu(iVar)); + if (config->GetTime_Domain() && !config->GetPseudoStatic()) { + SetHistoryOutputValue("SENS_RHO_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho(iVar)); + } + if (config->GetDeadLoad()) { + SetHistoryOutputValue("SENS_RHO_DL_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho_DL(iVar)); + } } - else { - for (unsigned short iVar = 0; iVar < config->GetnElasticityMod(); iVar++){ - Total_SensE += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_E(iVar),2); - Total_SensNu += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_Nu(iVar),2); + if (config->GetDE_Effects()) { + for (auto iVar = 0u; iVar < config->GetnElectric_Field(); iVar++) { + const auto iVarS = std::to_string(iVar); + SetHistoryOutputValue("SENS_EFIELD_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_EField(iVar)); } - Total_SensE = sqrt(Total_SensE); - Total_SensNu = sqrt(Total_SensNu); } - SetHistoryOutputValue("SENS_E", Total_SensE); - SetHistoryOutputValue("SENS_NU", Total_SensNu); SetHistoryOutputValue("LINSOL_ITER", solver[ADJFEA_SOL]->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL", log10(solver[ADJFEA_SOL]->GetResLinSolver())); diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index f11c9f580b6..34a19c89a3e 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -84,16 +84,7 @@ CDiscAdjFEASolver::CDiscAdjFEASolver(CGeometry *geometry, CConfig *config, CSolv /*--- Initialize vector structures for multiple material definition ---*/ - nMPROP = config->GetnElasticityMod(); - - /*--- For a material to be fully defined, we need to have the same number for all three parameters ---*/ - bool checkDef = ((config->GetnElasticityMod() == config->GetnPoissonRatio()) && - (config->GetnElasticityMod() == config->GetnMaterialDensity()) && - (config->GetnMaterialDensity() == config->GetnPoissonRatio())); - - if (!checkDef){ - SU2_MPI::Error("WARNING: For a material to be fully defined, E, Nu and Rho need to have the same dimensions.", CURRENT_FUNCTION); - } + nMPROP = config->GetnElasticityMat(); E.resize(nMPROP); Nu.resize(nMPROP); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index e4877179706..858257f094c 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -383,8 +383,8 @@ void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { /*--- Detect a wrong solution file ---*/ if (iElem_Global_Local != nElement) { - SU2_MPI::Error(string("The properties file ") + filename + string(" doesn't match with the mesh file!\n") + - string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); + SU2_MPI::Error("The properties file " + filename + " doesn't match with the mesh file!\n" + "It could be empty lines at the end of the file.", CURRENT_FUNCTION); } } diff --git a/SU2_PY/SU2/io/historyMap.py b/SU2_PY/SU2/io/historyMap.py index aca169155e2..a477478359c 100644 --- a/SU2_PY/SU2/io/historyMap.py +++ b/SU2_PY/SU2/io/historyMap.py @@ -1186,10 +1186,10 @@ "HEADER": "Sens_AoA", "TYPE": "COEFFICIENT", }, - "SENS_E": { + "SENS_E_0": { "DESCRIPTION": "d Objective / d Elasticity modulus", "GROUP": "SENSITIVITY", - "HEADER": "Sens[E]", + "HEADER": "Sens[E_0]", "TYPE": "DEFAULT", }, "SENS_GEO": { @@ -1207,10 +1207,10 @@ "HEADER": "Sens_Mach", "TYPE": "COEFFICIENT", }, - "SENS_NU": { + "SENS_NU_0": { "DESCRIPTION": "d Objective / d Poisson ratio", "GROUP": "SENSITIVITY", - "HEADER": "Sens[Nu]", + "HEADER": "Sens[Nu_0]", "TYPE": "DEFAULT", }, "SENS_PRESS": { diff --git a/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg b/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg index 6b58ea8e026..a484a26c79b 100755 --- a/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg +++ b/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg @@ -13,6 +13,6 @@ MESH_FILENAME= mesh.su2 OBJECTIVE_FUNCTION= CUSTOM_OBJFUNC -SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], LINSOL_RESIDUAL[0], SENS_E[1], SENS_NU[1] +SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], LINSOL_RESIDUAL[0], SENS_E_0[1], SENS_NU_0[1] %WRT_ZONE_CONV=YES diff --git a/TestCases/disc_adj_fsi/config.cfg b/TestCases/disc_adj_fsi/config.cfg index 35cf049125d..044f409302c 100644 --- a/TestCases/disc_adj_fsi/config.cfg +++ b/TestCases/disc_adj_fsi/config.cfg @@ -9,6 +9,6 @@ MESH_FILENAME= mesh.su2 OBJECTIVE_FUNCTION = REFERENCE_GEOMETRY -SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], RMS_ADJ_DENSITY[0], SENS_E[1], SENS_NU[1] +SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], RMS_ADJ_DENSITY[0], SENS_E_0[1], SENS_NU_0[1] %WRT_ZONE_CONV= YES diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py index 92e454083b3..cc85ca00239 100644 --- a/TestCases/py_wrapper/custom_load_fea/run_ad.py +++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py @@ -145,7 +145,10 @@ def RunPrimal(density, peak_load): def RunAdjoint(density, peak_load): - """Runs the adjoint solver and returns the sensitivity of the objective function to the peak load.""" + """ + Runs the adjoint solver and returns the sensitivity of the objective function to the peak + load and to the material density. + """ comm = MPI.COMM_WORLD with open('config_unsteady_ad.cfg', 'w') as f: @@ -187,10 +190,12 @@ def RunAdjoint(density, peak_load): driver.Monitor(time_iter) driver.Output(time_iter) + rho_sens = driver.GetOutputValue('SENS_RHO_0') + # Finalize the solver and exit cleanly. driver.Finalize() - return comm.allreduce(load_sens) + return comm.allreduce(load_sens), rho_sens def main(): @@ -204,10 +209,7 @@ def main(): sens_load_fd = (obj_pert_load - obj) / 0.002 sens_rho_fd = (obj_pert_rho - obj) / 0.0001 - sens_load = RunAdjoint(1, 2) - # Read the density sensitivity from the special output file. - with open('Results_Reverse_Adjoint.txt') as f: - sens_rho = float(f.readlines()[-1].strip().split('\t')[-1]) + sens_load, sens_rho = RunAdjoint(1, 2) if rank == 0: print(" Finite Differences\tDiscrete Adjoint") From 7548083521ea76395488974afe2ffb0b2b958f68 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 19 Apr 2023 21:55:51 -0700 Subject: [PATCH 17/20] fixes and update ref file --- SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp | 2 +- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 19 +++++++++++-------- .../disc_adj_fsi/dyn_fsi/grad_dv.opt.ref | 16 ++++++++-------- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp index 446306efdc3..2755731f5a7 100644 --- a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp @@ -92,7 +92,7 @@ class CDiscAdjFEASolver final : public CSolver { } } - void UpdateTotal() { + void Store() { /*--- Clears the old values such that on the next time step the total is * incremented instead of updated. ---*/ for (auto i = 0u; i < size; ++i) OldSens[i] = 0.0; diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 34a19c89a3e..5a03df9aff7 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -346,12 +346,17 @@ void CDiscAdjFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_cont void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*){ - E.UpdateTotal(); - Nu.UpdateTotal(); - Rho.UpdateTotal(); - Rho_DL.UpdateTotal(); - if (de_effects) EField.UpdateTotal(); - if (fea_dv) DV.UpdateTotal(); + const bool time_domain = config->GetTime_Domain(); + + /*--- Store the final material sensitivities for the time step to increment them on the next time step. ---*/ + if (time_domain) { + E.Store(); + Nu.Store(); + Rho.Store(); + Rho_DL.Store(); + if (de_effects) EField.Store(); + if (fea_dv) DV.Store(); + } /*--- Extract the topology optimization density sensitivities. ---*/ @@ -359,8 +364,6 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo /*--- Extract the geometric sensitivities ---*/ - const bool time_domain = config->GetTime_Domain(); - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { auto Coord = geometry->nodes->GetCoord(iPoint); diff --git a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref index e543b661a60..b2e428bed5a 100644 --- a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref +++ b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref @@ -1,9 +1,9 @@ INDEX GRAD -0 -3.461460667601001e-03 -1 -1.841786311588662e-03 -2 -7.915536257748965e-04 -3 -2.739622082729719e-04 -4 -2.734869133461108e-04 -5 -7.881162428890194e-04 -6 -1.828978290516676e-03 -7 -3.427219398258311e-03 +0 -3.461460667601000e-03 +1 -1.841786311588663e-03 +2 -7.915536257748967e-04 +3 -2.739622082729717e-04 +4 -2.734869133461104e-04 +5 -7.881162428890206e-04 +6 -1.828978290516677e-03 +7 -3.427219398258316e-03 From 1549d45be67b7056745fda07034215fffe4d70e3 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 22 Apr 2023 17:22:11 -0700 Subject: [PATCH 18/20] checkpoint fixing AD index issues in the flow tractions --- Common/include/option_structure.hpp | 2 +- SU2_CFD/include/drivers/CDriverBase.hpp | 9 +++ .../include/variables/CDiscAdjVariable.hpp | 1 + .../include/variables/CFEABoundVariable.hpp | 7 ++- SU2_CFD/include/variables/CVariable.hpp | 8 ++- SU2_CFD/src/output/CAdjElasticityOutput.cpp | 36 ++++++++++++ SU2_CFD/src/output/CFlowOutput.cpp | 4 +- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 16 ++++-- SU2_CFD/src/variables/CFEABoundVariable.cpp | 7 +++ TestCases/parallel_regression_AD.py | 4 +- .../py_wrapper/custom_load_fea/run_ad.py | 57 +++++++++++++------ 11 files changed, 121 insertions(+), 30 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index e661b451d1b..69acd1c0c67 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1439,7 +1439,7 @@ enum class STRUCT_TIME_INT { }; static const MapType Time_Int_Map_FEA = { MakePair("NEWMARK_IMPLICIT", STRUCT_TIME_INT::NEWMARK_IMPLICIT) - MakePair("GENERALIZED_ALPHA", STRUCT_TIME_INT::GENERALIZED_ALPHA) + // MakePair("GENERALIZED_ALPHA", STRUCT_TIME_INT::GENERALIZED_ALPHA) Not fully implemented. }; /*! diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp index 4e79537a829..d720246efb7 100644 --- a/SU2_CFD/include/drivers/CDriverBase.hpp +++ b/SU2_CFD/include/drivers/CDriverBase.hpp @@ -499,6 +499,15 @@ class CDriverBase { main_geometry->GetnVertex(iMarker), "MarkerPrimitives", false); } + /*! + * \brief Get a read-only view of the geometry sensitivity of a discrete adjoint solver. + */ + inline CPyWrapperMatrixView Sensitivity(unsigned short iSolver) { + auto* solver = GetSolverAndCheckMarker(iSolver); + auto& sensitivity = const_cast(solver->GetNodes()->GetSensitivity()); + return CPyWrapperMatrixView(sensitivity, "Sensitivity", true); + } + /*! * \brief Set the temperature of a vertex on a specified marker (MARKER_PYTHON_CUSTOM). * \note This can be the input of a heat or flow solver in a CHT setting. diff --git a/SU2_CFD/include/variables/CDiscAdjVariable.hpp b/SU2_CFD/include/variables/CDiscAdjVariable.hpp index 72ca67f9c24..bc72508de78 100644 --- a/SU2_CFD/include/variables/CDiscAdjVariable.hpp +++ b/SU2_CFD/include/variables/CDiscAdjVariable.hpp @@ -72,6 +72,7 @@ class CDiscAdjVariable : public CVariable { inline su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const final { return Sensitivity(iPoint,iDim); } + inline const MatrixType& GetSensitivity() const final { return Sensitivity; } /*! * \brief Set/store the dual time contributions to the adjoint variable. diff --git a/SU2_CFD/include/variables/CFEABoundVariable.hpp b/SU2_CFD/include/variables/CFEABoundVariable.hpp index e04d40467c1..950fef78ebe 100644 --- a/SU2_CFD/include/variables/CFEABoundVariable.hpp +++ b/SU2_CFD/include/variables/CFEABoundVariable.hpp @@ -159,10 +159,15 @@ class CFEABoundVariable final : public CFEAVariable { */ void RegisterFlowTraction() override; + /*! + * \brief Resets the AD indices of the flow tractions. + */ + void ResetInputFlowTraction() override; + /*! * \brief Extract the flow traction derivatives. */ - inline su2double ExtractFlowTraction_Sensitivity(unsigned long iPoint, unsigned long iDim) const override { + inline su2double ExtractFlowTractionSensitivity(unsigned long iPoint, unsigned long iDim) const override { if (!fsi_analysis) return 0.0; if (!VertexMap.GetVertexIndex(iPoint)) return 0.0; return SU2_TYPE::GetDerivative(FlowTraction(iPoint,iDim)); diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index da68a28f051..efdbc25c9ce 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2165,7 +2165,12 @@ class CVariable { /*! * \brief A virtual member. */ - inline virtual su2double ExtractFlowTraction_Sensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + inline virtual void ResetInputFlowTraction() { } + + /*! + * \brief A virtual member. + */ + inline virtual su2double ExtractFlowTractionSensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } /*! * \brief Register the variables in the solution array as input/output variable. @@ -2224,6 +2229,7 @@ class CVariable { * \return value of the Sensitivity */ inline virtual su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + inline virtual const MatrixType& GetSensitivity() const { AssertOverride(); return Solution; } inline virtual void SetTau_Wall(unsigned long iPoint, su2double tau_wall) {} diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp index f55e441c2d6..e30f2895f65 100644 --- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp +++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp @@ -193,6 +193,23 @@ void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, SetVolumeOutputValue("SENSITIVITY-Y", iPoint, Node_Struc->GetSensitivity(iPoint, 1)); if (nDim == 3) SetVolumeOutputValue("SENSITIVITY-Z", iPoint, Node_Struc->GetSensitivity(iPoint, 2)); + + if (!config->GetTime_Domain()) return; + + SetVolumeOutputValue("SENS_DISP-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, 0)); + SetVolumeOutputValue("SENS_DISP-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 1)); + if (nDim == 3) + SetVolumeOutputValue("SENS_DISP-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2)); + + SetVolumeOutputValue("SENS_VEL-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, nDim)); + SetVolumeOutputValue("SENS_VEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, nDim + 1)); + if (nDim == 3) + SetVolumeOutputValue("SENS_VEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 5)); + + SetVolumeOutputValue("SENS_ACCEL-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim)); + SetVolumeOutputValue("SENS_ACCEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim + 1)); + if (nDim == 3) + SetVolumeOutputValue("SENS_ACCEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 8)); } void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){ @@ -223,4 +240,23 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("SENSITIVITY-Z", "Sensitivity_z", "SENSITIVITY", "geometric sensitivity in the z direction"); /// END_GROUP + if (!config->GetTime_Domain()) return; + + /*--- Sensitivities with respect to initial conditions. ---*/ + + AddVolumeOutput("SENS_DISP-X", "SensInitialDisp_x", "SENSITIVITY_T0", "sensitivity to the initial x displacement"); + AddVolumeOutput("SENS_DISP-Y", "SensInitialDisp_y", "SENSITIVITY_T0", "sensitivity to the initial y displacement"); + if (nDim == 3) + AddVolumeOutput("SENS_DISP-Z", "SensInitialDisp_z", "SENSITIVITY_T0", "sensitivity to the initial z displacement"); + + AddVolumeOutput("SENS_VEL-X", "SensInitialVel_x", "SENSITIVITY_T0", "sensitivity to the initial x velocity"); + AddVolumeOutput("SENS_VEL-Y", "SensInitialVel_y", "SENSITIVITY_T0", "sensitivity to the initial y velocity"); + if (nDim == 3) + AddVolumeOutput("SENS_VEL-Z", "SensInitialVel_z", "SENSITIVITY_T0", "sensitivity to the initial z velocity"); + + AddVolumeOutput("SENS_ACCEL-X", "SensInitialAccel_x", "SENSITIVITY_T0", "sensitivity to the initial x acceleration"); + AddVolumeOutput("SENS_ACCEL-Y", "SensInitialAccel_y", "SENSITIVITY_T0", "sensitivity to the initial y acceleration"); + if (nDim == 3) + AddVolumeOutput("SENS_ACCEL-Z", "SensInitialAccel_z", "SENSITIVITY_T0", "sensitivity to the initial z acceleration"); + } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 3fc02d126e7..bff61c408c9 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -873,8 +873,8 @@ void CFlowOutput::SetCustomOutputs(const CSolver* const* solver, const CGeometry const auto varIdx = i % CustomOutput::MAX_VARS_PER_SOLVER; if (solIdx == FLOW_SOL) { return flowNodes->GetPrimitive(iPoint, varIdx); - } return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx); - + } + return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx); } else { return *output.otherOutputs[i - CustomOutput::NOT_A_VARIABLE]; } diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 5a03df9aff7..5b38d271497 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -125,14 +125,18 @@ CDiscAdjFEASolver::~CDiscAdjFEASolver() { delete nodes; } void CDiscAdjFEASolver::SetRecording(CGeometry* geometry, CConfig *config){ unsigned long iPoint; - unsigned short iVar; + unsigned short iDim, iVar; /*--- Reset the solution to the initial (converged) solution ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++) + for (iPoint = 0; iPoint < nPoint; iPoint++) { for (iVar = 0; iVar < nVar; iVar++) direct_solver->GetNodes()->SetSolution(iPoint, iVar, nodes->GetSolution_Direct(iPoint)[iVar]); + auto Coord = geometry->nodes->GetCoord(iPoint); + for (iDim = 0; iDim < nDim; iDim++) AD::ResetInput(Coord[iDim]); + } + /*--- Reset the input for time n ---*/ if (config->GetTime_Domain()) { @@ -198,6 +202,8 @@ void CDiscAdjFEASolver::RegisterVariables(CGeometry *geometry, CConfig *config, for (iVar = 0; iVar < nDV; iVar++) AD::ResetInput(DV[iVar]); } + direct_solver->GetNodes()->ResetInputFlowTraction(); + if (!reset) { E.Register(); Nu.Register(); @@ -294,7 +300,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *c if (config->GetnMarker_Fluid_Load() > 0) { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ for (unsigned short iDim = 0; iDim < nDim; iDim++){ - su2double val_sens = direct_solver->GetNodes()->ExtractFlowTraction_Sensitivity(iPoint,iDim); + su2double val_sens = direct_solver->GetNodes()->ExtractFlowTractionSensitivity(iPoint,iDim); nodes->SetFlowTractionSensitivity(iPoint, iDim, val_sens); } } @@ -366,12 +372,12 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - auto Coord = geometry->nodes->GetCoord(iPoint); + //auto Coord = geometry->nodes->GetCoord(iPoint); for (unsigned short iDim = 0; iDim < nDim; iDim++) { su2double Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); - AD::ResetInput(Coord[iDim]); + //AD::ResetInput(Coord[iDim]); if (!time_domain) { nodes->SetSensitivity(iPoint, iDim, Sensitivity); diff --git a/SU2_CFD/src/variables/CFEABoundVariable.cpp b/SU2_CFD/src/variables/CFEABoundVariable.cpp index d61f6b2e524..3124f0d93ca 100644 --- a/SU2_CFD/src/variables/CFEABoundVariable.cpp +++ b/SU2_CFD/src/variables/CFEABoundVariable.cpp @@ -75,3 +75,10 @@ void CFEABoundVariable::RegisterFlowTraction() { for (unsigned long iVar = 0; iVar < nVar; iVar++) AD::RegisterInput(FlowTraction(iVertex,iVar)); } + +void CFEABoundVariable::ResetInputFlowTraction() { + if (!fsi_analysis) return; + for (unsigned long iVertex = 0; iVertex < FlowTraction.rows(); iVertex++) + for (unsigned long iVar = 0; iVar < nVar; iVar++) + AD::ResetInput(FlowTraction(iVertex,iVar)); +} \ No newline at end of file diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 772b7012db7..e88a480b52d 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -276,7 +276,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-4.349377, 0.128475, -1.303589, 7.5407e-09, 2.3244] + discadj_fsi2.test_vals = [-4.349377, 0.189968, -1.303589, 7.5407e-09, 2.3244] discadj_fsi2.test_vals_aarch64 = [-3.479505, 0.127953, -1.303589, 7.5407e-09, 2.3244] discadj_fsi2.tol = 0.00001 test_list.append(discadj_fsi2) @@ -413,7 +413,7 @@ def main(): pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" pywrapper_Unst_FEA_AD.test_iter = 100 - pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.031988, 0.032015] + pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.031988, 0.032015, -18.44906, -18.45085] pywrapper_Unst_FEA_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py") pywrapper_Unst_FEA_AD.timeout = 1600 pywrapper_Unst_FEA_AD.tol = 0.00001 diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py index cc85ca00239..76a020631c6 100644 --- a/TestCases/py_wrapper/custom_load_fea/run_ad.py +++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py @@ -55,7 +55,7 @@ MESH_FORMAT= RECTANGLE MESH_BOX_SIZE= ( 17, 5, 0 ) -MESH_BOX_LENGTH= ( 0.5, 0.05, 0 ) +MESH_BOX_LENGTH= ( 0.5, __HEIGHT__, 0 ) OUTPUT_FILES= RESTART, PARAVIEW OUTPUT_WRT_FREQ= 1 @@ -84,6 +84,7 @@ UNST_ADJOINT_ITER= 21 ITER_AVERAGE_OBJ= 0 SOLUTION_FILENAME= restart.dat +VOLUME_OUTPUT= SENSITIVITY_T0 """ @@ -110,12 +111,16 @@ def ApplyLoad(driver, marker_id, peak_load): return derivatives -def RunPrimal(density, peak_load): - """Runs the primal solver for a given density and peak load and returns the time average objective function.""" +def RunPrimal(density, peak_load, height): + """ + Runs the primal solver for a given density, peak load, and beam height. + Returns the time average objective function. + """ comm = MPI.COMM_WORLD with open('config_unsteady.cfg', 'w') as f: - f.write(common_settings.replace('__DENSITY__', str(density)) + primal_settings) + f.write(common_settings.replace('__DENSITY__', str(density)).replace('__HEIGHT__', str(height)) + + primal_settings) # Initialize the primal driver of SU2, this includes solver preprocessing. try: @@ -144,15 +149,16 @@ def RunPrimal(density, peak_load): return tavg_stress_penalty -def RunAdjoint(density, peak_load): +def RunAdjoint(density, peak_load, height): """ Runs the adjoint solver and returns the sensitivity of the objective function to the peak - load and to the material density. + load, to the material density, and to the beam height. """ comm = MPI.COMM_WORLD with open('config_unsteady_ad.cfg', 'w') as f: - f.write(common_settings.replace('__DENSITY__', str(density)) + adjoint_settings) + f.write(common_settings.replace('__DENSITY__', str(density)).replace('__HEIGHT__', str(height)) + + adjoint_settings) # Initialize the adjoint driver of SU2, this includes solver preprocessing. try: @@ -182,7 +188,8 @@ def RunAdjoint(density, peak_load): driver.Postprocess() driver.Update() - # Accumulate sensitivies. + # Accumulate load sensitivies (the solver doesn't accumulate + # these for when they are used for FSI adjoints). for i_vertex in range(n_vertex): load_sens += derivatives[i_vertex] * driver.GetMarkerFEALoadSensitivity(marker_id, i_vertex)[1] @@ -192,10 +199,19 @@ def RunAdjoint(density, peak_load): rho_sens = driver.GetOutputValue('SENS_RHO_0') + height_sens = 0.0 + sensitivity = driver.Sensitivity(driver.GetSolverIndices()['ADJ.FEA']) + coords = driver.Coordinates() + + for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()): + y = coords(i_node, 1) + dy_dh = y / height + height_sens += dy_dh * sensitivity(i_node, 1) + # Finalize the solver and exit cleanly. driver.Finalize() - return comm.allreduce(load_sens), rho_sens + return comm.allreduce(load_sens), rho_sens, comm.allreduce(height_sens) def main(): @@ -203,26 +219,31 @@ def main(): rank = comm.Get_rank() # Run the primal with 2 loads to compute the sensitivity via finite differences. - obj_pert_load = RunPrimal(1, 2.002) - obj_pert_rho = RunPrimal(1.0001, 2) - obj = RunPrimal(1, 2) + obj_pert_height = RunPrimal(1, 2, 0.05001) + obj_pert_load = RunPrimal(1, 2.002, 0.05) + obj_pert_rho = RunPrimal(1.0001, 2, 0.05) + # Run the un-perturbed last to use the restarts for the adjoint. + obj = RunPrimal(1, 2, 0.05) + sens_height_fd = (obj_pert_height - obj) / 0.00001 sens_load_fd = (obj_pert_load - obj) / 0.002 sens_rho_fd = (obj_pert_rho - obj) / 0.0001 - sens_load, sens_rho = RunAdjoint(1, 2) + sens_load, sens_rho, sens_height = RunAdjoint(1, 2, 0.05) if rank == 0: - print(" Finite Differences\tDiscrete Adjoint") - print(f"Load {sens_load_fd}\t{sens_load}") - print(f"Rho {sens_rho_fd}\t{sens_rho}") + print(" Finite Differences\tDiscrete Adjoint") + print(f"Height {sens_height_fd}\t{sens_height}") + print(f"Load {sens_load_fd}\t{sens_load}") + print(f"Rho {sens_rho_fd}\t{sens_rho}") + assert abs(sens_height / sens_height_fd - 1) < 1e-4, "Error in geometric derivatives." assert abs(sens_load / sens_load_fd - 1) < 1e-4, "Error in load derivative." - assert abs(sens_rho / sens_rho_fd - 1) < 1e-3, "Error in density derivative." + assert abs(sens_rho / sens_rho_fd - 1) < 1e-3, "Error in material derivative." # Print results for the regression script to check. if rank == 0: print("\n") - print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho) + print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho, sens_height_fd, sens_height) if __name__ == '__main__': From 51a9cb948043481dbf4500711bf35474eeec7f0e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 22 Apr 2023 17:55:29 -0700 Subject: [PATCH 19/20] cleanup --- .../geometry/elements/CElementProperty.hpp | 8 +++-- .../include/variables/CFEABoundVariable.hpp | 7 +--- SU2_CFD/include/variables/CVariable.hpp | 7 +--- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 32 ++++++------------- SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 5 --- SU2_CFD/src/variables/CFEABoundVariable.cpp | 12 ++----- TestCases/parallel_regression_AD.py | 2 +- .../py_wrapper/custom_load_fea/run_ad.py | 2 +- 8 files changed, 23 insertions(+), 52 deletions(-) diff --git a/Common/include/geometry/elements/CElementProperty.hpp b/Common/include/geometry/elements/CElementProperty.hpp index fce7425694e..90efe7194b7 100644 --- a/Common/include/geometry/elements/CElementProperty.hpp +++ b/Common/include/geometry/elements/CElementProperty.hpp @@ -93,7 +93,7 @@ class CProperty { /*! * \brief Extract the derivative of the Design density. */ - inline virtual su2double GetAdjointDensity(void) const { return 0.0; } + inline virtual su2double GetAdjointDensity(void) { return 0.0; } /*! * \brief Register the Design density as an AD input variable. @@ -177,7 +177,11 @@ class CElementProperty final : public CProperty { /*! * \brief Extract the derivative of the Design density. */ - inline su2double GetAdjointDensity(void) const override { return SU2_TYPE::GetDerivative(design_rho); } + inline su2double GetAdjointDensity(void) override { + su2double der = SU2_TYPE::GetDerivative(design_rho); + AD::ResetInput(design_rho); + return der; + } /*! * \brief Register the Design density as an AD input variable. diff --git a/SU2_CFD/include/variables/CFEABoundVariable.hpp b/SU2_CFD/include/variables/CFEABoundVariable.hpp index 950fef78ebe..e0f8c273cc1 100644 --- a/SU2_CFD/include/variables/CFEABoundVariable.hpp +++ b/SU2_CFD/include/variables/CFEABoundVariable.hpp @@ -157,12 +157,7 @@ class CFEABoundVariable final : public CFEAVariable { /*! * \brief Register the flow tractions as input variable. */ - void RegisterFlowTraction() override; - - /*! - * \brief Resets the AD indices of the flow tractions. - */ - void ResetInputFlowTraction() override; + void RegisterFlowTraction(bool reset) override; /*! * \brief Extract the flow traction derivatives. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index efdbc25c9ce..aa75d07596d 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2160,12 +2160,7 @@ class CVariable { /*! * \brief A virtual member. */ - inline virtual void RegisterFlowTraction() { } - - /*! - * \brief A virtual member. - */ - inline virtual void ResetInputFlowTraction() { } + inline virtual void RegisterFlowTraction(bool reset) { } /*! * \brief A virtual member. diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 5b38d271497..b64d409e5cc 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -124,32 +124,21 @@ CDiscAdjFEASolver::~CDiscAdjFEASolver() { delete nodes; } void CDiscAdjFEASolver::SetRecording(CGeometry* geometry, CConfig *config){ - unsigned long iPoint; - unsigned short iDim, iVar; - /*--- Reset the solution to the initial (converged) solution ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto iVar = 0u; iVar < nVar; iVar++) direct_solver->GetNodes()->SetSolution(iPoint, iVar, nodes->GetSolution_Direct(iPoint)[iVar]); - - auto Coord = geometry->nodes->GetCoord(iPoint); - for (iDim = 0; iDim < nDim; iDim++) AD::ResetInput(Coord[iDim]); } /*--- Reset the input for time n ---*/ if (config->GetTime_Domain()) { - for (iPoint = 0; iPoint < nPoint; iPoint++) - for (iVar = 0; iVar < nVar; iVar++) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) + for (auto iVar = 0u; iVar < nVar; iVar++) AD::ResetInput(direct_solver->GetNodes()->GetSolution_time_n(iPoint)[iVar]); } - /*--- Set the Jacobian to zero since this is not done inside the meanflow iteration - * when running the discrete adjoint solver. ---*/ - - direct_solver->Jacobian.SetValZero(); - /*--- Set indices to zero ---*/ RegisterVariables(geometry, config, true); @@ -202,8 +191,6 @@ void CDiscAdjFEASolver::RegisterVariables(CGeometry *geometry, CConfig *config, for (iVar = 0; iVar < nDV; iVar++) AD::ResetInput(DV[iVar]); } - direct_solver->GetNodes()->ResetInputFlowTraction(); - if (!reset) { E.Register(); Nu.Register(); @@ -211,10 +198,11 @@ void CDiscAdjFEASolver::RegisterVariables(CGeometry *geometry, CConfig *config, Rho_DL.Register(); if (de_effects) EField.Register(); if (fea_dv) DV.Register(); + } - /*--- Register the flow tractions ---*/ - if (config->GetnMarker_Fluid_Load() > 0) - direct_solver->GetNodes()->RegisterFlowTraction(); + /*--- Register or reset the flow tractions ---*/ + if (config->GetnMarker_Fluid_Load() > 0) { + direct_solver->GetNodes()->RegisterFlowTraction(reset); } } @@ -372,12 +360,12 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - //auto Coord = geometry->nodes->GetCoord(iPoint); + auto Coord = geometry->nodes->GetCoord(iPoint); for (unsigned short iDim = 0; iDim < nDim; iDim++) { su2double Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); - //AD::ResetInput(Coord[iDim]); + AD::ResetInput(Coord[iDim]); if (!time_domain) { nodes->SetSensitivity(iPoint, iDim, Sensitivity); diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index d6d94cf0d58..e236b19114d 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -146,11 +146,6 @@ void CDiscAdjSolver::SetRecording(CGeometry* geometry, CConfig *config){ END_SU2_OMP_FOR } - /*--- Set the Jacobian to zero since this is not done inside the fluid iteration - * when running the discrete adjoint solver. ---*/ - - direct_solver->Jacobian.SetValZero(); - /*--- Set indices to zero ---*/ RegisterVariables(geometry, config, true); diff --git a/SU2_CFD/src/variables/CFEABoundVariable.cpp b/SU2_CFD/src/variables/CFEABoundVariable.cpp index 3124f0d93ca..012498896ad 100644 --- a/SU2_CFD/src/variables/CFEABoundVariable.cpp +++ b/SU2_CFD/src/variables/CFEABoundVariable.cpp @@ -69,16 +69,10 @@ void CFEABoundVariable::Clear_FlowTraction() { FlowTraction.setConstant(0.0); } void CFEABoundVariable::Clear_SurfaceLoad_Res() { Residual_Ext_Surf.setConstant(0.0); } -void CFEABoundVariable::RegisterFlowTraction() { +void CFEABoundVariable::RegisterFlowTraction(bool reset) { if (!fsi_analysis) return; for (unsigned long iVertex = 0; iVertex < FlowTraction.rows(); iVertex++) for (unsigned long iVar = 0; iVar < nVar; iVar++) - AD::RegisterInput(FlowTraction(iVertex,iVar)); + if (reset) AD::ResetInput(FlowTraction(iVertex,iVar)); + else AD::RegisterInput(FlowTraction(iVertex,iVar)); } - -void CFEABoundVariable::ResetInputFlowTraction() { - if (!fsi_analysis) return; - for (unsigned long iVertex = 0; iVertex < FlowTraction.rows(); iVertex++) - for (unsigned long iVar = 0; iVar < nVar; iVar++) - AD::ResetInput(FlowTraction(iVertex,iVar)); -} \ No newline at end of file diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index e88a480b52d..1f2d4297dd9 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -413,7 +413,7 @@ def main(): pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" pywrapper_Unst_FEA_AD.test_iter = 100 - pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.031988, 0.032015, -18.44906, -18.45085] + pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.319877, 0.320149, -0.184491, -0.184509] pywrapper_Unst_FEA_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py") pywrapper_Unst_FEA_AD.timeout = 1600 pywrapper_Unst_FEA_AD.tol = 0.00001 diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py index 76a020631c6..dd64982b51b 100644 --- a/TestCases/py_wrapper/custom_load_fea/run_ad.py +++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py @@ -243,7 +243,7 @@ def main(): # Print results for the regression script to check. if rank == 0: print("\n") - print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho, sens_height_fd, sens_height) + print(100, 100, sens_load_fd, sens_load, sens_rho_fd * 10, sens_rho * 10, sens_height_fd / 100, sens_height / 100) if __name__ == '__main__': From 2b7b16f456ea98b68cdd98dfd6819e0dd7294713 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 23 Apr 2023 10:53:45 -0700 Subject: [PATCH 20/20] update regression --- TestCases/parallel_regression_AD.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 1f2d4297dd9..61c79202fb8 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -276,7 +276,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-4.349377, 0.189968, -1.303589, 7.5407e-09, 2.3244] + discadj_fsi2.test_vals = [-4.349377, 0.192713, -1.303589, 7.5407e-09, 2.3244] discadj_fsi2.test_vals_aarch64 = [-3.479505, 0.127953, -1.303589, 7.5407e-09, 2.3244] discadj_fsi2.tol = 0.00001 test_list.append(discadj_fsi2)