Hybrid Parallel AD (Part 2/?) #1284

pcarruscag · 2021-05-09T09:53:32Z

Proposed Changes

1 - Local indices for solution variables (which are overwritten during recording), even in singlezone problems (it was already used in multizone), to make things consistent for the two drivers.

2 - SU2_TYPE::GetDerivative was based on a global counter and a shared vector of input indices, which meant derivatives had to be extracted in the same order they were registered, and made the whole process not thread-safe. This is gone, if an input variable is overwritten during recording its (initial) index needs to be managed explicitly.

3 - Simplify and improve the restart logic because of unsteady FSI issues (not an Hybrid AD topic, but I had enough branches).

Related Work

Continues #1214

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.

TobiKattmann · 2021-05-12T13:45:15Z

Do you plan to add any hybrid parallel AD stuff here? I would prefer if this was a follow-up to #1260 where the Single zone DA input index storing is changed to the multizone strategy with explicit index-containers owned by the solver. You know ... smaller PR, faster merging

I will also add a preliminary unsteady cht adjoint case myself for coverage. Then we could also move the !CrossTerm of ExtractAdjoint_Solution thing in this PR and see if it breaks stuff.

pcarruscag · 2021-05-12T13:55:13Z

To some extent this is a fix for hybrid parallel AD, because the single zone driver used the version of RegisterInput that pushed indices into that global "inputValues" vector, and that would not work with OpenMP.
This PR will not be very big anyway. You can test your testcase here, or in #1260, or in a different PR, and then we test the cross term thing.

…parallel_ad3

…arallel_ad3

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

…into hybrid_parallel_ad3

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

TobiKattmann · 2021-05-18T11:33:25Z

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

+  /*--- Residuals and time_n terms are not needed when evaluating multizone cross terms. ---*/
+  if (CrossTerm) return;
+
+  SetResToZero();


Could we have an unwanted addition to the residual because of moving the SetResToZero below the CrossTerm-return? The residual should compare the (full) adjoint solution so aren't we missing the External (crossterms + dualTime terms) here anyway. I guess just the dualTime Terms are constant so we can leave them out. Or is it correct now as CrossTerms are extracted first. I have to look that up

No, the residual should measure the norm of the residual of the (linear) adjoint equations, otherwise it would measure stagnation instead of convergence.
It just so happens that differences between consecutive fixed point iterations give you that residual (because the fixed point is the same as a Richardson iteration for the right-preconditioned adjoint equations).
The single zone uses differences of the total terms, the multizone uses differences of the product between total terms and the Jacobian of the iteration (it is easier to do it that way), this is equivalent since the right hand side is constant during inner iterations.

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

pcarruscag · 2021-05-25T10:23:58Z

SU2_CFD/src/drivers/CDriver.cpp

-  switch (config->GetKind_Solver()) {
-    case TEMPLATE_SOLVER: template_solver = true; break;
-    case EULER : case INC_EULER: euler = true; break;
-    case NEMO_EULER : NEMO_euler = true; break;
-    case NAVIER_STOKES: case INC_NAVIER_STOKES: ns = true; heat = config->GetWeakly_Coupled_Heat(); break;
-    case NEMO_NAVIER_STOKES: NEMO_ns = true; break;
-    case RANS : case INC_RANS: ns = true; turbulent = true; heat = config->GetWeakly_Coupled_Heat(); break;
-    case FEM_EULER : fem_euler = true; break;
-    case FEM_NAVIER_STOKES: fem_ns = true; break;
-    case FEM_RANS : fem_ns = true; break;
-    case FEM_LES : fem_ns = true; break;
-    case HEAT_EQUATION: heat = true; break;
-    case FEM_ELASTICITY: fem = true; break;
-    case ADJ_EULER : euler = true; adj_euler = true; break;
-    case ADJ_NAVIER_STOKES : ns = true; turbulent = (config->GetKind_Turb_Model() != NONE); adj_ns = true; break;
-    case ADJ_RANS : ns = true; turbulent = true; adj_ns = true; adj_turb = (!config->GetFrozen_Visc_Cont()); break;
-    case DISC_ADJ_EULER: case DISC_ADJ_INC_EULER: euler = true; disc_adj = true; break;
-    case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_INC_NAVIER_STOKES: ns = true; disc_adj = true; heat = config->GetWeakly_Coupled_Heat(); break;
-    case DISC_ADJ_RANS: case DISC_ADJ_INC_RANS: ns = true; turbulent = true; disc_adj = true; disc_adj_turb = (!config->GetFrozen_Visc_Disc()); heat = config->GetWeakly_Coupled_Heat(); break;
-    case DISC_ADJ_FEM_EULER: fem_euler = true; disc_adj = true; break;
-    case DISC_ADJ_FEM_NS: fem_ns = true; disc_adj = true; break;
-    case DISC_ADJ_FEM_RANS: fem_ns = true; turbulent = true; disc_adj = true; disc_adj_turb = (!config->GetFrozen_Visc_Disc()); break;
-    case DISC_ADJ_FEM: fem = true; disc_adj_fem = true; break;
-    case DISC_ADJ_HEAT: heat = true; disc_adj_heat = true; break;


🚀 🎆 ☠️

... this kind of cleanup will also makes the potential change to enum class Kind_Solver less tedious, so kudos from someone in the future and from me in the present :)

TobiKattmann

Thanks for the adjoint consolidation of index storing things and all the other stuff 💐 . Imo this helps quite a bit to understand the DA solver (single & multizone). I have mostly question for my own understanding below. (I also have to come back to some of your previous comments/asnwers which I have to understand ;) )

TobiKattmann · 2021-05-25T09:48:51Z

SU2_CFD/include/variables/CMeshVariable.hpp

+  inline void GetAdjoint_MeshCoord(unsigned long iPoint, su2double *adj_mesh) final {
+    for (unsigned long iDim = 0; iDim < nDim; iDim++) {
+      adj_mesh[iDim] = AD::GetDerivative(AD_InputIndex(iPoint,iDim));
+      AD::ResetInput(Mesh_Coord(iPoint,iDim));


A bit off-topic: I was digging a bit about AD::ResetInput and what it does. The code documentation says

* \brief Reset the variable (set index to zero). * \param[in] data - the variable to be unregistered from the tape.

and the code is the following

FORCEINLINE void ResetInput(su2double &data) {data.getGradientData() = su2double::GradientData();}

now data.getGradientData() returns some index on the tape and su2double::GradientData() is the default constructor (which I guess returns 0) of the datatype codi::RealForward::GradientData which matches the left sides datatype. Alright but why are we doing that? I am struggling a bit with that

We also AD::ResetInput in CDiscAdjSolver::SetRecording but only for Solution_n/1 not for Solution itself, i.e. in steady simulation AD::ResetInput is not called during the flow adjoint.

AD::ResetInput(direct_solver->GetNodes()->GetSolution_time_n(iPoint)[iVar]);

Thanks for any insight already

I think the trend is that if something is an input the index is not cleared by the dummy recording (since nothing is written to an input) and so you have to clear it explicitly.
Now, I guess this is only strictly required if sometimes (i.e. some recordings) you register that variable and sometimes you don't. If you always register it then the index is always valid, but if you leave it with "dangling" indices then you have a problem.
From my comment above it also seems that the location where this reset is performed matters... (but I might also have changed two things at the same time while testing).

SU2_CFD/include/variables/CVariable.hpp

TobiKattmann · 2021-05-25T10:01:35Z

SU2_CFD/src/drivers/CDriver.cpp

+    iteration_container[iZone]   = new CIteration*    [nInst[iZone]] ();
+    solver_container[iZone]      = new CSolver***     [nInst[iZone]] ();
+    integration_container[iZone] = new CIntegration** [nInst[iZone]] ();
+    numerics_container[iZone]    = new CNumerics****  [nInst[iZone]] ();


C++ understanding question: With () you call the default constructor, right? What is the benefit of that? Does it automatically initialize with nullptr as that set-ting was deleted below?

TobiKattmann · 2021-05-25T10:07:17Z

SU2_CFD/src/drivers/CDriver.cpp

+    for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+      auto sol = solver[MESH_0][iSol];
+      if (sol && !sol->GetAdjoint()) {
+        /*--- Note that the mesh solver always loads the most recent file (and not -2). ---*/
+        SU2_OMP_PARALLEL_(if(sol->GetHasHybridParallel()))
+        sol->LoadRestart(geometry, solver, config, val_iter + (iSol==MESH_SOL && dt_step_2nd), update_geo);
+        END_SU2_OMP_PARALLEL


Awesome cleanup of the restart logic🧹

Perhaps the only advantage of the solver container being an array.

TobiKattmann · 2021-05-25T10:17:47Z

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

-      direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution);
-    }
+    su2double Solution[MAXNVAR] = {0.0};
+    direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution);


Nice to have the input-index & output-index storing-strategies consolidated between single and multizone adjoint 👍

TestCases/hybrid_regression.py

TobiKattmann

Man this is good stuff 💐 Also without any problems (and modifications) with Register/Extract_Variable. I hope there is at least some Testcases that had e.g. AOA sensitivity in their screen output? (Did you maybe even check that)

And if I understood it correctly that this change revealed, that in two occasions

  FORCEINLINE passivedouble GetDerivative(const su2double& data) {
    return AD::getGlobalTape().getGradient(AD::inputValues[AD::adjointVectorPosition++]);
  }

fooled people with its fake input value.

Thanks again! As this is now at +500 -1000 and quite some changes some to the DA machine room I would personally prefer to have this merged without any major other updates. Also kudos for updating the PR explanation on top 👍

Common/include/basic_types/datatype_structure.hpp

Common/include/basic_types/ad_structure.hpp

TobiKattmann · 2021-05-26T07:33:48Z

SU2_CFD/include/variables/CVariable.hpp

        else AD::RegisterOutput(variable(iPoint,iVar));

-        AD::SetIndex(ad_index(iPoint,iVar), variable(iPoint,iVar));
+        if (ad_index) AD::SetIndex((*ad_index)(iPoint,iVar), variable(iPoint,iVar));


👍 working with additional index containers and without

TobiKattmann · 2021-05-26T07:36:51Z

SU2_CFD/include/variables/CVariable.hpp

-  su2matrix<int> AD_Time_n_InputIndex;  /*!< \brief Indices of Solution variables in the adjoint vector. */
-  su2matrix<int> AD_Time_n1_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */


now we are making explicit that solution_time_n/1 should not be changed in one dual-time step. I like that. The regression test didn't fail so I see this as good thing

TobiKattmann · 2021-05-26T07:50:00Z

UnitTests/Common/simple_ad_test.cpp

