Edge vector for boundary numerics

[maxaehle]

I wonder whether in the following code the arguments of SetCoord should be swapped. The code illustrates how the visc_numerics is initialized in CEulerSolver::BC_Far_Field:

iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();

V_domain = nodes->GetPrimitive(iPoint);
V_infty = ...

visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint),
                        geometry->nodes->GetCoord(Point_Normal));
visc_numerics->SetPrimitive(V_domain, V_infty);

So iPoint lies on the boundary. Point_Normal is its normal neighbour in the interior of the domain, and not some "dummy cell" outside. But the primitive variables at Point_Normal are set to the outside values.

As far as I can tell, the CNumerics::Coord_j member is only used in the child classes of CAvgGrad_Base, in order to compute

Edge_Vector[MAXNDIM] = {0.0},           /*!< \brief Vector from point i to point j. */

by

Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim];

which is then used to call CorrectGradient. It has the wrong sign if the above code is incorrect. (Apart from that, Coord_j is used to compute the distance in the NEMO viscous residuals.)

Note that similar calls to SetCoord can be found at many places, e.g. in CEulerSolver::BC_Riemann, CEulerSolver::BC_Giles, their siblings in CIncEulerSolver. and CTurbSASolver::BC_Inlet_Turbo.

[pcarruscag]

The gradient correction that involves the point coordinates is not active for the VISC_BOUND term of the numerics. And we flip the normal before passing it to the numerics... Can you check if the dot product of the normal with the edge vector is positive? Either the code has a silent bug, or it is correct but could use a comment about this :)

[maxaehle]

I tried it for the NACA0012 897x257 mesh with compressible RANS (in CAvgGrad_Flow::ComputeResidual).

The scalar product is negative at the far-field boundary. Coord_i is on the boundary, Coord_j is in the domain. Normal points outwards in the viscous numerics (but I don't know if this depends on the orientation of the boundary in the mesh file).

The scalar product is positive for the viscous fluxes between cells.

[pcarruscag]

It's a silent bug then. The direction of the vector does not matter for boundary terms but it should point vaguely in the direction of the normal.

[maxaehle]

Also I think that the _i and _j members of CNumerics should be used to store data that "belongs to one particular cell", respectively. But we cannot leave the second argument of SetCoord just empty (nullptr) because it is needed for dist_ij_2. We can add a comment saying that there is a bug without effect. We could also "extrapolate" coordinates like this:

Coord_i = (Coord of boundary point)
Coord_j = 2 * (Coord of boundary point) - (Coord of normal neighbour in the interior)

or we let the distance between the two cells be a member in CNumerics, so the BC residual functions can set only the distance (computed there with the normal neighbour) but not the coordinates.

[pcarruscag]

I would be fine with either swapping (if it is commented) or reflecting. I would prefer not to have extra checks for boundary logic in a class that operates 99% of the time on internal edges.