Fixing initial state values in MPhys

[lamkina]

Purpose

This PR fixes a bug where the initial ADflow states were being set to a vector of all ones when running with MPhys. This occurred because OpenMDAO defaults the initial value of output vectors to all ones.

To fix this issue, I overwrote the _set_vectors() method of the ADflowSolver component to update the initial value of the output vector.

This PR should be merged after PRs #307 and #308 because it includes changes from both to control the print outs in the _set_vectors() method and to ensure the adjoint works during mesh warping.

Expected time until merged

A week or two, depending on PRs #307 and #308

Type of change

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (non-backwards-compatible fix or feature)
• Code style update (formatting, renaming)
• Refactoring (no functional changes, no API changes)
• Documentation update
• Maintenance update
• Other (please describe)

Testing

We don't have testing for MPhys, but I tested this change on a handful of local test cases including BC and geometric design variables. Reviewers should test this with their own MPhys cases if possible.

Checklist

• I have run flake8 and black to make sure the Python code adheres to PEP-8 and is consistently formatted
• I have formatted the Fortran code with fprettify or C/C++ code with clang-format as applicable
• I have run unit and regression tests which pass locally with my changes
• I have added new tests that prove my fix is effective or that my feature works
• I have added necessary documentation

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 12 lines in your changes missing coverage. Please review.

Project coverage is 41.02%. Comparing base (6706920) to head (6c77d89).
Report is 3 commits behind head on main.

Files with missing lines	Patch %	Lines
adflow/mphys/mphys_adflow.py	0.00%	12 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #321      +/-   ##
==========================================
- Coverage   41.13%   41.02%   -0.11%     
==========================================
  Files          13       13              
  Lines        4108     4119      +11     
==========================================
  Hits         1690     1690              
- Misses       2418     2429      +11     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[lamkina]

bump @A-CGray @anilyil

[eirikurj]

@lamkina after merging #307, there seem to be still some items here that we resolved in #307. Please check.

[eirikurj]

Seems that this was retained after merging (same for above). I think you probably want to restore the previous state as done in PR #307. See here for original comment.

[lamkina]

Good catch. I'm not sure how I did this, but nonetheless I'll get it fixed.

[eirikurj]

Same here, this should not be here, like in PR #307. See here for original comment.

[A-CGray]

@lamkina do you have any suggestions on how I should test this?

[lamkina]

I would check out the main branch and run the tutorial wing case. Print the initial state vector before ADflow runs the first analysis and you'll see it's an array of all ones. Then switch to this branch and run the same case and you should see the initial state vector is the correct one from ADflow.

You can confirm the state vector is correct by running the case with MACH only and printing the initial state vector to compare.

[A-CGray]

I'm confused, if I add code printing out the state on the line in solve_nonlinear right before the ADflow solver is called then this change seems to be working, but the initial residual printed out by ADflow is the same in both cases:

Andrew's branch

self.solver.getStates()=array([1.00000000e+00, 9.46248398e-01, 0.00000000e+00, ...,
       2.47783863e-02, 2.94800000e+00, 1.86845485e-07])
-> Alpha... 1.500000 
#
# Grid 1: Performing 1000 iterations, unless converged earlier. Minimum required iteration before NK switch:      5. Switch to NK at totalR of:   8.73E+01
#
#------------------------------------------------------------------------------------------------------------------------------------------------------------
#  Grid  | Iter | Iter |  Iter  |   CFL   | Step | Lin  |        Res rho         |         C_lift         |        C_drag          |        totalRes        |
#  level |      | Tot  |  Type  |         |      | Res  |                        |                        |                        |                        |
#------------------------------------------------------------------------------------------------------------------------------------------------------------
      1       0      0     None     ----    ----   ----   4.7365425919241388E+03   9.5383731214195250E-04   1.2235096525765710E-01   8.7316352022811938E+06 
      1       1      3     *ANK   5.00E+00  0.15  0.001   4.1598265494843654E+03   1.7769816557941003E-02   1.4384704885726513E-01   7.7007228271569945E+06 
      1       2      8      ANK   5.00E+00  0.26  0.006   3.3403259619915884E+03   4.3599101742336291E-02   1.7417888366812920E-01   6.1704190701800054E+06 
      1       3     13      ANK   5.00E+00  0.32  0.027   2.5844088319778084E+03   6.7937514673665350E-02   2.0130232788327942E-01   4.7463852232853156E+06 
      1       4     19      ANK   5.00E+00  0.34  0.028   1.9900973496159793E+03   8.6145282454101885E-02   2.2036009068597789E-01   3.6286528693136014E+06

Main branch

self.solver.getStates()=array([1., 1., 1., ..., 1., 1., 1.])
-> Alpha... 1.500000 
#
# Grid 1: Performing 1000 iterations, unless converged earlier. Minimum required iteration before NK switch:      5. Switch to NK at totalR of:   8.73E+01
#
#------------------------------------------------------------------------------------------------------------------------------------------------------------
#  Grid  | Iter | Iter |  Iter  |   CFL   | Step | Lin  |        Res rho         |         C_lift         |        C_drag          |        totalRes        |
#  level |      | Tot  |  Type  |         |      | Res  |                        |                        |                        |                        |
#------------------------------------------------------------------------------------------------------------------------------------------------------------
      1       0      0     None     ----    ----   ----   4.7365425919241388E+03   9.5383731214195250E-04   1.2235096525765710E-01   8.7316352022811938E+06 
      1       1      3     *ANK   5.00E+00  0.15  0.001   4.1598265494843654E+03   1.7769816557941003E-02   1.4384704885726513E-01   7.7007228271569945E+06 
      1       2      8      ANK   5.00E+00  0.26  0.006   3.3403259619915884E+03   4.3599101742336291E-02   1.7417888366812920E-01   6.1704190701800054E+06 
      1       3     13      ANK   5.00E+00  0.32  0.027   2.5844088319778084E+03   6.7937514673665350E-02   2.0130232788327942E-01   4.7463852232853156E+06

~8e6 seems like a pretty reasonable initial residual, so my best guess is that there is something in ADflow's __call__ method that is overwriting the [1, 1, 1, ...] initial state anyway.

[lamkina]

Huh, I was not seeing this behavior in my initial testing. It's concerning that ADflow is overwriting a state vector internally that should be coming from OpenMDAO. This isn't happening other than the initialization, but still weird.

[lamkina]

@anilyil @eirikurj do either of you have any idea what could be happening here?

[A-CGray]

@lamkina I think we should just merge this PR, regardless of the reason that the initial state seemed to be being set correctly before by the ADflow, explicitly defining the actual correct initial states in OpenMDAO seems like a more correct approach.

[lamkina]

@anilyil can you take a quick look? Should be an easy review.