Time-accurate discrete adjoint gradients for unsteady problems with grid movement by cvencro · Pull Request #766 · su2code/SU2

cvencro · 2019-08-22T06:52:46Z

Proposed Changes

This PR, created together with @TobiKattmann , makes the necessary changes to obtain time-accurate discrete adjoint gradients for unsteady problems with grid movement. Main changes are for the handling of grid velocities of previous time steps.

Regression test case to follow.

Related Work

depends on #767 unsteady incompressible + ALE. Then a tested discrete adjoint for unsteady incompressible flow can be merged with this.
Edit @TobiKattmann: PR767 is merged so I remove the dependent label

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.

talbring · 2019-08-22T11:37:34Z

    GridVel  = new su2double[nDim];
+    GridVel_Old  = new su2double[nDim];
+    GridVel_n  = new su2double[nDim];
+    GridVel_n1  = new su2double[nDim];


Do we always need this arrays or just in case of the discrete adjoint ?

Wrapped the allocation it in a Discrete Adjoint conditional.

To give some background why: In the reverse sweep (say, with dual time stepping 2nd order) always the n-2 timestep is loaded from disk into memory. In that process, Coordinates and GridVels are stored. So compared to the corresponding primal timestep the adjoint had data from 2 timesteps before (...thinking in primal time) if these additional containers and corresponding data pushing wouldn't be there.
I thought it is more convenient than recomputing as all the data is already there in the restart files (and it was faster to code). I might also be beneficial for cases where GridVel is not a fixed function like in rigid motion but depends on the problem itself -> then it can be simply read from the restart file.

talbring · 2019-08-22T11:43:33Z

      Restart_Flow = false;

-      if (GetGrid_Movement()) {
-        SU2_MPI::Error("Dynamic mesh movement currently not supported for the discrete adjoint solver.", CURRENT_FUNCTION);


I guess it only works for rigid movement at the moment then, am I correct ?

Correct. Now checks Kind_GridMovement != RIGID_MOTION and throws an error.

Rigid motion cases have been used for the verification at the moment but GCL is included to handle changing volumes of a mesh in the primal solvers. Therefore the changes in this request should capture both rigid motion and deforming mesh problems for discrete adjoints as the grid velocities are now captured correctly. Please advise if I'm missing to account for something else here? If anyone can point me to a good single zone mesh deformation test case that already exists for the primal solver, I can investigate the discrete adjoint explicitly for a deforming mesh problem as well.

As of now, deformable meshes are mostly used in FSI (where the coupling terms should be accounted for) but outside of that, there are not many methods to run deformable, single-zone test cases available. I guess something could be done using SURFACE_MOVEMENT = DEFORMING and a pitching/plunging airfoil, but I don't believe the adjoint will work out of the box, as the Grid Velocities are recomputed several times in one round of deformation.

However, with the capabilities in #760, I believe it should be relatively easy to set up a case with time-dependent boundary displacements.

pcarruscag · 2019-08-22T13:10:47Z

 }

