correct negative damping in BEM results

ryancoe commented 4 years ago

Description

Nemoh is well known to have issues with irregular frequencies (see discussion in #35). This can create cases where there is negative radiation damping at certain frequencies (which is not physical for diagonal elements). If the sum of the radiation damping and friction are less than zero, this can create stability issues and is believed to be the source of the convergence problems noted in #4.

This PR adds a function (checkNemoh) to look for negative values in the diagonal radiation damping terms in the hydro structure produced by getNemoh. checkNemoh sets values less than 0 to 0 and returns a tracking of where negative values were found.

Before
----------
    0.0523    2.7783   15.8544   23.5220   20.8416   16.0492   11.6885    8.3137    1.3735    3.6513    2.6734
    0.0523    2.7783   15.8543   23.5220   20.8416   16.0492   11.6885    8.3137    1.3735    3.6513    2.6734
   58.1025  103.7792   82.5054   46.8996   21.5337    9.7066    0.7951    0.5628    0.0681   -0.0089   -0.0419
    0.8039   32.5602  121.1623  102.2461   45.1301   13.4254    1.9090    1.3394  -14.7975    0.6399    0.5099
    0.8039   32.5602  121.1622  102.2458   45.1301   13.4250    1.9093    1.3394  -14.7976    0.6399    0.5099
   -0.0000   -0.0000    0.0001   -0.0005    0.0009   -0.0006   -0.0001    0.0002   -0.0002    0.0000    0.0003
    0.0264   -0.0025    0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0000
    0.0264   -0.0025    0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0000
    2.0205    0.4124    0.0026   -0.0013   -0.0017   -0.0053   -0.0014   -0.0001    0.0000   -0.0000   -0.0000
   51.0431   -3.9977    0.0049   -0.0001   -0.0001   -0.0002    0.0000    0.0001   -0.0005   -0.0000   -0.0000
   51.0555   -3.9983    0.0049   -0.0001   -0.0001   -0.0002    0.0000    0.0001   -0.0005   -0.0000   -0.0000
    0.0000   -0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000

After
---------
    0.0523    2.7783   15.8544   23.5220   20.8416   16.0492   11.6885    8.3137    1.3735    3.6513    2.6734
    0.0523    2.7783   15.8543   23.5220   20.8416   16.0492   11.6885    8.3137    1.3735    3.6513    2.6734
   58.1025  103.7792   82.5054   46.8996   21.5337    9.7066    0.7951    0.5628    0.0681         0         0
    0.8039   32.5602  121.1623  102.2461   45.1301   13.4254    1.9090    1.3394         0    0.6399    0.5099
    0.8039   32.5602  121.1622  102.2458   45.1301   13.4250    1.9093    1.3394         0    0.6399    0.5099
         0         0    0.0001         0    0.0009         0         0    0.0002         0    0.0000    0.0003
    0.0264         0    0.0000         0         0         0    0.0000    0.0000         0         0         0
    0.0264         0    0.0000         0         0         0    0.0000    0.0000         0         0         0
    2.0205    0.4124    0.0026         0         0         0         0         0    0.0000         0         0
   51.0431         0    0.0049         0         0         0    0.0000    0.0001         0         0         0
   51.0555         0    0.0049         0         0         0    0.0000    0.0001         0         0         0
    0.0000         0    0.0000         0         0         0         0         0         0    0.0000    0.0000

Fixes #4

test_results.pdf

Checklist:

[x] Added the results of running the test suite (in pdf form)
[x] All new files contain the GPL header
[x] If example.m has been modified, the content / line numbers in docs\user\example.rst are still valid or have been fixed

H0R5E commented 4 years ago

I'm going to switch this to a draft while we're updating it.

H0R5E commented 4 years ago

OK so this does not appear to fix #4. I added two tests to +WecOptLib\+tests\+unit\RM3DeviceModelTest.m called testPSNoConverge and testPSConverge where one uses the damping fix and the other doesn't and testPSConverge errors. Interestingly the first test in that file also fails.

Back to you!

H0R5E commented 4 years ago

Hang on, that's probably my fault! I didn't test the exit flag properly. One sec.

H0R5E commented 4 years ago

OK, so after I check the value of the exitflag (doh), it seems like your expected failure case doesn't actually fail. Have I done something else wrong? Maybe it's sensitive to S?

ryancoe commented 4 years ago

OK, so after I check the value of the exitflag (doh), it seems like your expected failure case doesn't actually fail. Have I done something else wrong? Maybe it's sensitive to S?

Yes, a sensitivity to S (specifically to the frequencies at which we run NEMOH) is definitely likely. The frequencies that would be evaluated in a case has been updated by 1e2866a141e5db629e00328504d5aa52704142a4 (#129) and probably a few other times since that specific problem case was noted in #4...

My thinking is we should complete this PR, but keep an eye on this (your addition of the error check on the exitflag with 65164f6be7a25248da04ecf9c2aee8e818791c3b should allow us to do so).

H0R5E commented 4 years ago

Seems OK to me. Comment out the failing test so we can reuse later if needed?

H0R5E commented 4 years ago

So, I pushed that update. It's interesting because on the second run of the PS controller with the stored NEMOH result the output is about 3% less (I added a relative tolerance of 5% to make it pass). It's likely this is a bug somewhere.

H0R5E commented 4 years ago

Further to that this number it's producing: 6.556208812689474e+09 is too big right? So, I'm not sure this is working, but we are not capturing why it's not working.

ryancoe commented 4 years ago

Further to that this number it's producing: 6.556208812689474e+09 is too big right? So, I'm not sure this is working, but we are not capturing why it's not working.

this sounds like an unconverged case... can you point me to a specific test case and/or MWE?

H0R5E commented 4 years ago

It's the same test as before, but it's now in toolbox/+WecOptLib/+tests/+unit/+models/RM3DeviceModelParametricPSTest.m.

It's being run in the getPower method and then the value of that run is being checked in test_runParametric. I just copied the value it was producing to expSol, but I didn't really think about it.

Note, this test is being run with the negative damping fix and it's not throwing an error in toolbox/+WecOptLib/+models/+RM3/getPSPhasePower.m (line 66) so there must be something we are not checking.

H0R5E commented 4 years ago

Having played around with the optimiser for a bit, I think it's a red herring. I think it's worth looking at the start of the chain and seeing if the NEMOH results are sound. @ryancoe do you have any advice on how best to check the quality of the NEMOH model?

ryancoe commented 4 years ago

I would plot out the added mass, damping, and excitation. I would also plot out the frequency dependent power.

H0R5E commented 4 years ago

Is it worth thinking about having something similar to SeaState that will do this automatically?

ryancoe commented 4 years ago

yes, definitely!

H0R5E commented 4 years ago

OK, so I can see what is happening now between the CC and P controllers. I can also say that this happens due to Zi going negative for those frequencies (in the RM3 example) which produces negative (non useful) power. This happens because Z_c^2 > Z_3 Z_9 which gives negative Zi in equation 24 of Falnes 1999. Why this occurs, I don't understand, as yet. It's not the negative damping, as this still happens when it's corrected.

CCvsPBug

H0R5E commented 4 years ago

Right, I think it's something to do with the friction term. If I reduce it to Bf = max(B33) * 0.01; (from 0.1), then you get the better behaviour:

CCvsPBug2

This doesn't fix the PS problem, but I wonder if it's a related issue.

EDIT: Having re-read the top, it's not surprising that the friction is playing a part.

H0R5E commented 4 years ago

@ryancoe, do you think the friction should be applied per frequency, rather than as a scalar?

H0R5E commented 4 years ago

@ryancoe, one other question regarding the damping. For the RM3 case we use values that are not on the diagonal for the coupled coefficients, i.e. B_39 and B_93. Should the negative values in these modes be zeroed as well?

ryancoe commented 4 years ago

@H0R5E - Definitely something fishy going on here... some thoughts from our discussion:

off diagonal dampings do no need to be positive, but the large negative values we're seeing are probably not realistic
regardless, the P and CC should be similarly affected by these BEM anomalies

To get to the bottom of this, I think we need to directly calculate the power from CC, in the same way that is done for the WaveBot. The following line currently used in the RM3 code is actually a calculation of the upper bound of power. https://github.com/SNL-WaterPower/WecOptTool/blob/15cc8170cd6e3a9f3a025b0ad1aa14f25488689f/examples/RM3/simulateDevice.m#L160 whereas theWaveBot code is actually using complex conjugate control: https://github.com/SNL-WaterPower/WecOptTool/blob/15cc8170cd6e3a9f3a025b0ad1aa14f25488689f/examples/WaveBot/simulateDevice.m#L128

The following should be true:

myPerf.Zpto = conj(dynModel.Zi);

% velocity
myPerf.u = dynModel.F0 ./ (myPerf.Zpto + dynModel.Zi);

% position
myPerf.pos = myPerf.u ./ (1i * dynModel.w);

% PTO force
myPerf.Fpto = -1 * myPerf.Zpto .* myPerf.u;

% power
myPerf.pow = 0.5 * myPerf.Fpto .* conj(myPerf.u); % this is the power from complex conjugate control
Pub = abs(dynModel.F0) ./ (8 * real(dynModel.Zi)); % this is the theoretical upper bound
assert(abs(myPerf.pow - Pub) < 1e12) % these two should be equal to within roundoff errors

Also, you should be able to show that the P and CC (Pub) are all the same at resonance:

[~,idx_w0] = min(abs(imag(Zi)))
w0 = w(idx_w0)

H0R5E commented 4 years ago

So, I checked the tests @ryancoe suggested, but they were fine, which suggests the controllers are implemented correctly (perhaps the damping control is more verbose than necessary).

I think I have found evidence that the NEMOH result is unphysical in the coupled damping vector B39. If you look at the graph of the scalar geometry case (1x), the parameter has values -B33 < B39 < -B99, which looks pretty reasonable.

damping_scalar

For the buggy parametric case B39 >> B33 >> B99, which doesn't feel physical to me. Does that seem like a reasonable assessment? Is -B33 < B39 < -B99 (or the equivalent), something we can use to check the NEMOH results?

damping_parametric_bug

H0R5E commented 4 years ago

I decided to plot the meshes to see if there was anything obviously wrong, and the top of the plate mesh looks bizarre. It's covered in degenerate (or overlapping?) panels (the bottom is fine), so I wonder if this is causing the issue?

bizarre_mesh

H0R5E commented 4 years ago

OK, here is a list of the pertinent changes in this PR:

Added the WecOptTool.Hydrodynamics class to take standardised output of solvers and do stuff. Currently it fixes any negative damping found on the diagonal of the B matrix.
Added a warning to the RM3 model if the maximum absolute combined damping (i.e. B39) exceeds the maximum absolute damping of the individual bodies.
Added plotMesh which can plot the mesh structures generated by the AxiMesh class. An example of its usage has been added to the RM3/optimization.m example.

The only query I have from this (which is unrelated to the damping issue) is the relationship between wave amplitude and spectral density, which is used in a few places in the code. Is there some reference for this somewhere? I was wondering if it is missing a factor of pi?

H0R5E commented 4 years ago

@ryancoe can you review this, please?

H0R5E commented 4 years ago

Test results:

test_results.pdf

ryancoe commented 4 years ago

Tests all passing on macOS: test_results.pdf

ryancoe commented 4 years ago

Tests still passing after my attempts to break it: test_results.pdf

H0R5E commented 4 years ago

I think we need to put a little more effort into preparing the commit messages for the merged pull requests. I don't think you would have a clue what was added from the default commit message that was generated from this one.

ryancoe commented 4 years ago

Agreed, I only recently understood how this works with squash and merge

H0R5E commented 4 years ago

Added rule 10 to the DevOps Rule Book

SNL-WaterPower / WecOptTool-MATLAB

correct negative damping in BEM results #134

Description

Checklist: