diff --git a/Common/include/config_structure.hpp b/Common/include/config_structure.hpp index 364713ddde4..35dc671d0e2 100644 --- a/Common/include/config_structure.hpp +++ b/Common/include/config_structure.hpp @@ -58,7 +58,26 @@ using namespace std; * \version 4.1.2 "Cardinal" */ +struct StrainItem{ + + vector IndexCurr; + vector p1; + vector walldist; + vector IndexBndy; + vector strain_rate; + vector mu_t; + +}; + +struct TauItem{ + + vector IndexCurr; + vector TauTangent; + +}; + class CConfig { + private: unsigned short Kind_SU2; /*!< \brief Kind of SU2 software component.*/ unsigned short Ref_NonDim; /*!< \brief Kind of non dimensionalization.*/ @@ -730,7 +749,7 @@ class CConfig { bool ParMETIS; /*!< \brief Boolean for activating ParMETIS mode (while testing). */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */ - + double *betaArr; /*!< \brief Array for values of SA_Production_Factor >*/ /*--- all_options is a map containing all of the options. This is used during config file parsing to track the options which have not been set (so the default values can be used). Without this map @@ -1016,6 +1035,8 @@ class CConfig { public: + StrainItem StrainFile; + TauItem TauFile; vector fields; /*!< \brief Tags for the different fields in a restart file. */ /*! @@ -1469,6 +1490,8 @@ class CConfig { su2double GetSA_Production_Factor(void); + double GetbetaArr(int index); + /*! * \brief Get the value of the non-dimensionalized freestream pressure. * \return Non-dimensionalized freestream pressure. @@ -1862,6 +1885,8 @@ class CConfig { void SetSA_Production_Factor(su2double val_SA_Production_Factor); + void SetbetaArr(double *val_betaArr); + /*! * \brief Set the Froude number for free surface problems. * \return Value of the Froude number. diff --git a/Common/include/config_structure.inl b/Common/include/config_structure.inl index b6a41d30a72..ba193c2beea 100644 --- a/Common/include/config_structure.inl +++ b/Common/include/config_structure.inl @@ -206,6 +206,8 @@ inline su2double CConfig::GetPressure_FreeStream(void) { return Pressure_FreeStr inline su2double CConfig::GetSA_Production_Factor(void) { return SA_Production_Factor; } +inline double CConfig::GetbetaArr(int index) { return betaArr[index]; } + inline su2double CConfig::GetTemperature_ve_FreeStream(void) { return Temperature_ve_FreeStream; } inline su2double CConfig::GetPrandtl_Lam(void) { return Prandtl_Lam; } @@ -307,6 +309,8 @@ inline void CConfig::SetPressure_FreeStreamND(su2double val_pressure_freestreamn inline void CConfig::SetSA_Production_Factor(su2double val_SA_Production_Factor) { SA_Production_Factor = val_SA_Production_Factor; } +inline void CConfig::SetbetaArr(double *val_betaArr){ betaArr = val_betaArr; } + inline void CConfig::SetPressure_FreeStream(su2double val_pressure_freestream) { Pressure_FreeStream = val_pressure_freestream; } inline void CConfig::SetDensity_FreeStreamND(su2double val_density_freestreamnd) { Density_FreeStreamND = val_density_freestreamnd; } diff --git a/Common/include/dual_grid_structure.hpp b/Common/include/dual_grid_structure.hpp index 4165c4fc9b3..feb1d15b383 100644 --- a/Common/include/dual_grid_structure.hpp +++ b/Common/include/dual_grid_structure.hpp @@ -164,6 +164,7 @@ class CPoint : public CDualGrid { bool Move; /*!< \brief This flag indicates if the point is going to be move in the grid deformation process. */ unsigned short color; /*!< \brief Color of the point in the partitioning strategy. */ su2double Wall_Distance; /*!< \brief Distance to the nearest wall. */ + unsigned long Vertex_nearWall; /*!< \brief Global index of the wall vertex nearest to the cell >*/ su2double SharpEdge_Distance; /*!< \brief Distance to a sharp edge. */ su2double Curvature; /*!< \brief Value of the surface curvature (SU2_GEO). */ unsigned long GlobalIndex; /*!< \brief Global index in the parallel simulation. */ @@ -221,6 +222,9 @@ class CPoint : public CDualGrid { * \param[in] val_distance - Value of the distance. */ void SetWall_Distance(su2double val_distance); + + void SetVertex_nearWall(unsigned long val_Vertex_nearWall); + /*! * \brief Set the value of the distance to a sharp edge. @@ -233,7 +237,9 @@ class CPoint : public CDualGrid { * \return Value of the distance to the nearest wall. */ su2double GetWall_Distance(void); - + + unsigned long GetVertex_nearWall(void); + /*! * \brief Set the value of the curvature at a surface node. * \param[in] val_curvature - Value of the curvature. diff --git a/Common/include/dual_grid_structure.inl b/Common/include/dual_grid_structure.inl index 50814aa8d3c..bcd88c12b2c 100644 --- a/Common/include/dual_grid_structure.inl +++ b/Common/include/dual_grid_structure.inl @@ -215,12 +215,16 @@ inline bool CPoint::GetDomain(void) { return Domain; } inline void CPoint::SetWall_Distance(su2double val_distance) { Wall_Distance = val_distance; } +inline void CPoint::SetVertex_nearWall(unsigned long val_Vertex_nearWall) { Vertex_nearWall = val_Vertex_nearWall; } + inline void CPoint::SetCurvature(su2double val_curvature) { Curvature = val_curvature; } inline void CPoint::SetSharpEdge_Distance(su2double val_distance) { SharpEdge_Distance = val_distance; } inline su2double CPoint::GetWall_Distance(void) { return Wall_Distance; } +inline unsigned long CPoint::GetVertex_nearWall(void) { return Vertex_nearWall; } + inline su2double CPoint::GetCurvature(void) { return Curvature; } inline su2double CPoint::GetSharpEdge_Distance(void) { return SharpEdge_Distance; } diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp index 428821150b2..6af365e08cc 100644 --- a/Common/src/geometry_structure.cpp +++ b/Common/src/geometry_structure.cpp @@ -8677,8 +8677,9 @@ void CPhysicalGeometry::ComputeWall_Distance(CConfig *config) { /*--- Allocate an array to hold boundary node coordinates ---*/ - su2double **Coord_bound; + su2double **Coord_bound, **pNormal, *pind; Coord_bound = new su2double* [nVertex_SolidWall]; + pind = new su2double [nVertex_SolidWall]; for (iVertex = 0; iVertex < nVertex_SolidWall; iVertex++) Coord_bound[iVertex] = new su2double [nDim]; @@ -8692,6 +8693,7 @@ void CPhysicalGeometry::ComputeWall_Distance(CConfig *config) { iPoint = vertex[iMarker][iVertex]->GetNode(); for (iDim = 0; iDim < nDim; iDim++) Coord_bound[nVertex_SolidWall][iDim] = node[iPoint]->GetCoord(iDim); + pind[nVertex_SolidWall] = node[iPoint]->GetGlobalIndex(); nVertex_SolidWall++; } } @@ -8731,6 +8733,7 @@ void CPhysicalGeometry::ComputeWall_Distance(CConfig *config) { (coord[iDim] - Coord_bound[iVertex_nearestWall][iDim]); } node[iPoint]->SetWall_Distance(sqrt(dist)); + node[iPoint]->SetVertex_nearWall(pind[iVertex_nearestWall]); } } else { @@ -8776,9 +8779,12 @@ void CPhysicalGeometry::ComputeWall_Distance(CConfig *config) { /*--- Create and initialize to zero some buffers to hold the coordinates of the boundary nodes that are communicated from each partition (all-to-all). ---*/ - su2double *Buffer_Send_Coord = new su2double [MaxLocalVertex_NS*nDim]; - su2double *Buffer_Receive_Coord = new su2double [nProcessor*MaxLocalVertex_NS*nDim]; - unsigned long nBuffer = MaxLocalVertex_NS*nDim; + su2double *Buffer_Send_Coord = new su2double [MaxLocalVertex_NS*nDim]; + unsigned long *Buffer_Send_Index = new unsigned long [MaxLocalVertex_NS]; + su2double *Buffer_Receive_Coord = new su2double [nProcessor*MaxLocalVertex_NS*nDim]; + unsigned long *Buffer_Receive_Index = new unsigned long [nProcessor*MaxLocalVertex_NS]; + unsigned long nBuffer = MaxLocalVertex_NS*nDim; + unsigned long nBufferIndex = MaxLocalVertex_NS; for (iVertex = 0; iVertex < MaxLocalVertex_NS; iVertex++) for (iDim = 0; iDim < nDim; iDim++) @@ -8793,12 +8799,15 @@ void CPhysicalGeometry::ComputeWall_Distance(CConfig *config) { (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) ) for (iVertex = 0; iVertex < GetnVertex(iMarker); iVertex++) { iPoint = vertex[iMarker][iVertex]->GetNode(); - for (iDim = 0; iDim < nDim; iDim++) + for (iDim = 0; iDim < nDim; iDim++){ Buffer_Send_Coord[nVertex_SolidWall*nDim+iDim] = node[iPoint]->GetCoord(iDim); + } + Buffer_Send_Index[nVertex_SolidWall] = node[iPoint]->GetGlobalIndex(); nVertex_SolidWall++; } SU2_MPI::Allgather(Buffer_Send_Coord, nBuffer, MPI_DOUBLE, Buffer_Receive_Coord, nBuffer, MPI_DOUBLE, MPI_COMM_WORLD); + SU2_MPI::Allgather(Buffer_Send_Index, nBufferIndex, MPI_UNSIGNED_LONG, Buffer_Receive_Index, nBufferIndex, MPI_UNSIGNED_LONG, MPI_COMM_WORLD); /*--- Loop over all interior mesh nodes on the local partition and compute the distances to each of the no-slip boundary nodes in the entire mesh. @@ -8840,6 +8849,7 @@ void CPhysicalGeometry::ComputeWall_Distance(CConfig *config) { (coord[iDim] - Buffer_Receive_Coord[iVertex_nearestWall*nDim+iDim]); } node[iPoint]->SetWall_Distance(sqrt(dist)); + node[iPoint]->SetVertex_nearWall(Buffer_Receive_Index[iVertex_nearestWall]); } } else { @@ -8851,6 +8861,8 @@ void CPhysicalGeometry::ComputeWall_Distance(CConfig *config) { delete[] Buffer_Send_Coord; delete[] Buffer_Receive_Coord; + delete[] Buffer_Send_Index; + delete[] Buffer_Receive_Index; delete[] Buffer_Send_nVertex; delete[] Buffer_Receive_nVertex; diff --git a/SU2_AD/SU2_CFD/bin/SU2_CFD_AD b/SU2_AD/SU2_CFD/bin/SU2_CFD_AD index 07ab680308b..1c67fd87bf2 100755 Binary files a/SU2_AD/SU2_CFD/bin/SU2_CFD_AD and b/SU2_AD/SU2_CFD/bin/SU2_CFD_AD differ diff --git a/SU2_AD/SU2_DOT/bin/SU2_DOT_AD b/SU2_AD/SU2_DOT/bin/SU2_DOT_AD index 455baa88253..771746e0b66 100755 Binary files a/SU2_AD/SU2_DOT/bin/SU2_DOT_AD and b/SU2_AD/SU2_DOT/bin/SU2_DOT_AD differ diff --git a/SU2_BASE/SU2_CFD/bin/SU2_CFD b/SU2_BASE/SU2_CFD/bin/SU2_CFD index 26c5bf19149..2b6763bac02 100755 Binary files a/SU2_BASE/SU2_CFD/bin/SU2_CFD and b/SU2_BASE/SU2_CFD/bin/SU2_CFD differ diff --git a/SU2_BASE/SU2_DEF/bin/SU2_DEF b/SU2_BASE/SU2_DEF/bin/SU2_DEF index ee53be40fd2..45cf26225c5 100755 Binary files a/SU2_BASE/SU2_DEF/bin/SU2_DEF and b/SU2_BASE/SU2_DEF/bin/SU2_DEF differ diff --git a/SU2_BASE/SU2_DOT/bin/SU2_DOT b/SU2_BASE/SU2_DOT/bin/SU2_DOT index 7039ef8badd..a2201693ce5 100755 Binary files a/SU2_BASE/SU2_DOT/bin/SU2_DOT and b/SU2_BASE/SU2_DOT/bin/SU2_DOT differ diff --git a/SU2_BASE/SU2_GEO/bin/SU2_GEO b/SU2_BASE/SU2_GEO/bin/SU2_GEO index 3e598d188a1..71a3469e8f9 100755 Binary files a/SU2_BASE/SU2_GEO/bin/SU2_GEO and b/SU2_BASE/SU2_GEO/bin/SU2_GEO differ diff --git a/SU2_BASE/SU2_MSH/bin/SU2_MSH b/SU2_BASE/SU2_MSH/bin/SU2_MSH index 13478717982..a1530d9e43c 100755 Binary files a/SU2_BASE/SU2_MSH/bin/SU2_MSH and b/SU2_BASE/SU2_MSH/bin/SU2_MSH differ diff --git a/SU2_BASE/SU2_SOL/bin/SU2_SOL b/SU2_BASE/SU2_SOL/bin/SU2_SOL index 0cae08115a3..ecab036f65a 100755 Binary files a/SU2_BASE/SU2_SOL/bin/SU2_SOL and b/SU2_BASE/SU2_SOL/bin/SU2_SOL differ diff --git a/SU2_CFD/src/BAK.SU2_CFD.cpp b/SU2_CFD/src/BAK.SU2_CFD.cpp new file mode 100644 index 00000000000..221a1061844 --- /dev/null +++ b/SU2_CFD/src/BAK.SU2_CFD.cpp @@ -0,0 +1,848 @@ +/*! + * \file SU2_CFD.cpp + * \brief Main file of the Computational Fluid Dynamics code + * \author F. Palacios, T. Economon + * \version 4.1.2 "Cardinal" + * + * SU2 Lead Developers: Dr. Francisco Palacios (Francisco.D.Palacios@boeing.com). + * Dr. Thomas D. Economon (economon@stanford.edu). + * + * SU2 Developers: Prof. Juan J. Alonso's group at Stanford University. + * Prof. Piero Colonna's group at Delft University of Technology. + * Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. + * Prof. Alberto Guardone's group at Polytechnic University of Milan. + * Prof. Rafael Palacios' group at Imperial College London. + * + * Copyright (C) 2012-2016 SU2, the open-source CFD code. + * + * 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/SU2_CFD.hpp" + +using namespace std; + +int main(int argc, char *argv[]) { + + bool StopCalc = false; + su2double StartTime = 0.0, StopTime = 0.0, UsedTime = 0.0; + unsigned long ExtIter = 0; + unsigned short iMesh, iZone, nZone, nDim; + char config_file_name[MAX_STRING_SIZE]; + char runtime_file_name[MAX_STRING_SIZE]; + ofstream ConvHist_file; + int rank = MASTER_NODE; + int size = SINGLE_NODE; + + /*--- MPI initialization, and buffer setting ---*/ + +#ifdef HAVE_MPI + int *bptr, bl; + SU2_MPI::Init(&argc, &argv); + MPI_Buffer_attach( malloc(BUFSIZE), BUFSIZE ); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); +#endif + + /*--- Create pointers to all of the classes that may be used throughout + the SU2_CFD code. In general, the pointers are instantiated down a + heirarchy over all zones, multigrid levels, equation sets, and equation + terms as described in the comments below. ---*/ + + CDriver *driver = NULL; + CIteration **iteration_container = NULL; + COutput *output = NULL; + CIntegration ***integration_container = NULL; + CGeometry ***geometry_container = NULL; + CSolver ****solver_container = NULL; + CNumerics *****numerics_container = NULL; + CConfig **config_container = NULL; + CSurfaceMovement **surface_movement = NULL; + CVolumetricMovement **grid_movement = NULL; + CFreeFormDefBox*** FFDBox = NULL; + CInterpolator ***interpolator_container = NULL; + CTransfer ***transfer_container = NULL; + + /*--- 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 = NULL; + config = new CConfig(config_file_name, SU2_CFD); + + nZone = GetnZone(config->GetMesh_FileName(), config->GetMesh_FileFormat(), config); + nDim = GetnDim(config->GetMesh_FileName(), config->GetMesh_FileFormat()); + + delete config; + /*--- Definition and of the containers for all possible zones. ---*/ + + iteration_container = new CIteration*[nZone]; + solver_container = new CSolver***[nZone]; + integration_container = new CIntegration**[nZone]; + numerics_container = new CNumerics****[nZone]; + config_container = new CConfig*[nZone]; + geometry_container = new CGeometry**[nZone]; + surface_movement = new CSurfaceMovement*[nZone]; + grid_movement = new CVolumetricMovement*[nZone]; + FFDBox = new CFreeFormDefBox**[nZone]; + interpolator_container = new CInterpolator**[nZone]; + transfer_container = new CTransfer**[nZone]; + + for (iZone = 0; iZone < nZone; iZone++) { + solver_container[iZone] = NULL; + integration_container[iZone] = NULL; + numerics_container[iZone] = NULL; + config_container[iZone] = NULL; + geometry_container[iZone] = NULL; + surface_movement[iZone] = NULL; + grid_movement[iZone] = NULL; + FFDBox[iZone] = NULL; + interpolator_container[iZone] = NULL; + transfer_container[iZone] = NULL; + } + + /*--- 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. ---*/ + + config_container[iZone] = new CConfig(config_file_name, SU2_CFD, iZone, nZone, nDim, VERB_HIGH); + + /*--- Definition of the geometry class to store the primal grid in the + partitioning process. ---*/ + + CGeometry *geometry_aux = NULL; + + /*--- 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]); + + /*--- Allocate the memory of the current domain, and divide the grid + between the ranks. ---*/ + + geometry_container[iZone] = new CGeometry *[config_container[iZone]->GetnMGLevels()+1]; + + geometry_container[iZone][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][MESH_0]->SetSendReceive(config_container[iZone]); + + /*--- Add the Send/Receive boundaries ---*/ + + geometry_container[iZone][MESH_0]->SetBoundaries(config_container[iZone]); + + } + + ifstream infile("beta_for_su2.dat"); + unsigned long counter = geometry_container[ZONE_0][MESH_0]->GetGlobal_nPointDomain(); + double *val_betaArr = new double [geometry_container[ZONE_0][MESH_0]->GetGlobal_nPointDomain()]; + if (infile){ + for(int i=0; iGetGlobal_nPointDomain(); i++){ + int index; + infile>>index; + infile>>val_betaArr[index]; + } + } + else { + for(int i=0; iGetGlobal_nPointDomain(); i++){ + val_betaArr[i]=SU2_TYPE::GetValue(config_container[ZONE_0]->GetSA_Production_Factor()); + } + } + infile.close(); + config_container[ZONE_0]->SetbetaArr(val_betaArr); + + + if (rank == MASTER_NODE) + cout << endl <<"------------------------- Geometry Preprocessing ------------------------" << endl; + + /*--- Preprocessing of the geometry for all zones. In this routine, the edge- + based data structure is constructed, i.e. node and cell neighbors are + identified and linked, face areas and volumes of the dual mesh cells are + computed, and the multigrid levels are created using an agglomeration procedure. ---*/ + + Geometrical_Preprocessing(geometry_container, config_container, nZone); + + for (iZone = 0; iZone < nZone; iZone++) { + + /*--- Computation of wall distances for turbulence modeling ---*/ + + if (rank == MASTER_NODE) + cout << "Computing wall distances." << endl; + + if ((config_container[iZone]->GetKind_Solver() == RANS) || + (config_container[iZone]->GetKind_Solver() == ADJ_RANS) || + (config_container[iZone]->GetKind_Solver() == DISC_ADJ_RANS)) + geometry_container[iZone][MESH_0]->ComputeWall_Distance(config_container[iZone]); + + /*--- Computation of positive surface area in the z-plane which is used for + the calculation of force coefficient (non-dimensionalization). ---*/ + + geometry_container[iZone][MESH_0]->SetPositive_ZArea(config_container[iZone]); + + /*--- Set the near-field, interface and actuator disk boundary conditions, if necessary. ---*/ + + for (iMesh = 0; iMesh <= config_container[iZone]->GetnMGLevels(); iMesh++) { + geometry_container[iZone][iMesh]->MatchNearField(config_container[iZone]); + geometry_container[iZone][iMesh]->MatchInterface(config_container[iZone]); + geometry_container[iZone][iMesh]->MatchActuator_Disk(config_container[iZone]); + } + + } + + /*--- If activated by the compile directive, perform a partition analysis. ---*/ +#if PARTITION + Partition_Analysis(geometry_container[ZONE_0][MESH_0], config_container[ZONE_0]); +#endif + + if (rank == MASTER_NODE) + cout << endl <<"------------------------- Driver Preprocessing --------------------------" << endl; + + /*--- First, given the basic information about the number of zones and the + solver types from the config, instantiate the appropriate driver for the problem. ---*/ + + Driver_Preprocessing(&driver, iteration_container, solver_container, + geometry_container, integration_container, numerics_container, + interpolator_container, transfer_container, config_container, nZone, nDim); + + + /*--- Instantiate the geometry movement classes for the solution of unsteady + flows on dynamic meshes, including rigid mesh transformations, dynamically + deforming meshes, and time-spectral preprocessing. ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + + if (config_container[iZone]->GetGrid_Movement() || + (config_container[iZone]->GetDirectDiff() == D_DESIGN)) { + if (rank == MASTER_NODE) + cout << "Setting dynamic mesh structure." << endl; + grid_movement[iZone] = new CVolumetricMovement(geometry_container[iZone][MESH_0], config_container[iZone]); + FFDBox[iZone] = new CFreeFormDefBox*[MAX_NUMBER_FFD]; + surface_movement[iZone] = new CSurfaceMovement(); + surface_movement[iZone]->CopyBoundary(geometry_container[iZone][MESH_0], config_container[iZone]); + if (config_container[iZone]->GetUnsteady_Simulation() == TIME_SPECTRAL) + SetGrid_Movement(geometry_container[iZone], surface_movement[iZone], grid_movement[iZone], + FFDBox[iZone], solver_container[iZone], config_container[iZone], iZone, 0, 0); + } + + if (config_container[iZone]->GetDirectDiff() == D_DESIGN){ + if (rank == MASTER_NODE) + cout << "Setting surface/volume derivatives." << endl; + + /*--- Set the surface derivatives, i.e. the derivative of the surface mesh nodes with respect to the design variables ---*/ + + surface_movement[iZone]->SetSurface_Derivative(geometry_container[iZone][MESH_0],config_container[iZone]); + + /*--- Call the volume deformation routine with derivative mode enabled. + This computes the derivative of the volume mesh with respect to the surface nodes ---*/ + + grid_movement[iZone]->SetVolume_Deformation(geometry_container[iZone][MESH_0],config_container[iZone], true, true); + + /*--- Update the multi-grid structure to propagate the derivative information to the coarser levels ---*/ + + geometry_container[iZone][MESH_0]->UpdateGeometry(geometry_container[iZone],config_container[iZone]); + + /*--- Set the derivative of the wall-distance with respect to the surface nodes ---*/ + + if ( (config_container[iZone]->GetKind_Solver() == RANS) || + (config_container[iZone]->GetKind_Solver() == ADJ_RANS) || + (config_container[iZone]->GetKind_Solver() == DISC_ADJ_RANS)) + geometry_container[iZone][MESH_0]->ComputeWall_Distance(config_container[iZone]); + } + + + } + + /*--- Coupling between zones (limited to two zones at the moment) ---*/ + + bool fsi = config_container[ZONE_0]->GetFSI_Simulation(); + + if ((nZone == 2) && !(fsi)) { + if (rank == MASTER_NODE) + cout << endl <<"--------------------- Setting Coupling Between Zones --------------------" << endl; + geometry_container[ZONE_0][MESH_0]->MatchZone(config_container[ZONE_0], geometry_container[ZONE_1][MESH_0], + config_container[ZONE_1], ZONE_0, nZone); + geometry_container[ZONE_1][MESH_0]->MatchZone(config_container[ZONE_1], geometry_container[ZONE_0][MESH_0], + config_container[ZONE_0], ZONE_1, nZone); + } + + /*--- Definition of the output class (one for all zones). The output class + manages the writing of all restart, volume solution, surface solution, + surface comma-separated value, and convergence history files (both in serial + and in parallel). ---*/ + + output = new COutput(); + + /*--- Open the convergence history file ---*/ + + if (rank == MASTER_NODE) + output->SetConvHistory_Header(&ConvHist_file, config_container[ZONE_0]); + + /*--- Check for an unsteady restart. Update ExtIter if necessary. ---*/ + if (config_container[ZONE_0]->GetWrt_Unsteady() && config_container[ZONE_0]->GetRestart()) + ExtIter = config_container[ZONE_0]->GetUnst_RestartIter(); + + /*--- Check for a dynamic restart (structural analysis). Update ExtIter if necessary. ---*/ + if (config_container[ZONE_0]->GetKind_Solver() == FEM_ELASTICITY + && config_container[ZONE_0]->GetWrt_Dynamic() && config_container[ZONE_0]->GetRestart()) + ExtIter = config_container[ZONE_0]->GetDyn_RestartIter(); + + /*--- Initiate value at each interface for the mixing plane ---*/ + if(config_container[ZONE_0]->GetBoolMixingPlane()) + for (iZone = 0; iZone < nZone; iZone++) + iteration_container[iZone]->Preprocess(output, integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, iZone); + + + + + /*--- Main external loop of the solver. Within this loop, each iteration ---*/ + + if (rank == MASTER_NODE) + cout << endl <<"------------------------------ Begin Solver -----------------------------" << endl; + + /*--- Set up a timer for performance benchmarking (preprocessing time is not included) ---*/ + +#ifndef HAVE_MPI + StartTime = su2double(clock())/su2double(CLOCKS_PER_SEC); +#else + StartTime = MPI_Wtime(); +#endif + + /*--- This is temporal and just to check. It will have to be added to the regular history file ---*/ + + ofstream historyFile_FSI; + bool writeHistFSI = config_container[ZONE_0]->GetWrite_Conv_FSI(); + if (writeHistFSI && (rank == MASTER_NODE)){ + char cstrFSI[200]; + string filenameHistFSI = config_container[ZONE_0]->GetConv_FileName_FSI(); + strcpy (cstrFSI, filenameHistFSI.data()); + historyFile_FSI.open (cstrFSI); + historyFile_FSI << "Time,Iteration,Aitken,URes,logResidual,orderMagnResidual" << endl; + historyFile_FSI.close(); + } + + while (ExtIter < config_container[ZONE_0]->GetnExtIter()) { + + /*--- Set the value of the external iteration. ---*/ + for (iZone = 0; iZone < nZone; iZone++) config_container[iZone]->SetExtIter(ExtIter); + + + /*--- Read the target pressure ---*/ + + if (config_container[ZONE_0]->GetInvDesign_Cp() == YES) + output->SetCp_InverseDesign(solver_container[ZONE_0][MESH_0][FLOW_SOL], + geometry_container[ZONE_0][MESH_0], config_container[ZONE_0], ExtIter); + + /*--- Read the target heat flux ---*/ + + if (config_container[ZONE_0]->GetInvDesign_HeatFlux() == YES) + output->SetHeat_InverseDesign(solver_container[ZONE_0][MESH_0][FLOW_SOL], + geometry_container[ZONE_0][MESH_0], config_container[ZONE_0], ExtIter); + + /*--- Perform a single iteration of the chosen PDE solver. ---*/ + + /*--- Run a single iteration of the problem using the driver class. ---*/ + + driver->Run(iteration_container, output, integration_container, + geometry_container, solver_container, numerics_container, + config_container, surface_movement, grid_movement, FFDBox, + interpolator_container, transfer_container); + + + /*--- Synchronization point after a single solver iteration. Compute the + wall clock time required. ---*/ + +#ifndef HAVE_MPI + StopTime = su2double(clock())/su2double(CLOCKS_PER_SEC); +#else + StopTime = MPI_Wtime(); +#endif + + UsedTime = (StopTime - StartTime); + + /*--- For specific applications, evaluate and plot the equivalent area. ---*/ + + if (config_container[ZONE_0]->GetEquivArea() == YES) { + output->SetEquivalentArea(solver_container[ZONE_0][MESH_0][FLOW_SOL], + geometry_container[ZONE_0][MESH_0], config_container[ZONE_0], ExtIter); + } + + /*--- Check if there is any change in the runtime parameters ---*/ + + CConfig *runtime = NULL; + strcpy(runtime_file_name, "runtime.dat"); + runtime = new CConfig(runtime_file_name, config_container[ZONE_0]); + runtime->SetExtIter(ExtIter); + + /*--- Update the convergence history file (serial and parallel computations). ---*/ + + if (!fsi){ + + output->SetConvHistory_Body(&ConvHist_file, geometry_container, solver_container, + config_container, integration_container, false, UsedTime, ZONE_0); + + } + + + /*--- Evaluate the new CFL number (adaptive). ---*/ + + + if ( config_container[ZONE_0]->GetLocal_Relax_Factor() ) { + //output->SetLocalCFL_Number(solver_container, config_container, geometry_container, ZONE_0); + } + else if (config_container[ZONE_0]->GetCFL_Adapt() == YES) { + output->SetCFL_Number(solver_container, config_container, ZONE_0); + } + + + + + + /*--- Check whether the current simulation has reached the specified + convergence criteria, and set StopCalc to true, if so. ---*/ + + switch (config_container[ZONE_0]->GetKind_Solver()) { + case EULER: case NAVIER_STOKES: case RANS: + StopCalc = integration_container[ZONE_0][FLOW_SOL]->GetConvergence(); break; + case WAVE_EQUATION: + StopCalc = integration_container[ZONE_0][WAVE_SOL]->GetConvergence(); break; + case HEAT_EQUATION: + StopCalc = integration_container[ZONE_0][HEAT_SOL]->GetConvergence(); break; + case FEM_ELASTICITY: + StopCalc = integration_container[ZONE_0][FEA_SOL]->GetConvergence(); break; + case ADJ_EULER: case ADJ_NAVIER_STOKES: case ADJ_RANS: + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + StopCalc = integration_container[ZONE_0][ADJFLOW_SOL]->GetConvergence(); break; + } + + /*--- Solution output. Determine whether a solution needs to be written + after the current iteration, and if so, execute the output file writing + routines. ---*/ + + if ((ExtIter+1 >= config_container[ZONE_0]->GetnExtIter()) + + || + + ((ExtIter % config_container[ZONE_0]->GetWrt_Sol_Freq() == 0) && (ExtIter != 0) && + !((config_container[ZONE_0]->GetUnsteady_Simulation() == DT_STEPPING_1ST) || + (config_container[ZONE_0]->GetUnsteady_Simulation() == DT_STEPPING_2ND) || + (config_container[ZONE_0]->GetUnsteady_Simulation() == TIME_STEPPING))) + + || + + (StopCalc) + + || + + (((config_container[ZONE_0]->GetUnsteady_Simulation() == DT_STEPPING_1ST) || + (config_container[ZONE_0]->GetUnsteady_Simulation() == TIME_STEPPING)) && + ((ExtIter == 0) || (ExtIter % config_container[ZONE_0]->GetWrt_Sol_Freq_DualTime() == 0))) + + || + + ((config_container[ZONE_0]->GetUnsteady_Simulation() == DT_STEPPING_2ND) && (!fsi) && + ((ExtIter == 0) || ((ExtIter % config_container[ZONE_0]->GetWrt_Sol_Freq_DualTime() == 0) || + ((ExtIter-1) % config_container[ZONE_0]->GetWrt_Sol_Freq_DualTime() == 0)))) + + || + + ((config_container[ZONE_0]->GetUnsteady_Simulation() == DT_STEPPING_2ND) && (fsi) && + ((ExtIter == 0) || ((ExtIter % config_container[ZONE_0]->GetWrt_Sol_Freq_DualTime() == 0)))) + + || + + (((config_container[ZONE_0]->GetDynamic_Analysis() == DYNAMIC) && + ((ExtIter == 0) || (ExtIter % config_container[ZONE_0]->GetWrt_Sol_Freq_DualTime() == 0))))){ + + + /*--- Low-fidelity simulations (using a coarser multigrid level + approximation to the solution) require an interpolation back to the + finest grid. ---*/ + + if (config_container[ZONE_0]->GetLowFidelitySim()) { + integration_container[ZONE_0][FLOW_SOL]->SetProlongated_Solution(RUNTIME_FLOW_SYS, solver_container[ZONE_0][MESH_0][FLOW_SOL], solver_container[ZONE_0][MESH_1][FLOW_SOL], geometry_container[ZONE_0][MESH_0], geometry_container[ZONE_0][MESH_1], config_container[ZONE_0]); + integration_container[ZONE_0][FLOW_SOL]->Smooth_Solution(RUNTIME_FLOW_SYS, solver_container[ZONE_0][MESH_0][FLOW_SOL], geometry_container[ZONE_0][MESH_0], 3, 1.25, config_container[ZONE_0]); + solver_container[ZONE_0][MESH_0][config_container[ZONE_0]->GetContainerPosition(RUNTIME_FLOW_SYS)]->Set_MPI_Solution(geometry_container[ZONE_0][MESH_0], config_container[ZONE_0]); + solver_container[ZONE_0][MESH_0][config_container[ZONE_0]->GetContainerPosition(RUNTIME_FLOW_SYS)]->Preprocessing(geometry_container[ZONE_0][MESH_0], solver_container[ZONE_0][MESH_0], config_container[ZONE_0], MESH_0, 0, RUNTIME_FLOW_SYS, true); + } + + + if (rank == MASTER_NODE) cout << endl << "-------------------------- File Output Summary --------------------------"; + + /*--- Execute the routine for writing restart, volume solution, + surface solution, and surface comma-separated value files. ---*/ + + output->SetResult_Files(solver_container, geometry_container, config_container, ExtIter, nZone); + + /*--- Output a file with the forces breakdown. ---*/ + + output->SetForces_Breakdown(geometry_container, solver_container, + config_container, integration_container, ZONE_0); + + /*--- Compute the forces at different sections. ---*/ + + if (config_container[ZONE_0]->GetPlot_Section_Forces()) { + output->SetForceSections(solver_container[ZONE_0][MESH_0][FLOW_SOL], + geometry_container[ZONE_0][MESH_0], config_container[ZONE_0], ExtIter); + } + + if (rank == MASTER_NODE) cout << "-------------------------------------------------------------------------" << endl << endl; + + } + + /*--- If the convergence criteria has been met, terminate the simulation. ---*/ + + //if ( StopCalc || ExtIter == config_container[ZONE_0]->GetnExtIter()-1 ) { + // output->ComputeNozzleThrust(solver_container[ZONE_0][MESH_0][FLOW_SOL], + // geometry_container[ZONE_0][MESH_0], config_container[ZONE_0]); + //} + + + if ( StopCalc || ExtIter == config_container[ZONE_0]->GetnExtIter()-1 ) { + + /*--- Compute thrust ---*/ + if ( config_container[ZONE_0]->GetKind_ObjFunc() == THRUST_NOZZLE ) { + + //printf("COMPUTE NOZZLE THRUST\n"); + output->SetNozzleThrust(solver_container[ZONE_0][MESH_0][FLOW_SOL],geometry_container[ZONE_0][MESH_0], config_container[ZONE_0]); + + if (rank == MASTER_NODE){ + //char cstr[1024] = config_container[ZONE_0]->GetThrust_FileName(); + //std::sprintf(cstr,"thrust"); + ofstream file; + file.precision(20); + file.open(config_container[ZONE_0]->GetThrust_FileName().c_str(), ios::out); + file << solver_container[ZONE_0][MESH_0][FLOW_SOL]->GetThrust_Nozzle() << endl; + file.close(); + } + } + + } + + if (StopCalc) break; + + if(ExtIter>=(config_container[ZONE_0]->GetnExtIter()-1)){ + + int Strainsize = config->StrainFile.IndexCurr.size(); + int Tausize = config->TauFile.IndexCurr.size(); + int Ssize[size]; + int Tsize[size]; + int Scumlf[size]; + int Tcumlf[size]; + int Ssizef=0; + int Tsizef=0; + MPI_Allgather(&Strainsize,1,MPI_INT,Ssize,1,MPI_INT,MPI_COMM_WORLD); + MPI_Allgather(&Tausize,1,MPI_INT,Tsize,1,MPI_INT,MPI_COMM_WORLD); + + for(int i=0; iStrainFile.IndexCurr[j]; + Sp1p[j] = config->StrainFile.p1[j]; + Swalldistp[j] = config->StrainFile.walldist[j]; + SIndexBndyp[j] = config->StrainFile.IndexBndy[j]; + Sstrain_ratep[j] = config->StrainFile.strain_rate[j]; + Smu_tp[j] = config->StrainFile.mu_t[j]; + + } + + unsigned long TIndexCurrp[Tausize]; + double TTauTangentp[Tausize]; + + for(int j=0; jTauFile.IndexCurr[j]; + TTauTangentp[j]= config->TauFile.TauTangent[j]; + + } + + +// /////////////////////////////////////////////////////////////////////////////// + + + unsigned long* SIndexCurr; + double* Sp1; + double* Swalldist; + unsigned long* SIndexBndy; + double* Sstrain_rate; + double* Smu_t; + unsigned long* TIndexCurr; + double* TTauTangent; + + //if (rank==MASTER_NODE){ + + SIndexCurr = new unsigned long[Ssizef]; + Sp1 = new double[Ssizef]; + Swalldist = new double[Ssizef]; + SIndexBndy = new unsigned long[Ssizef]; + Sstrain_rate = new double[Ssizef]; + Smu_t = new double[Ssizef]; + TIndexCurr = new unsigned long[Tsizef]; + TTauTangent = new double[Tsizef]; + + //} + + MPI_Allgatherv( SIndexCurrp, Strainsize, MPI_UNSIGNED_LONG, SIndexCurr, Ssize, Scumlf, MPI_UNSIGNED_LONG, MPI_COMM_WORLD); + MPI_Allgatherv( Sp1p, Strainsize, MPI_DOUBLE, Sp1, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( Swalldistp, Strainsize, MPI_DOUBLE, Swalldist, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( SIndexBndyp, Strainsize, MPI_UNSIGNED_LONG, SIndexBndy, Ssize, Scumlf, MPI_UNSIGNED_LONG, MPI_COMM_WORLD); + MPI_Allgatherv( Sstrain_ratep, Strainsize, MPI_DOUBLE, Sstrain_rate, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( Smu_tp, Strainsize, MPI_DOUBLE, Smu_t, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( TIndexCurrp, Tausize, MPI_UNSIGNED_LONG, TIndexCurr, Tsize, Tcumlf, MPI_UNSIGNED_LONG, MPI_COMM_WORLD); + MPI_Allgatherv( TTauTangentp, Tausize, MPI_DOUBLE, TTauTangent, Tsize, Tcumlf, MPI_DOUBLE, MPI_COMM_WORLD); + + //if (rank==MASTER_NODE){ + + //int filenum=0; + + cout<<"Creating Features file for ML"<>PointID>>p1>>p3>>Neighbor>>strain_rate>>mu_t){ + p1f[PointID]=p1; + p3f[PointID]=p3; + srf[PointID]=strain_rate; + muf[PointID]=mu_t; + neighbf[PointID]=Neighbor; + } + } + else break; + infile.close(); + char command[50]; + sprintf(command, "rm Feature%d", filenum); + //stat = system(command); + + sprintf(Buffer_File, "Tau%d", filenum); + infile.open(Buffer_File); + double TauWall; + if (infile){ + while(infile>>PointID>>TauWall){ + tauwallf[PointID]=TauWall; + } + } + infile.close(); + sprintf(command, "rm Tau%d", filenum); + //stat = system(command); + + filenum++; + } + */ + + for(int i=0; iGetNonphysical_Points() > 0) + cout << "Warning: there are " << config_container[ZONE_0]->GetNonphysical_Points() << " non-physical points in the solution." << endl; + if (config_container[ZONE_0]->GetNonphysical_Reconstr() > 0) + cout << "Warning: " << config_container[ZONE_0]->GetNonphysical_Reconstr() << " reconstructed states for upwinding are non-physical." << endl; + + /*--- Close the convergence history file. ---*/ + + ConvHist_file.close(); + cout << "History file, closed." << endl; + } + + /*--- Deallocations: further work is needed, + * these routines can be used to check for memory leaks---*/ + /* + if (rank == MASTER_NODE) + cout << endl <<"------------------------ Driver Postprocessing ------------------------" << endl; + + driver->Postprocessing(iteration_container, solver_container, geometry_container, + integration_container, numerics_container, interpolator_container, + transfer_container, config_container, nZone); + + delete driver; + */ + + // Assemble all the vectors into these arrays + + + + /*--- Geometry class deallocation ---*/ + if (rank == MASTER_NODE) + cout << endl <<"------------------------ Geometry Postprocessing ------------------------" << endl; + for (iZone = 0; iZone < nZone; iZone++) { + if (geometry_container[iZone]!=NULL){ + for (unsigned short iMGlevel = 1; iMGlevel < config_container[iZone]->GetnMGLevels()+1; iMGlevel++){ + if (geometry_container[iZone][iMGlevel]!=NULL) delete geometry_container[iZone][iMGlevel]; + } + delete [] geometry_container[iZone]; + } + } + delete [] geometry_container; + + /*--- Free-form deformation class deallocation ---*/ + for (iZone = 0; iZone < nZone; iZone++) { + delete FFDBox[iZone]; + } + delete [] FFDBox; + + /*--- Grid movement and surface movement class deallocation ---*/ + delete [] surface_movement; + delete [] grid_movement; + + /*Deallocate config container*/ + if (rank == MASTER_NODE) + cout << endl <<"------------------------ Config Postprocessing ------------------------" << endl; + if (config_container!=NULL){ + for (iZone = 0; iZone < nZone; iZone++) { + if (config_container[iZone]!=NULL){ + delete config_container[iZone]; + } + } + delete [] config_container; + } + + /*--- Deallocate output container ---*/ + if (output!=NULL) delete output; + + /*--- Synchronization point after a single solver iteration. Compute the + wall clock time required. ---*/ + +#ifndef HAVE_MPI + StopTime = su2double(clock())/su2double(CLOCKS_PER_SEC); +#else + StopTime = MPI_Wtime(); +#endif + + /*--- 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_CFD) ------------------------" << endl << endl; + +#ifdef HAVE_MPI + /*--- Finalize MPI parallelization ---*/ + MPI_Buffer_detach(&bptr, &bl); + MPI_Finalize(); +#endif + + return EXIT_SUCCESS; + +} diff --git a/SU2_CFD/src/SU2_CFD.cpp b/SU2_CFD/src/SU2_CFD.cpp index 6322b50b7e1..98ad4737eb9 100644 --- a/SU2_CFD/src/SU2_CFD.cpp +++ b/SU2_CFD/src/SU2_CFD.cpp @@ -163,6 +163,25 @@ int main(int argc, char *argv[]) { geometry_container[iZone][MESH_0]->SetBoundaries(config_container[iZone]); } + + ifstream infile("beta_for_su2.dat"); + unsigned long counter = geometry_container[ZONE_0][MESH_0]->GetGlobal_nPointDomain(); + double *val_betaArr = new double [geometry_container[ZONE_0][MESH_0]->GetGlobal_nPointDomain()]; + if (infile){ + for(int i=0; iGetGlobal_nPointDomain(); i++){ + int index; + infile>>index; + infile>>val_betaArr[index]; + } + } + else { + for(int i=0; iGetGlobal_nPointDomain(); i++){ + val_betaArr[i]=SU2_TYPE::GetValue(config_container[ZONE_0]->GetSA_Production_Factor()); + } + } + infile.close(); + config_container[ZONE_0]->SetbetaArr(val_betaArr); + if (rank == MASTER_NODE) cout << endl <<"------------------------- Geometry Preprocessing ------------------------" << endl; @@ -568,6 +587,192 @@ int main(int argc, char *argv[]) { delete driver; */ + // Assemble all the vectors into these arrays + + int Strainsize = config->StrainFile.IndexCurr.size(); + int Tausize = config->TauFile.IndexCurr.size(); + int Ssize[size]; + int Tsize[size]; + int Scumlf[size]; + int Tcumlf[size]; + int Ssizef=0; + int Tsizef=0; + MPI_Allgather(&Strainsize,1,MPI_INT,Ssize,1,MPI_INT,MPI_COMM_WORLD); + MPI_Allgather(&Tausize,1,MPI_INT,Tsize,1,MPI_INT,MPI_COMM_WORLD); + + for(int i=0; iStrainFile.IndexCurr[j]; + Sp1p[j] = config->StrainFile.p1[j]; + Swalldistp[j] = config->StrainFile.walldist[j]; + SIndexBndyp[j] = config->StrainFile.IndexBndy[j]; + Sstrain_ratep[j] = config->StrainFile.strain_rate[j]; + Smu_tp[j] = config->StrainFile.mu_t[j]; + + } + + unsigned long TIndexCurrp[Tausize]; + double TTauTangentp[Tausize]; + + for(int j=0; jTauFile.IndexCurr[j]; + TTauTangentp[j]= config->TauFile.TauTangent[j]; + + } + + +///////////////////////////////////////////////////////////////////////////////// + + + unsigned long* SIndexCurr; + double* Sp1; + double* Swalldist; + unsigned long* SIndexBndy; + double* Sstrain_rate; + double* Smu_t; + unsigned long* TIndexCurr; + double* TTauTangent; + + //if (rank==MASTER_NODE){ + + SIndexCurr = new unsigned long[Ssizef]; + Sp1 = new double[Ssizef]; + Swalldist = new double[Ssizef]; + SIndexBndy = new unsigned long[Ssizef]; + Sstrain_rate = new double[Ssizef]; + Smu_t = new double[Ssizef]; + TIndexCurr = new unsigned long[Tsizef]; + TTauTangent = new double[Tsizef]; + + //} + + MPI_Allgatherv( SIndexCurrp, Strainsize, MPI_UNSIGNED_LONG, SIndexCurr, Ssize, Scumlf, MPI_UNSIGNED_LONG, MPI_COMM_WORLD); + MPI_Allgatherv( Sp1p, Strainsize, MPI_DOUBLE, Sp1, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( Swalldistp, Strainsize, MPI_DOUBLE, Swalldist, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( SIndexBndyp, Strainsize, MPI_UNSIGNED_LONG, SIndexBndy, Ssize, Scumlf, MPI_UNSIGNED_LONG, MPI_COMM_WORLD); + MPI_Allgatherv( Sstrain_ratep, Strainsize, MPI_DOUBLE, Sstrain_rate, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( Smu_tp, Strainsize, MPI_DOUBLE, Smu_t, Ssize, Scumlf, MPI_DOUBLE, MPI_COMM_WORLD); + MPI_Allgatherv( TIndexCurrp, Tausize, MPI_UNSIGNED_LONG, TIndexCurr, Tsize, Tcumlf, MPI_UNSIGNED_LONG, MPI_COMM_WORLD); + MPI_Allgatherv( TTauTangentp, Tausize, MPI_DOUBLE, TTauTangent, Tsize, Tcumlf, MPI_DOUBLE, MPI_COMM_WORLD); + + //if (rank==MASTER_NODE){ + + //int filenum=0; + + cout<<"Creating Features file for ML"<>PointID>>p1>>p3>>Neighbor>>strain_rate>>mu_t){ + p1f[PointID]=p1; + p3f[PointID]=p3; + srf[PointID]=strain_rate; + muf[PointID]=mu_t; + neighbf[PointID]=Neighbor; + } + } + else break; + infile.close(); + char command[50]; + sprintf(command, "rm Feature%d", filenum); + //stat = system(command); + + sprintf(Buffer_File, "Tau%d", filenum); + infile.open(Buffer_File); + double TauWall; + if (infile){ + while(infile>>PointID>>TauWall){ + tauwallf[PointID]=TauWall; + } + } + infile.close(); + sprintf(command, "rm Tau%d", filenum); + //stat = system(command); + + filenum++; + } + */ + + for(int i=0; iGetExtIter()==(config->GetnExtIter()-1)){ + //char buffer_file[20]; + //sprintf(buffer_file, "Tau%d", proc_rank); + //ofstream outfile(buffer_file, ofstream::app); + //outfile<node[iPoint]->GetGlobalIndex()<<'\t'<TauFile.IndexCurr.push_back(geometry->node[iPoint]->GetGlobalIndex()); + config->TauFile.TauTangent.push_back(sqrt(SU2_TYPE::GetValue(TauTangent[0]*TauTangent[0]+TauTangent[1]*TauTangent[1]))); + } + for (iDim = 0; iDim < nDim; iDim++) WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]); WallDistMod = 0.0; for (iDim = 0; iDim < nDim; iDim++) WallDistMod += WallDist[iDim]*WallDist[iDim]; WallDistMod = sqrt(WallDistMod); diff --git a/SU2_CFD/src/solver_direct_turbulent.cpp b/SU2_CFD/src/solver_direct_turbulent.cpp index cb1d1fe503b..dbbc0bfb539 100644 --- a/SU2_CFD/src/solver_direct_turbulent.cpp +++ b/SU2_CFD/src/solver_direct_turbulent.cpp @@ -1373,8 +1373,8 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai bool freesurface = (config->GetKind_Regime() == FREESURFACE); bool time_spectral = (config->GetUnsteady_Simulation() == TIME_SPECTRAL); bool transition = (config->GetKind_Trans_Model() == LM); - su2double epsilon = config->GetFreeSurface_Thickness(); - + su2double epsilon = config->GetFreeSurface_Thickness(); + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { /*--- Conservative variables w/o reconstruction ---*/ @@ -1384,6 +1384,38 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai /*--- Gradient of the primitive and conservative variables ---*/ numerics->SetPrimVarGradient(solver_container[FLOW_SOL]->node[iPoint]->GetGradient_Primitive(), NULL); + + //int proc_rank = MASTER_NODE; + //#ifdef HAVE_MPI + // MPI_Comm_rank(MPI_COMM_WORLD, &proc_rank); + //#endif + if(config->GetExtIter()==(config->GetnExtIter()-1)){ + //char buffer_file[20]; + //sprintf(buffer_file, "Feature%d", proc_rank); + //ofstream outfile(buffer_file, ofstream::app); + double dudx = SU2_TYPE::GetValue(solver_container[FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[1][0]); + double dudy = SU2_TYPE::GetValue(solver_container[FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[1][1]); + double dvdx = SU2_TYPE::GetValue(solver_container[FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[2][0]); + double dvdy = SU2_TYPE::GetValue(solver_container[FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[2][1]); + double sxy = 0.5*(dudy+dvdx); + double strain_rate = sqrt(2*(dudx*dudx+dvdy*dvdy+2.0*sxy*sxy)); + double visc = SU2_TYPE::GetValue(solver_container[FLOW_SOL]->node[iPoint]->GetLaminarViscosity()); + double dens = SU2_TYPE::GetValue(solver_container[FLOW_SOL]->node[iPoint]->GetDensity()); + double mu_t = SU2_TYPE::GetValue(solver_container[FLOW_SOL]->node[iPoint]->GetEddyViscosity()); + long IndexCurr = geometry->node[iPoint]->GetGlobalIndex(); + long IndexBndy = geometry->node[iPoint]->GetVertex_nearWall(); + double walldist = SU2_TYPE::GetValue(geometry->node[iPoint]->GetWall_Distance()); + double p1 = dens*strain_rate*walldist*walldist/visc; + + config->StrainFile.IndexCurr.push_back(IndexCurr); + config->StrainFile.p1.push_back(p1); + config->StrainFile.walldist.push_back(walldist); + config->StrainFile.IndexBndy.push_back(IndexBndy); + config->StrainFile.strain_rate.push_back(strain_rate); + config->StrainFile.mu_t.push_back(mu_t); + //outfile<SetDistance(geometry->node[iPoint]->GetWall_Distance(), 0.0); /*--- Compute the source term ---*/ + + config->SetSA_Production_Factor(config->GetbetaArr(geometry->node[iPoint]->GetGlobalIndex())); numerics->ComputeResidual(Residual, Jacobian_i, NULL, config);