Added Windowing Regularizers for time averaged outputs. by ScSteffen · Pull Request #836 · su2code/SU2

ScSteffen · 2019-12-11T09:58:20Z

Proposed Changes

Give a brief overview of your contribution here in a few sentences.
Based on the observations of Albring et al. in "Challenges in Sensitivity Computations for (D)DES and URANS" we implemented the windowing regularizers proposed in the work of Krakos et al. "Sensitivity analysis of limit cycle oscillations" and adapted it to the AD framework of SU2.
We call the feature in the following "windowing".
Windowing can be applied if one takes the (windowed) time average of a periodic function, i.e. in the case of an instationary flow field that exhibits periodic behaviour.

We changed the functionality of the output-field-prefix TAVG and D_TAVG from a simple time average to a time averaged weighted with a windowing function. For the direct run this is just a post-processing step and changes happen only in signal_processing_toolbox and Coutput.
We adapted the seeding for the adjoint run from a the classic time average to a windowed time average.
Tests have shown, that higher order windowing regularizers, tend to converge quickly with respect to the averaged time span. We have implemented a "timeConvergenceMonitor", that acts similarly to the convergence monitoring of the inner convergence. It applies a cauchy convergence criterion to the chosen windowed time average output coefficients.
The cauchy convergence criterion works for direct run, adjoint run and shape optimization scripts
In case of shape optimization, the script remembers the final time iteration of the direct run. If the final time step is smaller than the "original" final time step stated in the configuration file, it adapts the starting time of the adjoint run; and the time to seed the objective function, such that the seeding does not exceed the original time to end the seeding.

Related Work

Resolve any issues (bug fix or feature request), note any related PRs, or mention interactions with the work of others, if any.

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

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.

…d time average with window specified in WINDOW_FUNCTION.Successfully tested on local machine with 4 cores.

…test Adjoint seeding

…t for cluster

…ual 'old' on local machine, on cluster deviation below 1% still don't know why

…ith all dependencies; direct_differentiation.py: Cauchy Crit added, adjustments to windowed output made; discrete_adjoint: CauchyCrit addded. Code validated on local machine

pcarruscag

Hi @SteffenMath,

Welcome to SU2 (it seems it is your first PR), it looks like a good generalisation of the moving averages, I leave some small code comments below. As a general point does this averaging rely on fixed timesteps?

Cheers,
Pedro

pcarruscag · 2019-12-11T11:41:20Z

Common/include/config_structure.hpp

+  unsigned long GetStartWindowIteration(void);
+
+  /*!
+   * \brief Get Index of the window function used as weight in the cost functional
+   * \return
+   */
+  unsigned short GetWindowIdx(void);


These two methods should be const.
To be consistent with the prevailing naming convention GetWindowIdx should be called something like GetKindWindow (have a look at similar methods). Also, although not done in other methods, the return type for such functions should be the enumerator type (that can avoid some mix-ups).

pcarruscag · 2019-12-11T11:43:11Z

Common/include/option_structure.hpp

+/*!
+ * \brief types of window (weight) functions for cost functional
+ */
+enum WND_FUNCTION {


I would not abbreviate "window".

The enum is now called WINDOW_FUNCTION

pcarruscag · 2019-12-11T11:48:24Z

Common/include/toolboxes/signal_processing_toolbox.hpp

+    su2double GetWndWeight(int fctIdx, unsigned long CurTimeIdx, unsigned long endTimeIdx);
+
+protected:
+    const su2double PI_NUMBER = 4.0 * atan(1.0);	/*!< \brief Pi number. */


I think pi is already defined somewhere as a literal, or definitely in some library header file, it would be better to use one of those instead of storing it as a class member.

I'm using now the constant defined in "option_structure.hpp".

pcarruscag · 2019-12-11T11:50:50Z

Common/include/toolboxes/signal_processing_toolbox.hpp

+    WindowingTools(){}
+    ~WindowingTools(){}
+
+    /*! \brief Returns the value of a windowing function given by fctIdx at CurTimeIdx with given endTimeIdx (i.e. endTimeIdx=nTimeIter, if one starts  windowing at time t =0.) */
+    su2double GetWndWeight(int fctIdx, unsigned long CurTimeIdx, unsigned long endTimeIdx);


We do this in other places, but if the ctor and dtor are empty you are better off not defining them.
User defined destructors disable compiler generated move/copy constructors.

I've deleted both constructor and destructor.

pcarruscag · 2019-12-11T12:00:42Z

Common/src/toolboxes/signal_processing_toolbox.cpp

+su2double WindowingTools::GetWndWeight(int fctIdx, unsigned long CurTimeIdx, unsigned long endTimeIdx){
+    switch (fctIdx){
+      case 1: return HannWindow(CurTimeIdx, endTimeIdx);
+      case 2: return HannSquaredWindow(CurTimeIdx, endTimeIdx);
+      case 3: return BumpWindow(CurTimeIdx, endTimeIdx);
+      default:return 1.0;
+    }
+}


Here the type of fctIdx should be the enumerator, and the "cases" should be actual names not the numbers.
All these functions can be const or static if you use the globally defined pi.

GetWndWeight now uses the enum WINDOW_FUNCTION.

pcarruscag · 2019-12-11T12:01:10Z

Common/src/toolboxes/signal_processing_toolbox.cpp

+    }
+}
+su2double WindowingTools::HannWindow(unsigned long i, unsigned long endTimeIdx){
+  if(i==0) return 0; //catch div by zero error


There is no risk of division by 0

I refactored that

…lprocessingtoolbox to CWindowingTools; Implemented time convergence monitor in SU2_CFD for singlezone direct runs; Adapted python scripts such that shape_optimization.py, direct_differentiation.py and discrete_adjoint.py works with time convergence functionality; Adapted tavg ouptut in python scripts such that they can handle windowed averaged output safely.

…ols files

ScSteffen · 2019-12-17T12:39:39Z

Hi @SteffenMath,

Welcome to SU2 (it seems it is your first PR), it looks like a good generalisation of the moving averages, I leave some small code comments below. As a general point does this averaging rely on fixed timesteps?

Cheers,
Pedro

Hi @pcarruscag,

thanks for your quick reply and the usefull comments. I've pushed and updated version, that includes changes to your suggestions.
Regarding your general question: The windowing approach is in theory not dependent on a fixed time-stepsize, however the current implementation relies on a fixed stepsize. This can be changed quite easily in the CWindowingTools.cpp file and the corresponding class, should there be demand.

Best
Steffen

…ective function for sensitivity computation is defined

… feature_WND

pcarruscag

See my note on multizone cases but otherwise I think the c++ side is good to go.
It would be good to have a tutorial on the feature or to at least to update the list of options.
Someone more familiar than me with the python side should have a look. If no one has anything to say I'll approve next week, thank you again for the contribution!

pcarruscag · 2020-01-07T14:05:04Z

SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp

  if (time_stepping){
    if (TimeIter < IterAvg_Obj){
-      seeding = 1.0/((su2double)IterAvg_Obj);
+      seeding = windowEvaluator.GetWndWeight(config->GetKindWindow(),TimeIter, IterAvg_Obj-1)/ ((su2double)IterAvg_Obj);


Don't use "old style" casts (compiler warning) use static_cast or su2double(IterAvg_Obj)

Fixed it :)

pcarruscag · 2020-01-07T14:05:55Z

SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp

+
+bool CDiscAdjSinglezoneDriver::GetTimeConvergence() const{
+  return false;
+}


This could be inline in the hpp file.

Not addressed yet

Sry forgot about this one. Changed it in the last commit :)

pcarruscag · 2020-01-07T14:10:04Z

SU2_CFD/src/output/COutput.cpp

+  for (unsigned short iField = 0; iField < config->GetnWndConv_Field(); iField++){
+    wndConvFields.emplace_back(config->GetWndConv_Field(iField));
+  }


General note (not very important in this case) whenever you know how many things you are inserting into the vector you should reserve before.

Okay, I added it nevertheless. Thanks for the feedback!

pcarruscag · 2020-01-07T14:18:48Z

SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp

+   * \brief gets Convergence on physical time scale, (deactivated in adjoint case)
+   * \param none
+   */
+  virtual bool GetTimeConvergence() const;


The virtual specifier you should use here is override i.e. bool GetTimeConvergence() const override, this immediately says to anyone reading the code that you are overriding functionality inherited from a parent class.
virtual should be used on the first class in the hierarchy that defines the method.

Okay I'm going to fix that.

pcarruscag · 2020-01-07T14:20:17Z

SU2_CFD/include/drivers/CSinglezoneDriver.hpp

+   * \brief gets Convergence on physical time scale
+   * \param none
+   */
+  virtual bool GetTimeConvergence() const;


What are the implications for multizone problems?

At the moment, the convergence tools only work for single-zone problems, however it should be easy to extend it. I will speak about it with @talbring.

Ok. If you want to merge the PR before that is done please add an error message or warning somewhere if someone tries to use the feature in multizone problems.

I've added an error message for this case.

pcarruscag · 2020-01-07T14:30:26Z

TestCases/disc_adj_rans/cylinder/cylinder_Windowing.cfg

@@ -0,0 +1,416 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


As a general note on examples, I would try to clean the config as much as possible to really put the focus on the new options you added, for example if MG_LEVEL=0 all other MG options don't do anything, if you use the Roe scheme the JST coefficients are not important and so on.

I have removed default options, that are unnecessary for a meaningful test case.

talbring

Thanks @SteffenMath !

I was following up with him in person. So no further comments from my side, everything looks good.

@pcarruscag Small documentation/tutorials will be added soon.

…f one tries to use the windowing Convergence features in Multizone cases

pcarruscag

Lgtm. Thanks @SteffenMath, just the code factor issues left to address.

TobiKattmann · 2020-01-10T23:38:22Z

I'll read through the PR over the weekend, so I would appreciate this staying open until Monday.

TobiKattmann

Hi @SteffenMath ,

I have to read the mentioned papers and try it out to understand everything a bit better, but still it looks really good.
There are some points concerning comments / doxygen-documentation which can be improved, especially the hpp-file without doxygen-documentation. But I blame primarily your supervisor for that 😉

Cheers, Tobi

Common/include/config_structure.hpp

SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp

SU2_CFD/include/output/COutput.hpp

SU2_CFD/include/output/tools/CWindowingTools.hpp

SU2_CFD/src/drivers/CSinglezoneDriver.cpp

SU2_CFD/src/meson.build

SU2_CFD/src/output/tools/CWindowingTools.cpp

TobiKattmann

Hi Steffen,
just a few more things below.
If you address these then this can be merged from my point of view.

TobiKattmann · 2020-01-14T14:16:25Z

SU2_CFD/include/drivers/CSinglezoneDriver.hpp

+  /*!
+   * \brief gets Convergence on physical time scale
+   * \param none
+   */
+  inline virtual bool GetTimeConvergence() const{


Comment cleanup

TobiKattmann · 2020-01-14T14:45:14Z

SU2_CFD/include/output/tools/CWindowingTools.hpp

+  * \param endTimeIdx - Number of time steps to average over
+  * \return Value of the window-weight-function at time CurTimeIdx with end-time endTimeIdx
+  */
+  su2double HannWindow(unsigned long i, unsigned long endTimeIdx) const;


As far as I understand you are using unsigned long i CurTimeIdx and currentIter in this file which all have the same meaning : Current time iteration of the solver
For consistency it would be nice to choose one for your implementation, preferably CurTimeIdx or something new CurrTimeIter

I replaced the variable-names in questions with curTimeIter.

TobiKattmann · 2020-01-14T14:50:15Z

SU2_CFD/include/output/tools/CWindowingTools.hpp

+  /*! \brief Returns the value of a windowing function given by fctIdx at time-step CurTimeIdx with given time-frame endTimeIdx (i.e. endTimeIdx=nTimeIter, if one starts  windowing at time t =0.)
+   * \param windowId - enum specifing the used window
+   * \param CurTimeIdx - Current time iteration of the solver
+   * \param endTimeIdx - Number of time steps to average over
+   * \return Value of the window-weight-function at time CurTimeIdx with time-frame endTimeIdx
+   */
+  su2double GetWndWeight(WINDOW_FUNCTION windowId, unsigned long CurTimeIdx, unsigned long endTimeIdx) const;


here (in the \brief) it is probably windowId instead of fctIdx?

You are right.

ScSteffen added 10 commits October 23, 2019 10:02

initial commit, added only adoint wnd features

cea3d34

Added Windowing Functionality for direct runs, TAVG gives now windowe…

4060d8d

…d time average with window specified in WINDOW_FUNCTION.Successfully tested on local machine with 4 cores.

removed parts of signalprocessing, to debug false AD run

cf1ab30

refactoring of Signalprocessing, disabled all windowing functions to …

00a77ad

…test Adjoint seeding

minimal version for cluster

d6edb0c

Split windowing functionality between adoint run and post processing

c439158

inheritance of windowedAverages activated, COutput with windowing_Tes…

c8839ac

…t for cluster

Activated windowed direct output; Current Signalprocessing results eq…

acb757d

…ual 'old' on local machine, on cluster deviation below 1% still don't know why

shape_optimization.py: Cauchy Crit added (not yet added in SU2_CFD) w…

56ee75f

…ith all dependencies; direct_differentiation.py: Cauchy Crit added, adjustments to windowed output made; discrete_adjoint: CauchyCrit addded. Code validated on local machine

Merge successfull, scripts teste

714d7e9

pr-triage bot added the PR: draft label Dec 11, 2019

pcarruscag reviewed Dec 11, 2019

View reviewed changes

pr-triage bot removed the PR: draft label Dec 11, 2019

pr-triage bot added the PR: draft label Dec 17, 2019

removed useless comment in tools.py; added newly created CWindowingTo…

c8add4f

…ols files

pr-triage bot removed the PR: draft label Dec 17, 2019

re-added output functionality for time steady cases

145eaa0

pr-triage bot added the PR: draft label Dec 17, 2019

ScSteffen and others added 10 commits December 17, 2019 14:58

added error message in case steady outputfield was not found

3d444da

fixed failed test case, where restart file was not copied correctly

c1c6222

fixed issue, where parallel_computation fails, since naturally no obj…

dda46f7

…ective function for sensitivity computation is defined

adapted notation to SU2 python standard

7dcaa0e

small bug in CWindowingTools fixed; TestCases Added

315df14

updated test case name and reconfigured it

0957809

Merge branch 'develop' into feature_WND

0f18c53

configured test case

9c97e6f

Merge branch 'feature_WND' of https://github.com/SteffenMath/SU2 into…

320fadd

… feature_WND

further configuration of the test case

80c22e1

ScSteffen requested a review from pcarruscag January 6, 2020 15:34

small bug fix in output

d837e18

pcarruscag reviewed Jan 7, 2020

View reviewed changes

talbring approved these changes Jan 7, 2020

View reviewed changes

pr-triage bot added PR: reviewed-approved and removed PR: unreviewed labels Jan 7, 2020

removed old style casts, added override keyword to funcion

9091c9b

pr-triage bot added PR: unreviewed and removed PR: reviewed-approved labels Jan 8, 2020

ScSteffen requested a review from pcarruscag January 9, 2020 10:11

Added windowing support for Multizone Adjoint. Added Error Message, i…

25b1262

…f one tries to use the windowing Convergence features in Multizone cases

pcarruscag approved these changes Jan 10, 2020

View reviewed changes

pr-triage bot added PR: reviewed-approved and removed PR: unreviewed labels Jan 10, 2020

talbring added the changelog:feature label Jan 10, 2020

TobiKattmann requested changes Jan 12, 2020

View reviewed changes

pr-triage bot added PR: reviewed-changes-requested and removed PR: reviewed-approved labels Jan 12, 2020

Merge branch 'develop' into feature_WND

9d9e65b

pr-triage bot added PR: unreviewed and removed PR: reviewed-changes-requested labels Jan 14, 2020

added changes suggested in PR

527fe23

ScSteffen requested a review from TobiKattmann January 14, 2020 13:37

TobiKattmann approved these changes Jan 14, 2020

View reviewed changes

renaming of variables suggested in PR

e21fac2

talbring merged commit eb8584f into su2code:develop Jan 15, 2020

pr-triage bot added PR: merged and removed PR: unreviewed labels Jan 15, 2020

TobiKattmann mentioned this pull request Apr 11, 2022

Fix running averages for restarted simulations #1601

Merged

5 tasks

		@@ -0,0 +1,416 @@
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Conversation

ScSteffen commented Dec 11, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Related Work

PR Checklist

Uh oh!

pcarruscag 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!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ScSteffen Dec 17, 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!

Choose a reason for hiding this comment

Uh oh!

ScSteffen Dec 17, 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!

ScSteffen commented Dec 17, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pcarruscag 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!

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!

Choose a reason for hiding this comment

Uh oh!

talbring left a comment

Choose a reason for hiding this comment

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

TobiKattmann commented Jan 10, 2020

Uh oh!

ScSteffen commented Dec 11, 2019 •

edited

Loading

ScSteffen Dec 17, 2019 •

edited

Loading

ScSteffen Dec 17, 2019 •

edited

Loading

ScSteffen commented Dec 17, 2019 •

edited

Loading