Another charge against pointer to pointer #1312

pcarruscag · 2021-06-27T12:29:58Z

Proposed Changes

Introduce a "matrix view" type to avoid ** for example when passing gradients of primitives into the numerics.

PR Checklist

I am submitting my contribution to the develop branch.
My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags, or simply --warnlevel=2 when using meson).
My contribution is commented and consistent with SU2 style.
I have added a test case that demonstrates my contribution, if necessary.
I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp) , if necessary.

Common/include/containers/container_decorators.hpp

pcarruscag · 2021-06-27T12:37:21Z

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

-        visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), Grad_Reflected);
+        visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), CMatrixView<su2double>(Grad_Reflected));


Matrix types have an explicit conversion to this "matrix view", i.e. you can pass su2activematrix to CNumerics instead of a su2double**

pcarruscag · 2021-06-27T12:44:18Z

SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp

-    CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow)+1,Viscosity);
+    CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow,1), Viscosity);


Preferred way to get the gradient with an offset, although the "+" arithmetic still works.

pcarruscag · 2021-06-27T12:44:53Z

SU2_CFD/src/output/CFlowCompOutput.cpp

-    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(&(Node_Flow->GetGradient_Primitive(iPoint)[1])));
+    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1)));


But this (taking the address) will no longer work.

pcarruscag · 2021-06-27T12:47:18Z

Common/include/containers/container_decorators.hpp

- * \brief This contrived container is used to store small matrices in a contiguous manner
- *        but still present the "su2double**" interface to the outside world.
- *        The "interface" part should be replaced by something more efficient, e.g. a "matrix view".
- */
-class CVectorOfMatrix: public C3DDoubleMatrix {
-private:
-  su2matrix<Scalar*> interface;


The main motivation for this change is to store "vectors of matrices" e.g. gradients in a contiguous way without needing an extra matrix of pointers to provide the su2double** compatibility with other areas of the code.

SU2_CFD/include/numerics/CNumerics.hpp

pcarruscag · 2021-06-27T17:08:18Z

SU2_CFD/include/numerics/CNumerics.hpp

@@ -659,6 +627,58 @@ class CNumerics {
    }
  }

+  /*!
+   * \brief Project average gradient onto normal (with or w/o correction) for viscous fluxes of scalar quantities.


To reduce duplication

…radient_type

WallyMaier

@pcarruscag thanks for the work!

This looks good to me.

WallyMaier · 2021-06-28T17:42:00Z

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1053,8 +1053,7 @@ void CFVMFlowSolverBase<V, R>::BC_Sym_Plane(CGeometry* geometry, CSolver** solve
      Tangential[MAXNDIM] = {0.0}, GradNormVel[MAXNDIM] = {0.0}, GradTangVel[MAXNDIM] = {0.0};

  /*--- Allocation of primitive gradient arrays for viscous fluxes. ---*/
-  su2double** Grad_Reflected = new su2double*[nPrimVarGrad];
-  for (iVar = 0; iVar < nPrimVarGrad; iVar++) Grad_Reflected[iVar] = new su2double[nDim];
+  su2activematrix Grad_Reflected(nPrimVarGrad, nDim);


This looks much nicer, and less of a memory headache. 👍

WallyMaier · 2021-06-28T17:45:14Z

SU2_CFD/src/numerics/NEMO/convection/msw.cpp

@@ -198,7 +198,8 @@ CNumerics::ResidualType<> CUpwMSW_NEMO::ComputeResidual(const CConfig *config) {
                                                      epsilon*epsilon));

  /*--- Compute projected P, invP, and Lambda ---*/
-  CreateBasis(UnitNormal);
+  su2double l[MAXNDIM], m[MAXNDIM];


thanks for untangling some of this stuff. Its on the list to go through and completely abstract away....eventually.

WallyMaier · 2021-06-28T17:47:58Z

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -674,7 +674,7 @@ CNumerics::ResidualType<> CSourceWindGust::ComputeResidual(const CConfig* config
  smx = rho*(du_gust_dt + (u+u_gust)*du_gust_dx + (v+v_gust)*du_gust_dy);
  smy = rho*(dv_gust_dt + (u+u_gust)*dv_gust_dx + (v+v_gust)*dv_gust_dy);
  //smz = rho*(dw_gust_dt + (u+u_gust)*dw_gust_dx + (v+v_gust)*dw_gust_dy) + (w+w_gust)*dw_gust_dz;
-


WallyMaier · 2021-06-28T17:50:02Z

SU2_CFD/src/numerics/heat.cpp


-  SoundSpeed_i = sqrt(ProjVelocity_i*ProjVelocity_i + (BetaInc2_i/DensityInc_i)*Area*Area);
-  SoundSpeed_j = sqrt(ProjVelocity_j*ProjVelocity_j + (BetaInc2_j/DensityInc_j)*Area*Area);
+  su2double SoundSpeed_i = sqrt(pow(ProjVelocity_i,2) + (BetaInc2_i/DensityInc_i)*Area2);


Just for my own knowledge, is pow(x,2) superior to x*x?

I think it is more expressive, saves you the trouble of reading two variable names to make sure they are the same.
In terms of performance, the pow function is horrible, however, pow(x,2) always gets optimized to x*x.
Larger integer powers are only optimized in that way if some -ffast-math optimizations are allowed, since e.g. x*x*x is not a strict way to compute pow(x,3) (as per the floating-point standard).
That type of optimization is fairly innocuous in double-precision, the intel compilers actually do them by default.

pcarruscag · 2021-06-29T12:36:49Z

SU2_CFD/include/numerics/CNumerics.hpp

+                                                        const Vec2& var_i, const Vec2& var_j,
+                                                        su2double* projNormal,
+                                                        su2double* projCorrected) {
+    nDim = (nDim > 2)? 3 : 2;


This keeps the compiler from going crazy with loop unrolling thinking that nDim is large, and then creating unrolled code that is never used.

maybe it is time for

template<class nDim> class SU2_CFD { }

:) because this looks really weird, jk

btw I guess this loop unrolling thing is sth you checked with compiler explorer?

I should add an assert though, for 2 or 3.
Yes compiler explorer

Of course this is a valid check, but can you motivate that a bit why you think it is additionally necessary here? I guess if we would do this kind of sanity checks of sizes in more places, we would have quite a bit of them. At first glance that looks like a bit overkill but maybe i just dont understand it.

And might this already be enough for the compiler to not unroll further because it puts a hard constraint on the nDim value

Because the code is only right if called with nDim = 2 or 3, and it is good practice to assert that kind of assumption to facilitate debugging, keep in mind that asserts are no-ops in release builds!
This is also why they cannot be used to inform the compiler that we expect 2 or 3 -> https://gcc.godbolt.org/z/nbPv8KhPo

TobiKattmann

Thanks pedro for this 💐 which is a bit more of just a charge against pointer-to-pointer. This ComputeProjectedGradient func and the simplification of the viscous terms numerics for the heat solver. And the const-ing and cleanup.
I have a few things below, none of which are dealbreakers (mostly questions I guess)

TobiKattmann · 2021-07-15T09:46:50Z

Common/include/containers/container_decorators.hpp

+  const Scalar* operator[] (Index i) const noexcept { return &m_ptr[i*m_cols]; }
+  const Scalar& operator() (Index i, Index j) const noexcept { return m_ptr[i*m_cols + j]; }
+
+  template<class U = T, su2enable_if<!std::is_const<U>::value> = 0>
+  Scalar* operator[] (Index i) noexcept { return &m_ptr[i*m_cols]; }
+
+  template<class U = T, su2enable_if<!std::is_const<U>::value> = 0>
+  Scalar& operator() (Index i, Index j) noexcept { return m_ptr[i*m_cols + j]; }


