Linear solver changes to support hybrid parallel AD - #1228
Conversation
| /*--- Build preconditioner for the transposed Jacobian ---*/ | ||
|
|
||
| if (RequiresTranspose) Jacobian.TransposeInPlace(); | ||
|
|
There was a problem hiding this comment.
The transposition happens here.
Things that could go wrong with this approach:
- Solving the same system multiple times, without clearing the matrix or reverting the transposition in between. It does not happen atm.
I think it is worth the risk, since all preconditioners can be used now, and re-transposing the matrix if necessary is not very expensive (probably the same as applying the ILU once).
There was a problem hiding this comment.
Not that I used a preconditioner besides ILU for DA cases yet, but that is a nice thing 👍
TobiKattmann
left a comment
There was a problem hiding this comment.
Well, the code looks nice as always 💐 but I cannot really say much about the changes in general because I am not exactly familiar with that code :(
| LinearToleranceType tol_type = LinearToleranceType::RELATIVE; /*!< \brief How the linear solvers interpret the tolerance. */ | ||
| bool xIsZero = false; /*!< \brief If true assume the initial solution is always 0. */ | ||
| mutable LinearToleranceType tol_type = LinearToleranceType::ABSOLUTE; /*!< \brief How the linear solvers interpret the tolerance. */ | ||
| mutable bool xIsZero = false; /*!< \brief If true assume the initial solution is always 0. */ |
There was a problem hiding this comment.
Where there is const, there must be mutable somewhere
There was a problem hiding this comment.
Restarted FGMRES (which I put in a function) needs to force those values. But I'll do it some other way, mutable in this case is misleading as relevant state can be different after the solver is called.
I try to reserve mutable for variables that do not change the behavior of the object.
| /*--- To keep the behavior of SU2_DOT, but not strictly required since jacobian is symmetric(?). ---*/ | ||
| const bool RequiresTranspose = !mesh_deform || (config->GetKind_SU2() == SU2_DOT); |
There was a problem hiding this comment.
I guess SU2_DOT does a mesh deformation using the linear solvers so that is what the question mark refers to in the comment? Or am I missing sth here
There was a problem hiding this comment.
Well the Jacobian is not symmetric because of how we do the boundary conditions but conjugate gradient still works so... 🤷
| /*--- Build preconditioner for the transposed Jacobian ---*/ | ||
|
|
||
| if (RequiresTranspose) Jacobian.TransposeInPlace(); | ||
|
|
There was a problem hiding this comment.
Not that I used a preconditioner besides ILU for DA cases yet, but that is a nice thing 👍
| da_sp_pinArray_cht_2d_dp_hf.cfg_file = "DA_configMaster.cfg" | ||
| da_sp_pinArray_cht_2d_dp_hf.test_iter = 100 | ||
| da_sp_pinArray_cht_2d_dp_hf.test_vals = [-4.793283, -4.065832, -4.137121] #last 4 lines | ||
| da_sp_pinArray_cht_2d_dp_hf.test_vals = [-4.800583, -4.065533, -4.137185] #last 4 lines |
There was a problem hiding this comment.
Well the deviations are really tiny but both inc cases have INC_ENERGY_EQUATION= YES but discadj_cht is not here as the third case that has that... Not that I am worried but if you have an explanation why exactly those cases differ, I'd be interested?
There was a problem hiding this comment.
The common denominator seemed to be cases that had a significant number of linear iterations.
I ran them to convergence and there is 0 difference in final results.
I will still test these changes in my optimization cases before merging.
| enum ENUM_LINEAR_SOLVER { | ||
| STEEPEST_DESCENT = 1, /*!< \brief Steepest descent method for point inversion algoritm (Free-Form). */ | ||
| NEWTON = 2, /*!< \brief Newton method for point inversion algorithm (Free-Form). */ | ||
| QUASI_NEWTON = 3, /*!< \brief Quasi Newton method for point inversion algorithm (Free-Form). */ | ||
| CONJUGATE_GRADIENT = 4, /*!< \brief Preconditionated conjugate gradient method for grid deformation. */ | ||
| FGMRES = 5, /*!< \brief Flexible Generalized Minimal Residual method. */ | ||
| BCGSTAB = 6, /*!< \brief BCGSTAB - Biconjugate Gradient Stabilized Method (main solver). */ | ||
| RESTARTED_FGMRES = 7, /*!< \brief Flexible Generalized Minimal Residual method with restart. */ | ||
| SMOOTHER = 8, /*!< \brief Iterative smoother. */ | ||
| PASTIX_LDLT = 9, /*!< \brief PaStiX LDLT (complete) factorization. */ | ||
| PASTIX_LU = 10, /*!< \brief PaStiX LU (complete) factorization. */ | ||
| CONJUGATE_GRADIENT, /*!< \brief Preconditionated conjugate gradient method for grid deformation. */ | ||
| FGMRES, /*!< \brief Flexible Generalized Minimal Residual method. */ | ||
| BCGSTAB, /*!< \brief BCGSTAB - Biconjugate Gradient Stabilized Method (main solver). */ | ||
| RESTARTED_FGMRES, /*!< \brief Flexible Generalized Minimal Residual method with restart. */ | ||
| SMOOTHER, /*!< \brief Iterative smoother. */ | ||
| PASTIX_LDLT, /*!< \brief PaStiX LDLT (complete) factorization. */ | ||
| PASTIX_LU, /*!< \brief PaStiX LU (complete) factorization. */ |
| /*--- Reduce across all mpi ranks, only master thread communicates. ---*/ | ||
| SU2_OMP_BARRIER | ||
| SU2_OMP_MASTER { | ||
| sum = dotRes; | ||
| const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE; | ||
| SelectMPIWrapper<ScalarType>::W::Allreduce(&sum, &dotRes, 1, mpi_type, MPI_SUM, SU2_MPI::GetComm()); |
| sp_pinArray_cht_2d_dp_hf.test_iter = 100 | ||
| sp_pinArray_cht_2d_dp_hf.test_vals = [0.247022, -0.812199, -0.974877, -0.753315, 208.023676, 349.950000] #last 7 lines | ||
| sp_pinArray_cht_2d_dp_hf.test_vals = [0.247026, -0.811632, -0.982066, -0.753312, 208.023676, 350.180000] #last 7 lines | ||
| sp_pinArray_cht_2d_dp_hf.su2_exec = "mpirun -n 2 SU2_CFD" |
There was a problem hiding this comment.
I will redo the gradient validation on my side... and I like that commit message 👍
Proposed Changes
Instead of computing transposed matrix vector products, which are harder to parallelize with OpenMP, and having special handling for transposed preconditioners, the matrix is transposed in place, then any preconditioner can be applied.
Some tweaks to support better Krylov discrete adjoint.
Related Work
#1214
PR Checklist