WIP: Viscoelastic modeling and VIV for Lines

[RyanDavies19]

New Features

I am opening up this draft PR to get the conversation going on the best way to add new states. This is a WIP, the tests are currently failing. This work contains two new features:

1. The first is a vortex-induced-vibration model based on the work of Thorsen et al.. We have a paper submitted for this and when a preprint is released I will link it here. The semi-empirical model gives an oscillating lift force on each node. As part of this implementation, there is a synchronization model that looks at changes between the lift force phase and line motion phase. This model requires knowledge of the lines current acceleration (which we approximate by the previously calculated value) and an additional state variable that is integrated every time step. We validated this approach against the results in the original theory paper using rk4. We did not test this approach with implicit integration schemes, though my hope with the acceleration approximation was that it would work in that case. If, after we resolve the state approach, we can't get acceptable results then I will add a warning that requires users to use explicit schemes.
2. The second is a pair of viscoelastic models that have been implemented in MoorDyn-F to simulate non-linear rope behavior. They are based on the idea of dividing a segment into a spring-dashpot in parallel in series with a spring (a generalized Kelvin–Voigt model). A constant dynamic stiffness approach is explained in this paper. The mean-load dependent dynamic stiffness model builds off this work by changing the dynamic stiffness as a function of mean load. Test results for this are shown in the MoorDyn-F PR (MD: Adding Load dependent dynamic stiffness OpenFAST/openfast#2459). Both of these require an additional state to handle the length of one of the springs in the subdivided segment. Long term, it would be great to allow for n-number of spring damper parallel components, which would require n more states per node. This would give us access to the true generalized kelvin model. (A similar approach could be done with the generalized Maxwell viscoelastic model, also requiring n number of additional states per node).

The documentation updates reflect the input file changes that enable these new features. Other minor changes are:

1. # is now a comment character in the input file and all text on a line after # is ignored. This matches MD-F
2. Added space before messages output by WaterKin to make the terminal output more clear
3. The VIV force and the line curvature at each node are additional line output flags

State Approach

Prior to the improvements to states made by Alex things were working well, since then I have been fighting with failing tests in what seems to be allocation issues for this hackish misc state approach. Rather than continue to bang my head against a wall on something that will probably change, I figured it would be better to discuss approaches first with you all because you have a better sense of how this backend state stuff has been set up. Ideally we would include these two additional states in the line state structure, as they are unique to the lines. This is difficult to do without creating unused states, because the code is currently pretty rigidly structured around the pos, vel, and acc structures for each state (I.e. overloaded operators, implicit integration schemes, stationary solver, StateVar and StateVarDeriv classes). What we need to discuss is how we can adjust this so additional states can be added to objects without needing to completely rework the code each time. One approach I was considering is adjusting the vector lengths of the existing pos, velocity, and acc structures to accommodate these states. For example, in this case with lines. LineState.pos would be a vector of 5 element vectors, and DLineStateDt.vel would be a vector of 5 element vectors. We would then map the integrations of the first 3 components to the traditional approach of acc->vel->pos, and each time step also integrate the additional states held in DLineStateDt.vel. However, this approach isn’t foolproof because there is no check to ensure that the contents of DLineStateDt.vel and LineState.pos are in the right order. If we don’t want to leverage these existing structures then we will probably need to restructure a lot of the existing state and time code.

This addresses the conversation started in #261 and #269 (comment). @AlexWKinley and @sanguinariojoe let me know your thoughts.

TODO:

• Resolve the failing tests
• Develop test for VIV
• Develop test for viscoelastic models
• Clean up and adapt to new state approach
• Clean up commit history to remove duplicates
• Ensure performance is retained with new state approach
• Handle defaults for viscoelastic and VIV states when loading in exported states (from MoorPy, Map++, or MoorDyn)

[sanguinariojoe]

I was waiting for this! Niiiceee! I am checking out ASAP

…

On Tue, 14 Jan 2025, 18:59 Ryan Davies, ***@***.***> wrote: @RyanDavies19 <https://github.com/RyanDavies19> requested your review on: #290 <#290> WIP: Viscoelastic modeling and VIV for Lines. — Reply to this email directly, view it on GitHub <#290 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMXKKBUMFFVSK5KC6W57VL2KVGANAVCNFSM6AAAAABVFOTULSVHI2DSMVQWIX3LMV45UABCJFZXG5LFIV3GK3TUJZXXI2LGNFRWC5DJN5XDWMJVHEZDOOBUGU2TEOI> . You are receiving this because your review was requested.Message ID: ***@***.***>

[sanguinariojoe]

Hey guys!

We come with this solution:

The first significant change is that we are removing StateVarDeriv because:

• It seems that acceleration is becoming a state variable
• Some state variables, without associated derivative, are also landing on the system --e.g. the length of the segments--

The next significant change is that we are not specializing the state variables according to the kind of entity anymore. Instead of that the state variables will be typedefed Eigen matrices besides an indexer which is able to identify each entity with a submatrix. That should allow very idiomatic approaches to operate with the states. E.g. for a explicit 1st order Euler we would have the following pseudo-code:

State a;
(...)
a["r"] += dt * a["v"];

Finally, since we are setting the state variables as a mapped object, we can always install state variables on demand. That should give us functionality enough to have "n more states per node". We just need to install state variables with names "subspring_001", "subspring_002", ...

This is not requiring an unacceptable amount of work. It is not for free either... Please, give us some feedback before we start chewing keystrokes :-p

[RyanDavies19]

Hi @sanguinariojoe, I like this general approach, and I appreciate its flexibility in terms of future work adding new states to various objects. It's good to stay flexible because you never know what next crazy cool idea someone comes up with.

To address your two bullet points: I probably could've made it more clear but every state in the system is found by integrating its derivative in the previous step. Additionally, I'm not quite sure I would consider acceleration a state now as much as we are approximating it in the current step by using the value of the previous step (in the context of the VIV model).

Some clarifying questions so I can wrap my head around it:

1. This structure would give us a N length vector of matrices (NSTATES), where each matrix was N x len(StateVar) in size, and StateVar contains a vector of key indexed states (i.e. r, v ...)? So the only constraint to this approach then is that all nodes in a line need to have the same number of states, but every object can have an in theory abstract len(StateVar)?
2. How and where do we handle the derivatives of these states? For example in the viscoelastic approach we calculate the rate of change of the sub-segment ($\dot{l_1}$ in eqn 16 in the paper), and then integrate it to find the change of the subsegment to use for the tension force ($\Delta l_1$ in eqn 14 in the paper). The same goes for in the VIV approach the frequency being integrated to find the phase in the next step.

Let me know your thoughts on these comments. As always, impressive work from you all coming up with these nice structures.

[sanguinariojoe]

This structure would give us a N length vector of matrices (NSTATES), where each matrix was N x len(StateVar) in size, and StateVar contains a vector of key indexed states (i.e. r, v ...)? So the only constraint to this approach then is that all nodes in a line need to have the same number of states, but every object can have an in theory abstract len(StateVar)?

That's so far the idea. The hierarchy is not exactly that though.

The time integrator will own NSTATES States, and each State will own a map of StateVars, indexed by their names.

Each StateVar is actually a matrix, where the number of rows is invariably N = the number of DOFs (as 1 per point and body, and n_nodes per line and rod). Thus, we can model simple scalars as Eigen::Matrix<real, N, 1>, simple vectors as Eigen::Matrix<real, N, 3>, and for "more complex" types, like the subsprings on VIV, we can probably typedef something for them which gives the same flexibility, ending up with something like Eigen::Matrix<list, N, 1> (list is probably nothing else but another matrix, let's see...).

Besides that, the State will hold also a double-entry indexer which allows to associate instance -> rows and row -> instance (instance = a specific point/body/rod/line).

It will be hereby possible to access state variables to operate them indexing by the name, e.g. state["r"]. I also plan to add accessibility to specific instances, like state(line, "r") (I use parenthesis to get compatibility with C++<23).

Also the structure will make easy (but maybe not that efficient) to remove and add instances (e.g. line break).

The merit is that it will be idiomatically nice, and I suppose quite efficient. The drawback is that some StateVars will be useless for some entities (like the subsprings for points, bodies and rods). This will not impact the performance (I hope).

How and where do we handle the derivatives of these states? For example in the viscoelastic approach we calculate the rate of change of the sub-segment ( l 1 ˙ in eqn 16 in the paper), and then integrate it to find the change of the subsegment to use for the tension force ( Δ l 1 in eqn 14 in the paper). The same goes for in the VIV approach the frequency being integrated to find the phase in the next step.

I plan to get everything ready to can handle derivatives (which will be now StateVars anyway) either by external imposition (so far like it happens right now, with some kind of getStateDeriv) or by internal computation (e.g. from the beginning and ending values).

I have sort of set of methods in mind... But I am afraid it will probably change as I implement all the stuff.

[RyanDavies19]

Okay, I am getting a clearer picture of it now. It seems then that NSTATES in the integrator corresponds to the total number of objects (bodies, rods, points, lines) in the system. The hierarchy below that (State structure and StateVar structures) make sense now. I do like that this will make things easily adaptable, I was actually planning to go in and clean up the line failures feature soon too to bring it up to speed with what I have on the MD-F side of things.

The merit is that it will be idiomatically nice, and I suppose quite efficient. The drawback is that some StateVars will be useless for some entities (like the subsprings for points, bodies and rods). This will not impact the performance (I hope).

So does this mean the state structure will be the same regardless of the instance, but just some values will be empty? I.e. will a body have sub-spring and viv phase states that are just unused? My first impression here was that each state can have a custom set of StateVars.

I plan to get everything ready to can handle derivatives (which will be now StateVars anyway) either by external imposition (so far like it happens right now, with some kind of getStateDeriv) or by internal computation (e.g. from the beginning and ending values).

I'm still not sure how this will work out because every state needs a corresponding derivative, but I trust you have something good in mind.

I have sort of set of methods in mind... But I am afraid it will probably change as I implement all the stuff.

As they always do 🥲. This approach seems good to me if you have time to work on it. I am curious to hear if @AlexWKinley wants to chime in, especially with regards to #269.

[RyanDavies19]

Here's how I am interpreting things, based on your diagram and descriptions. Please let me know if things are off here. I can update accordingly.

[sanguinariojoe]

Hey guys!

Sorry for the delay, I got swamped with other stuff. I am addressing this task now, and let's where I am landing.

@RyanDavies19, on your diagram NSTATES is the number of states required by the Time integrator --so far as it happens now--