I guess this detects whether a var on the receiving side is const and then uses the latter versions of [] and () operators?
Edit: Or the matrix view - object has to be created with this property already and then it decides and what to do? I guess the latter

Hmmm no..
If the matrix view object is used in a const context then it returns const refs to su2double. Otherwise it will "try" to return non const ref.
However, if the view was returned by a const container, it cannot do that. Hence the use of sfinae to disable the non const version.

TobiKattmann · 2021-07-15T09:55:19Z

Common/include/containers/container_decorators.hpp

+  template<class U = T, su2enable_if<!std::is_const<U>::value> = 0>
+  Scalar& operator() (Index i, Index j) noexcept { return m_ptr[i*m_cols + j]; }
+
+  friend CMatrixView operator+ (CMatrixView mv, Index incr) { return CMatrixView(mv[incr], mv.m_cols); }


So this allows for the +1 logic which is equivalent to (i,j) logic to allow for similar use as with pointers. This of course might break less code of people and their code out there when pulling this, so 👍 ... but on the other hand might give some false security that the underlying datatype is still the same?

Well the type is not really important in this case, just the interface. And there is too much of this +stuff, so less work for me xD.

TobiKattmann · 2021-07-15T13:33:42Z

Common/include/containers/container_decorators.hpp

+  Matrix operator[] (Index i) noexcept { return Matrix(m_storage[i], m_innerSz); }
+  ConstMatrix operator[] (Index i) const noexcept { return ConstMatrix(m_storage[i], m_innerSz); }


So C3DContainerDecorator owns the data and can return MatrixView objects that provide access to the 2D structure for one entry after the first layer (or the 3rd dimension to stay in this dimension language)?

I think I was influenced by Eigen (inception style) for the "view"

SU2_CFD/include/numerics/CNumerics.hpp

TobiKattmann · 2021-07-15T13:49:22Z

SU2_CFD/include/numerics/CNumerics.hpp

+                                                        const Vec2& var_i, const Vec2& var_j,
+                                                        su2double* projNormal,
+                                                        su2double* projCorrected) {
+    nDim = (nDim > 2)? 3 : 2;


maybe it is time for

template<class nDim> class SU2_CFD { }

:) because this looks really weird, jk

btw I guess this loop unrolling thing is sth you checked with compiler explorer?

TobiKattmann · 2021-07-15T16:17:38Z

SU2_CFD/src/numerics/heat.cpp

-  /*--- Compute vector going from iPoint to jPoint ---*/
+  auto proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ConsVar_Grad_i,
+                                                 ConsVar_Grad_j, correct, &Temp_i, &Temp_j,
+                                                 NormalGrad, CorrectedGrad);


cut-paste-spezialize in action :)

TobiKattmann · 2021-07-15T16:23:49Z

SU2_CFD/src/solvers/CHeatSolver.cpp

+      const auto Temp_i_Grad = nodes->GetGradient(iPoint);
+      const auto Temp_j_Grad = nodes->GetGradient(jPoint);


now auto evaluates to the MatrixView type, right (that returns const values on its own)?

Correct. Before it was a const pointer that could return non const

SU2_CFD/src/solvers/CHeatSolver.cpp

TobiKattmann · 2021-07-15T16:30:57Z

TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref

@@ -1,2 +1,2 @@
 "VARIABLE"      , "AVG_DENSITY[0]", "AVG_ENTHALPY[0]", "AVG_NORMALVEL[0]", "DRAG[0]"       , "EFFICIENCY[0]" , "FORCE_X[0]"    , "FORCE_Y[0]"    , "FORCE_Z[0]"    , "LIFT[0]"       , "MOMENT_X[0]"   , "MOMENT_Y[0]"   , "MOMENT_Z[0]"   , "SIDEFORCE[0]"  , "SURFACE_MACH[0]", "SURFACE_MASSFLOW[0]", "SURFACE_MOM_DISTORTION[0]", "SURFACE_PRESSURE_DROP[0]", "SURFACE_SECONDARY[0]", "SURFACE_SECOND_OVER_UNIFORM[0]", "SURFACE_STATIC_PRESSURE[0]", "SURFACE_STATIC_TEMPERATURE[0]", "SURFACE_TOTAL_PRESSURE[0]", "SURFACE_TOTAL_TEMPERATURE[0]", "SURFACE_UNIFORMITY[0]", "AVG_TEMPERATURE[1]", "MAXIMUM_HEATFLUX[1]", "TOTAL_HEATFLUX[1]", "FINDIFF_STEP"  
-0               , 0.0             , -99999.96982514858, 2.2204999999731917e-08, 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , -0.04999999997368221, -1.1100000001884e-08, -2.069999999187999       , 0.0                     , 2.12000000054946    , 3.7100000016554446            , 330.00000030369847        , -39.999997625272954          , 315.0000011942211        , -40.000008993956726         , -1.400000004814217   , -149.9999996212864, 0.0                , -469.9999976764957, 1e-08           
+0               , 0.0             , -99999.96982514858, 3.330700000025998e-08, 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , 0.0             , -0.04999999997368221, -1.1100000001884e-08, -2.069999999187999       , 0.0                     , 2.12000000054946    , 3.7100000016554446            , 330.00000030369847        , -39.999997625272954          , 315.0000011942211        , -40.000008993956726         , -1.400000004814217   , -149.9999996212864, 0.0                , -469.9999976764957, 1e-08           


🔎 what's going on here? (breakfast 🐞 obviously, but jokes aside) was this a bug? I guess it is 1e-8 land ...

My different handling of division by zero maybe, or the optimizability of 2 or 3

Can you give me place to look for the different handling of division by zero 🔬 (so I dont have to think for myself)

TobiKattmann · 2021-07-15T16:33:03Z

TestCases/serial_regression.py

@@ -345,7 +345,7 @@ def main():
    turb_naca0012_sst.cfg_dir   = "rans/naca0012"
    turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
    turb_naca0012_sst.test_iter = 10
-    turb_naca0012_sst.test_vals = [-11.451011, -12.798258, -5.863895, 1.049989, 0.019163, -1.925028]
+    turb_naca0012_sst.test_vals = [-11.451010, -12.798258, -5.863895, 1.049989, 0.019163, -1.925018]


I have a bit of trouble in finding an explanation from the changes why these few have marginal differences ... because I would then expect more to fail ... but maybe you can shed some light on this if you have an idea

Co-authored-by: TobiKattmann <31306376+TobiKattmann@users.noreply.github.com>

SU2_CFD/include/numerics/CNumerics.hpp

pcarruscag · 2021-07-16T09:41:47Z

SU2_CFD/include/numerics/CNumerics.hpp

+      dist_ij_2 += pow(edgeVec[iDim], 2);
+      proj_vector_ij += edgeVec[iDim] * normal[iDim];
+    }
+    proj_vector_ij /= max(dist_ij_2,EPS);


create a type for "matrix views" (wraps a pointer) to avoid **

4f8fccc

pcarruscag added the changelog:chore label Jun 27, 2021

pr-triage bot added the PR: unreviewed label Jun 27, 2021

pcarruscag commented Jun 27, 2021

View reviewed changes

Common/include/containers/container_decorators.hpp Outdated Show resolved Hide resolved

small fix

6971ab0

pcarruscag commented Jun 27, 2021

View reviewed changes

const correctness and less duplication in CNumerics

31e94a0

pcarruscag commented Jun 27, 2021

View reviewed changes

SU2_CFD/include/numerics/CNumerics.hpp Outdated Show resolved Hide resolved

correct brief

7d1eb8c

pcarruscag commented Jun 27, 2021

View reviewed changes

pcarruscag added 3 commits June 27, 2021 22:22

more cleaning + update regressions

72e119a

Merge branch 'gradient_type' of https://github.com/su2code/SU2 into g…

973b331

…radient_type

not needed member variable

a488554

WallyMaier approved these changes Jun 28, 2021

