Skip to content

Commit 608cfb2

Browse files
thomasdickthomasdick
committed
Apply suggestions from code review.
1 parent db9d268 commit 608cfb2

2 files changed

Lines changed: 91 additions & 122 deletions

File tree

‎Common/include/option_structure.hpp‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1221,7 +1221,6 @@ enum BC_TYPE {
12211221
CHT_WALL_INTERFACE = 50, /*!< \brief Domain interface definition. */
12221222
SMOLUCHOWSKI_MAXWELL = 55, /*!< \brief Smoluchoski/Maxwell wall boundary condition. */
12231223
SEND_RECEIVE = 99, /*!< \brief Boundary send-receive definition. */
1224-
NOT_AVAILABLE = 1000, /*!< \brief Boundary not on this MPI rank. */
12251224
};
12261225

12271226
/*!

‎SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp‎

Lines changed: 91 additions & 121 deletions
Original file line numberDiff line numberDiff line change
@@ -25,9 +25,11 @@
2525
* License along with SU2. If not, see <http://www.gnu.org/licenses/>.
2626
*/
2727

28+
#include <algorithm>
29+
30+
#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
2831
#include "../../include/solvers/CGradientSmoothingSolver.hpp"
2932
#include "../../include/variables/CSobolevSmoothingVariable.hpp"
30-
#include <algorithm>
3133

3234
CGradientSmoothingSolver::CGradientSmoothingSolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(LINEAR_SOLVER_MODE::GRADIENT_MODE) {
3335
unsigned int iDim, marker_count = 0;
@@ -322,7 +324,7 @@ void CGradientSmoothingSolver::ApplyGradientSmoothingDV(CGeometry *geometry, CNu
322324
if (rank == MASTER_NODE) {
323325
/*--- For multizone append zone number to filename. ---*/
324326
string hess_filename = config->GetObjFunc_Hess_FileName();
325-
if (config->GetnZone() > 1) config->GetMultizone_FileName(hess_filename, config->GetiZone(), ".dat");
327+
if (config->GetMultizone_Problem()) config->GetMultizone_FileName(hess_filename, config->GetiZone(), ".dat");
326328
ofstream SysMatrix(hess_filename);
327329
SysMatrix.precision(config->GetOutput_Precision());
328330

@@ -416,62 +418,57 @@ void CGradientSmoothingSolver::Compute_Surface_StiffMatrix(CGeometry* geometry,
416418

417419
Jacobian_block.resize(nDim, nDim) = su2double(0.0);
418420

419-
/*--- Check if the current MPI rank has a part of the marker ---*/
420-
if (val_marker!=BC_TYPE::NOT_AVAILABLE) {
421-
422-
/*--- Loops over all the elements ---*/
423-
for (iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) {
421+
/*--- Loops over all the elements ---*/
422+
for (iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) {
424423

425-
/*--- Identify the kind of boundary element ---*/
426-
GetElemKindAndNumNodes(geometry->bound[val_marker][iElem]->GetVTK_Type(), EL_KIND, nNodes);
424+
/*--- Identify the kind of boundary element ---*/
425+
GetElemKindAndNumNodes(geometry->bound[val_marker][iElem]->GetVTK_Type(), EL_KIND, nNodes);
427426

428-
/*--- Retrieve the boundary reference and current coordinates ---*/
427+
/*--- Retrieve the boundary reference and current coordinates ---*/
429428

430-
for (iNode = 0; iNode < nNodes; iNode++) {
431-
indexNode[iNode] = geometry->bound[val_marker][iElem]->GetNode(iNode);
429+
for (iNode = 0; iNode < nNodes; iNode++) {
430+
indexNode[iNode] = geometry->bound[val_marker][iElem]->GetNode(iNode);
432431

433-
for (iDim = 0; iDim < nDim; iDim++) {
434-
val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
435-
element_container[GRAD_TERM][EL_KIND]->SetRef_Coord(iNode, iDim, val_Coord);
436-
}
432+
for (iDim = 0; iDim < nDim; iDim++) {
433+
val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
434+
element_container[GRAD_TERM][EL_KIND]->SetRef_Coord(iNode, iDim, val_Coord);
437435
}
436+
}
438437

439-
/*--- We need the indices of the vertices, which are "Dual Grid Info" ---*/
440-
for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
441-
iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
442-
for (iNode = 0; iNode < nNodes; iNode++) {
443-
if (iPoint == indexNode[iNode]) indexVertex[iNode] = iVertex;
444-
}
438+
/*--- We need the indices of the vertices, which are "Dual Grid Info" ---*/
439+
for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
440+
iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
441+
for (iNode = 0; iNode < nNodes; iNode++) {
442+
if (iPoint == indexNode[iNode]) indexVertex[iNode] = iVertex;
445443
}
444+
}
446445

447-
/*--- compute the contributions of the single elements inside the numerics container ---*/
448-
numerics->Compute_Tangent_Matrix(element_container[GRAD_TERM][EL_KIND], config);
446+
/*--- compute the contributions of the single elements inside the numerics container ---*/
447+
numerics->Compute_Tangent_Matrix(element_container[GRAD_TERM][EL_KIND], config);
449448

450-
NelNodes = element_container[GRAD_TERM][EL_KIND]->GetnNodes();
449+
NelNodes = element_container[GRAD_TERM][EL_KIND]->GetnNodes();
451450

452-
/*--- for all nodes add the contribution to the system Jacobian ---*/
451+
/*--- for all nodes add the contribution to the system Jacobian ---*/
453452

454-
for (iNode = 0; iNode < NelNodes; iNode++) {
455-
for (jNode = 0; jNode < NelNodes; jNode++) {
453+
for (iNode = 0; iNode < NelNodes; iNode++) {
454+
for (jNode = 0; jNode < NelNodes; jNode++) {
456455

457-
DHiDHj = element_container[GRAD_TERM][EL_KIND]->Get_DHiDHj(iNode, jNode);
458-
HiHj = element_container[GRAD_TERM][EL_KIND]->Get_HiHj(iNode, jNode);
456+
DHiDHj = element_container[GRAD_TERM][EL_KIND]->Get_DHiDHj(iNode, jNode);
457+
HiHj = element_container[GRAD_TERM][EL_KIND]->Get_HiHj(iNode, jNode);
459458

460-
/*--- Get the calculated stiffness blocks for the current element and add them to the block. ---*/
461-
for (iSurfDim=0; iSurfDim<nSurfDim; iSurfDim++) {
462-
Jacobian_block[iSurfDim][iSurfDim] = DHiDHj[0][0] + HiHj;
463-
}
464-
Jacobian.AddBlock(indexNode[iNode], indexNode[jNode], Jacobian_block);
459+
/*--- Get the calculated stiffness blocks for the current element and add them to the block. ---*/
460+
for (iSurfDim=0; iSurfDim<nSurfDim; iSurfDim++) {
461+
Jacobian_block[iSurfDim][iSurfDim] = DHiDHj[0][0] + HiHj;
462+
}
463+
Jacobian.AddBlock(indexNode[iNode], indexNode[jNode], Jacobian_block);
465464

466-
/*--- Store if we set values for a given node already ---*/
467-
if (visited[indexNode[iNode]]==false) { visited[indexNode[iNode]]=true; }
468-
if (visited[indexNode[jNode]]==false) { visited[indexNode[jNode]]=true; }
465+
/*--- Store if we set values for a given node already ---*/
466+
if (visited[indexNode[iNode]]==false) { visited[indexNode[iNode]]=true; }
467+
if (visited[indexNode[jNode]]==false) { visited[indexNode[jNode]]=true; }
469468

470-
}
471469
}
472470
}
473471
}
474-
475472
}
476473

477474
void CGradientSmoothingSolver::Compute_Residual(CGeometry* geometry, const CConfig* config) {
@@ -511,12 +508,12 @@ void CGradientSmoothingSolver::Compute_Residual(CGeometry* geometry, const CConf
511508

512509
for (unsigned int iNode = 0; iNode < nNodes; iNode++) {
513510
if (config->GetSmoothSepDim()) {
514-
Residual[GetCurrentDim()] += Weight * Jac_X * element_container[GRAD_TERM][EL_KIND]->GetNi(iNode,iGauss) * nodes->GetSensitivity(indexNode[iNode], GetCurrentDim());
511+
Residual[GetCurrentDim()] +=
512+
Weight * Jac_X * element_container[GRAD_TERM][EL_KIND]->GetNi(iNode, iGauss) * nodes->GetSensitivity(indexNode[iNode], GetCurrentDim());
515513
LinSysRes.AddBlock(indexNode[iNode], &Residual[GetCurrentDim()]);
516514
} else {
517515
for (iDim = 0; iDim < nDim; iDim++) {
518-
Residual[iDim] += Weight * Jac_X * element_container[GRAD_TERM][EL_KIND]->GetNi(iNode,iGauss) * nodes->GetSensitivity(indexNode[iNode], iDim);
519-
516+
Residual[iDim] += Weight * Jac_X * element_container[GRAD_TERM][EL_KIND]->GetNi(iNode, iGauss) * nodes->GetSensitivity(indexNode[iNode], iDim);
520517
}
521518
LinSysRes.AddBlock(indexNode[iNode], Residual);
522519
}
@@ -538,66 +535,55 @@ void CGradientSmoothingSolver::Compute_Surface_Residual(CGeometry* geometry, con
538535
std::array<unsigned long, MAXNNODE_2D> indexNode;
539536
std::array<unsigned long, MAXNNODE_2D> indexVertex;
540537
su2double Weight, Jac_X, normalSens = 0.0, norm, val_Coord;
541-
su2double* normal = NULL;
542-
543-
/*--- Check if the current MPI rank has a part of the marker ---*/
544-
if (val_marker!=BC_TYPE::NOT_AVAILABLE) {
545-
546-
for (iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) {
538+
su2double normal[MAXNDIM];
547539

548-
/*--- Identify the kind of boundary element ---*/
549-
GetElemKindAndNumNodes(geometry->bound[val_marker][iElem]->GetVTK_Type(), EL_KIND, nNodes);
540+
for (iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) {
541+
/*--- Identify the kind of boundary element ---*/
542+
GetElemKindAndNumNodes(geometry->bound[val_marker][iElem]->GetVTK_Type(), EL_KIND, nNodes);
550543

551-
/*--- Retrieve the boundary reference and current coordinates ---*/
552-
for (iNode = 0; iNode < nNodes; iNode++) {
553-
indexNode[iNode] = geometry->bound[val_marker][iElem]->GetNode(iNode);
544+
/*--- Retrieve the boundary reference and current coordinates ---*/
545+
for (iNode = 0; iNode < nNodes; iNode++) {
546+
indexNode[iNode] = geometry->bound[val_marker][iElem]->GetNode(iNode);
554547

555-
for (iDim = 0; iDim < nDim; iDim++) {
556-
val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
557-
element_container[GRAD_TERM][EL_KIND]->SetRef_Coord(iNode, iDim, val_Coord);
558-
}
548+
for (iDim = 0; iDim < nDim; iDim++) {
549+
val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
550+
element_container[GRAD_TERM][EL_KIND]->SetRef_Coord(iNode, iDim, val_Coord);
559551
}
552+
}
560553

561-
/*--- We need the indices of the vertices, which are "Dual Grid Info" ---*/
562-
for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
563-
iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
564-
for (iNode = 0; iNode < nNodes; iNode++) {
565-
if (iPoint == indexNode[iNode]) indexVertex[iNode] = iVertex;
566-
}
554+
/*--- We need the indices of the vertices, which are "Dual Grid Info" ---*/
555+
for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
556+
iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
557+
for (iNode = 0; iNode < nNodes; iNode++) {
558+
if (iPoint == indexNode[iNode]) indexVertex[iNode] = iVertex;
567559
}
560+
}
568561

569-
element_container[GRAD_TERM][EL_KIND]->ClearElement(); /*--- Restarts the element: avoids adding over previous results in other elements --*/
570-
element_container[GRAD_TERM][EL_KIND]->ComputeGrad_SurfaceEmbedded();
571-
unsigned int nGauss = element_container[GRAD_TERM][EL_KIND]->GetnGaussPoints();
572-
573-
for (unsigned int iGauss = 0; iGauss < nGauss; iGauss++) {
574-
575-
Weight = element_container[GRAD_TERM][EL_KIND]->GetWeight(iGauss);
576-
Jac_X = element_container[GRAD_TERM][EL_KIND]->GetJ_X(iGauss);
577-
578-
for (unsigned int iNode = 0; iNode < nNodes; iNode++) {
562+
element_container[GRAD_TERM][EL_KIND]
563+
->ClearElement(); /*--- Restarts the element: avoids adding over previous results in other elements --*/
564+
element_container[GRAD_TERM][EL_KIND]->ComputeGrad_SurfaceEmbedded();
565+
unsigned int nGauss = element_container[GRAD_TERM][EL_KIND]->GetnGaussPoints();
579566

580-
normal = geometry->vertex[val_marker][indexVertex[iNode]]->GetNormal();
581-
norm = 0.0;
582-
for (iDim = 0; iDim < nDim; iDim++) {
583-
norm += normal[iDim]*normal[iDim];
584-
}
585-
norm = sqrt(norm);
586-
for (iDim = 0; iDim < nDim; iDim++) {
587-
normal[iDim] = normal[iDim] / norm;
588-
}
567+
for (unsigned int iGauss = 0; iGauss < nGauss; iGauss++) {
568+
Weight = element_container[GRAD_TERM][EL_KIND]->GetWeight(iGauss);
569+
Jac_X = element_container[GRAD_TERM][EL_KIND]->GetJ_X(iGauss);
589570

590-
for (iDim = 0; iDim < nDim; iDim++) {
591-
normalSens += normal[iDim] * nodes->GetSensitivity(indexNode[iNode], iDim);
592-
}
571+
for (unsigned int iNode = 0; iNode < nNodes; iNode++) {
572+
geometry->vertex[val_marker][indexVertex[iNode]]->GetNormal(normal);
573+
norm = GeometryToolbox::Norm(nDim, normal);
574+
for (iDim = 0; iDim < nDim; iDim++) {
575+
normal[iDim] = normal[iDim] / norm;
576+
}
593577

594-
Residual[0] += Weight * Jac_X * element_container[GRAD_TERM][EL_KIND]->GetNi(iNode,iGauss) * normalSens;
595-
LinSysRes.AddBlock(indexNode[iNode], Residual);
578+
for (iDim = 0; iDim < nDim; iDim++) {
579+
normalSens += normal[iDim] * nodes->GetSensitivity(indexNode[iNode], iDim);
580+
}
596581

597-
Residual[0] = 0;
598-
normalSens = 0;
582+
Residual[0] += Weight * Jac_X * element_container[GRAD_TERM][EL_KIND]->GetNi(iNode, iGauss) * normalSens;
583+
LinSysRes.AddBlock(indexNode[iNode], Residual);
599584

600-
}
585+
Residual[0] = 0;
586+
normalSens = 0;
601587
}
602588
}
603589
}
@@ -638,18 +624,15 @@ void CGradientSmoothingSolver::BC_Surface_Dirichlet(const CGeometry* geometry, c
638624
unsigned long iPoint, iVertex;
639625
const su2double zeros[MAXNDIM] = {0.0};
640626

641-
/*--- Check if the current MPI rank has a part of the marker ---*/
642-
if (val_marker!=BC_TYPE::NOT_AVAILABLE) {
643-
for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
627+
for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
644628

645-
/*--- Get node index ---*/
646-
iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
629+
/*--- Get node index ---*/
630+
iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
647631

648-
if (nodes->GetIsBoundaryPoint(iPoint)) {
649-
LinSysSol.SetBlock(iPoint, zeros);
650-
LinSysRes.SetBlock(iPoint, zeros);
651-
Jacobian.EnforceSolutionAtNode(iPoint, zeros, LinSysRes);
652-
}
632+
if (nodes->GetIsBoundaryPoint(iPoint)) {
633+
LinSysSol.SetBlock(iPoint, zeros);
634+
LinSysRes.SetBlock(iPoint, zeros);
635+
Jacobian.EnforceSolutionAtNode(iPoint, zeros, LinSysRes);
653636
}
654637
}
655638
}
@@ -738,7 +721,7 @@ void CGradientSmoothingSolver::OutputDVGradient(const CConfig* config, string ou
738721
unsigned iDV;
739722
if (rank == MASTER_NODE) {
740723
/*--- For multizone append zone number to filename. ---*/
741-
if (config->GetnZone() > 1) config->GetMultizone_FileName(out_file, config->GetiZone(), ".dat");
724+
if (config->GetMultizone_Problem()) config->GetMultizone_FileName(out_file, config->GetiZone(), ".dat");
742725
ofstream delta_p (out_file);
743726
delta_p.precision(config->GetOutput_Precision());
744727

@@ -913,9 +896,9 @@ void CGradientSmoothingSolver::ReadSensFromGeometry(const CGeometry* geometry) {
913896
}
914897

915898
void CGradientSmoothingSolver::WriteSensitivity(CGeometry* geometry, const CConfig* config) {
916-
unsigned long iPoint, total_index;
899+
unsigned long iPoint;
917900
unsigned int iDim;
918-
su2double* normal;
901+
su2double normal[MAXNDIM];
919902
su2double norm;
920903

921904
/*--- Split between surface and volume first. ---*/
@@ -925,16 +908,10 @@ void CGradientSmoothingSolver::WriteSensitivity(CGeometry* geometry, const CConf
925908
if (config->GetMarker_All_DV(iMarker) == YES) {
926909
for (unsigned long iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
927910
iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
928-
normal = geometry->vertex[iMarker][iVertex]->GetNormal();
929-
norm = 0.0;
930-
for (iDim = 0; iDim < nDim; iDim++) {
931-
norm += normal[iDim] * normal[iDim];
932-
}
933-
norm = sqrt(norm);
911+
geometry->vertex[iMarker][iVertex]->GetNormal(normal);
912+
norm = GeometryToolbox::Norm(nDim, normal);
934913
for (iDim = 0; iDim < nDim; iDim++) {
935914
normal[iDim] = normal[iDim] / norm;
936-
}
937-
for (iDim = 0; iDim < nDim; iDim++) {
938915
nodes->SetSensitivity(iPoint, iDim, normal[iDim] * LinSysSol[iPoint]);
939916
}
940917
}
@@ -948,8 +925,7 @@ void CGradientSmoothingSolver::WriteSensitivity(CGeometry* geometry, const CConf
948925
} else {
949926
for (iPoint = 0; iPoint < nPoint; iPoint++) {
950927
for (iDim = 0; iDim < nDim; iDim++) {
951-
total_index = iPoint*nDim + iDim;
952-
nodes->SetSensitivity(iPoint, iDim, LinSysSol[total_index]);
928+
nodes->SetSensitivity(iPoint, iDim, LinSysSol(iPoint, iDim));
953929
}
954930
}
955931
}
@@ -1027,16 +1003,10 @@ void CGradientSmoothingSolver::Set_VertexEliminationSchedule(CGeometry* geometry
10271003

10281004
void CGradientSmoothingSolver::Complete_Surface_StiffMatrix(const CGeometry* geometry) {
10291005

1030-
su2activematrix mId_Aux;
1031-
mId_Aux.resize(nDim, nDim) = su2double(0.0);
1032-
for (unsigned iDim = 0; iDim < nDim; iDim++) {
1033-
mId_Aux[iDim][iDim] = su2double(1.0);
1034-
}
1035-
10361006
/*--- Assembling the stiffness matrix on the design surface means the Jacobian is the identity for nodes inside the domain. ---*/
10371007
for (unsigned long iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++){
10381008
if (visited[iPoint]==false) {
1039-
Jacobian.AddBlock(iPoint, iPoint, mId_Aux);
1009+
Jacobian.AddVal2Diag(iPoint, 1.0);
10401010
}
10411011
}
10421012

0 commit comments

Comments
 (0)