+inline void CPoint::SetGridVel_Old(su2double *val_gridvel) { 
+	for (unsigned short iDim = 0; iDim < nDim; iDim++)


Please change these tabs to spaces.

Yep, went over all 'Files changed' and fixed all indentation mismatches I saw.
Btw: In e.g. dual_grid_structure.hpp, most indentations are done using tabs and there are a lot of trailing whitespaces. Is it useful to change it in the whole file? Would be work of a few minutes but the latest git blame would be kindof misleading. I also noticed (which is probably) @WallyMaier in a branch cleanup just doing that cosmetic stuff. This is only a style(-guide) question.

Since we are currently restructuring the files anyway, we should do that clean up when creating/moving new files. This way we avoid duplicate conflicts.

Yep, please wait on the larger white space cleanup until we divide dual_grid_structure.hpp into files for each subclass.

Hi all, I did start to create a clean up branch, referencing the the style guide. I can continue with it (and skip the dual grid). But I think with new version coming in, we should clean up files that are going to be restructured.

economon

LGTM. Thanks @cvencro @TobiKattmann for getting this working (and finding/fixing a couple of bugs in the process).

economon · 2019-08-22T17:58:18Z

  unsigned short iMesh;
  int Direct_Iter;
  bool heat = config[val_iZone]->GetWeakly_Coupled_Heat();
+  bool grid_movement_bool = config[val_iZone]->GetGrid_Movement();


Can you please switch 'grid_movement_bool' to 'grid_movement' to keep with the convention throughout the rest of the code? Eventually, we could use something more expressive for the name like 'gridIsMoving' (find and replace everywhere), but for now, consistency is enough

Yep, but I would actually prefer bool gridIsMoving as grid_movement is used as a larger container for CVolumetricMovement ***grid_movement.

rsanfer

Hi @cvencro, @TobiKattmann ,
thank you so much for this PR, which addresses a certainly interesting capability to compute adjoints on moving meshes. However, I have a question and would like to open a discussion about the generalization of the feature.

I have seen that this PR allocates a lot of memory for GridVel_Old, GridVel_n, and GridVel_n1, and additionally Coord_Old. I see how this can work for rigid grid movement, as the Grid Velocities are then imposed (so, they are independent variables) and the grid is not really moving physically.

However, for deforming cases we might run into problems, as in that case the Grid Velocities are magnitudes that are derived from the grid positions. In fact, at time n, only the grid velocity u enters the dual-time stepping equation ALE-NS equation, but the dependency of u on previous time steps n, n-1 links directly to the mesh deformations z,
$\mathbf{u}_{grid}^{n+1}=f(\mathbf{z}^{n+1},\mathbf{z}^{n},\mathbf{z}^{n-1})$ ,
where at the same time we would need to differentiate with respect to the variables z for correct geometrical adjoint sensitivities.

Additionally, we already have those variables allocated as Coord_p1, Coord_n and Coord_n1 (in this mesh solver) or Solution, Solution_time_n and Solution_time_n1 in #760. So, really, obtaining the GridVel vector would just be an additional FD computation both for rigid movement and dynamic grids, and a lot of memory savings are possible.

Do you think it could be possible to adapt it for this PR? If not, it should be something to consider as further work (as I would expect there is quite a bit of interest on dynamic meshes), maybe using the basis of #760 (which should make it easier).

rsanfer · 2019-08-23T12:20:11Z

      Restart_Flow = false;

-      if (GetGrid_Movement()) {
-        SU2_MPI::Error("Dynamic mesh movement currently not supported for the discrete adjoint solver.", CURRENT_FUNCTION);


As of now, deformable meshes are mostly used in FSI (where the coupling terms should be accounted for) but outside of that, there are not many methods to run deformable, single-zone test cases available. I guess something could be done using SURFACE_MOVEMENT = DEFORMING and a pitching/plunging airfoil, but I don't believe the adjoint will work out of the box, as the Grid Velocities are recomputed several times in one round of deformation.

However, with the capabilities in #760, I believe it should be relatively easy to set up a case with time-dependent boundary displacements.

rsanfer · 2019-08-23T12:21:00Z

  unsigned short iMesh;
  int Direct_Iter;
  bool heat = config[val_iZone]->GetWeakly_Coupled_Heat();
+  bool grid_IsMoving = config[val_iZone]->GetGrid_Movement();


For consistency, could you rename grid_IsMoving to grid_movement?

@economon mentions that as well above :) but CVolumetricMovement ***grid_movement is an input to that function, so the compiler throws an error declaration of ‘bool grid_movement’ shadows a parameter

rsanfer · 2019-08-23T12:23:39Z

 	 * \return Grid velocity at the point.
 	 */	
 	su2double *GetGridVel(void);
+


Can you indent these functions to align with the rest of the file?

In the file mostly tabs are used for indentation, and github seems to display 4 spaces per tab. The added code follows https://su2code.github.io/docs/Style-Guide/ using 2 spaces as indentation.
I can of course change it to be consistent with the file using tabs if that is wanted. Above in pedros review that topic came already up.

Oh sorry, I didn't notice before. Let's leave it like that then.

rsanfer · 2019-08-23T12:23:49Z

-	 * \param[in] val_value - Value of the gradient.
-	 */
-	void SetGridVel_Grad(unsigned short val_var, unsigned short val_dim, su2double val_value);
+   * \brief Set the value of the grid velocity at the point.


rsanfer

I hadn't noticed there was an approval already, so I will just put the PR in standby temporarily to follow the discussion above.

TobiKattmann · 2019-08-23T15:29:17Z

So, really, obtaining the GridVel vector would just be an additional FD computation both for rigid movement and dynamic grids, and a lot of memory savings are possible.

Hey @rsanfer, I think that is ok for this PR. After a short search I guess such a routine does not exist yet. Or can you point me to the right direction?
Otherwise CVolumetricMovement::ComputeGridVel_from_Coord(..) would be a naming suggestion. What do think?

talbring · 2019-08-23T15:40:51Z

The routine you are looking for is probably this one:

void CPhysicalGeometry::SetGridVelocity(CConfig *config, unsigned long iter) geometry_structure.cpp:14760 (cannot paste a link, the file is too big)

cvencro · 2019-08-24T16:59:34Z

Hi @rsanfer and @TobiKattmann , I agree with the more generalised solution for rigid and deforming meshes too and this is excellent timing of #760 to combine the new approach.
I have FSI cases I can use, but thought it would be a good stepping stone to use fluid-only case to begin with to avoid the coupling terms for this verification. Using dynamic surface movement definitions for the single-zone pitching airfoil test which has previously been used for rigid motion, the primal solutions and gradients from finite differences for both are very similar as expected. Therefore this should be a good reference to check the discrete adjoints for deforming mesh when the adaptations for grid velocities are made using #760.

… Run()

…b.com/su2code/SU2 into feature_dynamic_mesh_discrete_adjoint

…current mesh solver for calculation of grid velocities

…sh_discrete_adjoint Conflicts: Common/src/config_structure.cpp

…sh_discrete_adjoint

TobiKattmann · 2019-09-16T08:57:32Z


-    /*--- Secondary sensitivities must be computed with a certain frequency. ---*/
-    /*--- It is also done at the beginning so all memory gets allocated.     ---*/
-    if ((Adjoint_Iter % config->GetWrt_Sol_Freq() == 0) && (SecondaryVariables != NONE)){


This SecondaryRecording is now removed entirely, as it introduces a bug for the unsteady adjoint. The Volume Mesh Sensitivity setting is a += in each physical timestep for unsteady cases and SecondaryRecording is also called in Postprocess() right after the Run(). I.e. if it is called multiple times in one physical timestep the sensitivity adds up. With this code modification SecondaryRecording gets only called once per timestep in Postprocess.

For steady cases the implications are: The memory for the secondary recording gets allocated after the 'flow adjoint iteration' ... so it could be really annoying if the simulation crashes due to lack of memory after waiting for convergence. More importantly: You cannot use an intermediate restart file for e.g. DOT_AD as it is does not hold any/correct sensitivities. So it might be better to only disable this conditional for unsteady computations as I see the advantage of having this for steady computations. What do you think @rsanfer

I honestly prefer it that way as it is much faster for steady-state computations in moving meshes (a lot less unnecessary recordings). In fact that was my first implementation, but then I added this SecondaryRecording after discussion with @talbring, because of the drawbacks you mention. No 100% what's the best approach...

My only concern is the memory thing, because usually the recording with coordinates as inputs takes more memory than the one with the solution variables. I am not so worried about the intermediate restart files, since you can just do a restart to get the correct sensitivity values.

Lets do it without the SecondaryRecording. If it really leads to problems, we can always change it back.

Ok, sounds good to me.

TobiKattmann · 2019-09-17T20:45:50Z

A little gradient validation using the added Testcase (comp euler, 150 timesteps, pitching airfoil, FFD points as DV, Efficiency obj function).

talbring · 2019-09-20T12:52:02Z

@TobiKattmann The validation looks perfect. Happy that it works now.

talbring · 2019-09-20T12:53:49Z

Remember to revert the Travis changes.

TobiKattmann · 2019-09-23T09:35:58Z

Alright, the su2code/Testcases repo is already merged, .travis.yml is reverted. If there are no further request than this can be merged from my standpoint.
Big thanks to @cvencro for the collaboration :)

talbring · 2019-09-23T10:12:55Z

Thanks. Merging it in. Just as a reminder: without a Tutorial or Documentation, there will be no advertising :/

cvencro added 4 commits August 21, 2019 22:42

tidy comments and indentation

9431716

correct column look-up in restart file

7caf5ea

consistent definition for adjoint boolean

7042b14

revert incompressible ale changes to prepare for PR

597c133

cvencro mentioned this pull request Aug 22, 2019

Unsteady incompressible flow simulation with grid movement #767

Merged

Revert excecutable file permission for cpp/hpp/inl files.

b805962

talbring reviewed Aug 22, 2019

View reviewed changes

pcarruscag reviewed Aug 22, 2019

View reviewed changes

Req. changes for PR766. GridVel allocation, indentation, error msg.

8a4c291

pr-triage Bot added the PR: unreviewed label Aug 22, 2019

economon approved these changes Aug 22, 2019

View reviewed changes

pr-triage Bot added PR: reviewed-approved and removed PR: unreviewed labels Aug 22, 2019

Renamed grid_movement_bool & fixed error catch for dynamic mesh+discAdj.

85047c4

pr-triage Bot added PR: unreviewed and removed PR: reviewed-approved labels Aug 22, 2019

rsanfer reviewed Aug 23, 2019

View reviewed changes

rsanfer requested changes Aug 23, 2019

View reviewed changes

pr-triage Bot added PR: reviewed-changes-requested and removed PR: unreviewed labels Aug 23, 2019

Wrapped GridVel allocation in discAdj conditional. PR766.

c9f002b

pr-triage Bot added PR: unreviewed and removed PR: reviewed-changes-requested labels Aug 23, 2019

cvencro and others added 4 commits August 24, 2019 18:25

Test without additional secondary and main recordings at the start of…

8e04358

… Run()

Merge branch 'feature_dynamic_mesh_discrete_adjoint' of https://githu…

272ee6a

…b.com/su2code/SU2 into feature_dynamic_mesh_discrete_adjoint

Comment out get, set of additional containers and initial test using …

585f702

…current mesh solver for calculation of grid velocities

Merge remote-tracking branch 'origin/develop' into feature_dynamic_me…

ae7460c

…sh_discrete_adjoint Conflicts: Common/src/config_structure.cpp

TobiKattmann added 4 commits September 7, 2019 00:21

Merge remote-tracking branch 'origin/develop' into feature_dynamic_me…

006e015

…sh_discrete_adjoint

RestartOld_Geometry conditional fix for FSI cases.

c3a6e7d

Re-add code deleted during merge conflicts.

2d8b0ea

Added serial+parallel reg test for pitching discAdj.

0d21ee0

TobiKattmann mentioned this pull request Sep 10, 2019

Feature dynamic mesh discrete adjoint su2code/TestCases#37

Merged

TobiKattmann removed the dependent label Sep 10, 2019

Removed doubled SecondaryRecording, which errors unsteady adjoints.

96dec9c

TobiKattmann reviewed Sep 16, 2019

View reviewed changes

TobiKattmann changed the title ~~WIP Feature dynamic mesh discrete adjoint~~ Feature dynamic mesh discrete adjoint Sep 17, 2019

dep Bot added the dependent label Sep 17, 2019

pr-triage Bot added the PR: unreviewed label Sep 17, 2019

Revert travis. Last little change for PR766 Regression test.

2d56130

rsanfer approved these changes Sep 23, 2019

View reviewed changes

pr-triage Bot added PR: reviewed-approved and removed PR: unreviewed labels Sep 23, 2019

Changed DV_* cfg lines such that parallel_computation.py works.

d1262f1

pr-triage Bot added PR: unreviewed and removed PR: reviewed-approved labels Sep 23, 2019

talbring merged commit 1ef7788 into develop Sep 23, 2019

talbring deleted the feature_dynamic_mesh_discrete_adjoint branch September 23, 2019 10:13

pr-triage Bot added PR: merged and removed PR: unreviewed labels Sep 23, 2019

talbring added the changelog:feature label Nov 6, 2019

talbring changed the title ~~Feature dynamic mesh discrete adjoint~~ Time-accurate discrete adjoint gradients for unsteady problems with grid movement Nov 8, 2019

cvencro mentioned this pull request Dec 8, 2019

Discrete adjoint for deforming meshes #833

Merged

5 tasks

Uh oh!

Conversation

cvencro commented Aug 22, 2019 • edited by TobiKattmann Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Related Work

PR Checklist

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

economon left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rsanfer left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

TobiKattmann Aug 23, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rsanfer left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

TobiKattmann commented Aug 23, 2019

Uh oh!

talbring commented Aug 23, 2019

Uh oh!

cvencro commented Aug 24, 2019

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

TobiKattmann commented Sep 17, 2019

Uh oh!

talbring commented Sep 20, 2019

Uh oh!

cvencro commented Aug 22, 2019 •

edited by TobiKattmann

Loading

TobiKattmann Aug 23, 2019 •

edited

Loading

rsanfer left a comment •

edited

Loading

TobiKattmann commented Sep 23, 2019 •

edited

Loading