View reviewed changes

pr-triage bot added PR: reviewed-approved and removed PR: unreviewed labels Jun 28, 2021

Merge branch 'develop' into gradient_type

936fe67

pr-triage bot added PR: unreviewed and removed PR: reviewed-approved labels Jun 29, 2021

pcarruscag commented Jun 29, 2021

View reviewed changes

pcarruscag mentioned this pull request Jul 12, 2021

Fix for axisymmetric terms in NEMO + general NEMO updates #1326

Merged

5 tasks

Merge branch 'develop' into gradient_type

e348065

TobiKattmann approved these changes Jul 15, 2021

View reviewed changes

Apply suggestions from code review

520b978

Co-authored-by: TobiKattmann <31306376+TobiKattmann@users.noreply.github.com>

pcarruscag commented Jul 15, 2021

View reviewed changes

SU2_CFD/include/numerics/CNumerics.hpp Show resolved Hide resolved

Update SU2_CFD/include/numerics/CNumerics.hpp

0ad7a02

pcarruscag commented Jul 16, 2021

View reviewed changes

pcarruscag merged commit 22bb669 into develop Jul 18, 2021

pcarruscag deleted the gradient_type branch July 18, 2021 21:45

pr-triage bot added PR: merged and removed PR: unreviewed labels Jul 18, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Another charge against pointer to pointer #1312

Another charge against pointer to pointer #1312

pcarruscag commented Jun 27, 2021

pcarruscag Jun 27, 2021

pcarruscag Jun 27, 2021

pcarruscag Jun 27, 2021

pcarruscag Jun 27, 2021 •

edited

Loading

pcarruscag Jun 27, 2021

WallyMaier left a comment

WallyMaier Jun 28, 2021

WallyMaier Jun 28, 2021

WallyMaier Jun 28, 2021

WallyMaier Jun 28, 2021

pcarruscag Jun 29, 2021 •

edited

Loading

pcarruscag Jun 29, 2021

TobiKattmann Jul 15, 2021

pcarruscag Jul 15, 2021

TobiKattmann Jul 16, 2021

pcarruscag Jul 16, 2021

TobiKattmann left a comment

TobiKattmann Jul 15, 2021

pcarruscag Jul 15, 2021

TobiKattmann Jul 15, 2021

pcarruscag Jul 15, 2021

TobiKattmann Jul 15, 2021

pcarruscag Jul 15, 2021 •

edited

Loading

TobiKattmann Jul 15, 2021

TobiKattmann Jul 15, 2021

TobiKattmann Jul 15, 2021

pcarruscag Jul 15, 2021

TobiKattmann Jul 15, 2021

pcarruscag Jul 15, 2021

TobiKattmann Jul 16, 2021

TobiKattmann Jul 15, 2021

pcarruscag Jul 16, 2021

		visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), Grad_Reflected);
		visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), CMatrixView<su2double>(Grad_Reflected));

		CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow)+1,Viscosity);
		CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow,1), Viscosity);

		SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(&(Node_Flow->GetGradient_Primitive(iPoint)[1])));
		SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1)));

		Matrix operator[] (Index i) noexcept { return Matrix(m_storage[i], m_innerSz); }
		ConstMatrix operator[] (Index i) const noexcept { return ConstMatrix(m_storage[i], m_innerSz); }

		const auto Temp_i_Grad = nodes->GetGradient(iPoint);
		const auto Temp_j_Grad = nodes->GetGradient(jPoint);

Another charge against pointer to pointer #1312

Another charge against pointer to pointer #1312

Conversation

pcarruscag commented Jun 27, 2021

Proposed Changes

PR Checklist

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pcarruscag Jun 27, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

WallyMaier left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pcarruscag Jun 29, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TobiKattmann left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pcarruscag Jul 15, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pcarruscag Jun 27, 2021 •

edited

Loading

pcarruscag Jun 29, 2021 •

edited

Loading

pcarruscag Jul 15, 2021 •

edited

Loading