@@ -56,5 +56,5 @@ TEST_CASE("Simple AD Test", "[AD tests]") {
  AD::ComputeAdjoint();

  CHECK(SU2_TYPE::GetValue(y) == Approx(64));
-  CHECK(SU2_TYPE::GetDerivative(y) == Approx(48));
+  CHECK(SU2_TYPE::GetDerivative(x) == Approx(48));


For the record: here GetDerivative(y) was not get-ting the derivative of y but out of the hidden datastructure inputValues which (is/)was filled with RegisterInput -> that's why it worked before

TobiKattmann · 2021-05-26T08:07:58Z

SU2_DOT/src/SU2_DOT.cpp

-
-    for (iDV_Value = 0; iDV_Value < nDV_Value; iDV_Value++){
+    for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++){

-      /*--- Initilization with su2double resets the index ---*/
+      config->SetDV_Value(iDV, iDV_Value, 0.0);

-      DV_Value = 0.0;
-
-      AD::RegisterInput(DV_Value);
-
-      config->SetDV_Value(iDV, iDV_Value, DV_Value);
+      AD::RegisterInput(config->GetDV_Value(iDV, iDV_Value));


I wondered what made this break before and I can only think that before the local DV_VALUE was registered and below we now(!) extracted what is in config->GetDV_Value(...). Before this PR it accessed the hidden inputValues so everything was good (no matter the input value to SU2_TYPE::GetDerivative(...)).

DV_Value = config->GetDV_Value(iDV, iDV_Value); my_Gradient = SU2_TYPE::GetDerivative(DV_Value);

Now you changed it to only work with the DV_Value that lives in the config and that is consistent with this change up here.

The problem before this very commit was probably index management optimization because config->SetDV_Value(iDV, iDV_Value, DV_Value); before was just an assignment. But that is sth I would ask @jblueh in the Dev-meeting as well

We were using always the same variable (declared at the top of the function).
On each loop iteration the =0 part would reset the previously registered index, then on extraction we were also writing over the same variable.

TobiKattmann · 2021-05-26T08:11:48Z

TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref

+0	-3.461460672772851e-03
+1	-1.841786315630031e-03


I know the changes are marginal and probably ok, but if this is not a fully converged adjoint and DOT step (not just 20 SU2_CFD_AD steps and then SU2_DOT_AD) I would prefer to have that checked on the 'converged' gradient. Maybe you did this already

I changed the settings of the case to make it converge better/faster.

Common/include/basic_types/ad_structure.hpp

…into hybrid_parallel_ad3

…on and MAXNVAR fix for flow.

pcarruscag · 2021-05-27T14:12:38Z

Back to WIP, sorry @TobiKattmann.
It may be worth git cherry-picking 19d3329 to see if it fixes the issues @oleburghardt found.

Common/include/code_config.hpp

TobiKattmann · 2021-05-29T15:05:04Z

Well thanks for reverting the latest two commits with RealReverseIndex and the CoDiPack update and merging this. I guess that work will go on in another PR

local indices to allow parallel registration

f7f1b96

pcarruscag added the changelog:feature label May 9, 2021

pr-triage bot added the PR: unreviewed label May 9, 2021

pcarruscag changed the title ~~Hybrid Parallel AD (Part 2/?)~~ [WIP] Hybrid Parallel AD (Part 2/?) May 9, 2021

pr-triage bot removed the PR: unreviewed label May 9, 2021

reset inputs as they are extracted

20cbe7a

This was referenced May 9, 2021

Discrete adjoint for dynamic FSI using multizone driver #1260

Merged

Incorrect use of short-circuit 'or'? #1285

Closed

pcarruscag added 2 commits May 10, 2021 11:07

fix #1285

e8522c6

address #1273

7f9f647

This was linked to issues May 10, 2021

Incorrect use of short-circuit 'or'? #1285

Closed

OpenMP version of the adjoint solver does not compile with the Intel compiler #1273

Closed

TobiKattmann mentioned this pull request May 12, 2021

Add unsteady cht adjoint testcase #1288

Merged

5 tasks

pcarruscag and others added 7 commits May 13, 2021 22:19

try to fix restarts

2144acd

Merge remote-tracking branch 'upstream/develop' into hybrid_parallel_ad3

7efd040

Merge remote-tracking branch 'upstream/feature_dynFSIAD' into hybrid_…

c4e3eb4

…parallel_ad3

update regressions

0d55909

Merge remote-tracking branch 'upstream/add_unstchtcase' into hybrid_p…

0f06d7a

…arallel_ad3

Merge branch 'feature_dynFSIAD' into hybrid_parallel_ad3

91503ea

Move \!Crossterm bool in front of loop for unsteady adjoint extraction.

8c6233a

TobiKattmann reviewed May 18, 2021

View reviewed changes

SU2_CFD/src/solvers/CDiscAdjSolver.cpp Outdated Show resolved Hide resolved

pcarruscag added 2 commits May 18, 2021 12:07

dont extract some terms during cross term evaluation

7422c00

Merge branch 'hybrid_parallel_ad3' of /~https://github.com/su2code/SU2 …

d598915

…into hybrid_parallel_ad3

TobiKattmann reviewed May 18, 2021

View reviewed changes

pcarruscag linked an issue May 18, 2021 that may be closed by this pull request

Unsteady FSI with SU2 Compressible Solver and SU2 FEM Solver Diverges when Restarted (v7.1.1) #1287

Closed

apply BGS relaxation also to velocities

d949cbe

TobiKattmann mentioned this pull request May 19, 2021

CFVMOutput & Streamwise+spanwise periodic #1290

Merged

5 tasks

cleanup overzealous SU2_OMP_MASTER

b7f5b05

pcarruscag commented May 25, 2021

View reviewed changes

TobiKattmann approved these changes May 25, 2021

View reviewed changes

pcarruscag added 2 commits May 25, 2021 16:16

local indices only, no more hidden global "adjointPosition"

d5715c4

fixes

567064a

TobiKattmann approved these changes May 26, 2021

View reviewed changes

proper-er way to reset input

ef4f233

pcarruscag changed the title ~~[WIP] Hybrid Parallel AD (Part 2/?)~~ Hybrid Parallel AD (Part 2/?) May 26, 2021

pr-triage bot added the PR: unreviewed label May 26, 2021

small things

9df5d02

TobiKattmann reviewed May 26, 2021

View reviewed changes

Common/include/basic_types/ad_structure.hpp Show resolved Hide resolved

TobiKattmann and others added 4 commits May 27, 2021 11:53

Merge branch 'hybrid_parallel_ad3' of /~https://github.com/su2code/SU2 …

9a6346f

…into hybrid_parallel_ad3

Remove errors from #1281. No auto for su2double with rhs math operati…

19d3329

…on and MAXNVAR fix for flow.

CoDiPack update.

e8cee0b

Testing RealReverseIndex.

be27e1d

pcarruscag changed the title ~~Hybrid Parallel AD (Part 2/?)~~ [WIP] Hybrid Parallel AD (Part 2/?) May 27, 2021

pr-triage bot removed the PR: unreviewed label May 27, 2021

TobiKattmann reviewed May 27, 2021

View reviewed changes

Common/include/code_config.hpp Outdated Show resolved Hide resolved

pcarruscag changed the title ~~[WIP] Hybrid Parallel AD (Part 2/?)~~ Hybrid Parallel AD (Part 2/?) May 29, 2021

pr-triage bot added the PR: unreviewed label May 29, 2021

Merge branch 'develop' into hybrid_parallel_ad3

4593083

pcarruscag merged commit dc2c9e1 into develop May 29, 2021

pcarruscag deleted the hybrid_parallel_ad3 branch May 29, 2021 14:27

pr-triage bot added PR: merged and removed PR: unreviewed labels May 29, 2021

jblueh mentioned this pull request Jun 1, 2021

Hybrid Parallel AD (Part 3/?) #1294

Merged

5 tasks

jblueh mentioned this pull request May 22, 2023

Hybrid Parallel AD Performance Improvements #2039

Merged

5 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Hybrid Parallel AD (Part 2/?) #1284

Hybrid Parallel AD (Part 2/?) #1284

pcarruscag commented May 9, 2021 •

edited

Loading

TobiKattmann commented May 12, 2021

pcarruscag commented May 12, 2021 •

edited

Loading

TobiKattmann May 18, 2021

pcarruscag May 18, 2021

pcarruscag May 25, 2021

TobiKattmann May 25, 2021 •

edited

Loading

TobiKattmann left a comment

TobiKattmann May 25, 2021

pcarruscag May 25, 2021

TobiKattmann May 25, 2021

pcarruscag May 25, 2021

TobiKattmann May 25, 2021

pcarruscag May 25, 2021

TobiKattmann May 25, 2021

TobiKattmann left a comment

TobiKattmann May 26, 2021

TobiKattmann May 26, 2021

TobiKattmann May 26, 2021 •

edited

Loading

TobiKattmann May 26, 2021

pcarruscag May 26, 2021

TobiKattmann May 26, 2021

pcarruscag May 26, 2021

pcarruscag commented May 27, 2021

TobiKattmann commented May 29, 2021

		su2matrix<int> AD_Time_n_InputIndex; /!< \brief Indices of Solution variables in the adjoint vector. /
		su2matrix<int> AD_Time_n1_InputIndex; /!< \brief Indices of Solution variables in the adjoint vector. /

Hybrid Parallel AD (Part 2/?) #1284

Hybrid Parallel AD (Part 2/?) #1284

Conversation

pcarruscag commented May 9, 2021 • edited Loading

Proposed Changes

Related Work

PR Checklist

TobiKattmann commented May 12, 2021

pcarruscag commented May 12, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TobiKattmann May 25, 2021 • edited Loading

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

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

TobiKattmann May 26, 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

pcarruscag commented May 27, 2021

TobiKattmann commented May 29, 2021

pcarruscag commented May 9, 2021 •

edited

Loading

pcarruscag commented May 12, 2021 •

edited

Loading

TobiKattmann May 25, 2021 •

edited

Loading

TobiKattmann May 26, 2021 •

edited

Loading