Skip to content

Commit

Permalink
Refactored CSolver
Browse files Browse the repository at this point in the history
  • Loading branch information
@maxi120
maxi120 committed Apr 29, 2024
1 parent ef9285b commit d94de4e
Show file tree
Hide file tree
Showing 3 changed files with 97 additions and 74 deletions.
20 changes: 17 additions & 3 deletions Common/include/toolboxes/geometry_toolbox.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -236,6 +236,18 @@ inline bool IntersectEdge(Int nDim, const T* p0, const T* n, const T* p1, const
else
return false;
}

/*!
* \brief Check if a point is inside a polygon
* \param[in] nDim - Number of dimensions of the particular problem.
* \param[in] vector - Vector to normalize.
*/
template <class T, class Int>
inline void NormalizeVector(Int nDim, T* vector) {
T norm = Norm(nDim, vector);
for (Int i = 0; i < nDim; i++) vector[i] /= norm;
}

/*!
* \brief Check if a point is inside a polygon
* \param[in] nDim - Number of dimensions of the particular problem.
Expand All @@ -257,11 +269,12 @@ inline bool PointInConvexPolygon(Int nDim, const Mat& pVert, const T* p0, int nV
if (nDim == 3) {
T plane_norm[3];
TriangleNormal(pVert, plane_norm);
auto normPlane_norm = Norm(3, plane_norm);
NormalizeVector(3, plane_norm);

for (unsigned short iDim = 0; iDim < nDim; iDim++) {
if (abs(plane_norm[iDim] / normPlane_norm) > max_proj) {
if (abs(plane_norm[iDim]) > max_proj) {
proj_index = iDim;
max_proj = abs(plane_norm[iDim] / normPlane_norm);
max_proj = abs(plane_norm[iDim]);
}
}

Expand Down Expand Up @@ -296,6 +309,7 @@ inline bool PointInConvexPolygon(Int nDim, const Mat& pVert, const T* p0, int nV

return nIntersections % 2 == 1;
}

// end added by max
/// @}
} // namespace GeometryToolbox
4 changes: 3 additions & 1 deletion SU2_CFD/include/solvers/CSolver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1583,7 +1583,9 @@ class CSolver {
//Added by max
void ReadVGConfigFile(CConfig* config);

//End added by max
void InitializeVGVariables(unsigned short nVgs_file, std::vector<std::string>& lines_configVg, CConfig* config);

// End added by max

/*!
* \brief A virtual member.
Expand Down
147 changes: 77 additions & 70 deletions SU2_CFD/src/solvers/CSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4316,29 +4316,16 @@ void CSolver::SavelibROM(CGeometry *geometry, CConfig *config, bool converged) {
//Added by max
void CSolver::ReadVGConfigFile(CConfig* config) {
string filename, line;
unsigned short nVgs_file = 0;
su2double **vg_b, **vg_t, **vg_n, **vg_uhat, ***vg_coord, *Svg, *betaVg;
su2double bNorm, nNorm, tNorm;
const unsigned short p1_idx = 4, l_idx = 0, h_idx = 1, nOpt = 13, un_idx = 7, u_idx = 10;
vector<string> lines_configVg;
int restart_iter;
unsigned short nVgs = 0;
vector<string> lines_vgConfig;
bool unsteady_restart = config->GetTime_Domain() && config->GetRestart();
bool ts_2nd = config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND;

/*--- Get the vortex generator filename. ---*/

if (config->GetTime_Domain() && config->GetRestart()) {

switch (config->GetTime_Marching())
{
case TIME_MARCHING::DT_STEPPING_1ST:
restart_iter = config->GetRestart_Iter() - 1;
break;
case TIME_MARCHING::DT_STEPPING_2ND:
restart_iter = config->GetRestart_Iter() - 2;
break;
default:
break;
}
if (unsteady_restart) {

int restart_iter = config->GetRestart_Iter()-1-ts_2nd;
filename = config->GetUnsteady_FileName(config->GetVGFileName(), restart_iter, ".cfg");

} else {
Expand All @@ -4351,87 +4338,107 @@ void CSolver::ReadVGConfigFile(CConfig* config) {

while (std::getline(file, line)) {
if (line.empty() || line.front() == '#') continue;
nVgs_file++;
lines_configVg.push_back(line);
nVgs++;
lines_vgConfig.push_back(line);
};
file.close();

InitializeVGVariables(nVgs, lines_vgConfig, config);
}


void CSolver::InitializeVGVariables(unsigned short nVgs, std::vector<std::string> &lines_vgConfig, CConfig * config)
{
su2double l,h1,h2,beta,p1[3],un[3],u[3],uc[3];

/*--- Allocate and compute the required varibales for the vortex generator model ---*/

vg_b = new su2double*[nVgs_file];
vg_n = new su2double*[nVgs_file];
vg_t = new su2double*[nVgs_file];
vg_uhat = new su2double*[nVgs_file];
vg_coord = new su2double**[nVgs_file];
Svg = new su2double[nVgs_file];
betaVg = new su2double[nVgs_file];

for (unsigned short iVG = 0; iVG < nVgs_file; iVG++) {
istringstream vg_line{lines_configVg[iVG]};
su2double opt[13];
su2double tmp;
unsigned short iOpt;

for (iOpt = 0; iOpt < nOpt; iOpt++) {
vg_line >> tmp;
opt[iOpt] = tmp;
}
su2double **vg_b = new su2double*[nVgs];
su2double **vg_n = new su2double*[nVgs];
su2double **vg_t = new su2double*[nVgs];
su2double **vg_uhat = new su2double*[nVgs];
su2double ***vg_coord = new su2double**[nVgs];
su2double *Svg = new su2double[nVgs];
su2double *betaVg = new su2double[nVgs];

for (unsigned short iVG = 0; iVG < nVgs; iVG++) {

/* --- Parse variables ---*/
istringstream vg_line{lines_vgConfig[iVG]};

betaVg[iVG] = opt[3];
vg_b[iVG] = new su2double[3]{opt[un_idx], opt[un_idx + 1], opt[un_idx + 2]};
vg_uhat[iVG] = new su2double[3]{opt[u_idx], opt[u_idx + 1], opt[u_idx + 2]};
vg_n[iVG] = new su2double[3]{0};
vg_t[iVG] = new su2double[3]{0};
vg_line >> l;
vg_line >> h1;
vg_line >> h2;
vg_line >> beta;

const auto beta = betaVg[iVG] * PI_NUMBER / 180.0;
for (unsigned short iDim = 0; iDim < 3; iDim++) {
vg_line >> p1[iDim];
}

su2double uc[3];
GeometryToolbox::CrossProduct(vg_b[iVG], vg_uhat[iVG], uc);
for (unsigned short iDim = 0; iDim < 3; iDim++) {
vg_line >> un[iDim];
}

for (unsigned short iDim = 0; iDim < 3; iDim++) {
vg_t[iVG][iDim] = vg_uhat[iVG][iDim] * cos(beta) + uc[iDim] * sin(beta);
vg_n[iVG][iDim] = vg_uhat[iVG][iDim] * sin(beta) + uc[iDim] * cos(beta);
vg_line >> u[iDim];
}

/*--- Allocate Variable ---*/

vg_b[iVG] = new su2double[3];
vg_uhat[iVG] = new su2double[3];
vg_n[iVG] = new su2double[3];
vg_t[iVG] = new su2double[3];

/*--- Compute the VG variables ---*/

bNorm = GeometryToolbox::Norm(nDim, vg_b[iVG]);
nNorm = GeometryToolbox::Norm(nDim, vg_n[iVG]);
tNorm = GeometryToolbox::Norm(nDim, vg_t[iVG]);
beta *= PI_NUMBER / 180.0;

GeometryToolbox::CrossProduct(un, u, uc);

for (unsigned short iDim = 0; iDim < 3; iDim++) {
vg_t[iVG][iDim] /= tNorm;
vg_b[iVG][iDim] /= bNorm;
vg_n[iVG][iDim] /= nNorm;
vg_t[iVG][iDim] = u[iDim] * cos(beta) + uc[iDim] * sin(beta);
vg_n[iVG][iDim] = u[iDim] * sin(beta) + uc[iDim] * cos(beta);
vg_b[iVG][iDim] = un[iDim];
vg_uhat[iVG][iDim] = u[iDim];
}
betaVg[iVG] = beta;

GeometryToolbox::NormalizeVector(nDim, vg_t[iVG]);
GeometryToolbox::NormalizeVector(nDim, vg_b[iVG]);
GeometryToolbox::NormalizeVector(nDim, vg_n[iVG]);

/*--- Set variables in CConfig ---*/

vg_coord[iVG] = new su2double*[4];

vg_coord[iVG][0] = new su2double[3]{opt[p1_idx], opt[p1_idx + 1], opt[p1_idx + 2]};
vg_coord[iVG][0] = new su2double[3]{p1[0], p1[1], p1[2]};

vg_coord[iVG][1] = new su2double[3]{p1[0] + vg_t[iVG][0] * l,
p1[1] + vg_t[iVG][1] * l,
p1[2] + vg_t[iVG][2] * l};

vg_coord[iVG][1] = new su2double[3]{opt[p1_idx] + vg_t[iVG][0] * opt[l_idx],
opt[p1_idx + 1] + vg_t[iVG][1] * opt[l_idx],
opt[p1_idx + 2] + vg_t[iVG][2] * opt[l_idx]};
su2double* p2=vg_coord[iVG][1];

vg_coord[iVG][2] = new su2double[3]{vg_coord[iVG][1][0] + vg_b[iVG][0] * opt[h_idx],
vg_coord[iVG][1][1] + vg_b[iVG][1] * opt[h_idx],
vg_coord[iVG][1][2] + vg_b[iVG][2] * opt[h_idx]};
vg_coord[iVG][2] = new su2double[3]{p2[0] + vg_b[iVG][0] * h2,
p2[1] + vg_b[iVG][1] * h2,
p2[2] + vg_b[iVG][2] * h2};

vg_coord[iVG][3] = new su2double[3]{opt[p1_idx] + vg_b[iVG][0] * opt[h_idx + 1],
opt[p1_idx + 1] + vg_b[iVG][1] * opt[h_idx + 1],
opt[p1_idx + 2] + vg_b[iVG][2] * opt[h_idx + 1]};
vg_coord[iVG][3] = new su2double[3]{p1[0] + vg_b[iVG][0] * h1,
p1[1] + vg_b[iVG][1] * h1,
p1[2] + vg_b[iVG][2] * h1};

Svg[iVG] = 0.5 * (opt[h_idx] + opt[h_idx + 1]) * opt[l_idx];
Svg[iVG] = 0.5 * (h1 + h2) * l;
};

/*--- Set the variables in CConfig ---*/

config->Set_nVG(vg_n);
config->Set_bVG(vg_b);
config->Set_tVG(vg_t);
config->SetVGCoord(vg_coord);
config->Set_Svg(Svg);
config->Set_nVGs(nVgs_file);
config->Set_nVGs(nVgs);
config->Set_uhatVg(vg_uhat);
config->Set_betaVg(betaVg);
}

//End added by max
}//End added by max

0 comments on commit d94de4e

Please sign in to comment.