Skip to content

Conversation

@pcarruscag
Copy link
Member

@pcarruscag pcarruscag commented Mar 20, 2020

Proposed Changes

Improvements and general cleanup / restructuring of CInterpolator and co. consisting of:

New features:

  • Control of the RBF interpolation matrix sparsity via option RADIAL_BASIS_FUNCTION_PRUNE_TOLERANCE to drop relatively small coefficients.
  • Generalize nearest neighbor interpolation to KNN + inverse distance weighting, "k" determined by option NUM_NEAREST_NEIGHBORS (because "why not").

Enhancements:

  • The setup of CRadialBasisFunction was optimized (more use of BLAS) and it is now parallel (MPI and threads) over multiple interface patches.
  • The setup of CIsoparametric was optimized, now it only evaluates isoparameters for a few nearby candidate donor elements (instead of all possible ones).
  • The setup of CMirror was optimized both for memory (by communicating only information as needed) and speed (by reducing algorithmic complexity from quadratic to linearithmic) this item is in Discrete adjoint for deforming meshes #833.
  • The setup of CIsoparametric, CNearestNeighbor, and CMirror are now hybrid parallel.

Maintenance:

  • Restructure (split) and cleanup the files.
  • Heavy refactoring of RBF interpolation.
  • Heavy refactoring of isoparametric interpolation, the method was incorrect for quadrilaterals (and the code too complex).
  • Heavy refactoring of mirror (conservative) interpolation, the algorithm is perhaps more complex but the code is cleaner and 10x faster (this is in Discrete adjoint for deforming meshes #833).
  • Medium cleanup of CInterpolator (const correctness, unnecessary mpi ifdefs, etc.).
  • Light cleanup of CDriver/CTransfer (deleting copy-pasted code)

Other enhancements:

  • Initialize the linear system solutions in CMeshSolver to speedup 3D FSI simulations, but especially to reduce overhead in FSI discrete adjoint cases.
  • Too improve the quality of deformed grids in FSI cases with large "out-of-plane" displacements, a thicker layer of stiff cells can be imposed around the deforming boundaries. The thickness of the layer is defined by option DEFORM_STIFF_LAYER_SIZE = [0] when DEFORM_STIFFNESS_TYPE = WALL_DISTANCE. This allows refined regions around a wing tip to be dragged instead of sheared (which happens with the simple 1/x variation).

Other fixes:

  • CFEASolver displacement boundary conditions (the weird bug detailed below).
  • CFEASolver load integration for FSI simplified to fix mpi-related bug.
  • Option CFL_REDUCTION_TURB having no effect.
  • Restarts with moving wall cases work by default.

I know, there's way too much stuff in this PR, but it was all required to make 3D FSI more reliable, and hey, it works!
3D_FSI

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).
  • 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.

AccessorImpl& operator= (AccessorImpl&& other) noexcept \
{ \
if(m_data!=nullptr) free(m_data); \
MemoryAllocation::aligned_free<Scalar_t>(m_data); \
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also fixes a potential memory leak for mac and windows.

Copy link
Member Author

@pcarruscag pcarruscag left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This PR turned into a bit of a bug fixing campaign, my solution for one of them is less straightforward than I would like. It would be great if someone who is very familiar with how markers are partitioned could have a look at this.

The problem

The displacement boundary conditions of CFEASolver (also used for mesh deformation in FSI problems) would sometimes break when running in parallel.
For example this would happen:
Selection_105
files.zip
For 3D FSI the mesh deformation is almost guaranteed to break due to every boundary being of this type, and we would get this type of artifact:
Selection_100

The cause (?)

For displacement boundary conditions we need to impose the solution of the linear system at some nodes, even on partitions where these nodes are halos otherwise matrix-vector products will be incorrect on that partition.
It seems that some partitions do not seem to receive enough of some markers to cover those halo points.

(A) solution

To get correct matrix vector products it is sufficient to communicate the linear system solution before solving the systems. This gives the halos the correct values, and since the matrix vector product (and all linear preconditioners) ignore the "halo rows" of the matrix, those values stay correct.

However, to maintain the symmetry of the stiffness matrix, we eliminate both the row and column associated with these nodes, while the row only needs to be eliminated on the partition for which the node is a domain point, the column needs to be eliminated on all partitions (thus updating the right hand side on any row with a non-zero block on the target column), this process was failing due to the "incomplete marker" problem.
Again, this is only a problem if we want to maintain strict symmetry which is required to use symmetric factorizations, but seemingly irrelevant for all our normal algorithms.
My current solution is to detect these very few problem nodes (10 or so on a case with 50k+ boundary nodes) and do extra row-column eliminations just before solving the linear system...
Ideally we would just have "enough marker" for these nodes to be covered during the normal boundary condition handling, but at the moment I have no idea how to go about fixing the root cause of the problem.

This does not seem to affect the FVM solvers, at least the results do not change with number of ranks...


}

void CFEASolver::Set_VertexEliminationSchedule(CGeometry *geometry, const vector<unsigned short>& markers) {
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the method to detect problem nodes.

Comment on lines +2735 to +2739
Jacobian.InitiateComms(LinSysSol, geometry, config, SOLUTION_MATRIX);
Jacobian.CompleteComms(LinSysSol, geometry, config, SOLUTION_MATRIX);

SU2_OMP(section)
for (auto iPoint = nPointDomain; iPoint < nPoint; iPoint++)
LinSysSol.SetBlock_Zero(iPoint);
for (auto iPoint : ExtraVerticesToEliminate) {
Jacobian.EnforceSolutionAtNode(iPoint, LinSysSol.GetBlock(iPoint), LinSysRes);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And this is what we need to do before solving the system.

@pcarruscag
Copy link
Member Author

I'll merge this in with #833...

@pcarruscag pcarruscag closed this Apr 7, 2020
@pcarruscag pcarruscag deleted the multiphysics_interpolation_improvements branch April 7, 2020 22:24
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants