diff --git a/Common/include/geometry/primal_grid/CPrimalGrid.hpp b/Common/include/geometry/primal_grid/CPrimalGrid.hpp index dfa2724c00d6..33f10270bf6f 100644 --- a/Common/include/geometry/primal_grid/CPrimalGrid.hpp +++ b/Common/include/geometry/primal_grid/CPrimalGrid.hpp @@ -30,6 +30,7 @@ #include #include +#include #include #include #include diff --git a/Common/include/grid_movement/CSurfaceMovement.hpp b/Common/include/grid_movement/CSurfaceMovement.hpp index 5558783a4e9a..564071e93383 100644 --- a/Common/include/grid_movement/CSurfaceMovement.hpp +++ b/Common/include/grid_movement/CSurfaceMovement.hpp @@ -443,7 +443,7 @@ class CSurfaceMovement : public CGridMovement { * \param[in] geometry - Geometrical definition of the problem. * \param[in] val_mesh_filename - Name of the grid output file. */ - void WriteFFDInfo(CSurfaceMovement **surface_movement, CGeometry **geometry, CConfig **config); + void WriteFFDInfo(CSurfaceMovement **surface_movement, CGeometry ****geometry, CConfig **config); /*! * \brief Get information about if there is a complete FFDBox definition, or it is necessary to diff --git a/Common/src/grid_movement/CSurfaceMovement.cpp b/Common/src/grid_movement/CSurfaceMovement.cpp index 308bd9a8b420..f4e7844ce32b 100644 --- a/Common/src/grid_movement/CSurfaceMovement.cpp +++ b/Common/src/grid_movement/CSurfaceMovement.cpp @@ -4672,7 +4672,7 @@ void CSurfaceMovement::MergeFFDInfo(CGeometry *geometry, CConfig *config) { } -void CSurfaceMovement::WriteFFDInfo(CSurfaceMovement** surface_movement, CGeometry **geometry, CConfig **config) { +void CSurfaceMovement::WriteFFDInfo(CSurfaceMovement** surface_movement, CGeometry ****geometry, CConfig **config) { unsigned short iOrder, jOrder, kOrder, iFFDBox, iCornerPoints, iParentFFDBox, iChildFFDBox, iZone; @@ -4683,13 +4683,13 @@ void CSurfaceMovement::WriteFFDInfo(CSurfaceMovement** surface_movement, CGeomet bool polar = (config[ZONE_0]->GetFFD_CoordSystem() == POLAR); - unsigned short nDim = geometry[ZONE_0]->GetnDim(); + unsigned short nDim = geometry[ZONE_0][INST_0][MESH_0]->GetnDim(); for (iZone = 0; iZone < config[ZONE_0]->GetnZone(); iZone++){ /*--- Merge the parallel FFD info ---*/ - surface_movement[iZone]->MergeFFDInfo(geometry[iZone], config[iZone]); + surface_movement[iZone]->MergeFFDInfo(geometry[iZone][INST_0][MESH_0], config[iZone]); if (iZone > 0){ diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 5481a8e9a977..66fb0665aa77 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -1,5 +1,5 @@ -/*! - * \file driver_structure.hpp +/*! + * \file CDriver.hpp * \brief Headers of the main subroutines for driving single or multi-zone problems. * The subroutines and functions are in the driver_structure.cpp file. * \author T. Economon, H. Kline, R. Sanchez @@ -28,13 +28,12 @@ #pragma once +#include "../../../Common/include/geometry/CGeometry.hpp" #include "../../../Common/include/parallelization/mpi_structure.hpp" - #include "../integration/CIntegration.hpp" -#include "../solvers/CSolver.hpp" #include "../interfaces/CInterface.hpp" - -#include "../../../Common/include/geometry/CGeometry.hpp" +#include "../solvers/CSolver.hpp" +#include "CDriverBase.hpp" using namespace std; @@ -49,69 +48,49 @@ class COutput; * \brief Parent class for driving an iteration of a single or multi-zone problem. * \author T. Economon */ -class CDriver { -protected: - int rank, /*!< \brief MPI Rank. */ - size; /*!< \brief MPI Size. */ - char* config_file_name; /*!< \brief Configuration file name of the problem.*/ - char runtime_file_name[MAX_STRING_SIZE]; - su2double StartTime, /*!< \brief Start point of the timer for performance benchmarking.*/ - StopTime, /*!< \brief Stop point of the timer for performance benchmarking.*/ - UsedTimePreproc, /*!< \brief Elapsed time between Start and Stop point of the timer for tracking preprocessing phase.*/ - UsedTimeCompute, /*!< \brief Elapsed time between Start and Stop point of the timer for tracking compute phase.*/ - UsedTimeOutput, /*!< \brief Elapsed time between Start and Stop point of the timer for tracking output phase.*/ - UsedTime; /*!< \brief Elapsed time between Start and Stop point of the timer.*/ - su2double BandwidthSum = 0.0; /*!< \brief Aggregate value of the bandwidth for writing restarts (to be average later).*/ - unsigned long IterCount, /*!< \brief Iteration count stored for performance benchmarking.*/ - OutputCount; /*!< \brief Output count stored for performance benchmarking.*/ - unsigned long DOFsPerPoint; /*!< \brief Number of unknowns at each vertex, i.e., number of equations solved. */ - su2double Mpoints; /*!< \brief Total number of grid points in millions in the calculation (including ghost points).*/ - su2double MpointsDomain; /*!< \brief Total number of grid points in millions in the calculation (excluding ghost points).*/ - su2double MDOFs; /*!< \brief Total number of DOFs in millions in the calculation (including ghost points).*/ - su2double MDOFsDomain; /*!< \brief Total number of DOFs in millions in the calculation (excluding ghost points).*/ - unsigned long TimeIter; /*!< \brief External iteration.*/ - ofstream **ConvHist_file; /*!< \brief Convergence history file.*/ - ofstream FSIHist_file; /*!< \brief FSI convergence history file.*/ - unsigned short iMesh, /*!< \brief Iterator on mesh levels.*/ - iZone, /*!< \brief Iterator on zones.*/ - nZone, /*!< \brief Total number of zones in the problem. */ - nDim, /*!< \brief Number of dimensions.*/ - iInst, /*!< \brief Iterator on instance levels.*/ - *nInst, /*!< \brief Total number of instances in the problem (per zone). */ - **interface_types; /*!< \brief Type of coupling between the distinct (physical) zones.*/ - bool StopCalc, /*!< \brief Stop computation flag.*/ - mixingplane, /*!< \brief mixing-plane simulation flag.*/ - fsi, /*!< \brief FSI simulation flag.*/ - fem_solver; /*!< \brief FEM fluid solver simulation flag. */ - CIteration ***iteration_container; /*!< \brief Container vector with all the iteration methods. */ - COutput **output_container; /*!< \brief Pointer to the COutput class. */ - CIntegration ****integration_container; /*!< \brief Container vector with all the integration methods. */ - CGeometry ****geometry_container; /*!< \brief Geometrical definition of the problem. */ - CSolver *****solver_container; /*!< \brief Container vector with all the solutions. */ - CNumerics ******numerics_container; /*!< \brief Description of the numerical method (the way in which the equations are solved). */ - CConfig **config_container; /*!< \brief Definition of the particular problem. */ - CConfig *driver_config; /*!< \brief Definition of the driver configuration. */ - COutput *driver_output; /*!< \brief Definition of the driver output. */ - CSurfaceMovement **surface_movement; /*!< \brief Surface movement classes of the problem. */ - CVolumetricMovement ***grid_movement; /*!< \brief Volume grid movement classes of the problem. */ - CFreeFormDefBox*** FFDBox; /*!< \brief FFD FFDBoxes of the problem. */ - vector > > - interpolator_container; /*!< \brief Definition of the interpolation method between non-matching discretizations of the interface. */ - CInterface ***interface_container; /*!< \brief Definition of the interface of information and physics. */ - bool dry_run; /*!< \brief Flag if SU2_CFD was started as dry-run via "SU2_CFD -d .cfg" */ - -public: +class CDriver : public CDriverBase { + protected: + char runtime_file_name[MAX_STRING_SIZE]; + su2double + UsedTimeOutput; /*!< \brief Elapsed time between Start and Stop point of the timer for tracking output phase.*/ + + su2double BandwidthSum = + 0.0; /*!< \brief Aggregate value of the bandwidth for writing restarts (to be average later).*/ + unsigned long IterCount, /*!< \brief Iteration count stored for performance benchmarking.*/ + OutputCount; /*!< \brief Output count stored for performance benchmarking.*/ + unsigned long DOFsPerPoint; /*!< \brief Number of unknowns at each vertex, i.e., number of equations solved. */ + su2double Mpoints; /*!< \brief Total number of grid points in millions in the calculation (including ghost points).*/ + su2double + MpointsDomain; /*!< \brief Total number of grid points in millions in the calculation (excluding ghost points).*/ + su2double MDOFs; /*!< \brief Total number of DOFs in millions in the calculation (including ghost points).*/ + su2double MDOFsDomain; /*!< \brief Total number of DOFs in millions in the calculation (excluding ghost points).*/ + + ofstream** ConvHist_file; /*!< \brief Convergence history file.*/ + ofstream FSIHist_file; /*!< \brief FSI convergence history file.*/ + + bool StopCalc, /*!< \brief Stop computation flag.*/ + mixingplane, /*!< \brief mixing-plane simulation flag.*/ + fsi, /*!< \brief FSI simulation flag.*/ + fem_solver; /*!< \brief FEM fluid solver simulation flag. */ + + CIteration*** iteration_container; /*!< \brief Container vector with all the iteration methods. */ + CIntegration**** integration_container; /*!< \brief Container vector with all the integration methods. */ + vector>> + interpolator_container; /*!< \brief Definition of the interpolation method between non-matching discretizations of + the interface. */ + CInterface*** interface_container; /*!< \brief Definition of the interface of information and physics. */ + bool dry_run; /*!< \brief Flag if SU2_CFD was started as dry-run via "SU2_CFD -d .cfg" */ + + public: /*! * \brief Constructor of the class. * \param[in] confFile - Configuration file name. * \param[in] val_nZone - Total number of zones. - * \param[in] val_nDim - Number of dimensions. * \param[in] MPICommunicator - MPI communicator for SU2. + * \param[in] dummy_geo - Dummy geometric definition of the problem. */ - CDriver(char* confFile, - unsigned short val_nZone, - SU2_Comm MPICommunicator, bool dummy_geo); + CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator, bool dummy_geo); /*! * \brief Destructor of the class. @@ -121,10 +100,9 @@ class CDriver { /*! * \brief A virtual member. */ - virtual void Run() { }; - -protected: + virtual void Run(){}; + protected: /*! * \brief Init_Containers */ @@ -132,145 +110,180 @@ class CDriver { /*! * \brief Read in the config and mesh files. + * \param[in] config - Definition of the particular problem. + * \param[in] driver_config - Definition of the driver configuration. */ - void Input_Preprocessing(CConfig **&config, CConfig *&driver_config); + void Input_Preprocessing(CConfig**& config, CConfig*& driver_config); /*! - * \brief Construction of the edge-based data structure and the multigrid structure. + * \brief Construction of the edge-based data structure and the multi-grid structure. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] dummy - Definition of the dummy driver. */ - void Geometrical_Preprocessing(CConfig *config, CGeometry **&geometry, bool dummy); + void Geometrical_Preprocessing(CConfig* config, CGeometry**& geometry, bool dummy); /*! * \brief Do the geometrical preprocessing for the DG FEM solver. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. */ - void Geometrical_Preprocessing_DGFEM(CConfig *config, CGeometry **&geometry); + void Geometrical_Preprocessing_DGFEM(CConfig* config, CGeometry**& geometry); /*! * \brief Geometrical_Preprocessing_FVM + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. */ - void Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geometry); + void Geometrical_Preprocessing_FVM(CConfig* config, CGeometry**& geometry); /*! * \brief Definition of the physics iteration class or within a single zone. - * \param[in] iteration_container - Pointer to the iteration container to be instantiated. * \param[in] config - Definition of the particular problem. - * \param[in] iZone - Index of the zone. + * \param[in] iteration - Pointer to the iteration container to be instantiated. */ - void Iteration_Preprocessing(CConfig *config, CIteration *&iteration) const; + void Iteration_Preprocessing(CConfig* config, CIteration*& iteration) const; /*! * \brief Definition and allocation of all solution classes. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. */ - void Solver_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***&solver); + void Solver_Preprocessing(CConfig* config, CGeometry** geometry, CSolver***& solver); /*! * \brief Restart of the solvers from the restart files. - * \param[in] solver_container - Container vector with all the solutions. + * \param[in] solver - Container vector with all the solutions. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] update_geo - Boolean to indicate if geometry should be updated. */ - void Solver_Restart(CSolver ***solver, CGeometry **geometry, CConfig *config, bool update_geo); + void Solver_Restart(CSolver*** solver, CGeometry** geometry, CConfig* config, bool update_geo); /*! * \brief Definition and allocation of all solution classes. - * \param[in] solver_container - Container vector with all the solutions. + * \param[in] solver - Container vector with all the solutions. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] val_iInst - Current solver instance. */ - void Solver_Postprocessing(CSolver ****solver, CGeometry **geometry, CConfig *config, unsigned short val_iInst); + void Solver_Postprocessing(CSolver**** solver, CGeometry** geometry, CConfig* config, unsigned short val_iInst); /*! * \brief Definition and allocation of all integration classes. * \param[in] config - Definition of the particular problem. * \param[in] solver - Container vector with all the solutions. - * \param[out] integration - Container vector with all the integration methods. + * \param[in] integration - Container vector with all the integration methods. */ - void Integration_Preprocessing(CConfig *config, CSolver **solver, CIntegration **&integration) const; + void Integration_Preprocessing(CConfig* config, CSolver** solver, CIntegration**& integration) const; /*! * \brief Definition and allocation of all integration classes. - * \param[in] integration_container - Container vector with all the integration methods. + * \param[in] integration - Container vector with all the integration methods. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] val_iInst - Current solver instance. */ - void Integration_Postprocessing(CIntegration ***integration, CGeometry **geometry, CConfig *config, unsigned short val_iInst); + void Integration_Postprocessing(CIntegration*** integration, CGeometry** geometry, CConfig* config, + unsigned short val_iInst); /*! * \brief Definition and allocation of all interface classes. + * \param[in] config - Definition of the particular problem. + * \param[in] solver - Container vector with all the solutions. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] interface_types - Type of coupling between the distinct (physical) zones. + * \param[in] interface - Class defining the physical transfer of information. + * \param[in] interpolation - Object defining the interpolation. */ - void Interface_Preprocessing(CConfig **config, CSolver *****solver, CGeometry ****geometry, - unsigned short **interface_types, CInterface ***interface, - vector > > &interpolation); + void Interface_Preprocessing(CConfig** config, CSolver***** solver, CGeometry**** geometry, + unsigned short** interface_types, CInterface*** interface, + vector>>& interpolation); /*! * \brief Definition and allocation of all solver classes. - * \param[in] numerics_container - Description of the numerical method (the way in which the equations are solved). - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method (the way in which the equations are solved). */ - void Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CNumerics ****&numerics) const; + void Numerics_Preprocessing(CConfig* config, CGeometry** geometry, CSolver*** solver, CNumerics****& numerics) const; /*! * \brief Helper to instantiate turbulence numerics specialized for different flow solvers. */ template - void InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, const CConfig *config, - const CSolver* turb_solver, CNumerics ****&numerics) const; + void InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, const CConfig* config, + const CSolver* turb_solver, CNumerics****& numerics) const; /*! * \brief Helper to instantiate transition numerics specialized for different flow solvers. */ template - void InstantiateTransitionNumerics(unsigned short nVar_Trans, int offset, const CConfig *config, - const CSolver* trans_solver, CNumerics ****&numerics) const; + void InstantiateTransitionNumerics(unsigned short nVar_Trans, int offset, const CConfig* config, + const CSolver* trans_solver, CNumerics****& numerics) const; /*! * \brief Helper to instantiate species transport numerics specialized for different flow solvers. */ template - void InstantiateSpeciesNumerics(unsigned short nVar_Species, int offset, const CConfig *config, - const CSolver* species_solver, CNumerics ****&numerics) const; + void InstantiateSpeciesNumerics(unsigned short nVar_Species, int offset, const CConfig* config, + const CSolver* species_solver, CNumerics****& numerics) const; /*! * \brief Definition and allocation of all solver classes. - * \param[in] numerics_container - Description of the numerical method (the way in which the equations are solved). - * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method (the way in which the equations are solved). + * \param[in] solver - Container vector with all the solutions. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] val_iInst - Current solver instance. */ - void Numerics_Postprocessing(CNumerics *****numerics, CSolver ***solver, CGeometry **geometry, CConfig *config, unsigned short val_iInst); + void Numerics_Postprocessing(CNumerics***** numerics, CSolver*** solver, CGeometry** geometry, CConfig* config, + unsigned short val_iInst); /*! * \brief GridMovement_Preprocessing - * \param config - * \param geometry - * \param solver - * \param iteration - * \param grid_movement - * \param surface_movement + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. + * \param[in] iteration - Container vector with all the iteration methods. + * \param[in] grid_movement - Volume grid movement classes of the problem. + * \param[in] surface_movement - Surface movement classes of the problem. */ - void DynamicMesh_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CIteration *iteration, CVolumetricMovement *&grid_movement, CSurfaceMovement *&surface_movement) const; + void DynamicMesh_Preprocessing(CConfig* config, CGeometry** geometry, CSolver*** solver, CIteration* iteration, + CVolumetricMovement*& grid_movement, CSurfaceMovement*& surface_movement) const; /*! * \brief Initialize Python interface functionalities + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. */ void PythonInterface_Preprocessing(CConfig** config, CGeometry**** geometry, CSolver***** solver); /*! * \brief Preprocess the output container. + * \param[in] config - Definition of the particular problem. + * \param[in] driver_config - Definition of the driver configuration. + * \param[in] output_container - Container vector with all the outputs. + * \param[in] driver_output - Definition of the driver output. */ - void Output_Preprocessing(CConfig **config, CConfig *driver_config, COutput **&output_container, COutput *&driver_output); + void Output_Preprocessing(CConfig** config, CConfig* driver_config, COutput**& output_container, + COutput*& driver_output); /*! * \brief Initiate value for static mesh movement such as the gridVel for the ROTATING frame. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. */ - void StaticMesh_Preprocessing(const CConfig *config, CGeometry **geometry); + void StaticMesh_Preprocessing(const CConfig* config, CGeometry** geometry); /*! * \brief Initiate value for static mesh movement such as the gridVel for the ROTATING frame. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. + * \param[in] interface - Class defining the physical transfer of information. */ void Turbomachinery_Preprocessing(CConfig** config, CGeometry**** geometry, CSolver***** solver, CInterface*** interface); @@ -320,14 +333,14 @@ class CDriver { virtual void Relaxation_Tractions(unsigned short donorZone, unsigned short targetZone, unsigned long iOuterIter) {} /*! - * \brief A virtual member to run a Block Gauss-Seidel iteration in multizone problems. + * \brief A virtual member to run a Block Gauss-Seidel iteration in multi-zone problems. */ - virtual void Run_GaussSeidel(){} + virtual void Run_GaussSeidel() {} /*! - * \brief A virtual member to run a Block-Jacobi iteration in multizone problems. + * \brief A virtual member to run a Block-Jacobi iteration in multi-zone problems. */ - virtual void Run_Jacobi(){} + virtual void Run_Jacobi() {} /*! * \brief A virtual member. @@ -340,8 +353,7 @@ class CDriver { */ void Print_DirectResidual(RECORDING kind_recording); -public: - + public: /*! * \brief Launch the computation for all zones and all physics. */ @@ -360,35 +372,37 @@ class CDriver { /*! * \brief Perform some pre-processing before an iteration of the physics. */ - virtual void Preprocess(unsigned long TimeIter){ } + virtual void Preprocess(unsigned long TimeIter) {} /*! * \brief Monitor the computation. */ - virtual bool Monitor(unsigned long TimeIter){ return false; } + virtual bool Monitor(unsigned long TimeIter) { return false; } /*! * \brief Output the solution in solution file. */ - virtual void Output(unsigned long TimeIter){ } + virtual void Output(unsigned long TimeIter) {} /*! - * \brief Perform a dynamic mesh deformation, including grid velocity computation and update of the multigrid structure. + * \brief Perform a dynamic mesh deformation, including grid velocity computation and update of the multi-grid + * structure. */ - virtual void DynamicMeshUpdate(unsigned long TimeIter) { } + virtual void DynamicMeshUpdate(unsigned long TimeIter) {} /*! - * \brief Perform a dynamic mesh deformation, including grid velocity computation and update of the multigrid structure. + * \brief Perform a dynamic mesh deformation, including grid velocity computation and update of the multi-grid + * structure. */ - virtual void DynamicMeshUpdate(unsigned short val_iZone, unsigned long TimeIter) { } + virtual void DynamicMeshUpdate(unsigned short val_iZone, unsigned long TimeIter) {} /*! * \brief Perform a mesh deformation as initial condition. */ - virtual void SetInitialMesh() { } + virtual void SetInitialMesh() {} /*! - * \brief Process the boundary conditions and update the multigrid structure. + * \brief Process the boundary conditions and update the multi-grid structure. */ void BoundaryConditionsUpdate(); @@ -434,28 +448,6 @@ class CDriver { */ passivedouble Get_LiftCoeff() const; - /*! - * \brief Get the number of vertices (halo nodes included) from a specified marker. - * \param[in] iMarker - Marker identifier. - * \return Number of vertices. - */ - unsigned long GetNumberVertices(unsigned short iMarker) const; - - /*! - * \brief Get the number of halo vertices from a specified marker. - * \param[in] iMarker - Marker identifier. - * \return Number of vertices. - */ - unsigned long GetNumberHaloVertices(unsigned short iMarker) const; - - /*! - * \brief Check if a vertex is physical or not (halo node) on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return True if the specified vertex is a halo node. - */ - bool IsAHaloNode(unsigned short iMarker, unsigned long iVertex) const; - /*! * \brief Get the number of external iterations. * \return Number of external iterations. @@ -480,22 +472,6 @@ class CDriver { */ string GetSurfaceFileName() const; - /*! - * \brief Get the global index of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return Vertex global index. - */ - unsigned long GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex) const; - - /*! - * \brief Get undeformed coordinates from the mesh solver. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return x,y,z coordinates of the vertex. - */ - vector GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) const; - /*! * \brief Get the temperature at a vertex on a specified marker. * \param[in] iMarker - Marker identifier. @@ -550,77 +526,19 @@ class CDriver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void Inlet_Preprocessing(CSolver ***solver, CGeometry **geometry, CConfig *config) const; - - /*! - * \brief Get the normal (vector) at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] unitNormal - Bool to normalise the vector. - * \return Normal (vector) at the vertex. - */ - vector GetVertexNormal(unsigned short iMarker, unsigned long iVertex, bool unitNormal = false) const; - - /*! - * \brief Get the unit normal (vector) at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return Unit normal (vector) at the vertex. - */ - inline vector GetVertexUnitNormal(unsigned short iMarker, unsigned long iVertex) const { - return GetVertexNormal(iMarker, iVertex, true); - } - - /*! - * \brief Get all the boundary markers tags. - * \return List of boundary markers tags. - */ - vector GetAllBoundaryMarkersTag() const; - - /*! - * \brief Get all the deformable boundary marker tags. - * \return List of deformable boundary markers tags. - */ - vector GetAllDeformMeshMarkersTag() const; + void Inlet_Preprocessing(CSolver*** solver, CGeometry** geometry, CConfig* config) const; /*! - * \brief Get all the heat transfer boundary markers tags. - * \return List of heat transfer boundary markers tags. + * \brief Get all the CHT boundary marker tags. + * \return List of CHT boundary markers tags. */ - vector GetAllCHTMarkersTag() const; + vector GetCHTMarkerTags() const; /*! - * \brief Get all the (subsonic) inlet boundary markers tags. + * \brief Get all the inlet boundary marker tags. * \return List of inlet boundary markers tags. */ - vector GetAllInletMarkersTag() const; - - /*! - * \brief Get all the boundary markers tags with their associated indices. - * \return List of boundary markers tags with their indices. - */ - map GetAllBoundaryMarkers() const; - - /*! - * \brief Get all the boundary markers tags with their associated types. - * \return List of boundary markers tags with their types. - */ - map GetAllBoundaryMarkersType() const; - - /*! - * \brief Set the mesh displacement for the elasticity mesh solver. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] DispX - Value of the mesh displacement in the direction X. - * \param[in] DispY - Value of the mesh displacement in the direction Y. - * \param[in] DispZ - Value of the mesh displacement in the direction Z. - */ - void SetMeshDisplacement(unsigned short iMarker, unsigned long iVertex, passivedouble DispX, passivedouble DispY, passivedouble DispZ); - - /*! - * \brief Communicate the boundary mesh displacements in a python call - */ - void CommunicateMeshDisplacement(void); + vector GetInletMarkerTags() const; /*! * \brief Return the sensitivities of the mesh boundary vertices. @@ -638,8 +556,8 @@ class CDriver { * \param[in] LoadX - Value of the load in the direction Y. * \param[in] LoadX - Value of the load in the direction Z. */ - void SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passivedouble LoadX, - passivedouble LoadY, passivedouble LoadZ); + void SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passivedouble LoadX, passivedouble LoadY, + passivedouble LoadZ); /*! * \brief Return the displacements from the FEA solver. @@ -693,7 +611,7 @@ class CDriver { * \param[in] val_AdjointZ - Value of the adjoint in the direction Z. */ void SetFlowLoad_Adjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, - passivedouble val_AdjointY, passivedouble val_AdjointZ); + passivedouble val_AdjointY, passivedouble val_AdjointZ); /*! * \brief Set the adjoint of the structural displacements (from an outside source) @@ -706,7 +624,7 @@ class CDriver { void SetSourceTerm_DispAdjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, passivedouble val_AdjointY, passivedouble val_AdjointZ); void SetSourceTerm_VelAdjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, - passivedouble val_AdjointY, passivedouble val_AdjointZ); + passivedouble val_AdjointY, passivedouble val_AdjointZ); /*! * \brief Set the position of the heat source. @@ -746,7 +664,7 @@ class CDriver { * \param[in] solution - Solution object with interface (iPoint,iVar). * \tparam Old - If true set "old solutions" instead. */ - template + template void SetAllSolutions(unsigned short iZone, bool adjoint, const Container& solution) { const auto nPoint = geometry_container[iZone][INST_0][MESH_0]->GetnPoint(); for (auto iSol = 0u, offset = 0u; iSol < MAX_SOLS; ++iSol) { @@ -754,8 +672,11 @@ class CDriver { if (!(solver && (solver->GetAdjoint() == adjoint))) continue; for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) - if (!Old) solver->GetNodes()->SetSolution(iPoint, iVar, solution(iPoint,offset+iVar)); - else solver->GetNodes()->SetSolution_Old(iPoint, iVar, solution(iPoint,offset+iVar)); + if (!Old) { + solver->GetNodes()->SetSolution(iPoint, iVar, solution(iPoint, offset + iVar)); + } else { + solver->GetNodes()->SetSolution_Old(iPoint, iVar, solution(iPoint, offset + iVar)); + } offset += solver->GetnVar(); } } @@ -763,9 +684,9 @@ class CDriver { /*! * \brief Set the "old solution" of all solvers (adjoint or primal) in a zone. */ - template + template void SetAllSolutionsOld(unsigned short iZone, bool adjoint, const Container& solution) { - SetAllSolutions(iZone, adjoint, solution); + SetAllSolutions(iZone, adjoint, solution); } /*! @@ -774,7 +695,7 @@ class CDriver { * \param[in] adjoint - True to consider adjoint solvers instead of primal. * \param[out] solution - Solution object with interface (iPoint,iVar). */ - template + template void GetAllSolutions(unsigned short iZone, bool adjoint, Container& solution) const { const auto nPoint = geometry_container[iZone][INST_0][MESH_0]->GetnPoint(); for (auto iSol = 0u, offset = 0u; iSol < MAX_SOLS; ++iSol) { @@ -783,11 +704,10 @@ class CDriver { const auto& sol = solver->GetNodes()->GetSolution(); for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) - solution(iPoint,offset+iVar) = SU2_TYPE::GetValue(sol(iPoint,iVar)); + solution(iPoint, offset + iVar) = SU2_TYPE::GetValue(sol(iPoint, iVar)); offset += solver->GetnVar(); } } - }; /*! @@ -797,24 +717,19 @@ class CDriver { * \author T. Economon, G. Gori */ class CFluidDriver : public CDriver { + protected: + unsigned long Max_Iter; -protected: - unsigned long Max_Iter; - -protected: - + protected: /*! * \brief Constructor of the class. * \param[in] confFile - Configuration file name. * \param[in] val_nZone - Total number of zones. - * \param[in] val_nDim - Number of dimensions. * \param[in] MPICommunicator - MPI communicator for SU2. */ - CFluidDriver(char* confFile, - unsigned short val_nZone, - SU2_Comm MPICommunicator); + CFluidDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator); -public: + public: /*! * \brief Destructor of the class. */ @@ -851,7 +766,8 @@ class CFluidDriver : public CDriver { void Preprocess(unsigned long Iter) override; /*! - * \brief Perform a dynamic mesh deformation, included grid velocity computation and the update of the multigrid structure (multiple zone). + * \brief Perform a dynamic mesh deformation, included grid velocity computation and the update of the multi-grid + * structure (multiple zone). */ void DynamicMeshUpdate(unsigned long TimeIter) override; @@ -859,10 +775,8 @@ class CFluidDriver : public CDriver { * \brief Transfer data among different zones (multiple zone). */ void Transfer_Data(unsigned short donorZone, unsigned short targetZone); - }; - /*! * \class CTurbomachineryDriver * \ingroup Drivers @@ -870,22 +784,17 @@ class CFluidDriver : public CDriver { * \author S. Vitale */ class CTurbomachineryDriver : public CFluidDriver { -private: + private: COutputLegacy* output_legacy; -public: - + public: /*! * \brief Constructor of the class. * \param[in] confFile - Configuration file name. * \param[in] val_nZone - Total number of zones. - * \param[in] val_nDim - Number of dimensions. - * \param[in] val_periodic - Bool for periodic BCs. * \param[in] MPICommunicator - MPI communicator for SU2. */ - CTurbomachineryDriver(char* confFile, - unsigned short val_nZone, - SU2_Comm MPICommunicator); + CTurbomachineryDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator); /*! * \brief Destructor of the class. @@ -912,9 +821,6 @@ class CTurbomachineryDriver : public CFluidDriver { * \brief Monitor the computation. */ bool Monitor(unsigned long TimeIter) override; - - - }; /*! @@ -924,24 +830,19 @@ class CTurbomachineryDriver : public CFluidDriver { * \author T. Economon */ class CHBDriver : public CFluidDriver { - -private: + private: COutputLegacy* output_legacy; unsigned short nInstHB; - su2double **D; /*!< \brief Harmonic Balance operator. */ - -public: + su2double** D; /*!< \brief Harmonic Balance operator. */ + public: /*! * \brief Constructor of the class. * \param[in] confFile - Configuration file name. * \param[in] val_nZone - Total number of zones. - * \param[in] val_nDim - Number of dimensions. * \param[in] MPICommunicator - MPI communicator for SU2. */ - CHBDriver(char* confFile, - unsigned short val_nZone, - SU2_Comm MPICommunicator); + CHBDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator); /*! * \brief Destructor of the class. diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp new file mode 100644 index 000000000000..132ee87ca717 --- /dev/null +++ b/SU2_CFD/include/drivers/CDriverBase.hpp @@ -0,0 +1,372 @@ +/*! + * \file CDriverBase.hpp + * \brief Base class for all drivers. + * \author H. Patel, A. Gastaldi + * \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 . + */ + +#pragma once + +#include "../../../Common/include/CConfig.hpp" +#include "../numerics/CNumerics.hpp" +#include "../output/COutput.hpp" +#include "../solvers/CSolver.hpp" + +/*! + * \class CDriverBase + * \ingroup Drivers + * \brief Base class for all drivers. + * \author H. Patel, A. Gastaldi + */ + +class CDriverBase { + protected: + int rank, /*!< \brief MPI Rank. */ + size; /*!< \brief MPI Size. */ + char* config_file_name; /*!< \brief Configuration file name of the problem. */ + + su2double StartTime, /*!< \brief Start point of the timer for performance benchmarking. */ + StopTime, /*!< \brief Stop point of the timer for performance benchmarking. */ + UsedTimePreproc, /*!< \brief Elapsed time between Start and Stop point of the timer for tracking preprocessing + phase. */ + UsedTimeCompute, /*!< \brief Elapsed time between Start and Stop point of the timer for tracking compute phase. */ + UsedTime; /*!< \brief Elapsed time between Start and Stop point of the timer. */ + + unsigned long TimeIter; + + unsigned short iMesh, /*!< \brief Iterator on mesh levels. */ + iZone, /*!< \brief Iterator on zones. */ + nZone, /*!< \brief Total number of zones in the problem. */ + nDim, /*!< \brief Number of dimensions. */ + iInst, /*!< \brief Iterator on instance levels. */ + *nInst, /*!< \brief Total number of instances in the problem (per zone). */ + **interface_types; /*!< \brief Type of coupling between the distinct (physical) zones. */ + + CConfig* driver_config; /*!< \brief Definition of the driver configuration. */ + COutput* driver_output; /*!< \brief Definition of the driver output. */ + + CConfig** config_container; /*!< \brief Definition of the particular problem. */ + COutput** output_container; /*!< \brief Pointer to the COutput class. */ + CGeometry**** geometry_container; /*!< \brief Geometrical definition of the problem. */ + CSolver***** solver_container; /*!< \brief Container vector with all the solutions. */ + CNumerics****** numerics_container; /*!< \brief Description of the numerical method (the way in which the equations + are solved). */ + CSurfaceMovement** surface_movement; /*!< \brief Surface movement classes of the problem. */ + CVolumetricMovement*** grid_movement; /*!< \brief Volume grid movement classes of the problem. */ + CFreeFormDefBox*** FFDBox; /*!< \brief FFD FFDBoxes of the problem. */ + + CConfig* main_config; /*!< \brief Reference to base (i.e. ZONE 0) configuration (used in driver API). */ + CGeometry* main_geometry; /*!< \brief Reference to base (i.e. ZONE, INST, MESH 0) geometry (used in driver API). */ + + public: + /*! + * \brief Constructor of the class. + * \param[in] confFile - Configuration file name. + * \param[in] val_nZone - Total number of zones. + * \param[in] MPICommunicator - MPI communicator for SU2. + */ + CDriverBase(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator); + + /*! + * \brief Destructor of the class. + */ + virtual ~CDriverBase(void); + + /*! + * \brief A virtual member. + */ + virtual void Preprocessing(){} + + /*! + * \brief A virtual member. + */ + virtual void Run(){} + + /*! + * \brief A virtual member. + */ + virtual void Update(){} + + /*! + * \brief A virtual member. + */ + virtual void Update_Legacy(){} + + /*! + * \brief A virtual member. + */ + virtual void Output(){} + + /*! + * \brief A virtual member. + */ + virtual void Postprocessing(){} + + /*! + * \brief Get the number of design variables. + * \return Number of design variables. + */ + unsigned short GetNumberDesignVariables() const; + + /*! + * \brief Get the number of FFD boxes. + * \return Number of FFD boxes. + */ + unsigned short GetNumberFFDBoxes() const; + + /*! + * \brief Get the number of dimensions of the mesh. + * \return Number of dimensions. + */ + unsigned long GetNumberDimensions() const; + + /*! + * \brief Get the number of elements in the mesh. + * \return Number of elements. + */ + unsigned long GetNumberElements() const; + + /*! + * \brief Get the global index of a mesh element. + * \param[in] iElem - Mesh element index. + * \return Global element index. + */ + unsigned long GetElementGlobalIndex(unsigned long iElem) const; + + /*! + * \brief Get the node indices of a mesh element. + * \param[in] iElem - Mesh element index. + * \return Element node indices (nNode). + */ + vector GetElementNodes(unsigned long iElem) const; + + /*! + * \brief Get the number of nodes in the mesh (including halos). + * \return Number of nodes. + */ + unsigned long GetNumberNodes() const; + + /*! + * \brief Get the number of halo nodes in the mesh. + * \return Number of halo nodes. + */ + unsigned long GetNumberHaloNodes() const; + + /*! + * \brief Get the global node index. + * \param[in] iPoint - Mesh node index. + * \return Global node index. + */ + unsigned long GetNodeGlobalIndex(unsigned long iPoint) const; + + /*! + * \brief Get the halo flag of a mesh node. + * \param[in] iPoint - Mesh node index. + * \return Node domain flag. + */ + bool GetNodeDomain(unsigned long iPoint) const; + + /*! + * \brief Get the initial (un-deformed) coordinates of a mesh node. + * \param[in] iPoint - Mesh node index. + * \return Initial node coordinates (nDim). + */ + vector GetInitialCoordinates(unsigned long iPoint) const; + + /*! + * \brief Get the current coordinates of a mesh node. + * \param[in] iPoint - Mesh node index. + * \return Node coordinates (nDim). + */ + vector GetCoordinates(unsigned long iPoint) const; + + /*! + * \brief Set the coordinates of a mesh node. + * \param[in] iPoint - Mesh node index. + * \param[in] values - Node coordinates (nDim). + */ + void SetCoordinates(unsigned long iPoint, vector values); + + /*! + * \brief Get the number of markers in the mesh. + * \return Number of markers. + */ + unsigned short GetNumberMarkers() const; + + /*! + * \brief Get all the boundary markers tags with their associated indices. + * \return Map of boundary markers tags to their indices. + */ + map GetMarkerIndices() const; + + /*! + * \brief Get all the boundary markers tags with their associated types. + * \return Map of boundary markers tags to their types. + */ + map GetMarkerTypes() const; + + /*! + * \brief Get all the boundary marker tags. + * \return List of boundary markers tags. + */ + vector GetMarkerTags() const; + + /*! + * \brief Get all the deformable boundary marker tags. + * \return List of deformable boundary markers tags. + */ + vector GetDeformableMarkerTags() const; + + /*! + * \brief Get the number of elements in the marker. + * \param[in] iMarker - Marker index. + * \return Number of elements. + */ + unsigned long GetNumberMarkerElements(unsigned short iMarker) const; + + /*! + * \brief Get the global index of a marker element. + * \param[in] iMarker - Marker index. + * \param[in] iElem - Marker element index. + * \return Global element index. + */ + unsigned long GetMarkerElementGlobalIndex(unsigned short iMarker, unsigned long iElem) const; + + /*! + * \brief Get the node indices of a marker element. + * \param[in] iMarker - Marker index. + * \param[in] iElem - Marker element index. + * \return Element node indices. + */ + vector GetMarkerElementNodes(unsigned short iMarker, unsigned long iElem) const; + + /*! + * \brief Get the number of nodes in the marker. + * \param[in] iMarker - Marker index. + * \return Number of nodes. + */ + unsigned long GetNumberMarkerNodes(unsigned short iMarker) const; + + /*! + * \brief Get the node index of a marker. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Marker vertex index. + * \return Marker vertex. + */ + unsigned long GetMarkerNode(unsigned short iMarker, unsigned long iVertex) const; + + /*! + * \brief Get the normal vector of a marker vertex. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Marker vertex index. + * \param[in] normalize - If true, the unit (i.e. normalized) normal vector is returned. + * \return Node normal vector (nDim). + */ + vector GetMarkerVertexNormals(unsigned short iMarker, unsigned long iVertex, + bool normalize = false) const; + + /*! + * \brief Get the displacements of a marker vertex. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Marker vertex index. + * \return Node displacements (nDim). + */ + vector GetMarkerDisplacements(unsigned short iMarker, unsigned long iVertex) const; + + /*! + * \brief Set the displacements of a marker vertex. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Marker vertex index. + * \param[in] values - Node displacements (nDim). + */ + void SetMarkerDisplacements(unsigned short iMarker, unsigned long iVertex, vector values); + + /*! + * \brief Get the velocities of a marker vertex. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Marker vertex index. + * \return Node velocities (nDim). + */ + vector GetMarkerVelocities(unsigned short iMarker, unsigned long iVertex) const; + + /*! + * \brief Set the velocities of a marker vertex. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Marker vertex index. + * \param[in] values - Node velocities (nDim). + */ + void SetMarkerVelocities(unsigned short iMarker, unsigned long iVertex, vector values); + + /*! + * \brief Communicate the boundary mesh displacements. + */ + void CommunicateMeshDisplacements(void); + + protected: + /*! + * \brief Initialize containers. + */ + void SetContainers_Null(); + + /*! + * \brief Read in the config and mesh files. + * \param[in] config - Definition of the particular problem. + * \param[in] driver_config - Definition of the driver configuration. + */ + void Input_Preprocessing(CConfig**& config, CConfig*& driver_config); + + /*! + * \brief Construction of the edge-based data structure and the multi-grid structure. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] dummy - Definition of the dummy driver. + */ + void Geometrical_Preprocessing(CConfig* config, CGeometry**& geometry, bool dummy); + + /*! + * \brief Definition and allocation of all solution classes. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. + */ + void Solver_Preprocessing(CConfig* config, CGeometry** geometry, CSolver***& solver); + + /*! + * \brief Definition and allocation of all solver classes. + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method (the way in which the equations are solved). + */ + void Numerics_Preprocessing(CConfig* config, CGeometry** geometry, CSolver*** solver, CNumerics****& numerics) const; + + /*! + * \brief Preprocess the output container. + * \param[in] config - Definition of the particular problem. + * \param[in] driver_config - Definition of the driver configuration. + * \param[in] output_container - Container vector with all the outputs. + * \param[in] driver_output - Definition of the driver output. + */ + void Output_Preprocessing(CConfig** config, CConfig* driver_config, COutput**& output_container, + COutput*& driver_output); +}; diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp index d03a841f5e38..79876e016c11 100644 --- a/SU2_CFD/include/solvers/CMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CMeshSolver.hpp @@ -125,6 +125,7 @@ class CMeshSolver final : public CFEASolver { /*! * \brief Grid deformation using the linear elasticity equations. * \param[in] geometry - Geometrical definition of the problem. + * \param[in] numerics - Numerics used in the solution. * \param[in] config - Definition of the particular problem. */ void DeformMesh(CGeometry **geometry, @@ -133,11 +134,10 @@ class CMeshSolver final : public CFEASolver { /*! * \brief Set the stiffness of the mesh. - * \param[in] geometry - Geometrical definition of the problem. + * \param[in] numerics - Numerics used in the solution. * \param[in] config - Definition of the particular problem. */ - void SetMesh_Stiffness(CGeometry **geometry, - CNumerics **numerics, + void SetMesh_Stiffness(CNumerics **numerics, CConfig *config) override; /*! diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index f87fcf58c44a..954b06ad75d4 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -4167,8 +4167,16 @@ class CSolver { * \param[in] config - Definition of the particular problem. * \param[in] referenceCoord - Determine if the mesh is deformed from the reference or from the current coordinates. */ - inline virtual void SetMesh_Stiffness(CGeometry **geometry, - CNumerics **numerics, + inline virtual void DeformMesh(CGeometry *geometry, + CNumerics **numerics, + CConfig *config) { } + + /*! + * \brief A virtual member. + * \param[in] config - Definition of the particular problem. + * \param[in] referenceCoord - Determine if the mesh is deformed from the reference or from the current coordinates. + */ + inline virtual void SetMesh_Stiffness(CNumerics **numerics, CConfig *config) { } /*! diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 121cd94e9a55..cf69fefeddde 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -107,10 +107,9 @@ #include CDriver::CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator, bool dummy_geo) : - config_file_name(confFile), StartTime(0.0), StopTime(0.0), UsedTime(0.0), - TimeIter(0), nZone(val_nZone), StopCalc(false), fsi(false), fem_solver(false), dry_run(dummy_geo) { +CDriverBase(confFile, val_nZone, MPICommunicator), StopCalc(false), fsi(false), fem_solver(false), dry_run(dummy_geo) { - /*--- Initialize Medipack (must also be here so it is initialized from python) ---*/ + /*--- Initialize Medipack (must also be here so it is initialized from python) ---*/ #ifdef HAVE_MPI #if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) SU2_MPI::Init_AMPI(); @@ -602,6 +601,10 @@ void CDriver::Input_Preprocessing(CConfig **&config, CConfig *&driver_config) { } } + /*--- Keep a reference to the main (ZONE 0) config. ---*/ + + main_config = config_container[ZONE_0]; + /*--- Determine whether or not the FEM solver is used, which decides the type of * geometry classes that are instantiated. Only adapted for single-zone problems ---*/ @@ -706,6 +709,9 @@ void CDriver::Geometrical_Preprocessing(CConfig* config, CGeometry **&geometry, } + /*--- Keep a reference to the main (ZONE_0, INST_0, MESH_0) geometry. ---*/ + + main_geometry = geometry_container[ZONE_0][INST_0][MESH_0]; } void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geometry) { diff --git a/SU2_CFD/src/drivers/CDriverBase.cpp b/SU2_CFD/src/drivers/CDriverBase.cpp new file mode 100644 index 000000000000..95d755c1e59d --- /dev/null +++ b/SU2_CFD/src/drivers/CDriverBase.cpp @@ -0,0 +1,377 @@ +/*! + * \file CDriverBase.hpp + * \brief Base class template for all drivers. + * \author H. Patel, A. Gastaldi + * \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 . + */ + +#include "../../include/drivers/CDriverBase.hpp" + +#include "../../../Common/include/geometry/CPhysicalGeometry.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + +using namespace std; + +CDriverBase::CDriverBase(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunicator) + : config_file_name(confFile), StartTime(0.0), StopTime(0.0), UsedTime(0.0), TimeIter(0), nZone(val_nZone) {} + +CDriverBase::~CDriverBase(void) {} + +void CDriverBase::SetContainers_Null() { + /*--- Create pointers to all the classes that may be used by drivers. In general, the pointers are instantiated + * down a hierarchy over all zones, multi-grid levels, equation sets, and equation terms as described in the comments + * below. ---*/ + + config_container = nullptr; + output_container = nullptr; + geometry_container = nullptr; + solver_container = nullptr; + numerics_container = nullptr; + + surface_movement = nullptr; + grid_movement = nullptr; + FFDBox = nullptr; + + config_container = new CConfig*[nZone](); + output_container = new COutput*[nZone](); + geometry_container = new CGeometry***[nZone](); + solver_container = new CSolver****[nZone](); + numerics_container = new CNumerics*****[nZone](); + surface_movement = new CSurfaceMovement*[nZone](); + grid_movement = new CVolumetricMovement**[nZone](); + FFDBox = new CFreeFormDefBox**[nZone](); + + nInst = new unsigned short[nZone]; + + for (iZone = 0; iZone < nZone; iZone++) { + nInst[iZone] = 1; + } + + driver_config = nullptr; + driver_output = nullptr; + + main_config = nullptr; + main_geometry = nullptr; +} + +unsigned short CDriverBase::GetNumberDesignVariables() const { return main_config->GetnDV(); } + +unsigned short CDriverBase::GetNumberFFDBoxes() const { return main_config->GetnFFDBox(); } + +unsigned long CDriverBase::GetNumberDimensions() const { return main_geometry->GetnDim(); } + +unsigned long CDriverBase::GetNumberElements() const { return main_geometry->GetnElem(); } + +unsigned long CDriverBase::GetElementGlobalIndex(unsigned long iElem) const { + if (iElem >= GetNumberElements()) { + SU2_MPI::Error("Element index exceeds size.", CURRENT_FUNCTION); + } + return main_geometry->elem[iElem]->GetGlobalIndex(); +} + +vector CDriverBase::GetElementNodes(unsigned long iElem) const { + if (iElem >= GetNumberElements()) { + SU2_MPI::Error("Element index exceeds size.", CURRENT_FUNCTION); + } + const auto nNode = main_geometry->elem[iElem]->GetnNodes(); + vector values(nNode); + + for (auto iNode = 0u; iNode < nNode; iNode++) { + values[iNode] = main_geometry->elem[iElem]->GetNode(iNode); + } + return values; +} + +unsigned long CDriverBase::GetNumberNodes() const { return main_geometry->GetnPoint(); } + +unsigned long CDriverBase::GetNumberHaloNodes() const { + return main_geometry->GetnPoint() - main_geometry->GetnPointDomain(); +} + +unsigned long CDriverBase::GetNodeGlobalIndex(unsigned long iPoint) const { + if (iPoint >= GetNumberNodes()) { + SU2_MPI::Error("Node index exceeds mesh size.", CURRENT_FUNCTION); + } + return main_geometry->nodes->GetGlobalIndex(iPoint); +} + +bool CDriverBase::GetNodeDomain(unsigned long iPoint) const { + if (iPoint >= GetNumberNodes()) { + SU2_MPI::Error("Node index exceeds mesh size.", CURRENT_FUNCTION); + } + return main_geometry->nodes->GetDomain(iPoint); +} + +vector CDriverBase::GetInitialCoordinates(unsigned long iPoint) const { + if (iPoint >= GetNumberNodes()) { + SU2_MPI::Error("Node index exceeds mesh size.", CURRENT_FUNCTION); + } + vector values(nDim, 0.0); + + if (main_config->GetDeform_Mesh()) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + const su2double value = solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord(iPoint, iDim); + values[iDim] = SU2_TYPE::GetValue(value); + } + } + return values; +} + +vector CDriverBase::GetCoordinates(unsigned long iPoint) const { + if (iPoint >= GetNumberNodes()) { + SU2_MPI::Error("Node index exceeds mesh size.", CURRENT_FUNCTION); + } + vector values; + + for (auto iDim = 0u; iDim < nDim; iDim++) { + const su2double value = main_geometry->nodes->GetCoord(iPoint, iDim); + values.push_back(SU2_TYPE::GetValue(value)); + } + return values; +} + +void CDriverBase::SetCoordinates(unsigned long iPoint, vector values) { + if (iPoint >= GetNumberNodes()) { + SU2_MPI::Error("Node index exceeds mesh size.", CURRENT_FUNCTION); + } + if (values.size() != nDim) { + SU2_MPI::Error("Invalid number of dimensions!", CURRENT_FUNCTION); + } + for (auto iDim = 0u; iDim < nDim; iDim++) { + main_geometry->nodes->SetCoord(iPoint, iDim, values[iDim]); + } +} + +unsigned short CDriverBase::GetNumberMarkers() const { return main_config->GetnMarker_All(); } + +map CDriverBase::GetMarkerIndices() const { + const auto nMarker = main_config->GetnMarker_All(); + map indexMap; + + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { + indexMap[main_config->GetMarker_All_TagBound(iMarker)] = iMarker; + } + return indexMap; +} + +map CDriverBase::GetMarkerTypes() const { + map typeMap; + string type; + + for (auto iMarker = 0u; iMarker < main_config->GetnMarker_All(); iMarker++) { + auto tag = main_config->GetMarker_All_TagBound(iMarker); + auto kindBC = main_config->GetMarker_All_KindBC(iMarker); + + switch (kindBC) { + case EULER_WALL: + type = "EULER_WALL"; + break; + case FAR_FIELD: + type = "FARFIELD"; + break; + case ISOTHERMAL: + type = "ISOTHERMAL"; + break; + case HEAT_FLUX: + type = "HEAT_FLUX"; + break; + case HEAT_TRANSFER: + type = "HEAT_TRANSFER"; + break; + case INLET_FLOW: + type = "INLET_FLOW"; + break; + case OUTLET_FLOW: + type = "OUTLET_FLOW"; + break; + case SUPERSONIC_INLET: + type = "SUPERSONIC_INLET"; + break; + case SUPERSONIC_OUTLET: + type = "SUPERSONIC_OUTLET"; + break; + case RIEMANN_BOUNDARY: + type = "RIEMANN"; + break; + case GILES_BOUNDARY: + type = "GILES"; + break; + case DISPLACEMENT_BOUNDARY: + type = "DISPLACEMENT"; + break; + case LOAD_BOUNDARY: + type = "LOAD"; + break; + case PERIODIC_BOUNDARY: + type = "PERIODIC"; + break; + case SYMMETRY_PLANE: + type = "SYMMETRY"; + break; + case SEND_RECEIVE: + type = "SEND_RECEIVE"; + break; + default: + type = "UNKNOWN_TYPE"; + } + typeMap[tag] = type; + } + + return typeMap; +} + +vector CDriverBase::GetMarkerTags() const { + const auto nMarker = main_config->GetnMarker_All(); + vector tags(nMarker); + + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { + tags[iMarker] = main_config->GetMarker_All_TagBound(iMarker); + } + return tags; +} + +vector CDriverBase::GetDeformableMarkerTags() const { + const auto nMarker = main_config->GetnMarker_Deform_Mesh(); + vector tags(nMarker); + + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { + tags[iMarker] = main_config->GetMarker_Deform_Mesh_TagBound(iMarker); + } + return tags; +} + +unsigned long CDriverBase::GetNumberMarkerElements(unsigned short iMarker) const { + if (iMarker >= GetNumberMarkers()) { + SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION); + } + return main_geometry->GetnElem_Bound(iMarker); +} + +unsigned long CDriverBase::GetMarkerElementGlobalIndex(unsigned short iMarker, unsigned long iElem) const { + if (iElem >= GetNumberMarkerElements(iMarker)) { + SU2_MPI::Error("Marker element index exceeds size.", CURRENT_FUNCTION); + } + return main_geometry->bound[iMarker][iElem]->GetGlobalIndex(); +} + +vector CDriverBase::GetMarkerElementNodes(unsigned short iMarker, unsigned long iElem) const { + if (iElem >= GetNumberMarkerElements(iMarker)) { + SU2_MPI::Error("Marker element index exceeds size.", CURRENT_FUNCTION); + } + unsigned short nNode = main_geometry->bound[iMarker][iElem]->GetnNodes(); + vector values(nNode); + + for (auto iNode = 0u; iNode < nNode; iNode++) { + values[iNode] = main_geometry->bound[iMarker][iElem]->GetNode(iNode); + } + return values; +} + +unsigned long CDriverBase::GetNumberMarkerNodes(unsigned short iMarker) const { + if (iMarker >= GetNumberMarkers()) { + SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION); + } + return main_geometry->GetnVertex(iMarker); +} + +unsigned long CDriverBase::GetMarkerNode(unsigned short iMarker, unsigned long iVertex) const { + if (iVertex >= GetNumberMarkerNodes(iMarker)) { + SU2_MPI::Error("Vertex index exceeds marker size.", CURRENT_FUNCTION); + } + return geometry_container[MESH_0][INST_0][ZONE_0]->vertex[iMarker][iVertex]->GetNode(); +} + +vector CDriverBase::GetMarkerVertexNormals(unsigned short iMarker, unsigned long iVertex, + bool normalize) const { + if (iVertex >= GetNumberMarkerNodes(iMarker)) { + SU2_MPI::Error("Vertex index exceeds marker size.", CURRENT_FUNCTION); + } + const auto* normal = main_geometry->vertex[iMarker][iVertex]->GetNormal(); + const su2double area = normalize ? GeometryToolbox::Norm(nDim, normal) : 1.0; + + vector values(nDim, 0.0); + + for (auto iDim = 0u; iDim < nDim; iDim++) { + values[iDim] = SU2_TYPE::GetValue(normal[iDim] / area); + } + return values; +} + +vector CDriverBase::GetMarkerDisplacements(unsigned short iMarker, unsigned long iVertex) const { + vector values(nDim, 0.0); + + if (main_config->GetDeform_Mesh()) { + const auto iPoint = GetMarkerNode(iMarker, iVertex); + for (auto iDim = 0u; iDim < nDim; iDim++) { + const su2double value = solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetBound_Disp(iPoint, iDim); + values[iDim] = SU2_TYPE::GetValue(value); + } + } + return values; +} + +void CDriverBase::SetMarkerDisplacements(unsigned short iMarker, unsigned long iVertex, vector values) { + if (!main_config->GetDeform_Mesh()) { + SU2_MPI::Error("Mesh solver is not defined!", CURRENT_FUNCTION); + } + if (values.size() != nDim) { + SU2_MPI::Error("Invalid number of dimensions!", CURRENT_FUNCTION); + } + const auto iPoint = GetMarkerNode(iMarker, iVertex); + + for (auto iDim = 0u; iDim < nDim; iDim++) { + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint, iDim, values[iDim]); + } +} + +vector CDriverBase::GetMarkerVelocities(unsigned short iMarker, unsigned long iVertex) const { + vector values(nDim, 0.0); + + if (main_config->GetDeform_Mesh()) { + const auto iPoint = GetMarkerNode(iMarker, iVertex); + for (auto iDim = 0u; iDim < nDim; iDim++) { + const su2double value = solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetBound_Vel(iPoint, iDim); + values[iDim] = SU2_TYPE::GetValue(value); + } + } + return values; +} + +void CDriverBase::SetMarkerVelocities(unsigned short iMarker, unsigned long iVertex, vector values) { + if (!main_config->GetDeform_Mesh()) { + SU2_MPI::Error("Mesh solver is not defined!", CURRENT_FUNCTION); + } + if (values.size() != nDim) { + SU2_MPI::Error("Invalid number of dimensions!", CURRENT_FUNCTION); + } + const auto iPoint = GetMarkerNode(iMarker, iVertex); + + for (auto iDim = 0u; iDim < nDim; iDim++) { + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Vel(iPoint, iDim, values[iDim]); + } +} + +void CDriverBase::CommunicateMeshDisplacements(void) { + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->InitiateComms(main_geometry, main_config, MESH_DISPLACEMENTS); + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->CompleteComms(main_geometry, main_config, MESH_DISPLACEMENTS); +} diff --git a/SU2_CFD/src/iteration/CIteration.cpp b/SU2_CFD/src/iteration/CIteration.cpp index 337e24edaea9..5b2eaa5cef16 100644 --- a/SU2_CFD/src/iteration/CIteration.cpp +++ b/SU2_CFD/src/iteration/CIteration.cpp @@ -163,7 +163,7 @@ void CIteration::SetMesh_Deformation(CGeometry** geometry, CSolver** solver, CNu /*--- Set the stiffness of each element mesh into the mesh numerics ---*/ - solver[MESH_SOL]->SetMesh_Stiffness(geometry, numerics[MESH_SOL], config); + solver[MESH_SOL]->SetMesh_Stiffness(numerics[MESH_SOL], config); /*--- Deform the volume grid around the new boundary locations ---*/ diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index fd5a5252d0e0..03fe477787df 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -162,7 +162,8 @@ su2_cfd_src += files(['drivers/CDriver.cpp', 'drivers/CSinglezoneDriver.cpp', 'drivers/CDiscAdjMultizoneDriver.cpp', 'drivers/CDiscAdjSinglezoneDriver.cpp', - 'drivers/CDummyDriver.cpp']) + 'drivers/CDummyDriver.cpp', + 'drivers/CDriverBase.cpp']) su2_cfd_src += files(['integration/CIntegration.cpp', 'integration/CIntegrationFactory.cpp', diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 32672e87eef3..95e8cd321f76 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -25,22 +25,18 @@ * License along with SU2. If not, see . */ - +#include "../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../include/drivers/CDriver.hpp" #include "../include/drivers/CSinglezoneDriver.hpp" -#include "../../Common/include/toolboxes/geometry_toolbox.hpp" - -void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geometry, CSolver *****solver){ +void CDriver::PythonInterface_Preprocessing(CConfig** config, CGeometry**** geometry, CSolver***** solver) { int rank = MASTER_NODE; SU2_MPI::Comm_rank(SU2_MPI::GetComm(), &rank); - /* --- Initialize boundary conditions customization, this is achieve through the Python wrapper --- */ - for(iZone=0; iZone < nZone; iZone++){ - - if (config[iZone]->GetnMarker_PyCustom() > 0){ - - if (rank == MASTER_NODE) cout << endl << "----------------- Python Interface Preprocessing ( Zone "<< iZone <<" ) -----------------" << endl; + /* --- Initialize boundary conditions customization, this is achieved through the Python wrapper. --- */ + for (iZone = 0; iZone < nZone; iZone++) { + if (config[iZone]->GetnMarker_PyCustom() > 0) { + if (rank == MASTER_NODE) cout << "----------------- Python Interface Preprocessing ( Zone " << iZone << " ) -----------------" << endl; if (rank == MASTER_NODE) cout << "Setting customized boundary conditions for zone " << iZone << endl; for (iMesh = 0; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) { @@ -50,21 +46,23 @@ void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geom if ((config[iZone]->GetKind_Solver() == MAIN_SOLVER::EULER) || (config[iZone]->GetKind_Solver() == MAIN_SOLVER::NAVIER_STOKES) || - (config[iZone]->GetKind_Solver() == MAIN_SOLVER::RANS)) { - + (config[iZone]->GetKind_Solver() == MAIN_SOLVER::RANS) || + (config[iZone]->GetKind_Solver() == MAIN_SOLVER::INC_EULER) || + (config[iZone]->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES) || + (config[iZone]->GetKind_Solver() == MAIN_SOLVER::INC_RANS) || + (config[iZone]->GetKind_Solver() == MAIN_SOLVER::NEMO_EULER) || + (config[iZone]->GetKind_Solver() == MAIN_SOLVER::NEMO_NAVIER_STOKES)) { solver[iZone][INST_0][MESH_0][FLOW_SOL]->UpdateCustomBoundaryConditions(geometry[iZone][INST_0], config[iZone]); } } } - } ///////////////////////////////////////////////////////////////////////////// -/* Functions related to the global performance indices (Lift, Drag, ecc..) */ +/* Functions related to the global performance indices (Lift, Drag, etc.) */ ///////////////////////////////////////////////////////////////////////////// passivedouble CDriver::Get_Drag() const { - unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); su2double CDrag, factor, val_Drag; @@ -73,13 +71,12 @@ passivedouble CDriver::Get_Drag() const { factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CDrag = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CD(); - val_Drag = CDrag*factor; + val_Drag = CDrag * factor; return SU2_TYPE::GetValue(val_Drag); } passivedouble CDriver::Get_Lift() const { - unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); su2double CLift, factor, val_Lift; @@ -88,13 +85,12 @@ passivedouble CDriver::Get_Lift() const { factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CLift = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CL(); - val_Lift = CLift*factor; + val_Lift = CLift * factor; return SU2_TYPE::GetValue(val_Lift); } passivedouble CDriver::Get_Mx() const { - unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); su2double CMx, RefLengthCoeff, factor, val_Mx; @@ -105,13 +101,12 @@ passivedouble CDriver::Get_Mx() const { factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CMx = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CMx(); - val_Mx = CMx*factor*RefLengthCoeff; + val_Mx = CMx * factor * RefLengthCoeff; return SU2_TYPE::GetValue(val_Mx); } passivedouble CDriver::Get_My() const { - unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); su2double CMy, RefLengthCoeff, factor, val_My; @@ -122,13 +117,12 @@ passivedouble CDriver::Get_My() const { factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CMy = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CMy(); - val_My = CMy*factor*RefLengthCoeff; + val_My = CMy * factor * RefLengthCoeff; return SU2_TYPE::GetValue(val_My); } passivedouble CDriver::Get_Mz() const { - unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); su2double CMz, RefLengthCoeff, factor, val_Mz; @@ -139,13 +133,12 @@ passivedouble CDriver::Get_Mz() const { factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CMz = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CMz(); - val_Mz = CMz*factor*RefLengthCoeff; + val_Mz = CMz * factor * RefLengthCoeff; return SU2_TYPE::GetValue(val_Mz); } passivedouble CDriver::Get_DragCoeff() const { - unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); su2double CDrag; @@ -156,7 +149,6 @@ passivedouble CDriver::Get_DragCoeff() const { } passivedouble CDriver::Get_LiftCoeff() const { - unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); su2double CLift; @@ -166,129 +158,25 @@ passivedouble CDriver::Get_LiftCoeff() const { return SU2_TYPE::GetValue(CLift); } -///////////////////////////////////////////////////////////////////////////// -/* Functions to obtain information from the geometry/mesh */ -///////////////////////////////////////////////////////////////////////////// - -unsigned long CDriver::GetNumberVertices(unsigned short iMarker) const { - - return geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; - -} - -unsigned long CDriver::GetNumberHaloVertices(unsigned short iMarker) const { - - unsigned long nHaloVertices, iVertex, iPoint; - - nHaloVertices = 0; - for(iVertex = 0; iVertex < geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; iVertex++){ - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - if(!(geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint))) nHaloVertices += 1; - } - - return nHaloVertices; - -} - -unsigned long CDriver::GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex) const { - - unsigned long iPoint, GlobalIndex; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - GlobalIndex = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetGlobalIndex(iPoint); - - return GlobalIndex; - -} - -bool CDriver::IsAHaloNode(unsigned short iMarker, unsigned long iVertex) const { - - unsigned long iPoint; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - if(geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint)) return false; - else return true; - -} - -vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) const { - - vector coord(3,0.0); - vector coord_passive(3, 0.0); - - auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - for (auto iDim = 0 ; iDim < nDim ; iDim++){ - coord[iDim] = solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord(iPoint,iDim); - } - - coord_passive[0] = SU2_TYPE::GetValue(coord[0]); - coord_passive[1] = SU2_TYPE::GetValue(coord[1]); - coord_passive[2] = SU2_TYPE::GetValue(coord[2]); - - return coord_passive; -} - -vector CDriver::GetVertexNormal(unsigned short iMarker, unsigned long iVertex, bool unitNormal) const { - - su2double *Normal; - su2double Area; - vector ret_Normal(3, 0.0); - vector ret_Normal_passive(3, 0.0); - - Normal = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNormal(); - - if (!unitNormal) { - - ret_Normal_passive[0] = SU2_TYPE::GetValue(Normal[0]); - ret_Normal_passive[1] = SU2_TYPE::GetValue(Normal[1]); - if(nDim>2) ret_Normal_passive[2] = SU2_TYPE::GetValue(Normal[2]); - - return ret_Normal_passive; - } - - Area = GeometryToolbox::Norm(nDim, Normal); - - ret_Normal[0] = Normal[0]/Area; - ret_Normal[1] = Normal[1]/Area; - if(nDim>2) ret_Normal[2] = Normal[2]/Area; - - ret_Normal_passive[0] = SU2_TYPE::GetValue(ret_Normal[0]); - ret_Normal_passive[1] = SU2_TYPE::GetValue(ret_Normal[1]); - ret_Normal_passive[2] = SU2_TYPE::GetValue(ret_Normal[2]); - - return ret_Normal_passive; -} - ////////////////////////////////////////////////////////////////////////////////// -/* Functions to obtain global parameters from SU2 (time steps, delta t, ecc...) */ +/* Functions to obtain global parameters from SU2 (time steps, delta t, etc.) */ ////////////////////////////////////////////////////////////////////////////////// -unsigned long CDriver::GetnTimeIter() const { - - return config_container[ZONE_0]->GetnTime_Iter(); -} +unsigned long CDriver::GetnTimeIter() const { return config_container[ZONE_0]->GetnTime_Iter(); } -unsigned long CDriver::GetTime_Iter() const{ - - return TimeIter; -} +unsigned long CDriver::GetTime_Iter() const { return TimeIter; } passivedouble CDriver::GetUnsteady_TimeStep() const { - return SU2_TYPE::GetValue(config_container[ZONE_0]->GetTime_Step()); } -string CDriver::GetSurfaceFileName() const { - - return config_container[ZONE_0]->GetSurfCoeff_FileName(); -} +string CDriver::GetSurfaceFileName() const { return config_container[ZONE_0]->GetSurfCoeff_FileName(); } /////////////////////////////////////////////////////////////////////////////// /* Functions related to CHT solver */ /////////////////////////////////////////////////////////////////////////////// passivedouble CDriver::GetVertexTemperature(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; su2double vertexWallTemp(0.0); @@ -296,43 +184,40 @@ passivedouble CDriver::GetVertexTemperature(unsigned short iMarker, unsigned lon iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - if(geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint) && compressible){ + if (geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint) && compressible) { vertexWallTemp = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetTemperature(iPoint); } return SU2_TYPE::GetValue(vertexWallTemp); - } -void CDriver::SetVertexTemperature(unsigned short iMarker, unsigned long iVertex, passivedouble val_WallTemp){ - +void CDriver::SetVertexTemperature(unsigned short iMarker, unsigned long iVertex, passivedouble val_WallTemp) { geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryTemperature(iMarker, iVertex, val_WallTemp); } vector CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; unsigned short iDim; - su2double Prandtl_Lam = config_container[ZONE_0]->GetPrandtl_Lam(); + su2double Prandtl_Lam = config_container[ZONE_0]->GetPrandtl_Lam(); su2double Gas_Constant = config_container[ZONE_0]->GetGas_ConstantND(); su2double Gamma = config_container[ZONE_0]->GetGamma(); su2double Gamma_Minus_One = Gamma - 1.0; su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; su2double laminar_viscosity, thermal_conductivity; - vector GradT (3,0.0); - vector HeatFlux (3,0.0); - vector HeatFluxPassive (3,0.0); + vector GradT(3, 0.0); + vector HeatFlux(3, 0.0); + vector HeatFluxPassive(3, 0.0); bool compressible = (config_container[ZONE_0]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE); iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - if(compressible){ - laminar_viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); - thermal_conductivity = Cp * (laminar_viscosity/Prandtl_Lam); - for(iDim=0; iDim < nDim; iDim++){ + if (compressible) { + laminar_viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); + thermal_conductivity = Cp * (laminar_viscosity / Prandtl_Lam); + for (iDim = 0; iDim < nDim; iDim++) { GradT[iDim] = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint, 0, iDim); - HeatFlux[iDim] = -thermal_conductivity*GradT[iDim]; + HeatFlux[iDim] = -thermal_conductivity * GradT[iDim]; } } @@ -343,19 +228,18 @@ vector CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsig return HeatFluxPassive; } -passivedouble CDriver::GetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex) const{ - +passivedouble CDriver::GetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; unsigned short iDim; su2double vertexWallHeatFlux; - su2double Prandtl_Lam = config_container[ZONE_0]->GetPrandtl_Lam(); + su2double Prandtl_Lam = config_container[ZONE_0]->GetPrandtl_Lam(); su2double Gas_Constant = config_container[ZONE_0]->GetGas_ConstantND(); su2double Gamma = config_container[ZONE_0]->GetGamma(); su2double Gamma_Minus_One = Gamma - 1.0; su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; su2double Area; su2double laminar_viscosity, thermal_conductivity, dTdn; - su2double *Normal, GradT[3] = {0.0,0.0,0.0}, UnitNormal[3] = {0.0,0.0,0.0}; + su2double *Normal, GradT[3] = {0.0, 0.0, 0.0}, UnitNormal[3] = {0.0, 0.0, 0.0}; bool compressible = (config_container[ZONE_0]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE); @@ -364,37 +248,34 @@ passivedouble CDriver::GetVertexNormalHeatFlux(unsigned short iMarker, unsigned iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - if(geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint) && compressible){ + if (geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint) && compressible) { Normal = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNormal(); Area = GeometryToolbox::Norm(nDim, Normal); - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = Normal[iDim]/Area; + for (iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area; - laminar_viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); - thermal_conductivity = Cp * (laminar_viscosity/Prandtl_Lam); + laminar_viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); + thermal_conductivity = Cp * (laminar_viscosity / Prandtl_Lam); /*Compute wall heat flux (normal to the wall) based on computed temperature gradient*/ - for(iDim=0; iDim < nDim; iDim++){ + for (iDim = 0; iDim < nDim; iDim++) { GradT[iDim] = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint, 0, iDim); - dTdn += GradT[iDim]*UnitNormal[iDim]; + dTdn += GradT[iDim] * UnitNormal[iDim]; } - vertexWallHeatFlux = -thermal_conductivity*dTdn; + vertexWallHeatFlux = -thermal_conductivity * dTdn; } return SU2_TYPE::GetValue(vertexWallHeatFlux); } -void CDriver::SetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex, passivedouble val_WallHeatFlux){ - +void CDriver::SetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex, passivedouble val_WallHeatFlux) { geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryHeatFlux(iMarker, iVertex, val_WallHeatFlux); } passivedouble CDriver::GetThermalConductivity(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; - su2double Prandtl_Lam = config_container[ZONE_0]->GetPrandtl_Lam(); + su2double Prandtl_Lam = config_container[ZONE_0]->GetPrandtl_Lam(); su2double Gas_Constant = config_container[ZONE_0]->GetGas_ConstantND(); su2double Gamma = config_container[ZONE_0]->GetGamma(); su2double Gamma_Minus_One = Gamma - 1.0; @@ -402,218 +283,124 @@ passivedouble CDriver::GetThermalConductivity(unsigned short iMarker, unsigned l su2double laminar_viscosity, thermal_conductivity; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - laminar_viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); - thermal_conductivity = Cp * (laminar_viscosity/Prandtl_Lam); + laminar_viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); + thermal_conductivity = Cp * (laminar_viscosity / Prandtl_Lam); return SU2_TYPE::GetValue(thermal_conductivity); - } //////////////////////////////////////////////////////////////////////////////// /* Functions related to the management of markers */ //////////////////////////////////////////////////////////////////////////////// -vector CDriver::GetAllBoundaryMarkersTag() const { - - vector boundariesTagList; - unsigned short iMarker,nBoundariesMarkers; - string Marker_Tag; - - nBoundariesMarkers = config_container[ZONE_0]->GetnMarker_All(); - boundariesTagList.resize(nBoundariesMarkers); - - for(iMarker=0; iMarker < nBoundariesMarkers; iMarker++){ - Marker_Tag = config_container[ZONE_0]->GetMarker_All_TagBound(iMarker); - boundariesTagList[iMarker] = Marker_Tag; - } - - return boundariesTagList; -} - -vector CDriver::GetAllDeformMeshMarkersTag() const { - - vector interfaceBoundariesTagList; - unsigned short iMarker, nBoundariesMarker; - string Marker_Tag; +vector CDriver::GetCHTMarkerTags() const { + vector tags; + const auto nMarker = config_container[ZONE_0]->GetnMarker_All(); - nBoundariesMarker = config_container[ZONE_0]->GetnMarker_Deform_Mesh(); - interfaceBoundariesTagList.resize(nBoundariesMarker); - - for(iMarker=0; iMarker < nBoundariesMarker; iMarker++){ - Marker_Tag = config_container[ZONE_0]->GetMarker_Deform_Mesh_TagBound(iMarker); - interfaceBoundariesTagList[iMarker] = Marker_Tag; - } - - return interfaceBoundariesTagList; -} - -vector CDriver::GetAllCHTMarkersTag() const { - - vector CHTBoundariesTagList; - unsigned short iMarker, nBoundariesMarker; - string Marker_Tag; - - nBoundariesMarker = config_container[ZONE_0]->GetnMarker_All(); - - //The CHT markers can be identified as the markers that are customizable with a BC type HEAT_FLUX or ISOTHERMAL. - for(iMarker=0; iMarkerGetMarker_All_KindBC(iMarker) == HEAT_FLUX || config_container[ZONE_0]->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) && config_container[ZONE_0]->GetMarker_All_PyCustom(iMarker)){ - Marker_Tag = config_container[ZONE_0]->GetMarker_All_TagBound(iMarker); - CHTBoundariesTagList.push_back(Marker_Tag); + // The CHT markers can be identified as the markers that are customizable with a BC type HEAT_FLUX or ISOTHERMAL. + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { + if ((config_container[ZONE_0]->GetMarker_All_KindBC(iMarker) == HEAT_FLUX || + config_container[ZONE_0]->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) && + config_container[ZONE_0]->GetMarker_All_PyCustom(iMarker)) { + tags.push_back(config_container[ZONE_0]->GetMarker_All_TagBound(iMarker)); } } - return CHTBoundariesTagList; + return tags; } -vector CDriver::GetAllInletMarkersTag() const { - - vector BoundariesTagList; - unsigned short iMarker, nBoundariesMarker; - string Marker_Tag; +vector CDriver::GetInletMarkerTags() const { + vector tags; + const auto nMarker = config_container[ZONE_0]->GetnMarker_All(); - nBoundariesMarker = config_container[ZONE_0]->GetnMarker_All(); - - for(iMarker=0; iMarkerGetMarker_All_PyCustom(iMarker); bool isInlet = (config_container[ZONE_0]->GetMarker_All_KindBC(iMarker) == INLET_FLOW); - if(isCustomizable && isInlet) { - Marker_Tag = config_container[ZONE_0]->GetMarker_All_TagBound(iMarker); - BoundariesTagList.push_back(Marker_Tag); - } - } - - return BoundariesTagList; -} - -map CDriver::GetAllBoundaryMarkers() const { - - map allBoundariesMap; - unsigned short iMarker, nBoundaryMarkers; - string Marker_Tag; - - nBoundaryMarkers = config_container[ZONE_0]->GetnMarker_All(); - - for(iMarker=0; iMarker < nBoundaryMarkers; iMarker++){ - Marker_Tag = config_container[ZONE_0]->GetMarker_All_TagBound(iMarker); - allBoundariesMap[Marker_Tag] = iMarker; - } - - return allBoundariesMap; -} - -map CDriver::GetAllBoundaryMarkersType() const { - map allBoundariesTypeMap; - unsigned short iMarker, KindBC; - string Marker_Tag, Marker_Type; - - for(iMarker=0; iMarker < config_container[ZONE_0]->GetnMarker_All(); iMarker++){ - Marker_Tag = config_container[ZONE_0]->GetMarker_All_TagBound(iMarker); - KindBC = config_container[ZONE_0]->GetMarker_All_KindBC(iMarker); - switch(KindBC){ - case EULER_WALL: - Marker_Type = "EULER_WALL"; - break; - case FAR_FIELD: - Marker_Type = "FARFIELD"; - break; - case ISOTHERMAL: - Marker_Type = "ISOTHERMAL"; - break; - case HEAT_FLUX: - Marker_Type = "HEATFLUX"; - break; - case INLET_FLOW: - Marker_Type = "INLET_FLOW"; - break; - case OUTLET_FLOW: - Marker_Type = "OUTLET_FLOW"; - break; - case SYMMETRY_PLANE: - Marker_Type = "SYMMETRY"; - break; - case SEND_RECEIVE: - Marker_Type = "SEND_RECEIVE"; - break; - default: - Marker_Type = "UNKNOWN_TYPE"; + if (isCustomizable && isInlet) { + tags.push_back(config_container[ZONE_0]->GetMarker_All_TagBound(iMarker)); } - allBoundariesTypeMap[Marker_Tag] = Marker_Type; } - return allBoundariesTypeMap; + return tags; } -void CDriver::SetHeatSource_Position(passivedouble alpha, passivedouble pos_x, passivedouble pos_y, passivedouble pos_z){ - - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][RAD_SOL]; +void CDriver::SetHeatSource_Position(passivedouble alpha, passivedouble pos_x, passivedouble pos_y, + passivedouble pos_z) { + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][RAD_SOL]; config_container[ZONE_0]->SetHeatSource_Rot_Z(alpha); config_container[ZONE_0]->SetHeatSource_Center(pos_x, pos_y, pos_z); solver->SetVolumetricHeatSource(geometry_container[ZONE_0][INST_0][MESH_0], config_container[ZONE_0]); - } -void CDriver::SetInlet_Angle(unsigned short iMarker, passivedouble alpha){ - - su2double alpha_rad = alpha * PI_NUMBER/180.0; +void CDriver::SetInlet_Angle(unsigned short iMarker, passivedouble alpha) { + su2double alpha_rad = alpha * PI_NUMBER / 180.0; unsigned long iVertex; - for (iVertex = 0; iVertex < geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; iVertex++){ + for (iVertex = 0; iVertex < geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; iVertex++) { solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 0, cos(alpha_rad)); solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 1, sin(alpha_rad)); } - } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// -/* Functions related to simulation control, high level functions (reset convergence, set initial mesh, ecc...) */ +/* Functions related to simulation control, high level functions (reset convergence, set initial mesh, etc.) */ ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// void CDriver::ResetConvergence() { + for (auto iZone = 0u; iZone < nZone; iZone++) { + switch (main_config->GetKind_Solver()) { + case MAIN_SOLVER::EULER: + case MAIN_SOLVER::NAVIER_STOKES: + case MAIN_SOLVER::RANS: + case MAIN_SOLVER::INC_EULER: + case MAIN_SOLVER::INC_NAVIER_STOKES: + case MAIN_SOLVER::INC_RANS: + case MAIN_SOLVER::NEMO_EULER: + case MAIN_SOLVER::NEMO_NAVIER_STOKES: + integration_container[iZone][INST_0][FLOW_SOL]->SetConvergence(false); + if (config_container[iZone]->GetKind_Solver() == MAIN_SOLVER::RANS) + integration_container[iZone][INST_0][TURB_SOL]->SetConvergence(false); + if (config_container[iZone]->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) + integration_container[iZone][INST_0][TRANS_SOL]->SetConvergence(false); + break; + + case MAIN_SOLVER::FEM_ELASTICITY: + integration_container[iZone][INST_0][FEA_SOL]->SetConvergence(false); + break; + + case MAIN_SOLVER::ADJ_EULER: + case MAIN_SOLVER::ADJ_NAVIER_STOKES: + case MAIN_SOLVER::ADJ_RANS: + case MAIN_SOLVER::DISC_ADJ_EULER: + case MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES: + case MAIN_SOLVER::DISC_ADJ_RANS: + case MAIN_SOLVER::DISC_ADJ_INC_EULER: + case MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES: + case MAIN_SOLVER::DISC_ADJ_INC_RANS: + integration_container[iZone][INST_0][ADJFLOW_SOL]->SetConvergence(false); + if ((config_container[iZone]->GetKind_Solver() == MAIN_SOLVER::ADJ_RANS) || + (config_container[iZone]->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_RANS)) + integration_container[iZone][INST_0][ADJTURB_SOL]->SetConvergence(false); + break; - for(iZone = 0; iZone < nZone; iZone++) { - switch (config_container[iZone]->GetKind_Solver()) { - - case MAIN_SOLVER::EULER: case MAIN_SOLVER::NAVIER_STOKES: case MAIN_SOLVER::RANS: - case MAIN_SOLVER::INC_EULER: case MAIN_SOLVER::INC_NAVIER_STOKES: case MAIN_SOLVER::INC_RANS: - integration_container[iZone][INST_0][FLOW_SOL]->SetConvergence(false); - if (config_container[iZone]->GetKind_Solver() == MAIN_SOLVER::RANS) integration_container[iZone][INST_0][TURB_SOL]->SetConvergence(false); - if(config_container[iZone]->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) integration_container[iZone][INST_0][TRANS_SOL]->SetConvergence(false); - break; - - case MAIN_SOLVER::FEM_ELASTICITY: - integration_container[iZone][INST_0][FEA_SOL]->SetConvergence(false); - break; - - case MAIN_SOLVER::ADJ_EULER: case MAIN_SOLVER::ADJ_NAVIER_STOKES: case MAIN_SOLVER::ADJ_RANS: case MAIN_SOLVER::DISC_ADJ_EULER: case MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES: case MAIN_SOLVER::DISC_ADJ_RANS: - case MAIN_SOLVER::DISC_ADJ_INC_EULER: case MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES: case MAIN_SOLVER::DISC_ADJ_INC_RANS: - integration_container[iZone][INST_0][ADJFLOW_SOL]->SetConvergence(false); - if( (config_container[iZone]->GetKind_Solver() == MAIN_SOLVER::ADJ_RANS) || (config_container[iZone]->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_RANS) ) - integration_container[iZone][INST_0][ADJTURB_SOL]->SetConvergence(false); - break; - - default: - break; + default: + break; } } - } void CSinglezoneDriver::SetInitialMesh() { - DynamicMeshUpdate(0); SU2_OMP_PARALLEL { - // Overwrite fictious velocities - for (iMesh = 0u; iMesh <= config_container[ZONE_0]->GetnMGLevels(); iMesh++) { - SU2_OMP_FOR_STAT(roundUpDiv(geometry_container[ZONE_0][INST_0][iMesh]->GetnPoint(),omp_get_max_threads())) - for (unsigned long iPoint = 0; iPoint < geometry_container[ZONE_0][INST_0][iMesh]->GetnPoint(); iPoint++) { - - /*--- Overwrite fictitious velocities ---*/ + for (iMesh = 0u; iMesh <= main_config->GetnMGLevels(); iMesh++) { + SU2_OMP_FOR_STAT(roundUpDiv(geometry_container[ZONE_0][INST_0][iMesh]->GetnPoint(), omp_get_max_threads())) + for (auto iPoint = 0ul; iPoint < geometry_container[ZONE_0][INST_0][iMesh]->GetnPoint(); iPoint++) { + /*--- Overwrite fictitious velocities. ---*/ su2double Grid_Vel[3] = {0.0, 0.0, 0.0}; /*--- Set the grid velocity for this coarse node. ---*/ @@ -624,23 +411,21 @@ void CSinglezoneDriver::SetInitialMesh() { geometry_container[ZONE_0][INST_0][iMesh]->nodes->SetVolume_n(); geometry_container[ZONE_0][INST_0][iMesh]->nodes->SetVolume_nM1(); } - /*--- Push back the solution so that there is no fictious velocity at the next step. ---*/ + /*--- Push back the solution so that there is no fictitious velocity at the next step. ---*/ solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n(); solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n1(); } END_SU2_OMP_PARALLEL } -void CDriver::BoundaryConditionsUpdate(){ - +void CDriver::BoundaryConditionsUpdate() { int rank = MASTER_NODE; - unsigned short iZone; SU2_MPI::Comm_rank(SU2_MPI::GetComm(), &rank); - if(rank == MASTER_NODE) cout << "Updating boundary conditions." << endl; - for(iZone = 0; iZone < nZone; iZone++){ - geometry_container[iZone][INST_0][MESH_0]->UpdateCustomBoundaryConditions(geometry_container[iZone][INST_0], config_container[iZone]); + if (rank == MASTER_NODE) cout << "Updating boundary conditions." << endl; + for (auto iZone = 0u; iZone < nZone; iZone++) { + geometry_container[iZone][INST_0][MESH_0]->UpdateCustomBoundaryConditions(geometry_container[iZone][INST_0],config_container[iZone]); } } @@ -648,29 +433,26 @@ void CDriver::BoundaryConditionsUpdate(){ /* Functions related to finite elements */ //////////////////////////////////////////////////////////////////////////////// -void CDriver::SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passivedouble LoadX, - passivedouble LoadY, passivedouble LoadZ) { - +void CDriver::SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passivedouble LoadX, passivedouble LoadY, + passivedouble LoadZ) { unsigned long iPoint; - su2double NodalForce[3] = {0.0,0.0,0.0}; + su2double NodalForce[3] = {0.0, 0.0, 0.0}; NodalForce[0] = LoadX; NodalForce[1] = LoadY; NodalForce[2] = LoadZ; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint,NodalForce); - + solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint, NodalForce); } vector CDriver::GetFEA_Displacements(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; vector Displacements(3, 0.0); vector Displacements_passive(3, 0.0); iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; Displacements[0] = solver->GetNodes()->GetSolution(iPoint, 0); Displacements[1] = solver->GetNodes()->GetSolution(iPoint, 1); @@ -686,18 +468,16 @@ vector CDriver::GetFEA_Displacements(unsigned short iMarker, unsi return Displacements_passive; } - vector CDriver::GetFEA_Velocity(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; vector Velocity(3, 0.0); - vector Velocity_passive(3,0.0); + vector Velocity_passive(3, 0.0); iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; - if (config_container[ZONE_0]->GetDynamic_Analysis() == DYNAMIC){ + if (config_container[ZONE_0]->GetDynamic_Analysis() == DYNAMIC) { Velocity[0] = solver->GetNodes()->GetSolution_Vel(iPoint, 0); Velocity[1] = solver->GetNodes()->GetSolution_Vel(iPoint, 1); if (geometry->GetnDim() == 3) @@ -714,16 +494,15 @@ vector CDriver::GetFEA_Velocity(unsigned short iMarker, unsigned } vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; vector Velocity_n(3, 0.0); vector Velocity_n_passive(3, 0.0); iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; - if (config_container[ZONE_0]->GetDynamic_Analysis() == DYNAMIC){ + if (config_container[ZONE_0]->GetDynamic_Analysis() == DYNAMIC) { Velocity_n[0] = solver->GetNodes()->GetSolution_Vel_time_n(iPoint, 0); Velocity_n[1] = solver->GetNodes()->GetSolution_Vel_time_n(iPoint, 1); if (geometry->GetnDim() == 3) @@ -737,7 +516,6 @@ vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigne Velocity_n_passive[2] = SU2_TYPE::GetValue(Velocity_n[2]); return Velocity_n_passive; - } //////////////////////////////////////////////////////////////////////////////// @@ -745,14 +523,13 @@ vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigne //////////////////////////////////////////////////////////////////////////////// vector CDriver::GetMeshDisp_Sensitivity(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; vector Disp_Sens(3, 0.0); vector Disp_Sens_passive(3, 0.0); iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][ADJMESH_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][ADJMESH_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; Disp_Sens[0] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 0); Disp_Sens[1] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 1); @@ -766,18 +543,16 @@ vector CDriver::GetMeshDisp_Sensitivity(unsigned short iMarker, u Disp_Sens_passive[2] = SU2_TYPE::GetValue(Disp_Sens[2]); return Disp_Sens_passive; - } vector CDriver::GetFlowLoad_Sensitivity(unsigned short iMarker, unsigned long iVertex) const { - unsigned long iPoint; vector FlowLoad_Sens(3, 0.0); vector FlowLoad_Sens_passive(3, 0.0); iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][ADJFEA_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][ADJFEA_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; FlowLoad_Sens[0] = solver->GetNodes()->GetFlowTractionSensitivity(iPoint, 0); FlowLoad_Sens[1] = solver->GetNodes()->GetFlowTractionSensitivity(iPoint, 1); @@ -791,91 +566,52 @@ vector CDriver::GetFlowLoad_Sensitivity(unsigned short iMarker, u FlowLoad_Sens_passive[2] = SU2_TYPE::GetValue(FlowLoad_Sens[2]); return FlowLoad_Sens_passive; - } void CDriver::SetFlowLoad_Adjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, passivedouble val_AdjointY, passivedouble val_AdjointZ) { - - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; solver->StoreVertexTractionsAdjoint(iMarker, iVertex, 0, val_AdjointX); solver->StoreVertexTractionsAdjoint(iMarker, iVertex, 1, val_AdjointY); - if (geometry->GetnDim() == 3) - solver->StoreVertexTractionsAdjoint(iMarker, iVertex, 2, val_AdjointZ); - + if (geometry->GetnDim() == 3) solver->StoreVertexTractionsAdjoint(iMarker, iVertex, 2, val_AdjointZ); } void CDriver::SetSourceTerm_DispAdjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, passivedouble val_AdjointY, passivedouble val_AdjointZ) { - unsigned long iPoint; - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][ADJFEA_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][ADJFEA_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); solver->GetNodes()->SetSourceTerm_DispAdjoint(iPoint, 0, val_AdjointX); solver->GetNodes()->SetSourceTerm_DispAdjoint(iPoint, 1, val_AdjointY); - if (geometry->GetnDim() == 3) - solver->GetNodes()->SetSourceTerm_DispAdjoint(iPoint, 2, val_AdjointZ); - + if (geometry->GetnDim() == 3) solver->GetNodes()->SetSourceTerm_DispAdjoint(iPoint, 2, val_AdjointZ); } void CDriver::SetSourceTerm_VelAdjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, - passivedouble val_AdjointY, passivedouble val_AdjointZ) { - - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][ADJFEA_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + passivedouble val_AdjointY, passivedouble val_AdjointZ) { + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][ADJFEA_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; const auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); solver->GetNodes()->SetSourceTerm_VelAdjoint(iPoint, 0, val_AdjointX); solver->GetNodes()->SetSourceTerm_VelAdjoint(iPoint, 1, val_AdjointY); - if (geometry->GetnDim() == 3) - solver->GetNodes()->SetSourceTerm_VelAdjoint(iPoint, 2, val_AdjointZ); - + if (geometry->GetnDim() == 3) solver->GetNodes()->SetSourceTerm_VelAdjoint(iPoint, 2, val_AdjointZ); } //////////////////////////////////////////////////////////////////////////////// -/* Functions related to mesh deformation */ -//////////////////////////////////////////////////////////////////////////////// - -void CDriver::SetMeshDisplacement(unsigned short iMarker, unsigned long iVertex, passivedouble DispX, passivedouble DispY, passivedouble DispZ) { - - unsigned long iPoint; - su2double MeshDispl[3] = {0.0,0.0,0.0}; - - MeshDispl[0] = DispX; - MeshDispl[1] = DispY; - MeshDispl[2] = DispZ; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint,MeshDispl); - -} - -void CDriver::CommunicateMeshDisplacement(void) { - - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->InitiateComms(geometry_container[ZONE_0][INST_0][MESH_0], - config_container[ZONE_0], MESH_DISPLACEMENTS); - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->CompleteComms(geometry_container[ZONE_0][INST_0][MESH_0], - config_container[ZONE_0], MESH_DISPLACEMENTS); - -} - -//////////////////////////////////////////////////////////////////////////////// -/* Functions related to flow loads */ +/* Functions related to flow loads */ //////////////////////////////////////////////////////////////////////////////// vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long iVertex) const { - vector FlowLoad(3, 0.0); vector FlowLoad_passive(3, 0.0); - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; if (config_container[ZONE_0]->GetSolid_Wall(iMarker)) { FlowLoad[0] = solver->GetVertexTractions(iMarker, iVertex, 0); @@ -891,5 +627,4 @@ vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long FlowLoad_passive[2] = SU2_TYPE::GetValue(FlowLoad[2]); return FlowLoad_passive; - } diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 77470e620fb6..57b4b25edb7e 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -415,7 +415,7 @@ void CMeshSolver::SetWallDistance(CGeometry *geometry, CConfig *config) { END_SU2_OMP_PARALLEL } -void CMeshSolver::SetMesh_Stiffness(CGeometry **geometry, CNumerics **numerics, CConfig *config){ +void CMeshSolver::SetMesh_Stiffness(CNumerics **numerics, CConfig *config){ if (stiffness_set) return; diff --git a/SU2_DEF/include/SU2_DEF.hpp b/SU2_DEF/include/SU2_DEF.hpp index 395fcc804e50..a0a9877b0897 100644 --- a/SU2_DEF/include/SU2_DEF.hpp +++ b/SU2_DEF/include/SU2_DEF.hpp @@ -1,7 +1,6 @@ /*! * \file SU2_DEF.hpp * \brief Headers of the main subroutines of the code SU2_DEF. - * The subroutines and functions are in the SU2_DEF.cpp file. * \author F. Palacios, T. Economon * \version 7.5.1 "Blackbird" * @@ -26,20 +25,16 @@ * License along with SU2. If not, see . */ - #pragma once -#include "../../Common/include/parallelization/mpi_structure.hpp" -#include "../../Common/include/parallelization/omp_structure.hpp" - +#include #include -#include #include -#include +#include -#include "../../SU2_CFD/include/solvers/CSolver.hpp" -#include "../../SU2_CFD/include/output/CMeshOutput.hpp" -#include "../../Common/include/geometry/CPhysicalGeometry.hpp" #include "../../Common/include/CConfig.hpp" +#include "../../Common/include/parallelization/mpi_structure.hpp" +#include "../../Common/include/parallelization/omp_structure.hpp" +#include "drivers/CDeformationDriver.hpp" using namespace std; diff --git a/SU2_DEF/include/drivers/CDeformationDriver.hpp b/SU2_DEF/include/drivers/CDeformationDriver.hpp new file mode 100644 index 000000000000..a562a07c1a5b --- /dev/null +++ b/SU2_DEF/include/drivers/CDeformationDriver.hpp @@ -0,0 +1,116 @@ +/*! + * \file CDeformationDriver.hpp + * \brief Headers of the main subroutines for driving the mesh deformation. + * \author A. Gastaldi, H. Patel + * \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 . + */ + +#pragma once + +#include "../../../Common/include/CConfig.hpp" +#include "../../../Common/include/geometry/CGeometry.hpp" +#include "../../../Common/include/grid_movement/CSurfaceMovement.hpp" +#include "../../../Common/include/grid_movement/CVolumetricMovement.hpp" +#include "../../../Common/include/parallelization/mpi_structure.hpp" +#include "../../../SU2_CFD/include/drivers/CDriverBase.hpp" +#include "../../../SU2_CFD/include/numerics/CNumerics.hpp" +#include "../../../SU2_CFD/include/output/COutput.hpp" + +class CDeformationDriver : public CDriverBase { + protected: + bool haveSurfaceDeformation = false; // flag used to determine whether surface deformation is available for output + + public: + /*! + * \brief Constructor of the class. + * \param[in] confFile - Configuration file name. + * \param[in] MPICommunicator - MPI communicator for SU2. + */ + CDeformationDriver(char* confFile, SU2_Comm MPICommunicator); + + /*! + * \brief Destructor of the class. + */ + ~CDeformationDriver(void); + + /*! + * \brief Preprocess the driver data. + */ + void Preprocess(); + + /*! + * \brief Launch the driver computation. + */ + void Run(); + + /*! + * \brief Output the mesh. + */ + void Output(); + + /*! + * \brief Deallocation routine. + */ + void Postprocessing(); + + /*! + * \brief Communicate boundary mesh displacements. + */ + void CommunicateMeshDisplacements(void); + + protected: + /*! + * \brief Read in the config and mesh files. + */ + void Input_Preprocessing(); + + /*! + * \brief Construction of the edge-based data structure. + */ + void Geometrical_Preprocessing(); + + /*! + * \brief Preprocess the output container. + */ + void Output_Preprocessing(); + + /*! + * \brief Preprocess the mesh solver container. + */ + void Solver_Preprocessing(); + + /*! + * \brief Preprocess the numerics container. + */ + void Numerics_Preprocessing(); + + /*! + * \brief Mesh deformation based on linear elasticity solver (CMeshSolver). + */ + void Update(); + + /*! + * \brief Mesh deformation based on legacy implementation. + */ + void Update_Legacy(); +}; diff --git a/SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp b/SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp new file mode 100644 index 000000000000..6b28f072b982 --- /dev/null +++ b/SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp @@ -0,0 +1,137 @@ +/*! + * \file CDiscAdjDeformationDriver.cpp + * \brief Headers of the main subroutines for driving the projection of sensitivities. + * \author T. Economon, H. Kline, R. Sanchez, A. Gastaldi, H. Patel + * \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 . + */ + +#pragma once + +#include "../../../Common/include/fem/fem_geometry_structure.hpp" +#include "../../../Common/include/parallelization/mpi_structure.hpp" +#include "../../../Common/include/parallelization/omp_structure.hpp" +#include "../../../SU2_CFD/include/drivers/CDriverBase.hpp" + +class CDiscAdjDeformationDriver : public CDriverBase { +protected: + su2double** Gradient; + ofstream Gradient_file; + +public: + /*! + * \brief Constructor of the class. + * \param[in] confFile - Configuration file name. + * \param[in] MPICommunicator - MPI communicator for SU2. + */ + CDiscAdjDeformationDriver(char* confFile, SU2_Comm MPICommunicator); + + /*! + * \brief Destructor of the class. + */ + ~CDiscAdjDeformationDriver(void); + + /*! + * \brief Preprocess the driver data (includes solution allocation and initialization). + */ + void Preprocess(); + + /*! + * \brief Launch the driver computation. + */ + void Run(); + + /*! + * \brief Deallocation routine. + */ + void Postprocessing(); + +protected: + /*! + * \brief Read in the config and mesh files. + */ + void Input_Preprocessing(); + + /*! + * \brief Construction of the edge-based data structure. + */ + void Geometrical_Preprocessing(); + + /*! + * \brief Preprocess the output container. + */ + void Output_Preprocessing(); + + /*! + * \brief Projection of the surface sensitivity using finite differences (FD). + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] surface_movement - Surface movement class of the problem. + * \param[in] Gradient_file - Output file to store the gradient data. + */ + void SetProjection_FD(CGeometry* geometry, CConfig* config, CSurfaceMovement* surface_movement, su2double** Gradient); + + /*! + * \brief Projection of the surface sensitivity using algorithmic differentiation (AD). + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] surface_movement - Surface movement class of the problem. + * \param[in] Gradient_file - Output file to store the gradient data. + */ + void SetProjection_AD(CGeometry* geometry, CConfig* config, CSurfaceMovement* surface_movement, su2double** Gradient); + + /*! + * \brief Prints the gradient information to a file. + * \param[in] Gradient - The gradient data. + * \param[in] config - Definition of the particular problem. + * \param[in] Gradient_file - Output file to store the gradient data. + */ + void OutputGradient(su2double** Gradient, CConfig* config, ofstream& Gradient_file); + + /*! + * \brief Write the sensitivity (including mesh sensitivity) computed with the discrete adjoint method + * on the surface and in the volume to a file. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] val_nZone - Number of Zones. + */ + void SetSensitivity_Files(CGeometry**** geometry, CConfig** config, unsigned short val_nZone); + + /*! + * \brief Treatment of derivatives with the Sobolev smoothing solver. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] grid_movement - Volumetric movement class of the problem. + */ + void DerivativeTreatment_MeshSensitivity(CGeometry* geometry, CConfig* config, CVolumetricMovement* grid_movement); + + /*! + * \brief Treatment of derivatives with the Sobolev smoothing solver. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] grid_movement - Volumetric movement class of the problem. + * \param[in] surface_movement - Surface movement class of the problem. + * \param[in] Gradient - Output array to store the gradient data. + */ + void DerivativeTreatment_Gradient(CGeometry* geometry, CConfig* config, CVolumetricMovement* grid_movement, + CSurfaceMovement* surface_movement, su2double** Gradient); +}; diff --git a/SU2_DEF/src/SU2_DEF.cpp b/SU2_DEF/src/SU2_DEF.cpp index 77bee3b5fab8..49cef08b52a1 100644 --- a/SU2_DEF/src/SU2_DEF.cpp +++ b/SU2_DEF/src/SU2_DEF.cpp @@ -25,17 +25,17 @@ * License along with SU2. If not, see . */ - #include "../include/SU2_DEF.hpp" + using namespace std; -int main(int argc, char *argv[]) { +int main(int argc, char* argv[]) { - unsigned short iZone, nZone = SINGLE_ZONE; - su2double StartTime = 0.0, StopTime = 0.0, UsedTime = 0.0; char config_file_name[MAX_STRING_SIZE]; - int rank, size; - string str; + + /*--- Create a pointer to the main SU2_DEF Driver ---*/ + + CDeformationDriver* driver = nullptr; /*--- MPI initialization ---*/ @@ -45,461 +45,38 @@ int main(int argc, char *argv[]) { #else SU2_MPI::Init(&argc, &argv); #endif - SU2_MPI::Comm MPICommunicator = SU2_MPI::GetComm(); - - rank = SU2_MPI::GetRank(); - size = SU2_MPI::GetSize(); - - /*--- Pointer to different structures that will be used throughout - the entire code ---*/ - - CConfig **config_container = nullptr; - CGeometry **geometry_container = nullptr; - CSurfaceMovement **surface_movement = nullptr; - CVolumetricMovement **grid_movement = nullptr; - COutput **output = nullptr; - CConfig *driver_config = nullptr; + SU2_MPI::Comm comm = SU2_MPI::GetComm(); /*--- Load in the number of zones and spatial dimensions in the mesh file - (if no config file is specified, default.cfg is used) ---*/ - - if (argc == 2) { strcpy(config_file_name, argv[1]); } - else { strcpy(config_file_name, "default.cfg"); } - - /*--- Read the name and format of the input mesh file to get from the mesh - file the number of zones and dimensions from the numerical grid (required - for variables allocation) ---*/ - - CConfig *config = nullptr; - config = new CConfig(config_file_name, SU2_COMPONENT::SU2_DEF); + (if no config file is specified, default.cfg is used). ---*/ - nZone = config->GetnZone(); - - /*--- Definition of the containers per zones ---*/ - - config_container = new CConfig*[nZone]; - geometry_container = new CGeometry*[nZone]; - surface_movement = new CSurfaceMovement*[nZone]; - grid_movement = new CVolumetricMovement*[nZone]; - output = new COutput*[nZone]; - - driver_config = nullptr; - - for (iZone = 0; iZone < nZone; iZone++) { - config_container[iZone] = nullptr; - geometry_container[iZone] = nullptr; - surface_movement[iZone] = nullptr; - grid_movement[iZone] = nullptr; - output[iZone] = nullptr; + if (argc == 2) { + strcpy(config_file_name, argv[1]); + } else { + strcpy(config_file_name, "default.cfg"); } - /*--- Initialize the configuration of the driver ---*/ - driver_config = new CConfig(config_file_name, SU2_COMPONENT::SU2_DEF, false); - - /*--- Initialize a char to store the zone filename ---*/ - char zone_file_name[MAX_STRING_SIZE]; - - /*--- Loop over all zones to initialize the various classes. In most - cases, nZone is equal to one. This represents the solution of a partial - differential equation on a single block, unstructured mesh. ---*/ - - for (iZone = 0; iZone < nZone; iZone++) { - - /*--- Definition of the configuration option class for all zones. In this - constructor, the input configuration file is parsed and all options are - read and stored. ---*/ - - if (driver_config->GetnConfigFiles() > 0){ - strcpy(zone_file_name, driver_config->GetConfigFilename(iZone).c_str()); - config_container[iZone] = new CConfig(driver_config, zone_file_name, SU2_COMPONENT::SU2_DEF, iZone, nZone, true); - } - else{ - config_container[iZone] = new CConfig(driver_config, config_file_name, SU2_COMPONENT::SU2_DEF, iZone, nZone, true); - } - config_container[iZone]->SetMPICommunicator(MPICommunicator); - } - - /*--- Set the multizone part of the problem. ---*/ - if (driver_config->GetMultizone_Problem()){ - for (iZone = 0; iZone < nZone; iZone++) { - /*--- Set the interface markers for multizone ---*/ - config_container[iZone]->SetMultizone(driver_config, config_container); - } - } - - for (iZone = 0; iZone < nZone; iZone++) { - - /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/ - - CGeometry *geometry_aux = nullptr; - - /*--- All ranks process the grid and call ParMETIS for partitioning ---*/ - - geometry_aux = new CPhysicalGeometry(config_container[iZone], iZone, nZone); - - /*--- Color the initial grid and set the send-receive domains (ParMETIS) ---*/ - - geometry_aux->SetColorGrid_Parallel(config_container[iZone]); - - /*--- Build the grid data structures using the ParMETIS coloring. ---*/ - - geometry_container[iZone] = new CPhysicalGeometry(geometry_aux, config_container[iZone]); - - /*--- Deallocate the memory of geometry_aux ---*/ - - delete geometry_aux; - - /*--- Add the Send/Receive boundaries ---*/ - - geometry_container[iZone]->SetSendReceive(config_container[iZone]); - - /*--- Add the Send/Receive boundaries ---*/ - - geometry_container[iZone]->SetBoundaries(config_container[iZone]); - - } - - /*--- Set up a timer for performance benchmarking (preprocessing time is included) ---*/ - - StartTime = SU2_MPI::Wtime(); - - for (iZone = 0; iZone < nZone; iZone++) { - - /*--- Computational grid preprocesing ---*/ - - if (rank == MASTER_NODE) cout << endl << "----------------------- Preprocessing computations ----------------------" << endl; - - /*--- Compute elements surrounding points, points surrounding points ---*/ - - if (rank == MASTER_NODE) cout << "Setting local point connectivity." <SetPoint_Connectivity(); - - /*--- Check the orientation before computing geometrical quantities ---*/ - - geometry_container[iZone]->SetBoundVolume(); - if (config_container[iZone]->GetReorientElements()) { - if (rank == MASTER_NODE) cout << "Checking the numerical grid orientation of the interior elements." <Check_IntElem_Orientation(config_container[iZone]); - geometry_container[iZone]->Check_BoundElem_Orientation(config_container[iZone]); - } - - /*--- Create the edge structure ---*/ - - if (rank == MASTER_NODE) cout << "Identify edges and vertices." <SetEdges(); - geometry_container[iZone]->SetVertex(config_container[iZone]); - - if (config_container[iZone]->GetDesign_Variable(0) != NO_DEFORMATION) { - - /*--- Create the dual control volume structures ---*/ - - if (rank == MASTER_NODE) cout << "Setting the bound control volume structure." << endl; - geometry_container[iZone]->SetControlVolume(config_container[iZone], ALLOCATE); - geometry_container[iZone]->SetBoundControlVolume(config_container[iZone], ALLOCATE); - - } - /*--- Create the point-to-point MPI communication structures. ---*/ - - geometry_container[iZone]->PreprocessP2PComms(geometry_container[iZone], config_container[iZone]); - - /*--- Allocate the mesh output ---*/ - - output[iZone] = new CMeshOutput(config_container[iZone], geometry_container[iZone]->GetnDim()); - - /*--- Preprocess the volume output ---*/ - - output[iZone]->PreprocessVolumeOutput(config_container[iZone]); - - /*--- Preprocess history --- */ - - output[iZone]->PreprocessHistoryOutput(config_container[iZone], false); - - } - - - /*--- Surface grid deformation using design variables ---*/ - - for (iZone = 0; iZone < nZone; iZone++){ - - if (config_container[iZone]->GetDesign_Variable(0) != NO_DEFORMATION) { - - /*--- Definition of the Class for grid movement ---*/ - grid_movement[iZone] = new CVolumetricMovement(geometry_container[iZone], config_container[iZone]); - - /*--- Save original coordinates to be reused in convexity checking procedure ---*/ - auto OriginalCoordinates = geometry_container[iZone]->nodes->GetCoord(); - - /*--- First check for volumetric grid deformation/transformations ---*/ - - if (config_container[iZone]->GetDesign_Variable(0) == SCALE_GRID) { - - if (rank == MASTER_NODE) - cout << endl << "--------------------- Volumetric grid scaling (ZONE " << iZone <<") ------------------" << endl; - grid_movement[iZone]->SetVolume_Scaling(geometry_container[iZone], config_container[iZone], false); - - } else if (config_container[iZone]->GetDesign_Variable(0) == TRANSLATE_GRID) { + /*--- Initialize the mesh deformation driver. ---*/ - if (rank == MASTER_NODE) - cout << endl << "------------------- Volumetric grid translation (ZONE " << iZone <<") ----------------" << endl; - grid_movement[iZone]->SetVolume_Translation(geometry_container[iZone], config_container[iZone], false); + driver = new CDeformationDriver(config_file_name, comm); - } else if (config_container[iZone]->GetDesign_Variable(0) == ROTATE_GRID) { + /*--- Preprocess the solver data. ---*/ - if (rank == MASTER_NODE) - cout << endl << "--------------------- Volumetric grid rotation (ZONE " << iZone <<") -----------------" << endl; - grid_movement[iZone]->SetVolume_Rotation(geometry_container[iZone], config_container[iZone], false); + driver->Preprocess(); - } else { + /*--- Launch the main external loop of the solver. ---*/ - /*--- If no volume-type deformations are requested, then this is a - surface-based deformation or FFD set up. ---*/ + driver->Run(); - if (rank == MASTER_NODE) - cout << endl << "--------------------- Surface grid deformation (ZONE " << iZone <<") -----------------" << endl; + /*--- Postprocess all the containers, close history file, and exit SU2. ---*/ - /*--- Definition and initialization of the surface deformation class ---*/ + driver->Postprocessing(); - surface_movement[iZone] = new CSurfaceMovement(); + delete driver; - /*--- Copy coordinates to the surface structure ---*/ + /*--- Finalize MPI parallelization. ---*/ - surface_movement[iZone]->CopyBoundary(geometry_container[iZone], config_container[iZone]); - - /*--- Surface grid deformation ---*/ - - if (rank == MASTER_NODE) cout << "Performing the deformation of the surface grid." << endl; - auto TotalDeformation = surface_movement[iZone]->SetSurface_Deformation(geometry_container[iZone], config_container[iZone]); - - if (config_container[iZone]->GetDesign_Variable(0) != FFD_SETTING) { - - if (rank == MASTER_NODE) - cout << endl << "------------------- Volumetric grid deformation (ZONE " << iZone <<") ----------------" << endl; - - if (rank == MASTER_NODE) - cout << "Performing the deformation of the volumetric grid." << endl; - grid_movement[iZone]->SetVolume_Deformation(geometry_container[iZone], config_container[iZone], false); - - /*--- Get parameters for convexity check ---*/ - bool ConvexityCheck; - unsigned short ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth; - - tie(ConvexityCheck, ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth) = config_container[iZone]->GetConvexityCheck(); - - /*--- Recursively change deformations if there are nonconvex elements. ---*/ - - if (ConvexityCheck && geometry_container[iZone]->GetnNonconvexElements() > 0) { - if (rank == MASTER_NODE) { - cout << "Nonconvex elements present after deformation. " << endl; - cout << "Recursively lowering deformation magnitude." << endl; - } - - /*--- Load initial deformation values ---*/ - auto InitialDeformation = TotalDeformation; - - unsigned short ConvexityCheckIter, RecursionDepth = 0; - su2double DeformationFactor = 1.0, DeformationDifference = 1.0; - for (ConvexityCheckIter = 1; ConvexityCheckIter <= ConvexityCheck_MaxIter; ConvexityCheckIter++) { - - /*--- Recursively change deformation magnitude: - decrease if there are nonconvex elements, increase otherwise ---*/ - DeformationDifference /= 2.0; - - if (geometry_container[iZone]->GetnNonconvexElements() > 0) { - DeformationFactor -= DeformationDifference; - } else { - RecursionDepth += 1; - - if (RecursionDepth == ConvexityCheck_MaxDepth) { - if (rank == MASTER_NODE) { - cout << "Maximum recursion depth reached." << endl; - cout << "Remaining amount of original deformation: "; - cout << DeformationFactor*100.0 << " percent. " << endl; - } - break; - } - - DeformationFactor += DeformationDifference; - } - - /*--- Load mesh to start every iteration with an undeformed grid ---*/ - for (auto iPoint = 0ul; iPoint < OriginalCoordinates.rows(); iPoint++) { - for (auto iDim = 0ul; iDim < OriginalCoordinates.cols(); iDim++) { - geometry_container[iZone]->nodes->SetCoord(iPoint, iDim, OriginalCoordinates(iPoint,iDim)); - } - } - - /*--- Set deformation magnitude as percentage of initial deformation ---*/ - for (auto iDV = 0u; iDV < config->GetnDV(); iDV++) { - for (auto iDV_Value = 0u; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) { - config_container[iZone]->SetDV_Value(iDV, iDV_Value, InitialDeformation[iDV][iDV_Value]*DeformationFactor); - } - } - - /*--- Surface grid deformation ---*/ - if (rank == MASTER_NODE) cout << "Performing the deformation of the surface grid." << endl; - - TotalDeformation = surface_movement[iZone]->SetSurface_Deformation(geometry_container[iZone], config_container[iZone]); - - if (rank == MASTER_NODE) - cout << endl << "------------------- Volumetric grid deformation (ZONE " << iZone <<") ----------------" << endl; - - if (rank == MASTER_NODE) - cout << "Performing the deformation of the volumetric grid." << endl; - grid_movement[iZone]->SetVolume_Deformation(geometry_container[iZone], config_container[iZone], false); - - if (rank == MASTER_NODE) { - cout << "Number of nonconvex elements for iteration " << ConvexityCheckIter << ": "; - cout << geometry_container[iZone]->GetnNonconvexElements() << endl; - cout << "Remaining amount of original deformation: "; - cout << DeformationFactor*100.0 << " percent. " << endl; - } - - } - - } - - } - - } - - } - - } - - /*--- Computational grid preprocesing ---*/ - - if (rank == MASTER_NODE) cout << endl << "----------------------- Write deformed grid files -----------------------" << endl; - - /*--- Output deformed grid for visualization, if requested (surface and volumetric), in parallel - requires to move all the data to the master node---*/ - - for (iZone = 0; iZone < nZone; iZone++){ - - /*--- Compute Mesh Quality if requested. Necessary geometry preprocessing re-done beforehand. ---*/ - - if (config_container[iZone]->GetWrt_MeshQuality() && !config->GetStructuralProblem()) { - - if (rank == MASTER_NODE) cout << "Recompute geometry properties necessary to evaluate mesh quality statistics.\n"; - - geometry_container[iZone]->SetPoint_Connectivity(); - geometry_container[iZone]->SetBoundVolume(); - geometry_container[iZone]->SetEdges(); - geometry_container[iZone]->SetVertex(config_container[iZone]); - geometry_container[iZone]->SetControlVolume(config_container[iZone], ALLOCATE); - geometry_container[iZone]->SetBoundControlVolume(config_container[iZone], ALLOCATE); - - if (rank == MASTER_NODE) cout << "Computing mesh quality statistics for the dual control volumes.\n"; - geometry_container[iZone]->ComputeMeshQualityStatistics(config_container[iZone]); - }// Mesh Quality Output - - /*--- Load the data --- */ - - output[iZone]->Load_Data(geometry_container[iZone], config_container[iZone], nullptr); - - output[iZone]->WriteToFile(config_container[iZone], geometry_container[iZone], OUTPUT_TYPE::MESH, config->GetMesh_Out_FileName()); - - /*--- Set the file names for the visualization files ---*/ - - output[iZone]->SetVolume_Filename("volume_deformed"); - output[iZone]->SetSurface_Filename("surface_deformed"); - - for (unsigned short iFile = 0; iFile < config_container[iZone]->GetnVolumeOutputFiles(); iFile++){ - auto FileFormat = config_container[iZone]->GetVolumeOutputFiles(); - if (FileFormat[iFile] != OUTPUT_TYPE::RESTART_ASCII && FileFormat[iFile] != OUTPUT_TYPE::RESTART_BINARY) - output[iZone]->WriteToFile(config_container[iZone], geometry_container[iZone], FileFormat[iFile]); - } - } - - - if ((config_container[ZONE_0]->GetDesign_Variable(0) != NO_DEFORMATION) && - (config_container[ZONE_0]->GetDesign_Variable(0) != SCALE_GRID) && - (config_container[ZONE_0]->GetDesign_Variable(0) != TRANSLATE_GRID) && - (config_container[ZONE_0]->GetDesign_Variable(0) != ROTATE_GRID)) { - - /*--- Write the the free-form deformation boxes after deformation. ---*/ - - if (rank == MASTER_NODE) cout << "Adding any FFD information to the SU2 file." << endl; - - surface_movement[ZONE_0]->WriteFFDInfo(surface_movement, geometry_container, config_container); - - } - - delete config; - config = nullptr; - if (rank == MASTER_NODE) - cout << endl <<"------------------------- Solver Postprocessing -------------------------" << endl; - - if (geometry_container != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (geometry_container[iZone] != nullptr) { - delete geometry_container[iZone]; - } - } - delete [] geometry_container; - } - if (rank == MASTER_NODE) cout << "Deleted CGeometry container." << endl; - - if (surface_movement != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (surface_movement[iZone] != nullptr) { - delete surface_movement[iZone]; - } - } - delete [] surface_movement; - } - if (rank == MASTER_NODE) cout << "Deleted CSurfaceMovement class." << endl; - - if (grid_movement != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (grid_movement[iZone] != nullptr) { - delete grid_movement[iZone]; - } - } - delete [] grid_movement; - } - if (rank == MASTER_NODE) cout << "Deleted CVolumetricMovement class." << endl; - - if (config_container != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (config_container[iZone] != nullptr) { - delete config_container[iZone]; - } - } - delete [] config_container; - } - if (output != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (output[iZone] != nullptr) { - delete output[iZone]; - } - } - delete [] output; - } - if (rank == MASTER_NODE) cout << "Deleted CConfig container." << endl; - - if (rank == MASTER_NODE) cout << "Deleted COutput class." << endl; - - /*--- Synchronization point after a single solver iteration. Compute the - wall clock time required. ---*/ - - StopTime = SU2_MPI::Wtime(); - - /*--- Compute/print the total time for performance benchmarking. ---*/ - - UsedTime = StopTime-StartTime; - if (rank == MASTER_NODE) { - cout << "\nCompleted in " << fixed << UsedTime << " seconds on "<< size; - if (size == 1) cout << " core." << endl; else cout << " cores." << endl; - } - - /*--- Exit the solver cleanly ---*/ - - if (rank == MASTER_NODE) - cout << endl << "------------------------- Exit Success (SU2_DEF) ------------------------" << endl << endl; - - /*--- Finalize MPI parallelization ---*/ SU2_MPI::Finalize(); return EXIT_SUCCESS; - } diff --git a/SU2_DEF/src/drivers/CDeformationDriver.cpp b/SU2_DEF/src/drivers/CDeformationDriver.cpp new file mode 100644 index 000000000000..856450217c25 --- /dev/null +++ b/SU2_DEF/src/drivers/CDeformationDriver.cpp @@ -0,0 +1,615 @@ +/*! + * \file CDeformationDriver.cpp + * \brief Main subroutines for driving the mesh deformation. + * \author A. Gastaldi, H. Patel + * \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 . + */ + +#include "../../include/drivers/CDeformationDriver.hpp" + +#include "../../../Common/include/geometry/CPhysicalGeometry.hpp" +#include "../../../SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp" +#include "../../../SU2_CFD/include/output/CMeshOutput.hpp" +#include "../../../SU2_CFD/include/solvers/CMeshSolver.hpp" + +using namespace std; + +CDeformationDriver::CDeformationDriver(char* confFile, SU2_Comm MPICommunicator) + : CDriverBase(confFile, 1, MPICommunicator) { + +/*--- Initialize MeDiPack (must also be here to initialize it from Python) ---*/ +#ifdef HAVE_MPI +#if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) + SU2_MPI::Init_AMPI(); +#endif +#endif + + SU2_MPI::SetComm(MPICommunicator); + + rank = SU2_MPI::GetRank(); + size = SU2_MPI::GetSize(); + + /*--- Initialize containers. --- */ + + SetContainers_Null(); + + /*--- Preprocessing of the config files. ---*/ + + Input_Preprocessing(); + + /*--- Set up a timer for performance benchmarking. ---*/ + + StartTime = SU2_MPI::Wtime(); + + /*--- Preprocessing of the geometry for all zones. ---*/ + + Geometrical_Preprocessing(); + + /*--- Preprocessing of the output for all zones. ---*/ + + Output_Preprocessing(); + + if (driver_config->GetDeform_Mesh()) { + /*--- Preprocessing of the mesh solver for all zones. ---*/ + + Solver_Preprocessing(); + + /*--- Preprocessing of the mesh solver for all zones. ---*/ + + Numerics_Preprocessing(); + + } + + /*--- Preprocessing time is reported now, but not included in the next compute portion. ---*/ + + StopTime = SU2_MPI::Wtime(); + + /*--- Compute the total time for performance benchmarking. ---*/ + + UsedTime = StopTime - StartTime; + UsedTimePreproc = UsedTime; + UsedTimeCompute = 0.0; +} + +CDeformationDriver::~CDeformationDriver(void) {} + +void CDeformationDriver::Input_Preprocessing() { + /*--- Initialize a char to store the zone filename. ---*/ + + char zone_file_name[MAX_STRING_SIZE]; + + /*--- Initialize the configuration of the driver. ---*/ + + driver_config = new CConfig(config_file_name, SU2_COMPONENT::SU2_DEF); + + nZone = driver_config->GetnZone(); + + /*--- Loop over all zones to initialize the various classes. In most + cases, nZone is equal to one. This represents the solution of a partial + differential equation on a single block, unstructured mesh. ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Definition of the configuration option class for all zones. In this + constructor, the input configuration file is parsed and all options are + read and stored. ---*/ + + if (driver_config->GetnConfigFiles() > 0) { + strcpy(zone_file_name, driver_config->GetConfigFilename(iZone).c_str()); + config_container[iZone] = new CConfig(driver_config, zone_file_name, SU2_COMPONENT::SU2_DEF, iZone, nZone, true); + } else { + config_container[iZone] = new CConfig(driver_config, config_file_name, SU2_COMPONENT::SU2_DEF, iZone, nZone, true); + } + + config_container[iZone]->SetMPICommunicator(SU2_MPI::GetComm()); + } + + /*--- Set the multi-zone part of the problem. ---*/ + + if (driver_config->GetMultizone_Problem()) { + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Set the interface markers for multi-zone. ---*/ + + config_container[iZone]->SetMultizone(driver_config, config_container); + } + } + + /*--- Keep a reference to the main (ZONE 0) config. ---*/ + + main_config = config_container[ZONE_0]; +} + +void CDeformationDriver::Geometrical_Preprocessing() { + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/ + + CGeometry* geometry_aux = nullptr; + + /*--- All ranks process the grid and call ParMETIS for partitioning. ---*/ + + geometry_aux = new CPhysicalGeometry(config_container[iZone], iZone, nZone); + + /*--- Color the initial grid and set the send-receive domains (ParMETIS). ---*/ + + geometry_aux->SetColorGrid_Parallel(config_container[iZone]); + + /*--- Build the grid data structures using the ParMETIS coloring. ---*/ + + unsigned short nInst_Zone = nInst[iZone]; + unsigned short nMesh = 1; + + geometry_container[iZone] = new CGeometry**[nInst_Zone](); + geometry_container[iZone][INST_0] = new CGeometry*[nMesh](); + geometry_container[iZone][INST_0][MESH_0] = new CPhysicalGeometry(geometry_aux, config_container[iZone]); + + /*--- Deallocate the memory of geometry_aux. ---*/ + + delete geometry_aux; + + /*--- Add the Send/Receive boundaries. ---*/ + + geometry_container[iZone][INST_0][MESH_0]->SetSendReceive(config_container[iZone]); + + /*--- Add the Send/Receive boundaries. ---*/ + + geometry_container[iZone][INST_0][MESH_0]->SetBoundaries(config_container[iZone]); + + /*--- Computational grid preprocessing. ---*/ + + if (rank == MASTER_NODE) cout << endl << "----------------------- Preprocessing computations ----------------------" << endl; + + /*--- Compute elements surrounding points, points surrounding points. ---*/ + + if (rank == MASTER_NODE) cout << "Setting local point connectivity." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetPoint_Connectivity(); + + /*--- Check the orientation before computing geometrical quantities. ---*/ + + geometry_container[iZone][INST_0][MESH_0]->SetBoundVolume(); + if (config_container[iZone]->GetReorientElements()) { + if (rank == MASTER_NODE) cout << "Checking the numerical grid orientation of the interior elements." << endl; + geometry_container[iZone][INST_0][MESH_0]->Check_IntElem_Orientation(config_container[iZone]); + geometry_container[iZone][INST_0][MESH_0]->Check_BoundElem_Orientation(config_container[iZone]); + } + + /*--- Create the edge structure. ---*/ + + if (rank == MASTER_NODE) cout << "Identify edges and vertices." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetEdges(); + geometry_container[iZone][INST_0][MESH_0]->SetVertex(config_container[iZone]); + + if (config_container[iZone]->GetDesign_Variable(0) != NO_DEFORMATION) { + /*--- Create the dual control volume structures. ---*/ + + if (rank == MASTER_NODE) cout << "Setting the bound control volume structure." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetControlVolume(config_container[iZone], ALLOCATE); + geometry_container[iZone][INST_0][MESH_0]->SetBoundControlVolume(config_container[iZone], ALLOCATE); + } + + /*--- Create the point-to-point MPI communication structures. ---*/ + + geometry_container[iZone][INST_0][MESH_0]->PreprocessP2PComms(geometry_container[iZone][INST_0][MESH_0],config_container[iZone]); + } + + /*--- Get the number of dimensions. ---*/ + + nDim = geometry_container[ZONE_0][INST_0][MESH_0]->GetnDim(); + + /*--- Keep a reference to the main (ZONE_0, INST_0, MESH_0) geometry. ---*/ + + main_geometry = geometry_container[ZONE_0][INST_0][MESH_0]; +} + +void CDeformationDriver::Output_Preprocessing() { + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Allocate the mesh output. ---*/ + + output_container[iZone] = new CMeshOutput(config_container[iZone], geometry_container[iZone][INST_0][MESH_0]->GetnDim()); + + /*--- Preprocess the volume output. ---*/ + + output_container[iZone]->PreprocessVolumeOutput(config_container[iZone]); + + /*--- Preprocess history. --- */ + + output_container[iZone]->PreprocessHistoryOutput(config_container[iZone], false); + } +} + +void CDeformationDriver::Solver_Preprocessing() { + for (iZone = 0; iZone < nZone; iZone++) { + unsigned short nInst_Zone = nInst[iZone]; + unsigned short nMesh = 1; + unsigned short nSols = MAX_SOLS; + + solver_container[iZone] = new CSolver***[nInst_Zone](); + solver_container[iZone][INST_0] = new CSolver**[nMesh](); + solver_container[iZone][INST_0][MESH_0] = new CSolver*[nSols](); + solver_container[iZone][INST_0][MESH_0][MESH_SOL] = + new CMeshSolver(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + } +} + +void CDeformationDriver::Numerics_Preprocessing() { + for (iZone = 0; iZone < nZone; iZone++) { + unsigned short nInst_Zone = nInst[iZone]; + unsigned short nMesh = 1; + unsigned short nSols = MAX_SOLS; + unsigned int nTerm = omp_get_num_threads() * MAX_TERMS; + + numerics_container[iZone] = new CNumerics****[nInst_Zone](); + numerics_container[iZone][INST_0] = new CNumerics***[nMesh](); + numerics_container[iZone][INST_0][MESH_0] = new CNumerics**[nSols](); + numerics_container[iZone][INST_0][MESH_0][MESH_SOL] = new CNumerics*[nTerm](); + + for (int thread = 0; thread < omp_get_max_threads(); ++thread) { + const int iTerm = FEA_TERM + thread * MAX_TERMS; + const int nDim = geometry_container[iZone][INST_0][MESH_0]->GetnDim(); + + numerics_container[iZone][INST_0][MESH_0][MESH_SOL][iTerm] = new CFEAMeshElasticity( + nDim, nDim, geometry_container[iZone][INST_0][MESH_0]->GetnElem(), config_container[iZone]); + } + } +} + +void CDeformationDriver::Preprocess() {} + +void CDeformationDriver::Run() { + /* --- Start measuring computation time. ---*/ + + StartTime = SU2_MPI::Wtime(); + + /*--- Surface grid deformation using design variables. ---*/ + + if (driver_config->GetDeform_Mesh()) { + Update(); + } else { + Update_Legacy(); + } + + /*--- Synchronization point after a single solver iteration. Compute the wall clock time required. ---*/ + + StopTime = SU2_MPI::Wtime(); + + UsedTimeCompute = StopTime - StartTime; + if (rank == MASTER_NODE) { + cout << "\nCompleted in " << fixed << UsedTimeCompute << " seconds on " << size; + + if (size == 1) cout << " core." << endl; else cout << " cores." << endl; + } + + /*--- Output the deformed mesh. ---*/ + + Output(); +} + +void CDeformationDriver::Update() { + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Set the stiffness of each element mesh into the mesh numerics. ---*/ + + solver_container[iZone][INST_0][MESH_0][MESH_SOL]->SetMesh_Stiffness(numerics_container[iZone][INST_0][MESH_0][MESH_SOL], config_container[iZone]); + + /*--- Deform the volume grid around the new boundary locations. ---*/ + /*--- Force the number of levels to be 0 because in this driver we do not build MG levels. ---*/ + const auto nMGLevels = config_container[iZone]->GetnMGLevels(); + config_container[iZone]->SetMGLevels(0); + solver_container[iZone][INST_0][MESH_0][MESH_SOL]->DeformMesh(geometry_container[iZone][INST_0], + numerics_container[iZone][INST_0][MESH_0][MESH_SOL], + config_container[iZone]); + config_container[iZone]->SetMGLevels(nMGLevels); + } +} + +void CDeformationDriver::Update_Legacy() { + for (iZone = 0; iZone < nZone; iZone++) { + if (config_container[iZone]->GetDesign_Variable(0) != NO_DEFORMATION) { + unsigned short nInst_Zone = nInst[iZone]; + + /*--- Definition of the Class for grid movement. ---*/ + + grid_movement[iZone] = new CVolumetricMovement*[nInst_Zone](); + grid_movement[iZone][INST_0] = new CVolumetricMovement(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + + /*--- Save original coordinates to be reused in convexity checking procedure. ---*/ + + auto OriginalCoordinates = geometry_container[iZone][INST_0][MESH_0]->nodes->GetCoord(); + + /*--- First check for volumetric grid deformation/transformations. ---*/ + + if (config_container[iZone]->GetDesign_Variable(0) == SCALE_GRID) { + if (rank == MASTER_NODE) + cout << endl << "--------------------- Volumetric grid scaling (ZONE " << iZone << ") ------------------" << endl; + grid_movement[iZone][INST_0]->SetVolume_Scaling(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], false); + + } else if (config_container[iZone]->GetDesign_Variable(0) == TRANSLATE_GRID) { + if (rank == MASTER_NODE) + cout << endl << "------------------- Volumetric grid translation (ZONE " << iZone << ") ----------------" << endl; + grid_movement[iZone][INST_0]->SetVolume_Translation(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], false); + + } else if (config_container[iZone]->GetDesign_Variable(0) == ROTATE_GRID) { + if (rank == MASTER_NODE) + cout << endl << "--------------------- Volumetric grid rotation (ZONE " << iZone << ") -----------------" << endl; + grid_movement[iZone][INST_0]->SetVolume_Rotation(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], false); + + } else { + /*--- If no volume-type deformations are requested, then this is a + * surface-based deformation or FFD set up. ---*/ + + if (rank == MASTER_NODE) + cout << endl << "--------------------- Surface grid deformation (ZONE " << iZone << ") -----------------" << endl; + + /*--- Definition and initialization of the surface deformation class. ---*/ + + surface_movement[iZone] = new CSurfaceMovement(); + haveSurfaceDeformation = true; + + /*--- Copy coordinates to the surface structure. ---*/ + + surface_movement[iZone]->CopyBoundary(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + + /*--- Surface grid deformation. ---*/ + + if (rank == MASTER_NODE) cout << "Performing the deformation of the surface grid." << endl; + auto TotalDeformation = surface_movement[iZone]->SetSurface_Deformation(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + + if (config_container[iZone]->GetDesign_Variable(0) != FFD_SETTING) { + if (rank == MASTER_NODE) + cout << endl << "------------------- Volumetric grid deformation (ZONE " << iZone << ") ----------------" << endl; + + if (rank == MASTER_NODE) cout << "Performing the deformation of the volumetric grid." << endl; + grid_movement[iZone][INST_0]->SetVolume_Deformation(geometry_container[iZone][INST_0][MESH_0],config_container[iZone], false); + + /*--- Get parameters for convexity check. ---*/ + bool ConvexityCheck; + unsigned short ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth; + + tie(ConvexityCheck, ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth) = config_container[iZone]->GetConvexityCheck(); + + /*--- Recursively change deformations if there are non-convex elements. ---*/ + + if (ConvexityCheck && geometry_container[iZone][INST_0][MESH_0]->GetnNonconvexElements() > 0) { + if (rank == MASTER_NODE) { + cout << "Non-convex elements present after deformation." << endl; + cout << "Recursively lowering deformation magnitude." << endl; + } + + /*--- Load initial deformation values. ---*/ + + auto InitialDeformation = TotalDeformation; + + unsigned short ConvexityCheckIter, RecursionDepth = 0; + su2double DeformationFactor = 1.0, DeformationDifference = 1.0; + for (ConvexityCheckIter = 1; ConvexityCheckIter <= ConvexityCheck_MaxIter; ConvexityCheckIter++) { + /*--- Recursively change deformation magnitude (decrease for non-convex elements, increase otherwise). ---*/ + + DeformationDifference /= 2.0; + + if (geometry_container[iZone][INST_0][MESH_0]->GetnNonconvexElements() > 0) { + DeformationFactor -= DeformationDifference; + } else { + RecursionDepth += 1; + + if (RecursionDepth == ConvexityCheck_MaxDepth) { + if (rank == MASTER_NODE) { + cout << "Maximum recursion depth reached." << endl; + cout << "Remaining amount of original deformation: "; + cout << DeformationFactor * 100.0 << " percent. " << endl; + } + break; + } + + DeformationFactor += DeformationDifference; + } + + /*--- Load mesh to start every iteration with an undeformed grid. ---*/ + + for (auto iPoint = 0ul; iPoint < OriginalCoordinates.rows(); iPoint++) { + for (auto iDim = 0ul; iDim < OriginalCoordinates.cols(); iDim++) { + geometry_container[iZone][INST_0][MESH_0]->nodes->SetCoord(iPoint, iDim, OriginalCoordinates(iPoint, iDim)); + } + } + + /*--- Set deformation magnitude as percentage of initial deformation. ---*/ + + for (auto iDV = 0u; iDV < driver_config->GetnDV(); iDV++) { + for (auto iDV_Value = 0u; iDV_Value < driver_config->GetnDV_Value(iDV); iDV_Value++) { + config_container[iZone]->SetDV_Value(iDV, iDV_Value,InitialDeformation[iDV][iDV_Value] * DeformationFactor); + } + } + + /*--- Surface grid deformation. ---*/ + + if (rank == MASTER_NODE) cout << "Performing the deformation of the surface grid." << endl; + + TotalDeformation = surface_movement[iZone]->SetSurface_Deformation(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + + if (rank == MASTER_NODE) + cout << endl << "------------------- Volumetric grid deformation (ZONE " << iZone << ") ----------------" << endl; + + if (rank == MASTER_NODE) cout << "Performing the deformation of the volumetric grid." << endl; + grid_movement[iZone][INST_0]->SetVolume_Deformation(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], false); + + if (rank == MASTER_NODE) { + cout << "Number of non-convex elements for iteration " << ConvexityCheckIter << ": "; + cout << geometry_container[iZone][INST_0][MESH_0]->GetnNonconvexElements() << endl; + cout << "Remaining amount of original deformation: "; + cout << DeformationFactor * 100.0 << " percent. " << endl; + } + } + } + } + } + } + } +} + +void CDeformationDriver::Output() { + /*--- Output deformed grid for visualization, if requested (surface and volumetric), in parallel + requires to move all the data to the master node. ---*/ + + if (rank == MASTER_NODE) + cout << endl << "----------------------- Write deformed grid files -----------------------" << endl; + + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Compute Mesh Quality if requested. Necessary geometry preprocessing re-done beforehand. ---*/ + + if (config_container[iZone]->GetWrt_MeshQuality() && !driver_config->GetStructuralProblem()) { + if (rank == MASTER_NODE) cout << "Recompute geometry properties necessary to evaluate mesh quality statistics.\n"; + + geometry_container[iZone][INST_0][MESH_0]->SetPoint_Connectivity(); + geometry_container[iZone][INST_0][MESH_0]->SetBoundVolume(); + geometry_container[iZone][INST_0][MESH_0]->SetEdges(); + geometry_container[iZone][INST_0][MESH_0]->SetVertex(config_container[iZone]); + geometry_container[iZone][INST_0][MESH_0]->SetControlVolume(config_container[iZone], ALLOCATE); + geometry_container[iZone][INST_0][MESH_0]->SetBoundControlVolume(config_container[iZone], ALLOCATE); + + if (rank == MASTER_NODE) cout << "Computing mesh quality statistics for the dual control volumes.\n"; + geometry_container[iZone][INST_0][MESH_0]->ComputeMeshQualityStatistics(config_container[iZone]); + } + + /*--- Load the data. --- */ + + output_container[iZone]->Load_Data(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], nullptr); + + output_container[iZone]->WriteToFile(config_container[iZone], geometry_container[iZone][INST_0][MESH_0],OUTPUT_TYPE::MESH, driver_config->GetMesh_Out_FileName()); + + /*--- Set the file names for the visualization files. ---*/ + + output_container[iZone]->SetVolume_Filename("volume_deformed"); + output_container[iZone]->SetSurface_Filename("surface_deformed"); + + for (unsigned short iFile = 0; iFile < config_container[iZone]->GetnVolumeOutputFiles(); iFile++) { + auto FileFormat = config_container[iZone]->GetVolumeOutputFiles(); + if (FileFormat[iFile] != OUTPUT_TYPE::RESTART_ASCII && FileFormat[iFile] != OUTPUT_TYPE::RESTART_BINARY && + FileFormat[iFile] != OUTPUT_TYPE::CSV) + output_container[iZone]->WriteToFile(config_container[iZone], geometry_container[iZone][INST_0][MESH_0],FileFormat[iFile]); + } + } + + if (!driver_config->GetDeform_Mesh()) { + if ((config_container[ZONE_0]->GetDesign_Variable(0) != NO_DEFORMATION) && + (config_container[ZONE_0]->GetDesign_Variable(0) != SCALE_GRID) && + (config_container[ZONE_0]->GetDesign_Variable(0) != TRANSLATE_GRID) && + (config_container[ZONE_0]->GetDesign_Variable(0) != ROTATE_GRID)) { + + /*--- Write the free form deformation boxes after deformation if defined. ---*/ + if (!haveSurfaceDeformation) { + if (rank == MASTER_NODE) cout << "No FFD information available." << endl; + } else { + if (rank == MASTER_NODE) cout << "Adding any FFD information to the SU2 file." << endl; + + surface_movement[ZONE_0]->WriteFFDInfo(surface_movement, geometry_container, config_container); + } + } + } +} + +void CDeformationDriver::Postprocessing() { + if (rank == MASTER_NODE) + cout << endl << "------------------------- Solver Postprocessing -------------------------" << endl; + + if (driver_config->GetDeform_Mesh()) { + for (iZone = 0; iZone < nZone; iZone++) { + if (numerics_container[iZone] != nullptr) { + for (unsigned int iTerm = 0; iTerm < MAX_TERMS * omp_get_max_threads(); iTerm++) { + delete numerics_container[iZone][INST_0][MESH_0][MESH_SOL][iTerm]; + delete[] numerics_container[iZone][INST_0][MESH_0][MESH_SOL]; + delete[] numerics_container[iZone][INST_0][MESH_0]; + delete[] numerics_container[iZone][INST_0]; + } + delete[] numerics_container[iZone]; + } + } + delete[] numerics_container; + if (rank == MASTER_NODE) cout << "Deleted CNumerics container." << endl; + + for (iZone = 0; iZone < nZone; iZone++) { + if (solver_container[iZone] != nullptr) { + delete solver_container[iZone][INST_0][MESH_0][MESH_SOL]; + delete[] solver_container[iZone][INST_0][MESH_0]; + delete[] solver_container[iZone][INST_0]; + delete[] solver_container[iZone]; + } + } + delete[] solver_container; + if (rank == MASTER_NODE) cout << "Deleted CSolver container." << endl; + } + + if (geometry_container != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete geometry_container[iZone][INST_0][MESH_0]; + delete[] geometry_container[iZone][INST_0]; + delete[] geometry_container[iZone]; + } + delete[] geometry_container; + } + if (rank == MASTER_NODE) cout << "Deleted CGeometry container." << endl; + + delete[] FFDBox; + if (rank == MASTER_NODE) cout << "Deleted CFreeFormDefBox class." << endl; + + if (surface_movement != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete surface_movement[iZone]; + } + delete[] surface_movement; + } + if (rank == MASTER_NODE) cout << "Deleted CSurfaceMovement class." << endl; + + if (grid_movement != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete grid_movement[iZone][INST_0]; + delete[] grid_movement[iZone]; + } + delete[] grid_movement; + } + if (rank == MASTER_NODE) cout << "Deleted CVolumetricMovement class." << endl; + + if (config_container != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete config_container[iZone]; + } + delete[] config_container; + } + delete driver_config; + if (rank == MASTER_NODE) cout << "Deleted CConfig container." << endl; + + if (output_container != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete output_container[iZone]; + } + delete[] output_container; + } + if (rank == MASTER_NODE) cout << "Deleted COutput class." << endl; + + if (nInst != nullptr) delete[] nInst; + + /*--- Exit the solver cleanly. ---*/ + + if (rank == MASTER_NODE) + cout << endl << "------------------------- Exit Success (SU2_DEF) ------------------------" << endl << endl; +} + +void CDeformationDriver::CommunicateMeshDisplacements(void) { + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->InitiateComms(geometry_container[ZONE_0][INST_0][MESH_0],config_container[ZONE_0], MESH_DISPLACEMENTS); + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->CompleteComms(geometry_container[ZONE_0][INST_0][MESH_0],config_container[ZONE_0], MESH_DISPLACEMENTS); +} diff --git a/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp b/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp new file mode 100644 index 000000000000..25ae73f9abd8 --- /dev/null +++ b/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp @@ -0,0 +1,1086 @@ +/*! + * \file CDiscAdjDeformationDriver.cpp + * \brief Main subroutines for driving the projection of sensitivities. + * \author T. Economon, H. Kline, R. Sanchez, A. Gastaldi, H. Patel + * \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 . + */ + +#define ENABLE_MAPS +#include "../../../Common/include/CConfig.hpp" +#undef ENABLE_MAPS + +#include "../../../Common/include/geometry/CPhysicalGeometry.hpp" +#include "../../../Common/include/grid_movement/CSurfaceMovement.hpp" +#include "../../../Common/include/grid_movement/CVolumetricMovement.hpp" +#include "../../../SU2_CFD/include/numerics/CGradSmoothing.hpp" +#include "../../../SU2_CFD/include/output/CBaselineOutput.hpp" +#include "../../../SU2_CFD/include/solvers/CBaselineSolver.hpp" +#include "../../../SU2_CFD/include/solvers/CGradientSmoothingSolver.hpp" +#include "../../include/drivers/CDiscAdjDeformationDriver.hpp" + +using namespace std; + +CDiscAdjDeformationDriver::CDiscAdjDeformationDriver(char* confFile, SU2_Comm MPICommunicator) + : CDriverBase(confFile, 1, MPICommunicator) { + +/*--- Initialize MeDiPack (must also be here to initialize it from Python). ---*/ +#ifdef HAVE_MPI +#if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) + SU2_MPI::Init_AMPI(); +#endif +#endif + + SU2_MPI::SetComm(MPICommunicator); + + rank = SU2_MPI::GetRank(); + size = SU2_MPI::GetSize(); + + /*--- Initialize containers. --- */ + + SetContainers_Null(); + + /*--- Preprocessing of the config files. ---*/ + + Input_Preprocessing(); + + /*--- Initialize structure to store the gradient. ---*/ + + unsigned short nDV = config_container[ZONE_0]->GetnDV(); + Gradient = new su2double*[nDV]; + + for (auto iDV = 0u; iDV < nDV; iDV++) { + /*--- Initialize to zero ---*/ + unsigned short nValue = config_container[ZONE_0]->GetnDV_Value(iDV); + + Gradient[iDV] = new su2double[nValue]; + } + + /*--- Set up a timer for performance benchmarking. ---*/ + + StartTime = SU2_MPI::Wtime(); + + /*--- Preprocessing of the geometry for all zones. ---*/ + + Geometrical_Preprocessing(); + + /*--- Preprocessing of the outputs for all zones. ---*/ + + Output_Preprocessing(); + + /*--- Preprocessing time is reported now, but not included in the next compute portion. ---*/ + + StopTime = SU2_MPI::Wtime(); + + /*--- Compute the total time for performance benchmarking. ---*/ + + UsedTime = StopTime - StartTime; + UsedTimePreproc = UsedTime; + UsedTimeCompute = 0.0; +} + +CDiscAdjDeformationDriver::~CDiscAdjDeformationDriver(void) {} + +void CDiscAdjDeformationDriver::Input_Preprocessing() { + /*--- Initialize a char to store the zone filename. ---*/ + + char zone_file_name[MAX_STRING_SIZE]; + + /*--- Initialize the configuration of the driver. ---*/ + + driver_config = new CConfig(config_file_name, SU2_COMPONENT::SU2_DEF); + + nZone = driver_config->GetnZone(); + + /*--- Loop over all zones to initialize the various classes. In most + * cases, nZone is equal to one. This represents the solution of a partial + * differential equation on a single block, unstructured mesh. ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Definition of the configuration option class for all zones. In this + * constructor, the input configuration file is parsed and all options are + * read and stored. ---*/ + + if (driver_config->GetnConfigFiles() > 0) { + strcpy(zone_file_name, driver_config->GetConfigFilename(iZone).c_str()); + config_container[iZone] = new CConfig(driver_config, zone_file_name, SU2_COMPONENT::SU2_DOT, iZone, nZone, true); + } else { + config_container[iZone] = new CConfig(driver_config, config_file_name, SU2_COMPONENT::SU2_DOT, iZone, nZone, true); + } + + config_container[iZone]->SetMPICommunicator(SU2_MPI::GetComm()); + + if (!config_container[iZone]->GetDiscrete_Adjoint() && !config_container[iZone]->GetContinuous_Adjoint()) { + SU2_MPI::Error("An adjoint solver (discrete or continuous) was not specified in the configuration file.", CURRENT_FUNCTION); + } + } + + /*--- Set the multi-zone part of the problem. ---*/ + + if (driver_config->GetMultizone_Problem()) { + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Set the interface markers for multi-zone. ---*/ + + config_container[iZone]->SetMultizone(driver_config, config_container); + } + } + + /*--- Keep a reference to the main (ZONE 0) config. ---*/ + + main_config = config_container[ZONE_0]; +} + +void CDiscAdjDeformationDriver::Geometrical_Preprocessing() { + /*--- Loop over all zones to initialize the various classes. In most + * cases, nZone is equal to one. This represents the solution of a partial + * differential equation on a single block, unstructured mesh. ---*/ + + unsigned short nMesh = 1; + + for (iZone = 0; iZone < nZone; iZone++) { + /*--- Determine if the FEM solver is used, which decides the type of geometry classes that are instantiated. ---*/ + + const bool fem_solver = config_container[iZone]->GetFEMSolver(); + + /*--- Read the number of instances for each zone. ---*/ + + nInst[iZone] = config_container[iZone]->GetnTimeInstances(); + + geometry_container[iZone] = new CGeometry**[nInst[iZone]]; + + for (iInst = 0; iInst < nInst[iZone]; iInst++) { + /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/ + + CGeometry* geometry_aux = nullptr; + + /*--- All ranks process the grid and call ParMETIS for partitioning. ---*/ + + geometry_aux = new CPhysicalGeometry(config_container[iZone], iZone, nZone); + + /*--- Color the initial grid and set the send-receive domains (ParMETIS). ---*/ + + if (fem_solver) { + geometry_aux->SetColorFEMGrid_Parallel(config_container[iZone]); + } else { + geometry_aux->SetColorGrid_Parallel(config_container[iZone]); + } + + /*--- Build the grid data structures using the ParMETIS coloring. ---*/ + + if (fem_solver) { + switch (config_container[iZone]->GetKind_FEM_Flow()) { + case DG: + geometry_container[iZone][iInst] = new CGeometry*[nMesh]; + geometry_container[iZone][iInst][MESH_0] = new CMeshFEM_DG(geometry_aux, config_container[iZone]); + break; + } + } else { + geometry_container[iZone][iInst] = new CGeometry*[nMesh]; + geometry_container[iZone][iInst][MESH_0] = new CPhysicalGeometry(geometry_aux, config_container[iZone]); + } + + /*--- Deallocate the memory of geometry_aux. ---*/ + + delete geometry_aux; + + /*--- Add the Send/Receive boundaries. ---*/ + + geometry_container[iZone][iInst][MESH_0]->SetSendReceive(config_container[iZone]); + + /*--- Add the Send/Receive boundaries. ---*/ + + geometry_container[iZone][iInst][MESH_0]->SetBoundaries(config_container[iZone]); + + /*--- Create the vertex structure (required for MPI). ---*/ + + if (rank == MASTER_NODE) cout << "Identify vertices." << endl; + geometry_container[iZone][iInst][MESH_0]->SetVertex(config_container[iZone]); + + /*--- Store the global to local mapping after preprocessing. ---*/ + + if (rank == MASTER_NODE) cout << "Storing a mapping from global to local point index." << endl; + geometry_container[iZone][iInst][MESH_0]->SetGlobal_to_Local_Point(); + + /*--- Test for a fem solver, because some more work must be done. ---*/ + + if (fem_solver) { + /*--- Carry out a dynamic cast to CMeshFEM_DG, such that it is not needed to + * define all virtual functions in the base class CGeometry. ---*/ + + CMeshFEM_DG* DGMesh = dynamic_cast(geometry_container[iZone][iInst][MESH_0]); + + /*--- Determine the standard elements for the volume elements. ---*/ + + if (rank == MASTER_NODE) cout << "Creating standard volume elements." << endl; + DGMesh->CreateStandardVolumeElements(config_container[iZone]); + + /*--- Create the face information needed to compute the contour integral + * for the elements in the Discontinuous Galerkin formulation. ---*/ + + if (rank == MASTER_NODE) cout << "Creating face information." << endl; + DGMesh->CreateFaces(config_container[iZone]); + } + } + + if (rank == MASTER_NODE) + cout << "\n----------------------- Preprocessing computations ----------------------" << endl; + + /*--- Compute elements surrounding points, points surrounding points. ---*/ + + if (rank == MASTER_NODE) cout << "Setting local point connectivity." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetPoint_Connectivity(); + + /*--- Check the orientation before computing geometrical quantities. ---*/ + + geometry_container[iZone][INST_0][MESH_0]->SetBoundVolume(); + if (config_container[iZone]->GetReorientElements()) { + if (rank == MASTER_NODE) cout << "Checking the numerical grid orientation of the elements." << endl; + geometry_container[iZone][INST_0][MESH_0]->Check_IntElem_Orientation(config_container[iZone]); + geometry_container[iZone][INST_0][MESH_0]->Check_BoundElem_Orientation(config_container[iZone]); + } + + /*--- Create the edge structure. ---*/ + + if (rank == MASTER_NODE) cout << "Identify edges and vertices." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetEdges(); + geometry_container[iZone][INST_0][MESH_0]->SetVertex(config_container[iZone]); + + /*--- Create the dual control volume structures. ---*/ + + if (rank == MASTER_NODE) cout << "Setting the bound control volume structure." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetBoundControlVolume(config_container[ZONE_0], ALLOCATE); + + /*--- Store the global to local mapping after preprocessing. ---*/ + + if (rank == MASTER_NODE) cout << "Storing a mapping from global to local point index." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetGlobal_to_Local_Point(); + + /*--- Create the point-to-point MPI communication structures. ---*/ + + geometry_container[iZone][INST_0][MESH_0]->PreprocessP2PComms(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone]); + } + + /*--- Keep a reference to the main (ZONE_0, INST_0, MESH_0) geometry. ---*/ + + main_geometry = geometry_container[ZONE_0][INST_0][MESH_0]; +} + +void CDiscAdjDeformationDriver::Output_Preprocessing() {} + +void CDiscAdjDeformationDriver::Preprocess() { + for (iZone = 0; iZone < nZone; iZone++) { + if (!config_container[iZone]->GetDiscrete_Adjoint()) { + if (rank == MASTER_NODE) cout << "Reading surface sensitivities at each node from file." << endl; + geometry_container[iZone][INST_0][MESH_0]->SetBoundSensitivity(config_container[iZone]); + + } else { + if (rank == MASTER_NODE) cout << "Reading volume sensitivities at each node from file." << endl; + unsigned short nInst_Zone = nInst[iZone]; + + grid_movement[iZone] = new CVolumetricMovement*[nInst_Zone](); + grid_movement[iZone][INST_0] = new CVolumetricMovement(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + + /*--- Read in sensitivities from file. ---*/ + + if (config_container[ZONE_0]->GetSensitivity_Format() == UNORDERED_ASCII) + geometry_container[iZone][INST_0][MESH_0]->ReadUnorderedSensitivity(config_container[iZone]); + else + geometry_container[iZone][INST_0][MESH_0]->SetSensitivity(config_container[iZone]); + } + } +} + +void CDiscAdjDeformationDriver::Run() { + for (iZone = 0; iZone < nZone; iZone++) { + if (!config_container[iZone]->GetDiscrete_Adjoint()) { + continue; + } else { + if (rank == MASTER_NODE) + cout << "\n---------------------- Mesh sensitivity computation ---------------------" << endl; + if (config_container[iZone]->GetDiscrete_Adjoint() && config_container[iZone]->GetSmoothGradient() && + config_container[iZone]->GetSobMode() == ENUM_SOBOLEV_MODUS::MESH_LEVEL) { + DerivativeTreatment_MeshSensitivity(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], + grid_movement[iZone][INST_0]); + } else { + grid_movement[iZone][INST_0]->SetVolume_Deformation(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone], false, true); + } + } + } + + if (config_container[ZONE_0]->GetDiscrete_Adjoint()) { + if (rank == MASTER_NODE) + cout << "\n------------------------ Mesh sensitivity Output ------------------------" << endl; + SetSensitivity_Files(geometry_container, config_container, nZone); + } + + Gradient_file.precision(driver_config->OptionIsSet("OUTPUT_PRECISION") ? driver_config->GetOutput_Precision() : 6); + + /*--- For multi-zone computations the gradient contributions are summed up and written into one file. ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + if ((config_container[iZone]->GetDesign_Variable(0) != NONE) && + (config_container[iZone]->GetDesign_Variable(0) != SURFACE_FILE)) { + if (rank == MASTER_NODE) + cout << "\n---------- Start gradient evaluation using sensitivity information ----------" << endl; + + /*--- Definition of the Class for surface deformation. ---*/ + + surface_movement[iZone] = new CSurfaceMovement(); + + /*--- Copy coordinates to the surface structure. ---*/ + + surface_movement[iZone]->CopyBoundary(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + + /*--- If AD mode is enabled we can use it to compute the projection, otherwise we use finite differences. ---*/ + + if (config_container[iZone]->GetAD_Mode()) { + if (config_container[iZone]->GetSmoothGradient()) { + DerivativeTreatment_Gradient(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], + grid_movement[iZone][INST_0], surface_movement[iZone], Gradient); + } else { + SetProjection_AD(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], surface_movement[iZone], Gradient); + } + } else { + SetProjection_FD(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], surface_movement[iZone], Gradient); + } + } + } + + /*--- Write the gradient to a file. ---*/ + + if (rank == MASTER_NODE) Gradient_file.open(config_container[ZONE_0]->GetObjFunc_Grad_FileName().c_str(), ios::out); + + /*--- Print gradients to screen and writes to file. ---*/ + + OutputGradient(Gradient, config_container[ZONE_0], Gradient_file); +} + +void CDiscAdjDeformationDriver::Postprocessing() { + for (auto iDV = 0u; iDV < config_container[ZONE_0]->GetnDV(); iDV++) { + delete[] Gradient[iDV]; + } + delete[] Gradient; + + delete driver_config; + driver_config = nullptr; + + if (rank == MASTER_NODE) + cout << "\n------------------------- Solver Postprocessing -------------------------" << endl; + + for (iZone = 0; iZone < nZone; iZone++) { + if (numerics_container[iZone] != nullptr) { + delete[] numerics_container[iZone]; + } + } + delete[] numerics_container; + if (rank == MASTER_NODE) cout << "Deleted CNumerics container." << endl; + + for (iZone = 0; iZone < nZone; iZone++) { + if (solver_container[iZone] != nullptr) { + delete[] solver_container[iZone]; + } + } + delete[] solver_container; + if (rank == MASTER_NODE) cout << "Deleted CSolver container." << endl; + + if (geometry_container != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + if (geometry_container[iZone] != nullptr) { + for (iInst = 0; iInst < nInst[iZone]; iInst++) { + delete geometry_container[iZone][iInst][MESH_0]; + delete[] geometry_container[iZone][iInst]; + } + delete[] geometry_container[iZone]; + } + } + delete[] geometry_container; + } + if (rank == MASTER_NODE) cout << "Deleted CGeometry container." << endl; + + for (iZone = 0; iZone < nZone; iZone++) { + delete[] FFDBox[iZone]; + } + delete[] FFDBox; + if (rank == MASTER_NODE) cout << "Deleted CFreeFormDefBox class." << endl; + + if (surface_movement != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete surface_movement[iZone]; + } + delete[] surface_movement; + } + if (rank == MASTER_NODE) cout << "Deleted CSurfaceMovement class." << endl; + + if (grid_movement != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + if (grid_movement[iZone] != nullptr) { + for (iInst = 0; iInst < nInst[iZone]; iInst++) { + delete grid_movement[iZone][iInst]; + } + delete[] grid_movement[iZone]; + } + } + delete[] grid_movement; + } + if (rank == MASTER_NODE) cout << "Deleted CVolumetricMovement class." << endl; + + if (config_container != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete config_container[iZone]; + } + delete[] config_container; + } + if (rank == MASTER_NODE) cout << "Deleted CConfig container." << endl; + + if (output_container != nullptr) { + for (iZone = 0; iZone < nZone; iZone++) { + delete output_container[iZone]; + } + delete[] output_container; + } + if (rank == MASTER_NODE) cout << "Deleted COutput class." << endl; + + if (nInst != nullptr) delete[] nInst; + + /*--- Exit the solver cleanly. ---*/ + + if (rank == MASTER_NODE) + cout << endl << "------------------------- Exit Success (SU2_DEF_AD) ------------------------" << endl << endl; +} + +void CDiscAdjDeformationDriver::SetProjection_FD(CGeometry* geometry, CConfig* config, + CSurfaceMovement* surface_movement, su2double** Gradient) { + unsigned short iDV, nDV, iFFDBox, nDV_Value, iMarker, iDim; + unsigned long iVertex, iPoint; + su2double delta_eps, my_Gradient, localGradient, *Normal, dS, *VarCoord, Sensitivity, dalpha[3], deps[3], dalpha_deps; + bool *UpdatePoint, MoveSurface, Local_MoveSurface; + CFreeFormDefBox** FFDBox; + + int rank = SU2_MPI::GetRank(); + + nDV = config->GetnDV(); + + /*--- Boolean controlling points to be updated. ---*/ + + UpdatePoint = new bool[geometry->GetnPoint()]; + + /*--- Definition of the FFD deformation class. ---*/ + + unsigned short nFFDBox = MAX_NUMBER_FFD; + FFDBox = new CFreeFormDefBox*[nFFDBox]; + for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) FFDBox[iFFDBox] = nullptr; + + for (iDV = 0; iDV < nDV; iDV++) { + nDV_Value = config->GetnDV_Value(iDV); + if (nDV_Value != 1) { + SU2_MPI::Error("The projection using finite differences currently only supports a fixed direction of movement for FFD points.", CURRENT_FUNCTION); + } + } + + /*--- Continuous adjoint gradient computation. ---*/ + + if (rank == MASTER_NODE) cout << "Evaluate functional gradient using Finite Differences." << endl; + + for (iDV = 0; iDV < nDV; iDV++) { + MoveSurface = true; + Local_MoveSurface = true; + + /*--- Free form deformation based. ---*/ + + if ((config->GetDesign_Variable(iDV) == FFD_CONTROL_POINT_2D) || + (config->GetDesign_Variable(iDV) == FFD_CAMBER_2D) || (config->GetDesign_Variable(iDV) == FFD_THICKNESS_2D) || + (config->GetDesign_Variable(iDV) == FFD_TWIST_2D) || (config->GetDesign_Variable(iDV) == FFD_CONTROL_POINT) || + (config->GetDesign_Variable(iDV) == FFD_NACELLE) || (config->GetDesign_Variable(iDV) == FFD_GULL) || + (config->GetDesign_Variable(iDV) == FFD_TWIST) || (config->GetDesign_Variable(iDV) == FFD_ROTATION) || + (config->GetDesign_Variable(iDV) == FFD_CAMBER) || (config->GetDesign_Variable(iDV) == FFD_THICKNESS) || + (config->GetDesign_Variable(iDV) == FFD_ANGLE_OF_ATTACK)) { + /*--- Read the FFD information in the first iteration. ---*/ + + if (iDV == 0) { + if (rank == MASTER_NODE) cout << "Read the FFD information from mesh file." << endl; + + /*--- Read the FFD information from the grid file. ---*/ + + surface_movement->ReadFFDInfo(geometry, config, FFDBox, config->GetMesh_FileName()); + + /*--- If the FFDBox was not defined in the input file. ---*/ + + if (!surface_movement->GetFFDBoxDefinition()) { + SU2_MPI::Error("The input grid doesn't have the entire FFD information!", CURRENT_FUNCTION); + } + + for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) { + if (rank == MASTER_NODE) cout << "Checking FFD box dimension." << endl; + surface_movement->CheckFFDDimension(geometry, config, FFDBox[iFFDBox], iFFDBox); + + if (rank == MASTER_NODE) cout << "Check the FFD box intersections with the solid surfaces." << endl; + surface_movement->CheckFFDIntersections(geometry, config, FFDBox[iFFDBox], iFFDBox); + } + + if (rank == MASTER_NODE) + cout << "-------------------------------------------------------------------------" << endl; + } + + if (rank == MASTER_NODE) { + cout << endl << "Design variable number " << iDV << "." << endl; + cout << "Performing 3D deformation of the surface." << endl; + } + + /*--- Apply the control point change. ---*/ + + MoveSurface = false; + + for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) { + /*--- Reset FFD box. ---*/ + + switch (config->GetDesign_Variable(iDV)) { + case FFD_CONTROL_POINT_2D: + Local_MoveSurface = + surface_movement->SetFFDCPChange_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_CAMBER_2D: + Local_MoveSurface = surface_movement->SetFFDCamber_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_THICKNESS_2D: + Local_MoveSurface = + surface_movement->SetFFDThickness_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_TWIST_2D: + Local_MoveSurface = surface_movement->SetFFDTwist_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_CONTROL_POINT: + Local_MoveSurface = surface_movement->SetFFDCPChange(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_NACELLE: + Local_MoveSurface = surface_movement->SetFFDNacelle(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_GULL: + Local_MoveSurface = surface_movement->SetFFDGull(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_TWIST: + Local_MoveSurface = surface_movement->SetFFDTwist(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_ROTATION: + Local_MoveSurface = surface_movement->SetFFDRotation(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_CAMBER: + Local_MoveSurface = surface_movement->SetFFDCamber(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_THICKNESS: + Local_MoveSurface = surface_movement->SetFFDThickness(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_CONTROL_SURFACE: + Local_MoveSurface = + surface_movement->SetFFDControl_Surface(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); + break; + case FFD_ANGLE_OF_ATTACK: + Gradient[iDV][0] = config->GetAoA_Sens(); + break; + } + + /*--- Recompute cartesian coordinates using the new control points position. ---*/ + + if (Local_MoveSurface) { + MoveSurface = true; + surface_movement->SetCartesianCoord(geometry, config, FFDBox[iFFDBox], iFFDBox, true); + } + } + } + + /*--- Hicks-Henne design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == HICKS_HENNE) { + surface_movement->SetHicksHenne(geometry, config, iDV, true); + } + + /*--- Surface bump design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == SURFACE_BUMP) { + surface_movement->SetSurface_Bump(geometry, config, iDV, true); + } + + /*--- Kulfan (CST) design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == CST) { + surface_movement->SetCST(geometry, config, iDV, true); + } + + /*--- Displacement design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == TRANSLATION) { + surface_movement->SetTranslation(geometry, config, iDV, true); + } + + /*--- Angle of Attack design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == ANGLE_OF_ATTACK) { + Gradient[iDV][0] = config->GetAoA_Sens(); + } + + /*--- Scale design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == SCALE) { + surface_movement->SetScale(geometry, config, iDV, true); + } + + /*--- Rotation design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == ROTATION) { + surface_movement->SetRotation(geometry, config, iDV, true); + } + + /*--- NACA_4Digits design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == NACA_4DIGITS) { + surface_movement->SetNACA_4Digits(geometry, config); + } + + /*--- Parabolic design variable. ---*/ + + else if (config->GetDesign_Variable(iDV) == PARABOLIC) { + surface_movement->SetParabolic(geometry, config); + } + + /*--- Design variable not implemented. ---*/ + + else { + if (rank == MASTER_NODE) cout << "Design Variable not implemented yet" << endl; + } + + /*--- Load the delta change in the design variable (finite difference step). ---*/ + + if ((config->GetDesign_Variable(iDV) != ANGLE_OF_ATTACK) && + (config->GetDesign_Variable(iDV) != FFD_ANGLE_OF_ATTACK)) { + /*--- If the Angle of attack is not involved, reset the value of the gradient. ---*/ + + my_Gradient = 0.0; + Gradient[iDV][0] = 0.0; + + if (MoveSurface) { + delta_eps = config->GetDV_Value(iDV); + + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) UpdatePoint[iPoint] = true; + + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if (config->GetMarker_All_DV(iMarker) == YES) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if ((iPoint < geometry->GetnPointDomain()) && UpdatePoint[iPoint]) { + Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); + Sensitivity = geometry->vertex[iMarker][iVertex]->GetAuxVar(); + + dS = 0.0; + for (iDim = 0; iDim < geometry->GetnDim(); iDim++) { + dS += Normal[iDim] * Normal[iDim]; + deps[iDim] = VarCoord[iDim] / delta_eps; + } + dS = sqrt(dS); + + dalpha_deps = 0.0; + for (iDim = 0; iDim < geometry->GetnDim(); iDim++) { + dalpha[iDim] = Normal[iDim] / dS; + dalpha_deps -= dalpha[iDim] * deps[iDim]; + } + + my_Gradient += Sensitivity * dalpha_deps; + UpdatePoint[iPoint] = false; + } + } + } + } + } + + SU2_MPI::Allreduce(&my_Gradient, &localGradient, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + Gradient[iDV][0] += localGradient; + } + } + + /*--- Delete memory for parameterization. ---*/ + + if (FFDBox != nullptr) { + for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) { + if (FFDBox[iFFDBox] != nullptr) { + delete FFDBox[iFFDBox]; + } + } + delete[] FFDBox; + } + + delete[] UpdatePoint; +} + +void CDiscAdjDeformationDriver::SetProjection_AD(CGeometry* geometry, CConfig* config, + CSurfaceMovement* surface_movement, su2double** Gradient) { + su2double *VarCoord = nullptr, Sensitivity, my_Gradient, localGradient, *Normal, Area = 0.0; + unsigned short iDV_Value = 0, iMarker, nMarker, iDim, nDim, iDV, nDV; + unsigned long iVertex, nVertex, iPoint; + + const int rank = SU2_MPI::GetRank(); + + nMarker = config->GetnMarker_All(); + nDim = geometry->GetnDim(); + nDV = config->GetnDV(); + + /*--- Discrete adjoint gradient computation. ---*/ + + if (rank == MASTER_NODE) + cout << endl + << "Evaluate functional gradient using Algorithmic Differentiation (ZONE " << config->GetiZone() << ")." + << endl; + + /*--- Start recording of operations. ---*/ + + AD::StartRecording(); + + /*--- Register design variables as input and set them to zero + * (since we want to have the derivative at alpha = 0, i.e. for the current design). ---*/ + + for (iDV = 0; iDV < nDV; iDV++) { + for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) { + config->SetDV_Value(iDV, iDV_Value, 0.0); + + AD::RegisterInput(config->GetDV_Value(iDV, iDV_Value)); + } + } + + /*--- Call the surface deformation routine. ---*/ + + surface_movement->SetSurface_Deformation(geometry, config); + + /*--- Stop the recording. --- */ + + AD::StopRecording(); + + /*--- Create a structure to identify points that have been already visited. + * We need that to make sure to set the sensitivity of surface points only once. + * Markers share points, so we would visit them more than once in the loop over the markers below). ---*/ + + vector visited(geometry->GetnPoint(), false); + + /*--- Initialize the derivatives of the output of the surface deformation routine + * with the discrete adjoints from the CFD solution. ---*/ + + for (iMarker = 0; iMarker < nMarker; iMarker++) { + if (config->GetMarker_All_DV(iMarker) == YES) { + nVertex = geometry->nVertex[iMarker]; + for (iVertex = 0; iVertex < nVertex; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if (!visited[iPoint]) { + VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); + Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + + Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Area += Normal[iDim] * Normal[iDim]; + } + Area = sqrt(Area); + + for (iDim = 0; iDim < nDim; iDim++) { + if (config->GetDiscrete_Adjoint()) { + Sensitivity = geometry->GetSensitivity(iPoint, iDim); + } else { + Sensitivity = -Normal[iDim] * geometry->vertex[iMarker][iVertex]->GetAuxVar() / Area; + } + SU2_TYPE::SetDerivative(VarCoord[iDim], SU2_TYPE::GetValue(Sensitivity)); + } + visited[iPoint] = true; + } + } + } + } + + /*--- Compute derivatives and extract gradient. ---*/ + + AD::ComputeAdjoint(); + + for (iDV = 0; iDV < nDV; iDV++) { + for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) { + my_Gradient = SU2_TYPE::GetDerivative(config->GetDV_Value(iDV, iDV_Value)); + + SU2_MPI::Allreduce(&my_Gradient, &localGradient, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + + /*--- Angle of Attack design variable (this is different, the value comes form the input file). ---*/ + + if ((config->GetDesign_Variable(iDV) == ANGLE_OF_ATTACK) || + (config->GetDesign_Variable(iDV) == FFD_ANGLE_OF_ATTACK)) { + Gradient[iDV][iDV_Value] = config->GetAoA_Sens(); + } + + Gradient[iDV][iDV_Value] += localGradient; + } + } + + AD::Reset(); +} + +void CDiscAdjDeformationDriver::OutputGradient(su2double** Gradient, CConfig* config, ofstream& Gradient_file) { + + unsigned short nDV, iDV, iDV_Value, nDV_Value; + + int rank = SU2_MPI::GetRank(); + + nDV = config->GetnDV(); + + /*--- Loop through all design variables and their gradients. ---*/ + + for (iDV = 0; iDV < nDV; iDV++) { + nDV_Value = config->GetnDV_Value(iDV); + if (rank == MASTER_NODE) { + /*--- Print the kind of design variable on screen. ---*/ + + cout << endl << "Design variable ("; + for (std::map::const_iterator it = Param_Map.begin(); it != Param_Map.end(); ++it) { + if (it->second == config->GetDesign_Variable(iDV)) { + cout << it->first << ") number " << iDV << "." << endl; + } + } + + /*--- Print the kind of objective function to screen. ---*/ + + for (std::map::const_iterator it = Objective_Map.begin(); it != Objective_Map.end(); ++it) { + if (it->second == config->GetKind_ObjFunc()) { + cout << it->first << " gradient : "; + if (iDV == 0) Gradient_file << it->first << " gradient " << endl; + } + } + + /*--- Print the gradient to file and screen. ---*/ + + for (iDV_Value = 0; iDV_Value < nDV_Value; iDV_Value++) { + cout << Gradient[iDV][iDV_Value]; + if (iDV_Value != nDV_Value - 1) { + cout << ", "; + } + Gradient_file << Gradient[iDV][iDV_Value] << endl; + } + cout << endl; + cout << "-------------------------------------------------------------------------" << endl; + } + } +} + +void CDiscAdjDeformationDriver::SetSensitivity_Files(CGeometry**** geometry, CConfig** config, + unsigned short val_nZone) { + unsigned short iMarker, iDim, nDim, nMarker, nVar; + unsigned long iVertex, iPoint, nPoint, nVertex; + su2double *Normal, Prod, Sens = 0.0, SensDim, Area; + + unsigned short iZone; + + CSolver* solver = nullptr; + COutput* output = nullptr; + + for (iZone = 0; iZone < val_nZone; iZone++) { + nPoint = geometry[iZone][INST_0][MESH_0]->GetnPoint(); + nDim = geometry[iZone][INST_0][MESH_0]->GetnDim(); + nMarker = config[iZone]->GetnMarker_All(); + nVar = nDim + 1; + + /*--- We create a baseline solver to easily merge the sensitivity information. ---*/ + + vector fieldnames; + fieldnames.push_back("\"Point\""); + fieldnames.push_back("\"x\""); + fieldnames.push_back("\"y\""); + if (nDim == 3) { + fieldnames.push_back("\"z\""); + } + fieldnames.push_back("\"Sensitivity_x\""); + fieldnames.push_back("\"Sensitivity_y\""); + if (nDim == 3) { + fieldnames.push_back("\"Sensitivity_z\""); + } + fieldnames.push_back("\"Surface_Sensitivity\""); + + solver = new CBaselineSolver(geometry[iZone][INST_0][MESH_0], config[iZone], nVar + nDim, fieldnames); + + for (iPoint = 0; iPoint < nPoint; iPoint++) { + for (iDim = 0; iDim < nDim; iDim++) { + solver->GetNodes()->SetSolution(iPoint, iDim, geometry[iZone][INST_0][MESH_0]->nodes->GetCoord(iPoint, iDim)); + solver->GetNodes()->SetSolution(iPoint, iDim + nDim, geometry[iZone][INST_0][MESH_0]->GetSensitivity(iPoint, iDim)); + } + } + + /*--- Compute the sensitivity in normal direction. ---*/ + + for (iMarker = 0; iMarker < nMarker; iMarker++) { + if (config[iZone]->GetSolid_Wall(iMarker) || (config[iZone]->GetMarker_All_DV(iMarker) == YES)) { + nVertex = geometry[iZone][INST_0][MESH_0]->GetnVertex(iMarker); + + for (iVertex = 0; iVertex < nVertex; iVertex++) { + iPoint = geometry[iZone][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + Normal = geometry[iZone][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNormal(); + Prod = 0.0; + Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + /*--- Retrieve the gradient calculated with discrete adjoint method. ---*/ + + SensDim = geometry[iZone][INST_0][MESH_0]->GetSensitivity(iPoint, iDim); + + /*--- Calculate scalar product for projection onto the normal vector. ---*/ + + Prod += Normal[iDim] * SensDim; + + Area += Normal[iDim] * Normal[iDim]; + } + + Area = sqrt(Area); + + /*--- Projection of the gradient onto the normal vector of the surface. ---*/ + + Sens = Prod / Area; + + solver->GetNodes()->SetSolution(iPoint, 2 * nDim, Sens); + } + } + } + + output = new CBaselineOutput(config[iZone], nDim, solver); + output->PreprocessVolumeOutput(config[iZone]); + output->PreprocessHistoryOutput(config[iZone], false); + + /*--- Load the data. --- */ + + output->Load_Data(geometry[iZone][INST_0][MESH_0], config[iZone], &solver); + + /*--- Set the surface filename. ---*/ + + output->SetSurface_Filename(config[iZone]->GetSurfSens_FileName()); + + /*--- Set the volume filename. ---*/ + + output->SetVolume_Filename(config[iZone]->GetVolSens_FileName()); + + /*--- Write to file. ---*/ + + for (unsigned short iFile = 0; iFile < config[iZone]->GetnVolumeOutputFiles(); iFile++) { + auto FileFormat = config[iZone]->GetVolumeOutputFiles(); + if (FileFormat[iFile] != OUTPUT_TYPE::RESTART_ASCII && FileFormat[iFile] != OUTPUT_TYPE::RESTART_BINARY && + FileFormat[iFile] != OUTPUT_TYPE::CSV) + output->WriteToFile(config[iZone], geometry[iZone][INST_0][MESH_0], FileFormat[iFile]); + } + + /*--- Free memory. ---*/ + + delete output; + delete solver; + } +} + +void CDiscAdjDeformationDriver::DerivativeTreatment_MeshSensitivity(CGeometry* geometry, CConfig* config, + CVolumetricMovement* grid_movement) { + int rank = SU2_MPI::GetRank(); + + /*--- Warning if chosen smoothing mode is unsupported. + * This is the default option if the user has not specified a mode in the config file. ---*/ + + if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::NONE) { + SU2_MPI::Error("Unsupported operation modus for the Sobolev Smoothing Solver.", CURRENT_FUNCTION); + } + + /*-- Construct the smoothing solver and numerics. ---*/ + + std::unique_ptr solver(new CGradientSmoothingSolver(geometry, config)); + unsigned dim = (config->GetSmoothOnSurface() ? geometry->GetnDim() - 1 : geometry->GetnDim()); + std::unique_ptr numerics(new CGradSmoothing(dim, config)); + + if (rank == MASTER_NODE) cout << "Sobolev Smoothing of derivatives is active." << endl; + + /*--- Apply the smoothing procedure on the mesh level. ---*/ + + if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::MESH_LEVEL) { + if (rank == MASTER_NODE) cout << " working on mesh level" << endl; + + /*--- Work with the surface derivatives. ---*/ + if (config->GetSmoothOnSurface()) { + /*--- Project to surface and then smooth on surface. ---*/ + + grid_movement->SetVolume_Deformation(geometry, config, false, true); + + /*--- Get the sensitivities from the geometry class to work with. ---*/ + + solver->ReadSensFromGeometry(geometry); + + /*--- Perform the smoothing procedure on all boundaries marked as DV marker. ---*/ + + solver->ApplyGradientSmoothingSurface(geometry, numerics.get(), config); + + /*--- After applying the solver write the results back. ---*/ + + solver->WriteSensToGeometry(geometry); + + /*--- Work with the volume derivatives. ---*/ + } else { + /*--- Get the sensitivities from the geometry class to work with. ---*/ + + solver->ReadSensFromGeometry(geometry); + + solver->ApplyGradientSmoothingVolume(geometry, numerics.get(), config); + + /*--- After applying the solver write the results back. ---*/ + + solver->WriteSensToGeometry(geometry); + + /*--- Projection is applied after smoothing. ---*/ + + grid_movement->SetVolume_Deformation(geometry, config, false, true); + } + } +} + +void CDiscAdjDeformationDriver::DerivativeTreatment_Gradient(CGeometry* geometry, CConfig* config, + CVolumetricMovement* grid_movement, + CSurfaceMovement* surface_movement, su2double** Gradient) { + int rank = SU2_MPI::GetRank(); + + /*--- Error if chosen smoothing mode is unsupported. + * This is the default option if the user has not specified a mode in the config file. ---*/ + + if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::NONE) { + SU2_MPI::Error("Unsupported operation modus for the Sobolev Smoothing Solver.", CURRENT_FUNCTION); + } + + /*-- Construct the smoothing solver and numerics. ---*/ + + std::unique_ptr solver(new CGradientSmoothingSolver(geometry, config)); + unsigned dim = (config->GetSmoothOnSurface() ? geometry->GetnDim() - 1 : geometry->GetnDim()); + std::unique_ptr numerics(new CGradSmoothing(dim, config)); + + if (rank == MASTER_NODE) cout << "Sobolev Smoothing of derivatives is active." << endl; + + /*--- Get the sensitivities from the geometry class to work with. ---*/ + + solver->ReadSensFromGeometry(geometry); + + /*--- Apply the smoothing procedure on the DV level. ---*/ + if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::PARAM_LEVEL_COMPLETE) { + solver->ApplyGradientSmoothingDV(geometry, numerics.get(), surface_movement, grid_movement, config, Gradient); + + /*--- If smoothing already took place on the mesh level, or none is requested, just do standard projection. ---*/ + } else if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::ONLY_GRAD || + config->GetSobMode() == ENUM_SOBOLEV_MODUS::MESH_LEVEL) { + solver->RecordTapeAndCalculateOriginalGradient(geometry, surface_movement, grid_movement, config, Gradient); + } +} diff --git a/SU2_DEF/src/meson.build b/SU2_DEF/src/meson.build index e4f3705111dd..3fcca7809d79 100644 --- a/SU2_DEF/src/meson.build +++ b/SU2_DEF/src/meson.build @@ -1,4 +1,9 @@ -su2_def_src = ['SU2_DEF.cpp'] +su2_def_include = include_directories('./') +su2_def_src = files([ + 'SU2_DEF.cpp', + 'drivers/CDeformationDriver.cpp', + 'drivers/CDiscAdjDeformationDriver.cpp' +]) if get_option('enable-normal') su2_def = executable('SU2_DEF', @@ -6,4 +11,27 @@ if get_option('enable-normal') install: true, dependencies: [su2_deps, common_dep, su2_cfd_dep], cpp_args :[default_warning_flags, su2_cpp_args]) + + su2_def_lib = static_library('SU2core_DEF', + su2_def_src, + install : false, + dependencies : [su2_deps, common_dep, su2_cfd_dep], + cpp_args: [default_warning_flags, su2_cpp_args]) + + su2_def_dep = declare_dependency(link_with: su2_def_lib, + include_directories: su2_def_include) + +endif + +if get_option('enable-autodiff') + + su2_def_lib_ad = static_library('SU2core_DEF_AD', + su2_def_src, + install : false, + dependencies : [su2_deps, codi_dep, commonAD_dep, su2_cfd_dep_ad], + cpp_args: [default_warning_flags, su2_cpp_args, codi_rev_args]) + + su2_def_dep_ad = declare_dependency(link_with: su2_def_lib_ad, + include_directories: su2_def_include) + endif diff --git a/SU2_DOT/include/SU2_DOT.hpp b/SU2_DOT/include/SU2_DOT.hpp index 0a0d9b2482bd..f6b6d39e5083 100644 --- a/SU2_DOT/include/SU2_DOT.hpp +++ b/SU2_DOT/include/SU2_DOT.hpp @@ -1,7 +1,6 @@ /*! * \file SU2_DOT.hpp * \brief Headers of the main subroutines of the code SU2_DOT. - * The subroutines and functions are in the SU2_DOT.cpp file. * \author F. Palacios, T. Economon * \version 7.5.1 "Blackbird" * @@ -26,87 +25,16 @@ * License along with SU2. If not, see . */ - #pragma once -#define ENABLE_MAPS -#include "../../Common/include/CConfig.hpp" -#undef ENABLE_MAPS - -#include "../../Common/include/parallelization/mpi_structure.hpp" -#include "../../Common/include/parallelization/omp_structure.hpp" - +#include #include -#include #include -#include - -#include "../../Common/include/geometry/CPhysicalGeometry.hpp" -#include "../../Common/include/fem/fem_geometry_structure.hpp" -#include "../../SU2_CFD/include/output/CBaselineOutput.hpp" -#include "../../SU2_CFD/include/solvers/CBaselineSolver.hpp" +#include -#include "../../SU2_CFD/include/solvers/CGradientSmoothingSolver.hpp" -#include "../../SU2_CFD/include/numerics/CGradSmoothing.hpp" +#include "../../Common/include/CConfig.hpp" +#include "../../Common/include/parallelization/mpi_structure.hpp" +#include "../../Common/include/parallelization/omp_structure.hpp" +#include "../../SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp" using namespace std; - - -/*! - * \brief Projection of the surface sensitivity using finite differences (FD). - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] surface_movement - Surface movement class of the problem. - * \param[in] Gradient_file - Output file to store the gradient data. - */ - -void SetProjection_FD(CGeometry *geometry, CConfig *config, CSurfaceMovement *surface_movement, su2double **Gradient); - -/*! - * \brief Projection of the surface sensitivity using algorithmic differentiation (AD). - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] surface_movement - Surface movement class of the problem. - * \param[in] Gradient_file - Output file to store the gradient data. - */ - -void SetProjection_AD(CGeometry *geometry, CConfig *config, CSurfaceMovement *surface_movement, su2double **Gradient); - -/*! - * \brief Prints the gradient information to a file. - * \param[in] Gradient - The gradient data. - * \param[in] config - Definition of the particular problem. - * \param[in] Gradient_file - Output file to store the gradient data. - */ - -void OutputGradient(su2double** Gradient, CConfig* config, ofstream& Gradient_file); - -/*! - * \brief Write the sensitivity (including mesh sensitivity) computed with the discrete adjoint method - * on the surface and in the volume to a file. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] val_nZone - Number of Zones. - */ - -void SetSensitivity_Files(CGeometry ***geometry, CConfig **config, unsigned short val_nZone); - -/*! - * \brief Treatment of derivatives with the Sobolev smoothing solver. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] grid_movement - Volumetric movement class of the problem. - */ - -void DerivativeTreatment_MeshSensitivity(CGeometry *geometry, CConfig *config, CVolumetricMovement *grid_movement); - -/*! - * \brief Treatment of derivatives with the Sobolev smoothing solver. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] grid_movement - Volumetric movement class of the problem. - * \param[in] surface_movement - Surface movement class of the problem. - * \param[in] Gradient - Output array to store the gradient data. - */ - -void DerivativeTreatment_Gradient(CGeometry *geometry, CConfig *config, CVolumetricMovement *grid_movement, CSurfaceMovement *surface_movement, su2double **Gradient); diff --git a/SU2_DOT/src/SU2_DOT.cpp b/SU2_DOT/src/SU2_DOT.cpp index 75e057d319dd..6c269be249c8 100644 --- a/SU2_DOT/src/SU2_DOT.cpp +++ b/SU2_DOT/src/SU2_DOT.cpp @@ -25,22 +25,19 @@ * License along with SU2. If not, see . */ - #include "../include/SU2_DOT.hpp" -using namespace std; -int main(int argc, char *argv[]) { +using namespace std; - unsigned short iZone, iInst; - su2double StartTime = 0.0, StopTime = 0.0, UsedTime = 0.0; +int main(int argc, char* argv[]) { char config_file_name[MAX_STRING_SIZE]; - /*--- OpenMP initialization ---*/ + /*--- Create a pointer to the main SU2_DEF Driver. ---*/ - omp_initialize(); + CDiscAdjDeformationDriver* driver = nullptr; - /*--- MPI initialization, and buffer setting ---*/ + /*--- MPI initialization. ---*/ #if defined(HAVE_OMP) && defined(HAVE_MPI) int provided; @@ -48,1009 +45,38 @@ int main(int argc, char *argv[]) { #else SU2_MPI::Init(&argc, &argv); #endif - SU2_MPI::Comm MPICommunicator = SU2_MPI::GetComm(); - - const int rank = SU2_MPI::GetRank(); - const int size = SU2_MPI::GetSize(); - - /*--- AD initialization ---*/ -#ifdef HAVE_OPDI - AD::getGlobalTape().initialize(); -#endif - - /*--- Pointer to different structures that will be used throughout the entire code ---*/ - - CConfig **config_container = nullptr; - CConfig *driver_config = nullptr; - CGeometry ***geometry_container = nullptr; - CSurfaceMovement **surface_movement = nullptr; - CVolumetricMovement **grid_movement = nullptr; - unsigned short *nInst = nullptr; - - /*--- Load in the number of zones and spatial dimensions in the mesh file (if no config - file is specified, default.cfg is used) ---*/ - - if (argc == 2) { strcpy(config_file_name,argv[1]); } - else { strcpy(config_file_name, "default.cfg"); } - - /*--- Read the name and format of the input mesh file to get from the mesh - file the number of zones and dimensions from the numerical grid (required - for variables allocation) ---*/ - - CConfig *config = nullptr; - config = new CConfig(config_file_name, SU2_COMPONENT::SU2_DOT); - - const auto nZone = config->GetnZone(); - - /*--- Definition of the containers per zones ---*/ - - config_container = new CConfig*[nZone] (); - geometry_container = new CGeometry**[nZone] (); - surface_movement = new CSurfaceMovement*[nZone] (); - grid_movement = new CVolumetricMovement*[nZone] (); - - nInst = new unsigned short[nZone]; - driver_config = nullptr; - - for (iZone = 0; iZone < nZone; iZone++) { - nInst[iZone] = 1; - } - - /*--- Initialize the configuration of the driver ---*/ - driver_config = new CConfig(config_file_name, SU2_COMPONENT::SU2_DOT, false); - - /*--- Initialize a char to store the zone filename ---*/ - char zone_file_name[MAX_STRING_SIZE]; - - /*--- Loop over all zones to initialize the various classes. In most - cases, nZone is equal to one. This represents the solution of a partial - differential equation on a single block, unstructured mesh. ---*/ - - for (iZone = 0; iZone < nZone; iZone++) { - - /*--- Definition of the configuration option class for all zones. In this - constructor, the input configuration file is parsed and all options are - read and stored. ---*/ - - if (driver_config->GetnConfigFiles() > 0){ - strcpy(zone_file_name, driver_config->GetConfigFilename(iZone).c_str()); - config_container[iZone] = new CConfig(driver_config, zone_file_name, SU2_COMPONENT::SU2_DOT, iZone, nZone, true); - } - else{ - config_container[iZone] = new CConfig(driver_config, config_file_name, SU2_COMPONENT::SU2_DOT, iZone, nZone, true); - } - config_container[iZone]->SetMPICommunicator(MPICommunicator); - - } - - /*--- Set the multizone part of the problem. ---*/ - if (driver_config->GetMultizone_Problem()){ - for (iZone = 0; iZone < nZone; iZone++) { - /*--- Set the interface markers for multizone ---*/ - config_container[iZone]->SetMultizone(driver_config, config_container); - } - } - - /*--- Loop over all zones to initialize the various classes. In most - cases, nZone is equal to one. This represents the solution of a partial - differential equation on a single block, unstructured mesh. ---*/ - - for (iZone = 0; iZone < nZone; iZone++) { - - /*--- Determine whether or not the FEM solver is used, which decides the - type of geometry classes that are instantiated. ---*/ - const bool fem_solver = config_container[iZone]->GetFEMSolver(); - - /*--- Read the number of instances for each zone ---*/ - - nInst[iZone] = config_container[iZone]->GetnTimeInstances(); - - geometry_container[iZone] = new CGeometry*[nInst[iZone]]; - - for (iInst = 0; iInst < nInst[iZone]; iInst++){ - - /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/ - - CGeometry *geometry_aux = nullptr; - - /*--- All ranks process the grid and call ParMETIS for partitioning ---*/ - - geometry_aux = new CPhysicalGeometry(config_container[iZone], iZone, nZone); - - /*--- Color the initial grid and set the send-receive domains (ParMETIS) ---*/ - - if ( fem_solver ) geometry_aux->SetColorFEMGrid_Parallel(config_container[iZone]); - else geometry_aux->SetColorGrid_Parallel(config_container[iZone]); - - /*--- Build the grid data structures using the ParMETIS coloring. ---*/ - - if( fem_solver ) { - switch( config_container[iZone]->GetKind_FEM_Flow() ) { - case DG: { - geometry_container[iZone][iInst] = new CMeshFEM_DG(geometry_aux, config_container[iZone]); - break; - } - } - } - else { - geometry_container[iZone][iInst] = new CPhysicalGeometry(geometry_aux, config_container[iZone]); - } - - /*--- Deallocate the memory of geometry_aux ---*/ - - delete geometry_aux; - - /*--- Add the Send/Receive boundaries ---*/ - - geometry_container[iZone][iInst]->SetSendReceive(config_container[iZone]); - - /*--- Add the Send/Receive boundaries ---*/ - - geometry_container[iZone][iInst]->SetBoundaries(config_container[iZone]); - - /*--- Create the vertex structure (required for MPI) ---*/ - - if (rank == MASTER_NODE) cout << "Identify vertices." << endl; - geometry_container[iZone][iInst]->SetVertex(config_container[iZone]); - - /*--- Store the global to local mapping after preprocessing. ---*/ - - if (rank == MASTER_NODE) cout << "Storing a mapping from global to local point index." << endl; - geometry_container[iZone][iInst]->SetGlobal_to_Local_Point(); - - /* Test for a fem solver, because some more work must be done. */ - - if (fem_solver) { - - /*--- Carry out a dynamic cast to CMeshFEM_DG, such that it is not needed to - define all virtual functions in the base class CGeometry. ---*/ - CMeshFEM_DG *DGMesh = dynamic_cast(geometry_container[iZone][iInst]); - - /*--- Determine the standard elements for the volume elements. ---*/ - if (rank == MASTER_NODE) cout << "Creating standard volume elements." << endl; - DGMesh->CreateStandardVolumeElements(config_container[iZone]); - - /*--- Create the face information needed to compute the contour integral - for the elements in the Discontinuous Galerkin formulation. ---*/ - if (rank == MASTER_NODE) cout << "Creating face information." << endl; - DGMesh->CreateFaces(config_container[iZone]); - } - } - } - - /*--- Set up a timer for performance benchmarking (preprocessing time is included) ---*/ - - StartTime = SU2_MPI::Wtime(); - - for (iZone = 0; iZone < nZone; iZone++) { - - if (rank == MASTER_NODE) - cout << "\n----------------------- Preprocessing computations ----------------------" << endl; - - /*--- Compute elements surrounding points, points surrounding points ---*/ - - if (rank == MASTER_NODE) cout << "Setting local point connectivity." << endl; - geometry_container[iZone][INST_0]->SetPoint_Connectivity(); - - /*--- Check the orientation before computing geometrical quantities ---*/ - - geometry_container[iZone][INST_0]->SetBoundVolume(); - if (config_container[iZone]->GetReorientElements()) { - if (rank == MASTER_NODE) cout << "Checking the numerical grid orientation of the elements." << endl; - geometry_container[iZone][INST_0]->Check_IntElem_Orientation(config_container[iZone]); - geometry_container[iZone][INST_0]->Check_BoundElem_Orientation(config_container[iZone]); - } - - /*--- Create the edge structure ---*/ - - if (rank == MASTER_NODE) cout << "Identify edges and vertices." << endl; - geometry_container[iZone][INST_0]->SetEdges(); geometry_container[iZone][INST_0]->SetVertex(config_container[iZone]); - - /*--- Create the dual control volume structures ---*/ - - if (rank == MASTER_NODE) cout << "Setting the bound control volume structure." << endl; - geometry_container[iZone][INST_0]->SetBoundControlVolume(config_container[ZONE_0], ALLOCATE); - - /*--- Store the global to local mapping after preprocessing. ---*/ - - if (rank == MASTER_NODE) cout << "Storing a mapping from global to local point index." << endl; - geometry_container[iZone][INST_0]->SetGlobal_to_Local_Point(); - - /*--- Create the point-to-point MPI communication structures. ---*/ - - geometry_container[iZone][INST_0]->PreprocessP2PComms(geometry_container[iZone][INST_0], config_container[iZone]); - - /*--- Load the surface sensitivities from file. This is done only - once: if this is an unsteady problem, a time-average of the surface - sensitivities at each node is taken within this routine. ---*/ - - if (!config_container[iZone]->GetDiscrete_Adjoint()) { - if (rank == MASTER_NODE) cout << "Reading surface sensitivities at each node from file." << endl; - geometry_container[iZone][INST_0]->SetBoundSensitivity(config_container[iZone]); - } - else { - if (rank == MASTER_NODE) cout << "Reading volume sensitivities at each node from file." << endl; - grid_movement[iZone] = new CVolumetricMovement(geometry_container[iZone][INST_0], config_container[iZone]); - - /*--- Read in sensitivities from file. ---*/ - if (config_container[ZONE_0]->GetSensitivity_Format() == UNORDERED_ASCII) - geometry_container[iZone][INST_0]->ReadUnorderedSensitivity(config_container[iZone]); - else - geometry_container[iZone][INST_0]->SetSensitivity(config_container[iZone]); - - if (rank == MASTER_NODE) - cout << "\n---------------------- Mesh sensitivity computation ---------------------" << endl; - - if(config_container[iZone]->GetDiscrete_Adjoint() && config_container[iZone]->GetSmoothGradient() && - config_container[iZone]->GetSobMode() == ENUM_SOBOLEV_MODUS::MESH_LEVEL) { - DerivativeTreatment_MeshSensitivity(geometry_container[iZone][INST_0], config_container[iZone], grid_movement[iZone]); - } else { - grid_movement[iZone]->SetVolume_Deformation(geometry_container[iZone][INST_0], config_container[iZone], false, true); - } + SU2_MPI::Comm comm = SU2_MPI::GetComm(); - } - } - - - if (config_container[ZONE_0]->GetDiscrete_Adjoint()) { - if (rank == MASTER_NODE) - cout << "\n------------------------ Mesh sensitivity Output ------------------------" << endl; - SetSensitivity_Files(geometry_container, config_container, nZone); - } - - /*--- Initialize structure to store the gradient ---*/ - su2double** Gradient = new su2double*[config_container[ZONE_0]->GetnDV()]; - - for (auto iDV = 0u; iDV < config_container[ZONE_0]->GetnDV(); iDV++) { - /*--- Initialize to zero ---*/ - Gradient[iDV] = new su2double[config_container[ZONE_0]->GetnDV_Value(iDV)](); - } - - ofstream Gradient_file; - Gradient_file.precision(config->OptionIsSet("OUTPUT_PRECISION") ? config->GetOutput_Precision() : 6); - - /*--- For multizone computations the gradient contributions are summed up and written into one file. ---*/ - for (iZone = 0; iZone < nZone; iZone++){ - if ((config_container[iZone]->GetDesign_Variable(0) != NONE) && - (config_container[iZone]->GetDesign_Variable(0) != SURFACE_FILE)) { - - if (rank == MASTER_NODE) - cout << "\n---------- Start gradient evaluation using sensitivity information ----------" << endl; - - /*--- Definition of the Class for surface deformation ---*/ - - surface_movement[iZone] = new CSurfaceMovement(); - - /*--- Copy coordinates to the surface structure ---*/ - - surface_movement[iZone]->CopyBoundary(geometry_container[iZone][INST_0], config_container[iZone]); - - /*--- If AD mode is enabled we can use it to compute the projection, - * otherwise we use finite differences. ---*/ - - if (config_container[iZone]->GetAD_Mode()) { - if (config_container[iZone]->GetSmoothGradient()) { - DerivativeTreatment_Gradient(geometry_container[iZone][INST_0], config_container[iZone], grid_movement[iZone], surface_movement[iZone] , Gradient); - } else { - SetProjection_AD(geometry_container[iZone][INST_0], config_container[iZone], surface_movement[iZone] , Gradient); - } - } else { - SetProjection_FD(geometry_container[iZone][INST_0], config_container[iZone], surface_movement[iZone] , Gradient); - } - } - } // for iZone - - /*--- Write the gradient to a file ---*/ - - if (rank == MASTER_NODE) - Gradient_file.open(config_container[ZONE_0]->GetObjFunc_Grad_FileName().c_str(), ios::out); - - /*--- Print gradients to screen and writes to file ---*/ - - OutputGradient(Gradient, config_container[ZONE_0], Gradient_file); + /*--- Load in the number of zones and spatial dimensions in the mesh file + (if no config file is specified, default.cfg is used). ---*/ - for (auto iDV = 0u; iDV < config_container[ZONE_0]->GetnDV(); iDV++){ - delete [] Gradient[iDV]; + if (argc == 2) { + strcpy(config_file_name, argv[1]); + } else { + strcpy(config_file_name, "default.cfg"); } - delete [] Gradient; - delete config; - config = nullptr; + /*--- Initialize the mesh deformation driver. ---*/ - if (rank == MASTER_NODE) - cout << "\n------------------------- Solver Postprocessing -------------------------" << endl; + driver = new CDiscAdjDeformationDriver(config_file_name, comm); - if (geometry_container != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (geometry_container[iZone] != nullptr) { - for (iInst = 0; iInst < nInst[iZone]; iInst++){ - if (geometry_container[iZone][iInst] != nullptr) { - delete geometry_container[iZone][iInst]; - } - } - delete geometry_container[iZone]; - } - } - delete [] geometry_container; - } - if (rank == MASTER_NODE) cout << "Deleted CGeometry container." << endl; - - if (surface_movement != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (surface_movement[iZone] != nullptr) { - delete surface_movement[iZone]; - } - } - delete [] surface_movement; - } - if (rank == MASTER_NODE) cout << "Deleted CSurfaceMovement class." << endl; - - if (grid_movement != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (grid_movement[iZone] != nullptr) { - delete grid_movement[iZone]; - } - } - delete [] grid_movement; - } - if (rank == MASTER_NODE) cout << "Deleted CVolumetricMovement class." << endl; + /*--- Preprocess the solver data. ---*/ - delete config; - config = nullptr; - if (config_container != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (config_container[iZone] != nullptr) { - delete config_container[iZone]; - } - } - delete [] config_container; - } - if (rank == MASTER_NODE) cout << "Deleted CConfig container." << endl; - - /*--- Synchronization point after a single solver iteration. Compute the - wall clock time required. ---*/ + driver->Preprocess(); - StopTime = SU2_MPI::Wtime(); + /*--- Launch the main external loop of the solver. ---*/ - /*--- Compute/print the total time for performance benchmarking. ---*/ + driver->Run(); - UsedTime = StopTime-StartTime; - if (rank == MASTER_NODE) { - cout << "\nCompleted in " << fixed << UsedTime << " seconds on "<< size; - if (size == 1) cout << " core." << endl; else cout << " cores." << endl; - } + /*--- Postprocess all the containers, close history file, and exit SU2. ---*/ - /*--- Exit the solver cleanly ---*/ + driver->Postprocessing(); - if (rank == MASTER_NODE) - cout << "\n------------------------- Exit Success (SU2_DOT) ------------------------\n" << endl; - - /*--- Finalize AD, if necessary. ---*/ -#ifdef HAVE_OPDI - AD::getGlobalTape().finalize(); -#endif + delete driver; /*--- Finalize MPI parallelization. ---*/ - SU2_MPI::Finalize(); - /*--- Finalize OpenMP. ---*/ - omp_finalize(); + SU2_MPI::Finalize(); return EXIT_SUCCESS; - -} - -void SetProjection_FD(CGeometry *geometry, CConfig *config, CSurfaceMovement *surface_movement, su2double** Gradient){ - - unsigned short iDV, nDV, iFFDBox, nDV_Value, iMarker, iDim; - unsigned long iVertex, iPoint; - su2double delta_eps, my_Gradient, localGradient, *Normal, dS, *VarCoord, Sensitivity, - dalpha[3], deps[3], dalpha_deps; - bool *UpdatePoint, MoveSurface, Local_MoveSurface; - CFreeFormDefBox **FFDBox; - - int rank = SU2_MPI::GetRank(); - - nDV = config->GetnDV(); - - /*--- Boolean controlling points to be updated ---*/ - - UpdatePoint = new bool[geometry->GetnPoint()]; - - /*--- Definition of the FFD deformation class ---*/ - - unsigned short nFFDBox = MAX_NUMBER_FFD; - FFDBox = new CFreeFormDefBox*[nFFDBox]; - for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) FFDBox[iFFDBox] = nullptr; - - for (iDV = 0; iDV < nDV; iDV++){ - nDV_Value = config->GetnDV_Value(iDV); - if (nDV_Value != 1){ - SU2_MPI::Error("The projection using finite differences currently only supports a fixed direction of movement for FFD points.", CURRENT_FUNCTION); - } - } - - /*--- Continuous adjoint gradient computation ---*/ - - if (rank == MASTER_NODE) - cout << "Evaluate functional gradient using Finite Differences." << endl; - - for (iDV = 0; iDV < nDV; iDV++) { - - MoveSurface = true; - Local_MoveSurface = true; - - /*--- Free Form deformation based ---*/ - - if ((config->GetDesign_Variable(iDV) == FFD_CONTROL_POINT_2D) || - (config->GetDesign_Variable(iDV) == FFD_CAMBER_2D) || - (config->GetDesign_Variable(iDV) == FFD_THICKNESS_2D) || - (config->GetDesign_Variable(iDV) == FFD_TWIST_2D) || - (config->GetDesign_Variable(iDV) == FFD_CONTROL_POINT) || - (config->GetDesign_Variable(iDV) == FFD_NACELLE) || - (config->GetDesign_Variable(iDV) == FFD_GULL) || - (config->GetDesign_Variable(iDV) == FFD_TWIST) || - (config->GetDesign_Variable(iDV) == FFD_ROTATION) || - (config->GetDesign_Variable(iDV) == FFD_CAMBER) || - (config->GetDesign_Variable(iDV) == FFD_THICKNESS) || - (config->GetDesign_Variable(iDV) == FFD_ANGLE_OF_ATTACK)) { - - /*--- Read the FFD information in the first iteration ---*/ - - if (iDV == 0) { - - if (rank == MASTER_NODE) - cout << "Read the FFD information from mesh file." << endl; - - /*--- Read the FFD information from the grid file ---*/ - - surface_movement->ReadFFDInfo(geometry, config, FFDBox, config->GetMesh_FileName()); - - /*--- If the FFDBox was not defined in the input file ---*/ - if (!surface_movement->GetFFDBoxDefinition()) { - SU2_MPI::Error("The input grid doesn't have the entire FFD information!", CURRENT_FUNCTION); - } - - for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) { - - if (rank == MASTER_NODE) cout << "Checking FFD box dimension." << endl; - surface_movement->CheckFFDDimension(geometry, config, FFDBox[iFFDBox], iFFDBox); - - if (rank == MASTER_NODE) cout << "Check the FFD box intersections with the solid surfaces." << endl; - surface_movement->CheckFFDIntersections(geometry, config, FFDBox[iFFDBox], iFFDBox); - - } - - if (rank == MASTER_NODE) - cout <<"-------------------------------------------------------------------------" << endl; - - } - - if (rank == MASTER_NODE) { - cout << endl << "Design variable number "<< iDV <<"." << endl; - cout << "Performing 3D deformation of the surface." << endl; - } - - /*--- Apply the control point change ---*/ - - MoveSurface = false; - - for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) { - - /*--- Reset FFD box ---*/ - - switch (config->GetDesign_Variable(iDV) ) { - case FFD_CONTROL_POINT_2D : Local_MoveSurface = surface_movement->SetFFDCPChange_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_CAMBER_2D : Local_MoveSurface = surface_movement->SetFFDCamber_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_THICKNESS_2D : Local_MoveSurface = surface_movement->SetFFDThickness_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_TWIST_2D : Local_MoveSurface = surface_movement->SetFFDTwist_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_CONTROL_POINT : Local_MoveSurface = surface_movement->SetFFDCPChange(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_NACELLE : Local_MoveSurface = surface_movement->SetFFDNacelle(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_GULL : Local_MoveSurface = surface_movement->SetFFDGull(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_TWIST : Local_MoveSurface = surface_movement->SetFFDTwist(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_ROTATION : Local_MoveSurface = surface_movement->SetFFDRotation(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_CAMBER : Local_MoveSurface = surface_movement->SetFFDCamber(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_THICKNESS : Local_MoveSurface = surface_movement->SetFFDThickness(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_CONTROL_SURFACE : Local_MoveSurface = surface_movement->SetFFDControl_Surface(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, true); break; - case FFD_ANGLE_OF_ATTACK : Gradient[iDV][0] = config->GetAoA_Sens(); break; - } - - /*--- Recompute cartesian coordinates using the new control points position ---*/ - - if (Local_MoveSurface) { - MoveSurface = true; - surface_movement->SetCartesianCoord(geometry, config, FFDBox[iFFDBox], iFFDBox, true); - } - - } - - } - - /*--- Hicks Henne design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == HICKS_HENNE) { - surface_movement->SetHicksHenne(geometry, config, iDV, true); - } - - /*--- Surface bump design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == SURFACE_BUMP) { - surface_movement->SetSurface_Bump(geometry, config, iDV, true); - } - - /*--- Kulfan (CST) design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == CST) { - surface_movement->SetCST(geometry, config, iDV, true); - } - - /*--- Displacement design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == TRANSLATION) { - surface_movement->SetTranslation(geometry, config, iDV, true); - } - - /*--- Angle of Attack design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == ANGLE_OF_ATTACK) { - Gradient[iDV][0] = config->GetAoA_Sens(); - } - - /*--- Scale design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == SCALE) { - surface_movement->SetScale(geometry, config, iDV, true); - } - - /*--- Rotation design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == ROTATION) { - surface_movement->SetRotation(geometry, config, iDV, true); - } - - /*--- NACA_4Digits design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == NACA_4DIGITS) { - surface_movement->SetNACA_4Digits(geometry, config); - } - - /*--- Parabolic design variable ---*/ - - else if (config->GetDesign_Variable(iDV) == PARABOLIC) { - surface_movement->SetParabolic(geometry, config); - } - - /*--- Design variable not implemented ---*/ - - else { - if (rank == MASTER_NODE) - cout << "Design Variable not implemented yet" << endl; - } - - /*--- Load the delta change in the design variable (finite difference step). ---*/ - - if ((config->GetDesign_Variable(iDV) != ANGLE_OF_ATTACK) && - (config->GetDesign_Variable(iDV) != FFD_ANGLE_OF_ATTACK)) { - - /*--- If the Angle of attack is not involved, reset the value of the gradient ---*/ - - my_Gradient = 0.0; Gradient[iDV][0] = 0.0; - - if (MoveSurface) { - - delta_eps = config->GetDV_Value(iDV); - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) - UpdatePoint[iPoint] = true; - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if (config->GetMarker_All_DV(iMarker) == YES) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - if ((iPoint < geometry->GetnPointDomain()) && UpdatePoint[iPoint]) { - - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); - Sensitivity = geometry->vertex[iMarker][iVertex]->GetAuxVar(); - - dS = 0.0; - for (iDim = 0; iDim < geometry->GetnDim(); iDim++) { - dS += Normal[iDim]*Normal[iDim]; - deps[iDim] = VarCoord[iDim] / delta_eps; - } - dS = sqrt(dS); - - dalpha_deps = 0.0; - for (iDim = 0; iDim < geometry->GetnDim(); iDim++) { - dalpha[iDim] = Normal[iDim] / dS; - dalpha_deps -= dalpha[iDim]*deps[iDim]; - } - - my_Gradient += Sensitivity*dalpha_deps; - UpdatePoint[iPoint] = false; - } - } - } - } - } - - SU2_MPI::Allreduce(&my_Gradient, &localGradient, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - Gradient[iDV][0] += localGradient; - } - } - - /*--- Delete memory for parameterization. ---*/ - - if (FFDBox != nullptr) { - for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) { - if (FFDBox[iFFDBox] != nullptr) { - delete FFDBox[iFFDBox]; - } - } - delete [] FFDBox; - } - - delete [] UpdatePoint; - -} - - -void SetProjection_AD(CGeometry *geometry, CConfig *config, CSurfaceMovement *surface_movement, su2double** Gradient){ - - su2double *VarCoord = nullptr, Sensitivity, my_Gradient, localGradient, *Normal, Area = 0.0; - unsigned short iDV_Value = 0, iMarker, nMarker, iDim, nDim, iDV, nDV; - unsigned long iVertex, nVertex, iPoint; - - const int rank = SU2_MPI::GetRank(); - - nMarker = config->GetnMarker_All(); - nDim = geometry->GetnDim(); - nDV = config->GetnDV(); - - /*--- Discrete adjoint gradient computation ---*/ - - if (rank == MASTER_NODE) - cout << endl << "Evaluate functional gradient using Algorithmic Differentiation (ZONE " << config->GetiZone() << ")." << endl; - - /*--- Start recording of operations ---*/ - - AD::StartRecording(); - - /*--- Register design variables as input and set them to zero - * (since we want to have the derivative at alpha = 0, i.e. for the current design) ---*/ - - for (iDV = 0; iDV < nDV; iDV++){ - - for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++){ - - config->SetDV_Value(iDV, iDV_Value, 0.0); - - AD::RegisterInput(config->GetDV_Value(iDV, iDV_Value)); - } - } - - /*--- Call the surface deformation routine ---*/ - - surface_movement->SetSurface_Deformation(geometry, config); - - /*--- Stop the recording --- */ - - AD::StopRecording(); - - /*--- Create a structure to identify points that have been already visited. - * We need that to make sure to set the sensitivity of surface points only once - * (Markers share points, so we would visit them more than once in the loop over the markers below) ---*/ - - vector visited(geometry->GetnPoint(), false); - - /*--- Initialize the derivatives of the output of the surface deformation routine - * with the discrete adjoints from the CFD solution ---*/ - - for (iMarker = 0; iMarker < nMarker; iMarker++) { - if (config->GetMarker_All_DV(iMarker) == YES) { - nVertex = geometry->nVertex[iMarker]; - for (iVertex = 0; iVertex vertex[iMarker][iVertex]->GetNode(); - if (!visited[iPoint]){ - VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++){ - Area += Normal[iDim]*Normal[iDim]; - } - Area = sqrt(Area); - - for (iDim = 0; iDim < nDim; iDim++){ - if (config->GetDiscrete_Adjoint()){ - Sensitivity = geometry->GetSensitivity(iPoint, iDim); - } else { - Sensitivity = -Normal[iDim]*geometry->vertex[iMarker][iVertex]->GetAuxVar()/Area; - } - SU2_TYPE::SetDerivative(VarCoord[iDim], SU2_TYPE::GetValue(Sensitivity)); - } - visited[iPoint] = true; - } - } - } - } - - /*--- Compute derivatives and extract gradient ---*/ - - AD::ComputeAdjoint(); - - for (iDV = 0; iDV < nDV; iDV++){ - - for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++){ - - my_Gradient = SU2_TYPE::GetDerivative(config->GetDV_Value(iDV, iDV_Value)); - - SU2_MPI::Allreduce(&my_Gradient, &localGradient, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - - /*--- Angle of Attack design variable (this is different, - the value comes form the input file) ---*/ - - if ((config->GetDesign_Variable(iDV) == ANGLE_OF_ATTACK) || - (config->GetDesign_Variable(iDV) == FFD_ANGLE_OF_ATTACK)) { - Gradient[iDV][iDV_Value] = config->GetAoA_Sens(); - } - - Gradient[iDV][iDV_Value] += localGradient; - } - } - - AD::Reset(); - -} - -void OutputGradient(su2double** Gradient, CConfig* config, ofstream& Gradient_file){ - - unsigned short nDV, iDV, iDV_Value, nDV_Value; - - int rank = SU2_MPI::GetRank(); - - nDV = config->GetnDV(); - - /*--- Loop through all design variables and their gradients ---*/ - - for (iDV = 0; iDV < nDV; iDV++){ - nDV_Value = config->GetnDV_Value(iDV); - if (rank == MASTER_NODE){ - - /*--- Print the kind of design variable on screen ---*/ - - cout << endl << "Design variable ("; - for (std::map::const_iterator it = Param_Map.begin(); it != Param_Map.end(); ++it ){ - if (it->second == config->GetDesign_Variable(iDV)){ - cout << it->first << ") number "<< iDV << "." << endl; - } - } - - /*--- Print the kind of objective function to screen ---*/ - - for (std::map::const_iterator it = Objective_Map.begin(); it != Objective_Map.end(); ++it ){ - if (it->second == config->GetKind_ObjFunc()){ - cout << it->first << " gradient : "; - if (iDV == 0) Gradient_file << it->first << " gradient " << endl; - } - } - - /*--- Print the gradient to file and screen ---*/ - - for (iDV_Value = 0; iDV_Value < nDV_Value; iDV_Value++){ - cout << Gradient[iDV][iDV_Value]; - if (iDV_Value != nDV_Value-1 ){ - cout << ", "; - } - Gradient_file << Gradient[iDV][iDV_Value] << endl; - } - cout << endl; - cout <<"-------------------------------------------------------------------------" << endl; - } - } -} - - -void SetSensitivity_Files(CGeometry ***geometry, CConfig **config, unsigned short val_nZone) { - - unsigned short iMarker,iDim, nDim, nMarker, nVar; - unsigned long iVertex, iPoint, nPoint, nVertex; - su2double *Normal, Prod, Sens = 0.0, SensDim, Area; - - unsigned short iZone; - - CSolver *solver = nullptr; - COutput *output = nullptr; - - for (iZone = 0; iZone < val_nZone; iZone++) { - - nPoint = geometry[iZone][INST_0]->GetnPoint(); - nDim = geometry[iZone][INST_0]->GetnDim(); - nMarker = config[iZone]->GetnMarker_All(); - nVar = nDim + 1; - - /*--- We create a baseline solver to easily merge the sensitivity information ---*/ - - vector fieldnames; - fieldnames.push_back("\"Point\""); - fieldnames.push_back("\"x\""); - fieldnames.push_back("\"y\""); - if (nDim == 3) { - fieldnames.push_back("\"z\""); - } - fieldnames.push_back("\"Sensitivity_x\""); - fieldnames.push_back("\"Sensitivity_y\""); - if (nDim == 3) { - fieldnames.push_back("\"Sensitivity_z\""); - } - fieldnames.push_back("\"Surface_Sensitivity\""); - - solver = new CBaselineSolver(geometry[iZone][INST_0], config[iZone], nVar+nDim, fieldnames); - - for (iPoint = 0; iPoint < nPoint; iPoint++) { - for (iDim = 0; iDim < nDim; iDim++) { - solver->GetNodes()->SetSolution(iPoint, iDim, geometry[iZone][INST_0]->nodes->GetCoord(iPoint, iDim)); - solver->GetNodes()->SetSolution(iPoint, iDim+nDim, geometry[iZone][INST_0]->GetSensitivity(iPoint, iDim)); - } - } - - /*--- Compute the sensitivity in normal direction ---*/ - - for (iMarker = 0; iMarker < nMarker; iMarker++) { - - if(config[iZone]->GetSolid_Wall(iMarker) || (config[iZone]->GetMarker_All_DV(iMarker) == YES )) { - - nVertex = geometry[iZone][INST_0]->GetnVertex(iMarker); - - for (iVertex = 0; iVertex < nVertex; iVertex++) { - iPoint = geometry[iZone][INST_0]->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry[iZone][INST_0]->vertex[iMarker][iVertex]->GetNormal(); - Prod = 0.0; - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - - /*--- Retrieve the gradient calculated with discrete adjoint method ---*/ - - SensDim = geometry[iZone][INST_0]->GetSensitivity(iPoint, iDim); - - /*--- Calculate scalar product for projection onto the normal vector ---*/ - - Prod += Normal[iDim]*SensDim; - - Area += Normal[iDim]*Normal[iDim]; - } - - Area = sqrt(Area); - - /*--- Projection of the gradient onto the normal vector of the surface ---*/ - - Sens = Prod/Area; - - solver->GetNodes()->SetSolution(iPoint, 2*nDim, Sens); - - } - } - } - - output = new CBaselineOutput(config[iZone], nDim, solver); - output->PreprocessVolumeOutput(config[iZone]); - output->PreprocessHistoryOutput(config[iZone], false); - - /*--- Load the data --- */ - - output->Load_Data(geometry[iZone][INST_0], config[iZone], &solver); - - /*--- Set the surface filename ---*/ - - output->SetSurface_Filename(config[iZone]->GetSurfSens_FileName()); - - /*--- Set the volume filename ---*/ - - output->SetVolume_Filename(config[iZone]->GetVolSens_FileName()); - - /*--- Write to file ---*/ - - for (unsigned short iFile = 0; iFile < config[iZone]->GetnVolumeOutputFiles(); iFile++){ - auto FileFormat = config[iZone]->GetVolumeOutputFiles(); - if (FileFormat[iFile] != OUTPUT_TYPE::RESTART_ASCII && - FileFormat[iFile] != OUTPUT_TYPE::RESTART_BINARY && - FileFormat[iFile] != OUTPUT_TYPE::CSV) - output->WriteToFile(config[iZone], geometry[iZone][INST_0], FileFormat[iFile]); - } - - /*--- Free memory ---*/ - - delete output; - delete solver; - - } - -} - -void DerivativeTreatment_MeshSensitivity(CGeometry *geometry, CConfig *config, CVolumetricMovement *grid_movement) { - - int rank = SU2_MPI::GetRank(); - - /*--- Warning if choosen smoothing mode is unsupported. - * This is the default option if the user has not specified a mode in the config file. ---*/ - if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::NONE) { - SU2_MPI::Error("Unsupported operation modus for the Sobolev Smoothing Solver.", CURRENT_FUNCTION); - } - - /*-- Construct the smoothing solver and numerics ---*/ - std::unique_ptr solver(new CGradientSmoothingSolver(geometry, config)); - unsigned dim = (config->GetSmoothOnSurface() ? geometry->GetnDim() - 1 : geometry->GetnDim()); - std::unique_ptr numerics(new CGradSmoothing(dim, config)); - - if (rank == MASTER_NODE) cout << "Sobolev Smoothing of derivatives is active." << endl; - - /*--- Apply the smoothing procedure on the mesh level. ---*/ - if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::MESH_LEVEL) { - if (rank == MASTER_NODE) cout << " working on mesh level" << endl; - - /*--- Work with the surface derivatives. ---*/ - if (config->GetSmoothOnSurface()) { - - /*--- Project to surface and then smooth on surface. ---*/ - grid_movement->SetVolume_Deformation(geometry, config, false, true); - - /*--- Get the sensitivities from the geometry class to work with. ---*/ - solver->ReadSensFromGeometry(geometry); - - /*--- Perform the smoothing procedure on all boundaries marked as DV marker. ---*/ - solver->ApplyGradientSmoothingSurface(geometry, numerics.get(), config); - - /*--- After appling the solver write the results back ---*/ - solver->WriteSensToGeometry(geometry); - - /*--- Work with the volume derivatives. ---*/ - } else { - - /*--- Get the sensitivities from the geometry class to work with. ---*/ - solver->ReadSensFromGeometry(geometry); - - solver->ApplyGradientSmoothingVolume(geometry, numerics.get(), config); - - /*--- After appling the solver write the results back ---*/ - solver->WriteSensToGeometry(geometry); - - /*--- Projection is applied after smoothing. ---*/ - grid_movement->SetVolume_Deformation(geometry, config, false, true); - } - - } - -} - -void DerivativeTreatment_Gradient(CGeometry *geometry, CConfig *config, CVolumetricMovement* grid_movement, CSurfaceMovement *surface_movement, su2double** Gradient) { - - int rank = SU2_MPI::GetRank(); - - /*--- Error if choosen smoothing mode is unsupported. - * This is the default option if the user has not specified a mode in the config file. ---*/ - if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::NONE) { - SU2_MPI::Error("Unsupported operation modus for the Sobolev Smoothing Solver.", CURRENT_FUNCTION); - } - - /*-- Construct the smoothing solver and numerics ---*/ - std::unique_ptr solver(new CGradientSmoothingSolver(geometry, config)); - unsigned dim = (config->GetSmoothOnSurface() ? geometry->GetnDim() - 1 : geometry->GetnDim()); - std::unique_ptr numerics(new CGradSmoothing(dim, config)); - - if (rank == MASTER_NODE) cout << "Sobolev Smoothing of derivatives is active." << endl; - - /*--- Get the sensitivities from the geometry class to work with. ---*/ - solver->ReadSensFromGeometry(geometry); - - /*--- Apply the smoothing procedure on the DV level. ---*/ - if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::PARAM_LEVEL_COMPLETE) { - solver->ApplyGradientSmoothingDV(geometry, numerics.get(), surface_movement, grid_movement, config, Gradient); - - /*--- If smoothing already took place on the mesh level, or none is requested, just do standard projection. ---*/ - } else if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::ONLY_GRAD || - config->GetSobMode() == ENUM_SOBOLEV_MODUS::MESH_LEVEL) { - solver->RecordTapeAndCalculateOriginalGradient(geometry, surface_movement, grid_movement, config, Gradient); - } - } diff --git a/SU2_DOT/src/meson.build b/SU2_DOT/src/meson.build index e6a92e38a647..3cda6c961ed9 100644 --- a/SU2_DOT/src/meson.build +++ b/SU2_DOT/src/meson.build @@ -1,11 +1,11 @@ -su2_dot_src = ['SU2_DOT.cpp'] +su2_dot_src = 'SU2_DOT.cpp' if get_option('enable-normal') su2_dot = executable('SU2_DOT', su2_dot_src, install: true, - dependencies: [su2_deps, common_dep, su2_cfd_dep], + dependencies: [su2_deps, common_dep, su2_cfd_dep, su2_def_dep], cpp_args :[default_warning_flags, su2_cpp_args]) endif @@ -15,7 +15,7 @@ if get_option('enable-autodiff') su2_dot_ad = executable('SU2_DOT_AD', su2_dot_src, install: true, - dependencies: [su2_deps, codi_dep, commonAD_dep, su2_cfd_dep_ad], - cpp_args : [default_warning_flags, su2_cpp_args, codi_rev_args]) + dependencies: [su2_deps, codi_dep, commonAD_dep, su2_cfd_dep_ad, su2_def_dep_ad], + cpp_args : [default_warning_flags, su2_cpp_args, codi_rev_args]) endif diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py index 7e2f19403e2e..99bab9afe9a3 100644 --- a/SU2_PY/FSI_tools/FSIInterface.py +++ b/SU2_PY/FSI_tools/FSIInterface.py @@ -227,15 +227,15 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): if FluidSolver is not None: print('Fluid solver is initialized on process {}'.format(myid)) self.haveFluidSolver = True - allMovingMarkersTags = FluidSolver.GetAllDeformMeshMarkersTag() - allMarkersID = FluidSolver.GetAllBoundaryMarkers() + allMovingMarkersTags = FluidSolver.GetDeformableMarkerTags() + allMarkersID = FluidSolver.GetMarkerTags() if not allMovingMarkersTags: raise Exception('No interface for FSI was defined.') else: if allMovingMarkersTags[0] in allMarkersID.keys(): self.fluidInterfaceIdentifier = allMarkersID[allMovingMarkersTags[0]] if self.fluidInterfaceIdentifier is not None: - self.nLocalFluidInterfaceNodes = FluidSolver.GetNumberVertices(self.fluidInterfaceIdentifier) + self.nLocalFluidInterfaceNodes = FluidSolver.GetNumberMarkerNodes(self.fluidInterfaceIdentifier) if self.nLocalFluidInterfaceNodes != 0: self.haveFluidInterface = True print('Number of interface fluid nodes (halo nodes included) on proccess {} : {}'.format(myid,self.nLocalFluidInterfaceNodes)) @@ -300,8 +300,9 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): # Calculate the number of halo nodes on each partition self.nLocalFluidInterfaceHaloNode = 0 for iVertex in range(self.nLocalFluidInterfaceNodes): - if FluidSolver.IsAHaloNode(self.fluidInterfaceIdentifier, iVertex): - GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) + iPoint = FluidSolver.GetMarkerNode(self.fluidInterfaceIdentifier, iVertex) + if not FluidSolver.GetNodeDomain(iPoint): + GlobalIndex = FluidSolver.GetNodeGlobalIndex(iPoint) self.FluidHaloNodeList[GlobalIndex] = iVertex self.nLocalFluidInterfaceHaloNode += 1 # Calculate the number of physical (= not halo) nodes on each partition @@ -565,8 +566,13 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): # Note that the fluid solver is separated in more processors outside the python script # thus when, from a core, we request for the vertices on the interface, we only obtain # those in that core - GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) - posx, posy, posz = FluidSolver.GetInitialMeshCoord(self.fluidInterfaceIdentifier, iVertex) + iPoint = FluidSolver.GetMarkerNode(self.fluidInterfaceIdentifier, iVertex) + GlobalIndex = FluidSolver.GetNodeGlobalIndex(iPoint) + if self.nDim == 2: + posx, posy = FluidSolver.GetInitialCoordinates(iPoint) + posz = 0 + else: + posx, posy, posz = FluidSolver.GetInitialCoordinates(iPoint) if GlobalIndex not in self.FluidHaloNodeList[myid].keys(): fluidIndexing_temp[GlobalIndex] = self.__getGlobalIndex('fluid', myid, localIndex) self.localFluidInterface_array_X_init[localIndex] = posx @@ -1456,7 +1462,7 @@ def getFluidInterfaceNodalForce(self, FSI_config, FluidSolver): # --- Get the fluid interface loads from the fluid solver and directly fill the corresponding PETSc vector --- for iVertex in range(self.nLocalFluidInterfaceNodes): - GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) + GlobalIndex = FluidSolver.GetNodeGlobalIndex(FluidSolver.GetMarkerNode(self.fluidInterfaceIdentifier, iVertex)) if GlobalIndex not in self.FluidHaloNodeList[myid].keys(): loadX, loadY, loadZ = FluidSolver.GetFlowLoad(self.fluidInterfaceIdentifier, iVertex) iGlobalVertex = self.__getGlobalIndex('fluid', myid, localIndex) @@ -1486,15 +1492,15 @@ def setFluidInterfaceVarCoord(self, FluidSolver): # --- Send the new fluid interface position to the fluid solver (on each partition, halo nodes included) --- localIndex = 0 for iVertex in range(self.nLocalFluidInterfaceNodes): - GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) + GlobalIndex = FluidSolver.GetNodeGlobalIndex(FluidSolver.GetMarkerNode(self.fluidInterfaceIdentifier, iVertex)) if GlobalIndex in self.FluidHaloNodeList[myid].keys(): DispX, DispY, DispZ = self.haloNodesDisplacements[GlobalIndex] - FluidSolver.SetMeshDisplacement(self.fluidInterfaceIdentifier, int(iVertex), DispX, DispY, DispZ) + FluidSolver.SetMarkerDisplacements(self.fluidInterfaceIdentifier, int(iVertex), np.array([DispX, DispY, DispZ])) else: DispX = self.localFluidInterface_array_DispX[localIndex] DispY = self.localFluidInterface_array_DispY[localIndex] DispZ = self.localFluidInterface_array_DispZ[localIndex] - FluidSolver.SetMeshDisplacement(self.fluidInterfaceIdentifier, int(iVertex), DispX, DispY, DispZ) + FluidSolver.SetMarkerDisplacements(self.fluidInterfaceIdentifier, int(iVertex), np.array([DispX, DispY, DispZ])) localIndex += 1 @@ -2142,8 +2148,12 @@ def MapModes(self, FSI_config, FluidSolver, SolidSolver): nodeNormals = {} for iVertex in range(self.nLocalFluidInterfaceNodes): - nx, ny, nz = FluidSolver.GetVertexNormal(self.fluidInterfaceIdentifier, iVertex, False) - GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) + if self.nDim == 2: + nx, ny = FluidSolver.GetMarkerVertexNormals(self.fluidInterfaceIdentifier, iVertex, False) + nz = 0 + else: + nx, ny, nz = FluidSolver.GetMarkerVertexNormals(self.fluidInterfaceIdentifier, iVertex, False) + GlobalIndex = FluidSolver.GetNodeGlobalIndex(FluidSolver.GetMarkerNode(self.fluidInterfaceIdentifier, iVertex)) nodeNormals[GlobalIndex] = [nx, ny, nz] nodeNormals = self.comm.gather(nodeNormals, root=self.rootProcess) diff --git a/SU2_PY/pySU2/Makefile.am b/SU2_PY/pySU2/Makefile.am index febedbc98b67..a5185f54e3f2 100644 --- a/SU2_PY/pySU2/Makefile.am +++ b/SU2_PY/pySU2/Makefile.am @@ -132,4 +132,3 @@ clean: install: all ${INSTALL} ${PYLIB} ${bindir} ${INSTALL} ${SOLIB} ${bindir} - diff --git a/SU2_PY/pySU2/meson.build b/SU2_PY/pySU2/meson.build index a71e6e307b9e..3d9abc747ec5 100644 --- a/SU2_PY/pySU2/meson.build +++ b/SU2_PY/pySU2/meson.build @@ -37,7 +37,10 @@ if get_option('enable-normal') '_pysu2', cpp_source, dependencies: [wrapper_deps, common_dep, su2_deps], - objects: su2_cfd_lib.extract_all_objects(), + objects: [ + su2_cfd_lib.extract_all_objects(), + su2_def_lib.extract_objects('drivers/CDeformationDriver.cpp') + ], install: true, include_directories : mpi4py_include, cpp_args : [default_warning_flags,su2_cpp_args], @@ -53,7 +56,10 @@ if get_option('enable-autodiff') '_pysu2ad', cpp_source, dependencies: [wrapper_deps, commonAD_dep, su2_deps, codi_dep], - objects: su2_cfd_lib_ad.extract_all_objects(), + objects: [ + su2_cfd_lib_ad.extract_all_objects(), + su2_def_lib_ad.extract_objects('drivers/CDiscAdjDeformationDriver.cpp') + ], install: true, include_directories : mpi4py_include, cpp_args : [default_warning_flags, su2_cpp_args, codi_rev_args], diff --git a/SU2_PY/pySU2/pySU2.i b/SU2_PY/pySU2/pySU2.i index 69384d70f3db..a05049798b7b 100644 --- a/SU2_PY/pySU2/pySU2.i +++ b/SU2_PY/pySU2/pySU2.i @@ -37,11 +37,12 @@ directors="1", threads="1" ) pysu2 %{ - +#include "../../SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp" #include "../../SU2_CFD/include/drivers/CDriver.hpp" -#include "../../SU2_CFD/include/drivers/CSinglezoneDriver.hpp" +#include "../../SU2_CFD/include/drivers/CDriverBase.hpp" #include "../../SU2_CFD/include/drivers/CMultizoneDriver.hpp" -#include "../../SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp" +#include "../../SU2_CFD/include/drivers/CSinglezoneDriver.hpp" +#include "../../SU2_DEF/include/drivers/CDeformationDriver.hpp" %} @@ -49,6 +50,7 @@ threads="1" %import "../../Common/include/code_config.hpp" %import "../../Common/include/basic_types/datatype_structure.hpp" %import "../../Common/include/parallelization/mpi_structure.hpp" + %include "std_string.i" %include "std_vector.i" %include "std_map.i" @@ -60,10 +62,14 @@ threads="1" #endif namespace std { - %template() vector; + %template() vector; + %template() vector; + %template() vector; + %template() vector>; %template() vector; + %template() vector>; %template() vector; - %template() map; + %template() map; %template() map; } @@ -86,8 +92,9 @@ const unsigned int MESH_1 = 1; /*!< \brief Definition of the finest grid level. const unsigned int ZONE_0 = 0; /*!< \brief Definition of the first grid domain. */ const unsigned int ZONE_1 = 1; /*!< \brief Definition of the first grid domain. */ -// CDriver class +%include "../../SU2_CFD/include/drivers/CDriverBase.hpp" %include "../../SU2_CFD/include/drivers/CDriver.hpp" %include "../../SU2_CFD/include/drivers/CSinglezoneDriver.hpp" %include "../../SU2_CFD/include/drivers/CMultizoneDriver.hpp" %include "../../SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp" +%include "../../SU2_DEF/include/drivers/CDeformationDriver.hpp" diff --git a/SU2_PY/pySU2/pySU2ad.i b/SU2_PY/pySU2/pySU2ad.i index fe329e7ab8be..cba62127e92c 100644 --- a/SU2_PY/pySU2/pySU2ad.i +++ b/SU2_PY/pySU2/pySU2ad.i @@ -37,11 +37,12 @@ directors="1", threads="1" ) pysu2ad %{ - +#include "../../SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp" #include "../../SU2_CFD/include/drivers/CDriver.hpp" -#include "../../SU2_CFD/include/drivers/CSinglezoneDriver.hpp" +#include "../../SU2_CFD/include/drivers/CDriverBase.hpp" #include "../../SU2_CFD/include/drivers/CMultizoneDriver.hpp" -#include "../../SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp" +#include "../../SU2_CFD/include/drivers/CSinglezoneDriver.hpp" +#include "../../SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp" %} @@ -49,6 +50,7 @@ threads="1" %import "../../Common/include/code_config.hpp" %import "../../Common/include/basic_types/datatype_structure.hpp" %import "../../Common/include/parallelization/mpi_structure.hpp" + %include "std_string.i" %include "std_vector.i" %include "std_map.i" @@ -60,10 +62,13 @@ threads="1" #endif namespace std { - %template() vector; + %template() vector; + %template() vector; + %template() vector>; %template() vector; + %template() vector>; %template() vector; - %template() map; + %template() map; %template() map; } @@ -86,8 +91,9 @@ const unsigned int MESH_1 = 1; /*!< \brief Definition of the finest grid level. const unsigned int ZONE_0 = 0; /*!< \brief Definition of the first grid domain. */ const unsigned int ZONE_1 = 1; /*!< \brief Definition of the first grid domain. */ -// CDriver class +%include "../../SU2_CFD/include/drivers/CDriverBase.hpp" %include "../../SU2_CFD/include/drivers/CDriver.hpp" %include "../../SU2_CFD/include/drivers/CSinglezoneDriver.hpp" %include "../../SU2_CFD/include/drivers/CMultizoneDriver.hpp" %include "../../SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp" +%include "../../SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp" diff --git a/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/run_adjoint.py b/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/run_adjoint.py index 49d0194410ae..3b33bd46d538 100755 --- a/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/run_adjoint.py +++ b/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/run_adjoint.py @@ -68,10 +68,10 @@ def main(): MarkerName = 'RightBeamS' # Specified by the user # Get all the boundary tags - MarkerList = SU2Driver.GetAllBoundaryMarkersTag() + MarkerList = SU2Driver.GetMarkerTags() # Get all the markers defined on this rank and their associated indices. - allMarkerIDs = SU2Driver.GetAllBoundaryMarkers() + allMarkerIDs = SU2Driver.GetMarkerIndices() #Check if the specified marker exists and if it belongs to this rank. if MarkerName in MarkerList and MarkerName in allMarkerIDs.keys(): diff --git a/TestCases/py_wrapper/disc_adj_flow/mesh_disp_sens/run_adjoint.py b/TestCases/py_wrapper/disc_adj_flow/mesh_disp_sens/run_adjoint.py index 8a236c4ba6a1..99d321264147 100755 --- a/TestCases/py_wrapper/disc_adj_flow/mesh_disp_sens/run_adjoint.py +++ b/TestCases/py_wrapper/disc_adj_flow/mesh_disp_sens/run_adjoint.py @@ -67,12 +67,12 @@ def main(): MarkerName = 'wallF' # Specified by the user # Get all the boundary tags - MarkerList = SU2Driver.GetAllBoundaryMarkersTag() + MarkerList = SU2Driver.GetMarkerTags() # Get all the markers defined on this rank and their associated indices. - allMarkerIDs = SU2Driver.GetAllBoundaryMarkers() + allMarkerIDs = SU2Driver.GetMarkerIndices() - #Check if the specified marker exists and if it belongs to this rank. + # Check if the specified marker exists and if it belongs to this rank. if MarkerName in MarkerList and MarkerName in allMarkerIDs.keys(): MarkerID = allMarkerIDs[MarkerName] @@ -80,7 +80,7 @@ def main(): nVertex_Marker = 0 #total number of vertices (physical + halo) if MarkerID != None: - nVertex_Marker = SU2Driver.GetNumberVertices(MarkerID) + nVertex_Marker = SU2Driver.GetNumberMarkerNodes(MarkerID) # Time loop is defined in Python so that we have acces to SU2 functionalities at each time step if rank == 0: diff --git a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py index da2ecb284337..37e7a1127f6d 100755 --- a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py +++ b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py @@ -76,10 +76,10 @@ def main(): MovingMarker = 'plate' #specified by the user # Get all the tags with the moving option - MovingMarkerList = SU2Driver.GetAllDeformMeshMarkersTag() + MovingMarkerList = SU2Driver.GetMarkerTags() # Get all the markers defined on this rank and their associated indices. - allMarkerIDs = SU2Driver.GetAllBoundaryMarkers() + allMarkerIDs = SU2Driver.GetMarkerIndices() # Check if the specified marker has a moving option and if it exists on this rank. if MovingMarker in MovingMarkerList and MovingMarker in allMarkerIDs.keys(): @@ -87,13 +87,9 @@ def main(): # Number of vertices on the specified marker (per rank) nVertex_MovingMarker = 0 #total number of vertices (physical + halo) - nVertex_MovingMarker_HALO = 0 #number of halo vertices - nVertex_MovingMarker_PHYS = 0 #number of physical vertices if MovingMarkerID != None: - nVertex_MovingMarker = SU2Driver.GetNumberVertices(MovingMarkerID) - nVertex_MovingMarker_HALO = SU2Driver.GetNumberHaloVertices(MovingMarkerID) - nVertex_MovingMarker_PHYS = nVertex_MovingMarker - nVertex_MovingMarker_HALO + nVertex_MovingMarker = SU2Driver.GetNumberMarkerNodes(MovingMarkerID)\ # Retrieve some control parameters from the driver deltaT = SU2Driver.GetUnsteady_TimeStep() @@ -104,9 +100,9 @@ def main(): # Extract the initial position of each node on the moving marker CoordX = np.zeros(nVertex_MovingMarker) CoordY = np.zeros(nVertex_MovingMarker) - CoordZ = np.zeros(nVertex_MovingMarker) for iVertex in range(nVertex_MovingMarker): - CoordX[iVertex], CoordY[iVertex], CoordZ[iVertex] = SU2Driver.GetInitialMeshCoord(MovingMarkerID, iVertex) + iPoint = SU2Driver.GetMarkerNode(MovingMarkerID, iVertex) + CoordX[iVertex], CoordY[iVertex] = SU2Driver.GetInitialCoordinates(iPoint) # Time loop is defined in Python so that we have acces to SU2 functionalities at each time step if rank == 0: @@ -117,9 +113,11 @@ def main(): while (TimeIter < nTimeIter): # Define the rigid body displacement and set the new coords of each node on the marker - d_y = 0.0175*sin(2*pi*time) + value = 0.0, 0.0175*sin(2*pi*time) + for iVertex in range(nVertex_MovingMarker): - SU2Driver.SetMeshDisplacement(MovingMarkerID, int(iVertex), 0.0, d_y, 0.0) + SU2Driver.SetMarkerDisplacements(MovingMarkerID, int(iVertex), value) + # Time iteration preprocessing SU2Driver.Preprocess(TimeIter) # Run one time iteration (e.g. dual-time) diff --git a/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py b/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py index 34007a95eb8b..6a613b1fa180 100755 --- a/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py +++ b/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py @@ -75,24 +75,20 @@ def main(): CHTMarker = 'plate' # Specified by the user # Get all the tags with the CHT option - CHTMarkerList = SU2Driver.GetAllCHTMarkersTag() + CHTMarkerList = SU2Driver.GetCHTMarkerTags() # Get all the markers defined on this rank and their associated indices. - allMarkerIDs = SU2Driver.GetAllBoundaryMarkers() + allMarkerIDs = SU2Driver.GetMarkerIndices() #Check if the specified marker has a CHT option and if it exists on this rank. if CHTMarker in CHTMarkerList and CHTMarker in allMarkerIDs.keys(): CHTMarkerID = allMarkerIDs[CHTMarker] # Number of vertices on the specified marker (per rank) - nVertex_CHTMarker = 0 #total number of vertices (physical + halo) - nVertex_CHTMarker_HALO = 0 #number of halo vertices - nVertex_CHTMarker_PHYS = 0 #number of physical vertices + nVertex_CHTMarker = 0 # total number of vertices (physical + halo) if CHTMarkerID != None: - nVertex_CHTMarker = SU2Driver.GetNumberVertices(CHTMarkerID) - nVertex_CHTMarker_HALO = SU2Driver.GetNumberHaloVertices(CHTMarkerID) - nVertex_CHTMarker_PHYS = nVertex_CHTMarker - nVertex_CHTMarker_HALO + nVertex_CHTMarker = SU2Driver.GetNumberMarkerNodes(CHTMarkerID) # Retrieve some control parameters from the driver deltaT = SU2Driver.GetUnsteady_TimeStep() @@ -115,6 +111,7 @@ def main(): # Set this temperature to all the vertices on the specified CHT marker for iVertex in range(nVertex_CHTMarker): SU2Driver.SetVertexTemperature(CHTMarkerID, iVertex, WallTemp) + # Tell the SU2 drive to update the boundary conditions SU2Driver.BoundaryConditionsUpdate() # Run one time iteration (e.g. dual-time) diff --git a/TestCases/py_wrapper/translating_NACA0012/run_su2.py b/TestCases/py_wrapper/translating_NACA0012/run_su2.py index 731c109184d5..7d638c92b08c 100644 --- a/TestCases/py_wrapper/translating_NACA0012/run_su2.py +++ b/TestCases/py_wrapper/translating_NACA0012/run_su2.py @@ -30,26 +30,28 @@ def run_solver(self): self.FluidSolver.Run() self.FluidSolver.Postprocess() # write outputs - self.FluidSolver.Monitor(0) + self.FluidSolver.Monitor(0) self.FluidSolver.Output(0) self.comm.barrier() def save_forces(self): - solver_all_moving_markers = np.array(self.FluidSolver.GetAllDeformMeshMarkersTag()) - solver_marker_ids = self.FluidSolver.GetAllBoundaryMarkers() + solver_all_moving_markers = np.array(self.FluidSolver.GetDeformableMarkerTags()) + solver_markers = self.FluidSolver.GetMarkerTags() + solver_marker_ids = self.FluidSolver.GetMarkerIndices() # The surface marker and the partitioning of the solver usually don't agree. # Thus, it is necessary to figure out if the partition of the current mpi process has # a node that belongs to a moving surface marker. - has_moving_marker = [marker in solver_marker_ids.keys() for marker in solver_all_moving_markers] + has_moving_marker = [marker in solver_markers for marker in solver_all_moving_markers] f = open('forces_'+str(self.myid)+'.csv','w') for marker in solver_all_moving_markers[has_moving_marker]: solver_marker_id = solver_marker_ids[marker] - n_vertices = self.FluidSolver.GetNumberVertices(solver_marker_id) + n_vertices = self.FluidSolver.GetNumberMarkerNodes(solver_marker_id) for i_vertex in range(n_vertices): fxyz = self.FluidSolver.GetFlowLoad(solver_marker_id, i_vertex) - GlobalIndex = self.FluidSolver.GetVertexGlobalIndex(solver_marker_id, i_vertex) + iPoint = self.FluidSolver.GetMarkerNode(solver_marker_id, i_vertex) + GlobalIndex = self.FluidSolver.GetNodeGlobalIndex(iPoint) f.write('{}, {:.2f}, {:.2f}, {:.2f}\n'.format(GlobalIndex, fxyz[0], fxyz[1], fxyz[2])) f.close() diff --git a/meson.build b/meson.build index ca86992ea88b..84999a5476b9 100644 --- a/meson.build +++ b/meson.build @@ -268,10 +268,10 @@ endif subdir('Common/src') # compile SU2_CFD executable subdir('SU2_CFD/src') -# compile SU2_DOT executable -subdir('SU2_DOT/src') # compile SU2_DEF executable subdir('SU2_DEF/src') +# compile SU2_DOT executable +subdir('SU2_DOT/src') # compile SU2_GEO executable subdir('SU2_GEO/src') # compile SU2_SOL executable