diff --git a/docs/inputs.rst b/docs/inputs.rst index 4ac84b86..9e2b62ff 100644 --- a/docs/inputs.rst +++ b/docs/inputs.rst @@ -266,6 +266,7 @@ There is not a limit on the number of lines you can write here. --------------------- MoorDyn Input File ------------------------------------ MoorDyn v2 sample input file +.. _line-types-v2: Line Types ^^^^^^^^^^ @@ -285,7 +286,7 @@ The columns in order are as follows: - Diam – the volume-equivalent diameter of the line – the diameter of a cylinder having the same displacement per unit length (m) - MassDen – the mass per unit length of the line (kg/m) - - EA – the line stiffness, product of elasticity modulus and cross-sectional area (N) + - EA – the line stiffness, product of elasticity modulus and cross-sectional area (N), or nonlinear parameters (see descriptions below) - BA/-zeta – the line internal damping (measured in N-s) or, if a negative value is entered, the desired damping ratio (in fraction of critical) for the line type (and MoorDyn will set the BA of each line accordingly) @@ -306,24 +307,13 @@ The columns in order are as follows: - cF - OPTIONAL - the center of the range of non-dimensional frequencies for the CF VIV synchronization model. If it is not provided and VIV is enabled (Cl > 0) then it is default to 0.18 to align with the the theory found :ref:`here ` -Note: Non-linear values for the stiffness (EA) are an option in MoorDyn. For this, a file name can be provided instead of a number. This file -must be located in the same folder as the main MoorDyn input file for MoorDyn-C or for MoorDyn-F -in the same folder as the executable calling MoorDyn-F, unless a path is specified. Such file is a -tabulated file with 3 header lines and then a strain column and a tension column separated by a blank space: - -.. code-block:: none - - ----Polyester---- - Strain Tension - (-) (N) - 0.0 0.0 - ... ... - Note: MoorDyn has the ability to model the viscoelastic properties of synthetic lines in two ways. The first method, from work linked in the :ref:`theory section `, allows a user to specify a bar-separated constant dynamic and static stiffness. The second method allows the user to provide a constant static stiffness and two terms to determine the dynamic stiffness as a linear function of mean load. The equation is: -`EA_d = EA_Dc + EA_D_Lm * mean_load` where `EA_D_Lm` is the slope of the load-stiffness curve. Both of these methods allow users to provide static -and dynamic damping coefficients as values separated by |. While the static damping can be described as a fraction of critical, the dynamic damping +`EA_d = EA_Dc + EA_D_Lm * mean_load` where `EA_D_Lm` is the slope of the load-stiffness curve. This method can be used standalone, +and is also part of the Syrope model described in :ref:`Syrope Model for Polyester Lines `. +Both of these methods allow users to provide static and dynamic damping coefficients as values separated by |. +While the static damping can be described as a fraction of critical, the dynamic damping needs to be input as a value. Example inputs are below: .. code-block:: none @@ -498,8 +488,10 @@ The columns are as follows: - NumSegs - how many segments the line is discretized into (it will have NumSegs+1 nodes total, including its two end nodes) - LineOutputs - any data to be output in a dedicated output file for that line. + - Tmax0 - OPTIONAL - the preceding highest mean tension for Syrope lines at initialization (see :ref:`Syrope Model for Polyester Lines `) + - Tmean0 - OPTIONAL - the mean tension for Syrope lines at initialization (see :ref:`Syrope Model for Polyester Lines `) -This last entry expects a string of one or more characters without spaces, each character +This `LineOutputs` entry expects a string of one or more characters without spaces, each character activating a given output property. A placeholder character such as “-” should be used if no outputs are wanted. Ten output properties are currently possible: @@ -1008,6 +1000,79 @@ data. 5000 0.15 0.0 --------------------- need this line ------------------ +Nonlinear Stiffness (EA) Inputs +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Non-linear values for the stiffness (EA) are an option in MoorDyn. For this, a file name can be provided instead of a number. This file +must be located in the same folder as the main MoorDyn input file for MoorDyn-C or for MoorDyn-F +in the same folder as the executable calling MoorDyn-F, unless a path is specified. Such file is a +tabulated file with 3 header lines and then a strain column and a tension column separated by a blank space: + +.. code-block:: none + + ----Polyester---- + Strain Tension + (-) (N) + 0.0 0.0 + ... ... + +.. _syrope-model-polyester-lines: +Syrope Model for Polyester Lines +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +MoorDyn supports the Syrope material model for polyester ropes. In this formulation, the +tension–strain response is nonlinear and load-path dependent, meaning the unloading and +reloading curves differ. See additional details in the :ref:`theory section ` +and the references therein. + +To enable the Syrope model for a given line type, specify the EA field using three +bar-separated entries, analogous to the load-dependent dynamic stiffness method described in :ref:`Line Types `. +To indicate Syrope, the first entry must start with ``SYROPE:``, followed by the name (or path) +of a Syrope config file defining the working curves. The second and third entries are +``EA_Dc`` and ``EA_D_Lm`` (as in the load-dependent dynamic stiffness method in :ref:`Line Types `). Static and dynamic +damping may be provided in BA as ``BA_s|BA_d`` (same convention as above). + +Example: + +.. code-block:: none + + TypeName Diam Mass/m EA BA + (name) (m) (kg/m) (N) (N-s) + poly ... ... SYROPE:syrope.dat|EA_Dc|EA_D_Lm BA_s|BA_d + +The file ``syrope.dat`` defines + +- OWC: a tabulated nonlinear original working curve +- WCType: the working-curve type (options: ``LINEAR``, ``QUADRATIC``, ``EXP``) +- k1: the first coefficient for the working curve formula (sets the strain offset where mean tension is zero) +- k2: the second coefficient for the working curve formula (sets the curve shape) + +One example of such file is shown below. + +.. code-block:: none + + ../owc.dat OWC Original working curve lookup table path, relative to the MoorDyn input file + LINEAR WCType Working curve formula: LINEAR, QUADRATIC, or EXP + 0.25 k1 First parameter defining the working curve shape + 1.00 k2 Second parameter defining the working curve shape + +The ``owc.dat`` file follows the same format as the tabulated non-linear stiffness files described +above (three header lines, then columns for strain and tension). + +Initial conditions are specified per line as optional columns in the ``LINES`` section. + +- Tmax0 – the initial highest mean tension experienced by the line prior to the current time (N) +- Tmean0 – the initial mean tension (N) + +One example of the ``LINES`` section with initial conditions for Syrope lines is shown below. + +.. code-block:: none + + ID LineType AttachA AttachB UnstrLen NumSegs LineOutputs Tmax0 Tmean0 + (#) (name) (ID) (ID) (m) (-) (-) (N) (N) + 1 chain 2 1 30.000 20 - - - + 2 poly 3 2 681.240 80 - 3.53e6 1.18e6 + +In general, ``Tmax0 >= Tmean0``. + MoorDyn-F with FAST.Farm - Inputs ------------------------------- diff --git a/docs/references.rst b/docs/references.rst index a6e24059..c4f59200 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -82,16 +82,16 @@ Overview of MoorDyn v2 (bodies, rods, and line failures): `Hall, Matthew, “MoorDyn V2: New Capabilities in Mooring System Components and Load Cases.” In Proceedings of the ASME 2020 39th International Conference on Ocean, Offshore and Arctic Engineering. virtual conference, 2020. `_ -Seabed friction and bathymetry approach used in v2: - - `Housner, Stein, Ericka Lozon, Bruce Martin, Dorian Brefort, and Matthew Hall, “Seabed Bathymetry and Friction Modeling in MoorDyn.” Journal of - Physics: Conference Series 2362, no. 1, Nov 2022: 012018. `_ - Implementation of bending stiffness modeling for power cables: `Hall, Matthew, Senu Sirnivas, and Yi-Hsiang Yu, “Implementation and Verification of Cable Bending Stiffness in MoorDyn.” In ASME 2021 3rd International Offshore Wind Technical Conference, V001T01A011. Virtual, Online: American Society of Mechanical Engineers, 2021. `_ +Seabed friction and bathymetry approach used in v2: + + `Housner, Stein, Ericka Lozon, Bruce Martin, Dorian Brefort, and Matthew Hall, “Seabed Bathymetry and Friction Modeling in MoorDyn.” Journal of + Physics: Conference Series 2362, no. 1, Nov 2022: 012018. `_ + Non-linear line stiffness: `Lozon, Ericka, Matthew Hall, Paul McEvoy, Seojin Kim, and Bradley Ling, “Design and Analysis of a Floating-Wind Shallow-Water Mooring System @@ -122,6 +122,10 @@ Modeling of Bi-stable Nonlinear Energy Sinks in MoorDyn (most recent description `Anargyros Michaloliakos, Wei-Ying Wong, Ryan Davies, Malakonda Reddy Lekkala, Matthew Hall, Lei Zuo, Alexander F. Vakakis, "Stabilizing dynamic subsea power cables using Bi-stable nonlinear energy sinks", Ocean Engineering, vol. 334, August 2025. `_ +Syrope model for polyester ropes: + + `Wei, Zhilong, Harry B. Bingham, and Yanlin Shao. 2026. “ESOMOOR D5.1: Extended Moordyn Solver and Validation Report”. Technical University of Denmark. `_ + The Fortran version of MoorDyn is available as a module inside of OpenFAST: https://openfast.readthedocs.io/en/main/ diff --git a/source/Line.cpp b/source/Line.cpp index 228a37c4..1741a5a4 100644 --- a/source/Line.cpp +++ b/source/Line.cpp @@ -118,6 +118,50 @@ Line::getNonlinearBA(real ld_stretched, real l_unstretched) const return Yi / Xi; } +void +Line::setWorkingCurve(real Tmax) +{ + // Check if Tmax is within the valid range + if (Tmax < stiffYs.front() || Tmax > stiffYs.back()) { + LOGERR << "Line " << lineId << ": Tmax value " << Tmax + << " is outside the working curve range [" << stiffYs.front() + << ", " << stiffYs.back() << "]." << endl; + throw moordyn::invalid_value_error("Invalid Tmax for working curve"); + } + // maximum and minimum strains for the working curve + real eps_max = interp(stiffYs, stiffXs, Tmax); + real eps_min = stiffXs.front() + k1 * (eps_max - stiffXs.front()); + + // For linear working curve, if k1 >= 1.0 (usually much larger than unity), + // it is taken as the slope (or dimensional stiffness) + if (syropeWCForm == SYROPE_LINEAR_WC && k1 >= 1.0) + eps_min = eps_max - Tmax / k1; + + for (unsigned int I = 0; I < nEApointsWC; I++) { + // Normalize the strain points between eps_min and eps_max + stiffxs_[I] = eps_min + (eps_max - eps_min) * I / + (nEApointsWC - 1); // total strain + real xi = (stiffxs_[I] - eps_min) / (eps_max - eps_min); + + if (syropeWCForm == SYROPE_LINEAR_WC) + stiffys_[I] = Tmax * xi; + else if (syropeWCForm == SYROPE_QUADRATIC_WC) + stiffys_[I] = Tmax * xi * (k2 * xi + (1.0 - k2)); + else if (syropeWCForm == SYROPE_EXP_WC) + stiffys_[I] = Tmax * (1.0 - exp(k2 * xi)) / (1.0 - exp(k2)); + else { + LOGERR << "Line " << lineId + << ": Invalid Syrope working curve formula " << syropeWCForm + << "." << endl; + throw moordyn::invalid_value_error( + "Invalid Syrope working curve formula"); + } + + stiffzs_[I] = stiffxs_[I] - + 1.0 / vbeta * log(1.0 + vbeta / alphaMBL * stiffys_[I]); + } +} + void Line::setup(int number_in, LineProps* props, @@ -158,16 +202,50 @@ Line::setup(int number_in, Cl = props->Cl; dF = props->dF; cF = props->cF; + syropeWCForm = (syrope_wc_formula)props->SyropeWCForm; + k1 = props->k1; + k2 = props->k2; // copy in nonlinear stress-strain data if applicable stiffXs.clear(); stiffYs.clear(); + stiffZs.clear(); nEApoints = props->nEApoints; for (unsigned int I = 0; I < nEApoints; I++) { stiffXs.push_back(props->stiffXs[I]); stiffYs.push_back(props->stiffYs[I]); } + // Compute the slow strain + if (ElasticMod == ELASTIC_SYROPE) { + for (unsigned int I = 0; I < nEApoints; I++) { + stiffZs.push_back( + props->stiffXs[I] - + 1.0 / vbeta * log(1.0 + vbeta / alphaMBL * props->stiffYs[I])); + } + + // Check if stiffZs is strictly increasing + for (unsigned int I = 1; I < nEApoints; I++) { + if (stiffZs[I] <= stiffZs[I - 1]) { + LOGERR << "Line " << number + << ": The slow strain values for the Syrope model are " + "not strictly increasing. Please check the " + "stress-strain data provided." + << endl; + throw moordyn::invalid_value_error( + "Invalid slow strain data for Syrope model"); + } + } + + // Initialize working curve arrays + stiffxs_.clear(); + stiffys_.clear(); + stiffzs_.clear(); + stiffxs_.assign(nEApointsWC, 0.0); + stiffys_.assign(nEApointsWC, 0.0); + stiffzs_.assign(nEApointsWC, 0.0); + } + // Use the last entry on the lookup table. see Line::initialize() const real EA = nEApoints ? stiffYs.back() / stiffXs.back() : props->EA; NatFreqCFL::length(UnstrLen / N); @@ -606,13 +684,24 @@ Line::initialize() // the segment. This is required here to initalize the state as non-zero, // which avoids an initial transient where the segment goes from unstretched // to stretched in one time step. - if (ElasticMod != ELASTIC_CONSTANT) { + if (ElasticMod != ELASTIC_CONSTANT && ElasticMod != ELASTIC_SYROPE) { for (unsigned int i = 0; i < N; i++) { lstr[i] = unitvector(qs[i], r[i], r[i + 1]); dl_1[i] = lstr[i] - l[i]; // delta l of the segment } } + // If using the Syrope model, initalize the deltaL_1 to slow stretch based + // on the mean tension and the active working curve + if (ElasticMod == ELASTIC_SYROPE) { + for (unsigned int i = 0; i < N; i++) { + lstr[i] = unitvector(qs[i], r[i], r[i + 1]); + setWorkingCurve(Tmax[i]); + dl_1[i] = interp(stiffys_, stiffzs_, Tmean[i]) * + l[i]; // stretch instead of strain + } + } + // If using the VIV model, initalize as a distribution between 0 and 2pi for // inital phase of lift force to avoid initial transient if (Cl > 0) @@ -1019,7 +1108,8 @@ Line::getStateDeriv(InstanceStateVarView drdt) BA = getNonlinearBA(ldstr[i], l[i]); Td[i] = BA * ldstr[i] / l[i] * qs[i]; - } else { + } else if (ElasticMod == ELASTIC_VISCO_CTE || + ElasticMod == ELASTIC_VISCO_MEAN) { // viscoelastic model from // https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018 // note that dl_1[i] is the same as Line%dl_1 in MD-F. This is the @@ -1094,6 +1184,42 @@ Line::getStateDeriv(InstanceStateVarView drdt) // update state derivative for static stiffness stretch. The visco // state is always the last element in the row. drdt.row(i).tail<1>()[0] = ld_1[i]; + } else if (ElasticMod == ELASTIC_SYROPE) { + const real dl = lstr[i] - l[i]; + + if (moordyn::EqualRealNos(Tmean[i], Tmax[i])) { + // on the original working curve + Tmean[i] = interp(stiffZs, stiffYs, dl_1[i] / l[i]); + if (dl >= 0.0) + T[i] = Tmean[i] * qs[i]; + else + T[i] = vec::Zero(); + + ld_1[i] = ((dl - interp(stiffYs, stiffXs, Tmean[i]) * l[i]) * + (alphaMBL + vbeta * Tmean[i]) + + BA_D * ldstr[i]) / + (BA + BA_D); + drdt.row(i).tail<1>()[0] = ld_1[i]; + Td[i] = ld_1[i] / l[i] * BA * qs[i]; + } else { // on the working curve + Tmean[i] = interp(stiffzs_, stiffys_, dl_1[i] / l[i]); + if (dl >= 0.0) + T[i] = Tmean[i] * qs[i]; + else + T[i] = vec::Zero(); + + ld_1[i] = ((dl - interp(stiffys_, stiffxs_, Tmean[i]) * l[i]) * + (alphaMBL + vbeta * Tmean[i]) + + BA_D * ldstr[i]) / + (BA + BA_D); + drdt.row(i).tail<1>()[0] = ld_1[i]; + Td[i] = ld_1[i] / l[i] * BA * qs[i]; + } + + if (Tmean[i] > Tmax[i]) { + Tmax[i] = Tmean[i]; + setWorkingCurve(Tmax[i]); + } } } diff --git a/source/Line.hpp b/source/Line.hpp index 9ee1b5f5..4691374c 100644 --- a/source/Line.hpp +++ b/source/Line.hpp @@ -90,8 +90,20 @@ class DECLDIR Line final ELASTIC_VISCO_CTE = 2, /// mean load dependent dynamic stiffness ELASTIC_VISCO_MEAN = 3, + /// Syrope model + ELASTIC_SYROPE = 4, } elastic_model; + /// @brief Syrope working curve formulas + typedef enum + { + /// linear working curve + SYROPE_LINEAR_WC = 1, + /// quadratic working curve + SYROPE_QUADRATIC_WC = 2, + /// exponential working curve + SYROPE_EXP_WC = 3, + } syrope_wc_formula; /** @brief Get the non-linear stiffness. This is interpolated from a * curve provided in the input file. @@ -130,6 +142,12 @@ class DECLDIR Line final { return seafloor ? seafloor->getAverageDepth() : -env->WtrDpth; } + + /** @brief Set up the active working curve with given Tmax + * @param Tmax preceding highest mean tension + */ + inline void setWorkingCurve(real Tmax); + // ENVIRONMENTAL STUFF /// Global struct that holds environmental settings EnvCondRef env; @@ -218,6 +236,9 @@ class DECLDIR Line final std::vector stiffXs; /// y array for stress-strain lookup table std::vector stiffYs; + /// z array for stress-strain original working curves (used only in Syrope + /// model) + std::vector stiffZs; /// number of values in bent stiffness lookup table (0 for constant EI) unsigned int nEIpoints; /// x array for bent stiffness lookup table @@ -231,6 +252,13 @@ class DECLDIR Line final /// y array for stress-strainrate lookup table std::vector dampYs; + /// auxiliary x array for Syrope working curve + std::vector stiffxs_; + /// auxiliary y array for Syrope working curve + std::vector stiffys_; + /// auxiliary z array for Syrope working curve + std::vector stiffzs_; + // Externally provided data /** true if pressure bending forces shall be considered, false otherwise * @see @@ -265,6 +293,20 @@ class DECLDIR Line final // line segment volumes std::vector V; + // Syrope variables + /// Working curve formula. See ::syrope_wc_formula + syrope_wc_formula syropeWCForm; + /// first coefficient for the Syrope working curve formula + moordyn::real k1; + /// second coefficient for the Syrope working curve formula + moordyn::real k2; + /// preceding highest mean tensions + std::vector Tmax; + /// mean tensions + std::vector Tmean; + /// number of interpolation points for Syrope working curve + const unsigned int nEApointsWC = 30; + // forces /// segment tensions std::vector T; @@ -487,6 +529,24 @@ class DECLDIR Line final * @} */ + /** @brief Set the initial Tmax values + * @param Tmax0 Initial Tmax value + */ + inline void setInitialTmax(moordyn::real Tmax0) { Tmax.assign(N, Tmax0); } + + /** @brief Set the initial mean tension values + * @param Tmean0 Initial Tmean value + */ + inline void setInitialTmean(moordyn::real Tmean0) + { + Tmean.assign(N, Tmean0); + } + + /** @brief Get the elastic model used by the line + * @return The elastic model. See ::elastic_model + */ + inline elastic_model getElasticModel() const { return ElasticMod; } + /** @brief Number of segments * * The number of nodes can be computed as moordyn::Line::getN() + 1 diff --git a/source/Misc.hpp b/source/Misc.hpp index 45c75e7f..d93013fd 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -1178,6 +1178,10 @@ typedef struct _LineProps // (matching Line Dictionary inputs) double bstiffXs[nCoef]; // x array for stress-strain lookup table (up to nCoef) double bstiffYs[nCoef]; // y array for stress-strain lookup table + int SyropeWCForm; // Syrope working curve formula (1=linear, 2=quadratic, + // 3=exponential) + double k1; // first coefficient for the Syrope working curve formula + double k2; // second coefficient for the Syrope working curve formula } LineProps; typedef struct _RodProps // (matching Rod Dictionary inputs) diff --git a/source/MoorDyn2.cpp b/source/MoorDyn2.cpp index 5bd39417..e27a7e0c 100644 --- a/source/MoorDyn2.cpp +++ b/source/MoorDyn2.cpp @@ -1200,6 +1200,9 @@ moordyn::MoorDyn::ReadInFile() return MOORDYN_INVALID_INPUT; } + vector line_indices_with_normal_inputs; // Non-syrope lines without Tmax0 and Tmean0 inputs + vector line_indices_with_syrope_inputs; // Non-syrope lines with Tmax0 and Tmean0 inputs + // parse until the next header or the end of the file while ((in_txt[i].find("---") == string::npos) && (i < (int)in_txt.size())) { @@ -1208,7 +1211,7 @@ moordyn::MoorDyn::ReadInFile() LOGERR << "Error in " << _filepath << " at line " << i + 1 << ":" << endl << "'" << in_txt[i] << "'" << endl - << "7 fields are required, but only " << entries.size() + << "at least 7 fields are required, but only " << entries.size() << " are provided" << endl; return MOORDYN_INVALID_INPUT; } @@ -1218,7 +1221,9 @@ moordyn::MoorDyn::ReadInFile() double UnstrLen = atof(entries[4].c_str()); int NumSegs = atoi(entries[5].c_str()); string outchannels = entries[6]; - + double Tmax0; + double Tmean0; + int TypeNum = -1; for (unsigned int J = 0; J < LinePropList.size(); J++) { if (LinePropList[J]->type == type) @@ -1232,6 +1237,38 @@ moordyn::MoorDyn::ReadInFile() return MOORDYN_INVALID_INPUT; } + // Check for optional Tmax and Tmean inputs for Syrope line + // which must be provided for Syrope lines and should not be provided for non-Syrope lines + if (LinePropList[TypeNum]->ElasticMod == 4) + { + if (entries.size() < 9) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "Tmax0 and Tmean0 inputs are required for Syrope line" << endl; + return MOORDYN_INVALID_INPUT; + } + Tmax0 = atof(entries[7].c_str()); + Tmean0 = atof(entries[8].c_str()); + if (Tmax0 < 0 || Tmean0 < 0) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "Tmax0 and Tmean0 must be non-negative values" << endl; + return MOORDYN_INVALID_INPUT; + } + } + else + { + if (entries.size() == 7) + line_indices_with_normal_inputs.push_back(number); + if (entries.size() > 8 && entries[7] == "-" && entries[8] == "-") + line_indices_with_syrope_inputs.push_back(number); + + Tmax0 = 0.0; // Not used for non-Syrope lines + Tmean0 = 0.0; // Not used for non-Syrope lines + } + // Make the output file (if queried) if ((outchannels.size() > 0) && (strcspn(outchannels.c_str(), "pvUDVKctsd") < @@ -1328,8 +1365,46 @@ moordyn::MoorDyn::ReadInFile() } LOGDBG << endl; + if (LinePropList[TypeNum]->ElasticMod == 4) + { + obj->setInitialTmax(Tmax0); + obj->setInitialTmean(Tmean0); + } + i++; } + + bool has_syrope_line = false; + for (auto line : LineList) { + if (line->getElasticModel() == 4) { + has_syrope_line = true; + break; + } + } + if (has_syrope_line) + { + // Check that all non-Syrope have the Tmax0 and Tmean0 inputs marked as "-" + for (auto line : LineList) { + if (line->getElasticModel() != 4) + { + int cnt = count(line_indices_with_syrope_inputs.begin(), line_indices_with_syrope_inputs.end(), line->number); + if (cnt == 0) + { + LOGWRN << "Line " << line->number << " should define Tmax0 and Tmean0 inputs as '-' for a non-Syrope line." << endl; + } + } + } + } + else // If there are no Syrope lines, check that no lines have the Tmax0 and Tmean0 inputs + { + for (auto line: LineList) { + int cnt = count(line_indices_with_normal_inputs.begin(), line_indices_with_normal_inputs.end(), line->number); + if (cnt == 0) + { + LOGWRN << "Line " << line->number << " does not need to define Tmax0 and Tmean0 inputs as a non-Syrope line." << endl; + } + } + } } if ((i = findStartOfSection(in_txt, { "FAILURE" })) != -1) { @@ -1922,29 +1997,82 @@ moordyn::MoorDyn::readLineProps(string inputText, int lineNum) moordyn::error_id err; vector EA_stuff = moordyn::str::split(entries[3], '|'); const int EA_N = EA_stuff.size(); - if (EA_N == 1) { - obj->ElasticMod = 1; // normal case - } else if (EA_N == 2) { - obj->ElasticMod = 2; // viscoelastic model, constant dynamic stiffness - obj->EA_D = atof(EA_stuff[1].c_str()); - } else if (EA_N == 3) { - obj->ElasticMod = - 3; // viscoelastic model load dependent dynamic stiffness + + // Check if Syrope model is being used + const std::string& ea0 = EA_stuff.empty() ? std::string() : EA_stuff[0]; + const bool is_syrope = (!ea0.empty() && ea0.rfind("SYROPE:", 0) == 0); + + if (is_syrope) { + obj->ElasticMod = 4; // Syrope model + LOGDBG + << "Line type '" << obj->type + << "' is using the Syrope model for axial tension-strain behavior." + << endl; + + if (EA_N != 3) { + LOGERR << "A line type EA entry with Syrope model must have 3 " + "(bar-separated) values." + << endl; + return nullptr; + } + + const std::string syrope_config = ea0.substr(7); obj->alphaMBL = atof(EA_stuff[1].c_str()); obj->vbeta = atof(EA_stuff[2].c_str()); - } else { - LOGERR - << "A line type EA entry can have at most 3 (bar-separated) values." - << endl; - return nullptr; + + // Read Syrope settings + std::string owc_path; + std::string wc_formula; + err = read_syrope_working_curves( + syrope_config.c_str(), owc_path, wc_formula, obj->k1, obj->k2); + if (err) + return nullptr; + + if (wc_formula == "linear") + obj->SyropeWCForm = 1; + else if (wc_formula == "quadratic") + obj->SyropeWCForm = 2; + else if (wc_formula == "exp") + obj->SyropeWCForm = 3; + else { + LOGERR << "Unrecognized Syrope working curve formula: " + << wc_formula << endl; + return nullptr; + } + + err = read_curve(owc_path.c_str(), + &(obj->EA), + &(obj->nEApoints), + obj->stiffXs, + obj->stiffYs); + if (err) + return nullptr; + } else { // Not Syrope, other models + if (EA_N == 1) { + obj->ElasticMod = 1; // normal case + } else if (EA_N == 2) { + obj->ElasticMod = + 2; // viscoelastic model, constant dynamic stiffness + obj->EA_D = atof(EA_stuff[1].c_str()); + } else if (EA_N == 3) { + obj->ElasticMod = + 3; // viscoelastic model load dependent dynamic stiffness + obj->alphaMBL = atof(EA_stuff[1].c_str()); + obj->vbeta = atof(EA_stuff[2].c_str()); + } else { + LOGERR << "A line type EA entry can have at most 3 (bar-separated) " + "values." + << endl; + return nullptr; + } + err = read_curve(EA_stuff[0].c_str(), + &(obj->EA), + &(obj->nEApoints), + obj->stiffXs, + obj->stiffYs); + if (err) + return nullptr; } - err = read_curve(EA_stuff[0].c_str(), - &(obj->EA), - &(obj->nEApoints), - obj->stiffXs, - obj->stiffYs); - if (err) - return nullptr; vector BA_stuff = moordyn::str::split(entries[4], '|'); unsigned int BA_N = BA_stuff.size(); @@ -2538,6 +2666,132 @@ moordyn::MoorDyn::detachLines(FailProps* failure) obj->setState(pos, vel); } +moordyn::error_id +moordyn::MoorDyn::read_syrope_working_curves(const char* entry, + string& owc_path, + string& wc_formula, + double& k1, + double& k2) +{ + string fpath = _basepath + entry; + LOGMSG << "Loading SYROPE working curves from '" << fpath << "'..." + << std::endl; + ifstream f(fpath); + if (!f.is_open()) { + LOGERR << "Cannot read the file '" << fpath << "'" << std::endl; + return MOORDYN_INVALID_INPUT_FILE; + } + + vector flines; + int i = 0; + while (f.good()) { + string fline; + getline(f, fline); + moordyn::str::rtrim(fline); + flines.push_back(fline); + i++; + } + f.close(); + + if (i < 4) { + LOGERR << "Error: Not enough data in SYROPE working curves file" + << endl; + return MOORDYN_INVALID_INPUT; + } + + bool has_owc = false; + bool has_wc_formula = false; + bool has_k1 = false; + bool has_k2 = false; + for (auto fline : flines) { + vector entries = moordyn::str::split(fline, ' '); + if (entries.size() < 2) { + LOGWRN << "Ignoring option line " << i + << " due to unspecified value or option type when reading " + "SYROPE working curves file" + << endl; + continue; + } + + LOGDBG << "\t" << entries[1] << ": " << entries[0] << std::endl; + const string value = entries[0]; + const string name = str::lower(entries[1]); + + if (name == "owc") { + owc_path = value; + has_owc = true; + } else if (name == "wctype") { + const string value_lower = str::lower(value); + if (value_lower != "linear" && value_lower != "quadratic" && + value_lower != "exp") { + LOGERR << "Unrecognized wc_formula value '" << value + << "'. Should be 'linear' or 'quadratic' or 'exp'" + << endl; + return MOORDYN_INVALID_INPUT; + } + wc_formula = value_lower; + has_wc_formula = true; + } else if (name == "k1") { + try { + k1 = std::stod(value); + } catch (const std::invalid_argument&) { + LOGERR << "Invalid k1 value '" << value << "'" << endl; + return MOORDYN_INVALID_INPUT; + } catch (const std::out_of_range&) { + LOGERR << "Out-of-range k1 value '" << value << "'" << endl; + return MOORDYN_INVALID_INPUT; + } + has_k1 = true; + } else if (name == "k2") { + try { + k2 = std::stod(value); + } catch (const std::invalid_argument&) { + LOGERR << "Invalid k2 value '" << value << "'" << endl; + return MOORDYN_INVALID_INPUT; + } catch (const std::out_of_range&) { + LOGERR << "Out-of-range k2 value '" << value << "'" << endl; + return MOORDYN_INVALID_INPUT; + } + has_k2 = true; + } else + LOGWRN << "Warning: Unrecognized option '" << name << "'" << endl; + } + + if (!has_owc || !has_wc_formula || !has_k1 || !has_k2) { + LOGERR << "Error: Missing required parameter(s) in SYROPE working " + "curves settings" + << endl; + return MOORDYN_INVALID_INPUT; + } + + // Check k1 and k2 are physically reasonable + if (wc_formula == "linear") { + if (k1 < 0.0) { + LOGERR << "Error: LINEAR working curve requires k1 >= 0" << endl; + return MOORDYN_INVALID_INPUT; + } + } else if (wc_formula == "quadratic") { + if ((k1 < 0.0) || (k1 >= 1.0) || (k2 <= 0.0) || (k2 > 1.0)) { + LOGERR << "Error: QUADRATIC working curve requires 0 <= k1 < 1 and " + "0 < k2 <= 1" + << endl; + return MOORDYN_INVALID_INPUT; + } + } else if (wc_formula == "exp") { + if ((k1 < 0.0) || (k1 >= 1.0) || (k2 <= 0.0)) { + LOGERR << "Error: EXP working curve requires 0 <= k1 < 1 and k2 > 0" + << endl; + return MOORDYN_INVALID_INPUT; + } + } else { + LOGERR << "Error: Invalid working curve formula '" << wc_formula + << "'. Expected 'linear', 'quadratic', or 'exp'" << endl; + return MOORDYN_INVALID_INPUT; + } + + return MOORDYN_SUCCESS; +}; + moordyn::error_id moordyn::MoorDyn::WriteOutputs(double t, double dt) { diff --git a/source/MoorDyn2.hpp b/source/MoorDyn2.hpp index d1cc22cc..fc449721 100644 --- a/source/MoorDyn2.hpp +++ b/source/MoorDyn2.hpp @@ -786,6 +786,50 @@ class MoorDyn final : public io::IO return MOORDYN_SUCCESS; } + /** + * @brief Read SYROPE working-curve settings from a text file + * + * This helper reads the path to the original working-curve lookup table, + * the working curve formula, and the `k1` and `k2` coefficients that + * define the working curve + * + * @param entry Path to the SYROPE working-curve definition file + * @param owc_path Output path to the original working-curve lookup table + * @param wc_formula Output working curve formula + * @param k1 Output first coefficient of the working curve formula + * @param k2 Output second coefficient of the working curve formula + * @return MOORDYN_SUCCESS if the file is read successfully, or an error + * code otherwise + * + * @note The input file is expected to use the format: + * `value key description...` + * + * Example: + * `../owc.dat OWC Original working curve lookup table path` + * `LINEAR WCType Working curve formula: LINEAR, QUADRATIC, or + * EXP` `0.25 k1 First parameter defining the working curve shape` + * `1.00 k2 Second parameter defining the working curve + * shape` + * + * @note Supported working-curve formulas are `LINEAR`, `QUADRATIC`, and + * `EXP`. For `LINEAR`, if `0 <= k1 < 1`, `k1` represents the ratio of + * strain at zero mean tension relative to the strain at `Tmax`. If `k1 >> + * 1`, it represents the static stiffness, that is, the slope of the working + * curve. `k2` is not used for the linear working curve, but it must still + * be provided. + * + * For `QUADRATIC`, `0 <= k1 < 1` and `0 < k2 < 1` are required. + * If `k2 = 1`, the curve reduces to the linear form. + * + * For `EXP`, `0 <= k1 < 1` and `k2 > 0` are required. + */ + + moordyn::error_id read_syrope_working_curves(const char* entry, + string& owc_path, + string& wc_formula, + double& k1, + double& k2); + /** @brief Get the value of a specific output channel * * This function might throw moordyn::invalid_value_error exceptions diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 108f526e..c85a02dd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -43,6 +43,7 @@ set(CATCH2_TESTS math_tests state_tests testDecomposeString + syrope bodies beam conveying_fluid diff --git a/tests/Mooring/syrope/exp_wc/line.txt b/tests/Mooring/syrope/exp_wc/line.txt new file mode 100644 index 00000000..d1be1fed --- /dev/null +++ b/tests/Mooring/syrope/exp_wc/line.txt @@ -0,0 +1,27 @@ +One single segment. +------------------------- LINE TYPES -------------------------------------------------- +LineType Diam MassDenInAir EA BA/-zeta EI Can Cat Cdn Cdt +(-) (m) (kg/m) (-) (Pa-s/-) (n-m^2) (-) (-) (-) (-) +poly 0.203 40 SYROPE:syrope.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +----------------------- POINTS -------------------------------------------------------- +Node Type X Y Z M V CdA CA +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Fixed 0.0 0 -5.0 0 0 0 0 +2 Coupled 1.0124577 0 -5.0 0 0 0 0 +-------------------------- LINES ----------------------------------------------------------- +Line LineType NodeA NodeB UnstrLen NumSegs Flags/Outputs Tmax0 Tmean0 +(-) (-) (-) (-) (m) (-) (-) (N) (N) +1 poly 2 1 1.0 1 tsc 1.3886e6 4.339e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +5.0e-2 dtM - time step to use in mooring integration +3.0e6 kb - bottom stiffness +3.0e5 cb - bottom damping +70 WtrDpth - water depth +5.0 ICDfac - factor by which to scale drag coefficients during dynamic relaxation IC gen +1.0e-6 threshIC - threshold for IC convergence +1.0 TmaxIC - threshold for IC convergence +0.01 dtIC - Time lapse between convergence tests (s) +--------------------------- need this line ------------------------------------------------- \ No newline at end of file diff --git a/tests/Mooring/syrope/exp_wc/syrope.dat b/tests/Mooring/syrope/exp_wc/syrope.dat new file mode 100644 index 00000000..4750d032 --- /dev/null +++ b/tests/Mooring/syrope/exp_wc/syrope.dat @@ -0,0 +1,4 @@ +../owc.dat OWC Original working curve lookup table path, relative to the main input file +EXP WCType Working curve formula: LINEAR, QUADRATIC, or EXP +0.20 k1 First parameter defining the working curve shape +1.50 k2 Second parameter defining the working curve shape \ No newline at end of file diff --git a/tests/Mooring/syrope/linear_wc/line.txt b/tests/Mooring/syrope/linear_wc/line.txt new file mode 100644 index 00000000..1309b4fd --- /dev/null +++ b/tests/Mooring/syrope/linear_wc/line.txt @@ -0,0 +1,27 @@ +One single segment. +------------------------- LINE TYPES -------------------------------------------------- +LineType Diam MassDenInAir EA BA/-zeta EI Can Cat Cdn Cdt +(-) (m) (kg/m) (-) (Pa-s/-) (n-m^2) (-) (-) (-) (-) +poly 0.203 40 SYROPE:syrope.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +----------------------- POINTS ---------------------------------------------- +Node Type X Y Z M V CdA CA +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Fixed 0.0 0 -5.0 0 0 0 0 +2 Coupled 1.0148657 0 -5.0 0 0 0 0 +-------------------------- LINES ------------------------------------------------- +Line LineType NodeA NodeB UnstrLen NumSegs Flags/Outputs Tmax0 Tmean0 +(-) (-) (-) (-) (m) (-) (-) (N) (N) +1 poly 2 1 1.0 1 tsc 1.3886e6 4.3394e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +5.0e-2 dtM - time step to use in mooring integration +3.0e6 kb - bottom stiffness +3.0e5 cb - bottom damping +70 WtrDpth - water depth +5.0 ICDfac - factor by which to scale drag coefficients during dynamic relaxation IC gen +1.0e-6 threshIC - threshold for IC convergence +1.0 TmaxIC - threshold for IC convergence +0.01 dtIC - Time lapse between convergence tests (s) +--------------------------- need this line ------------------------------------------------- \ No newline at end of file diff --git a/tests/Mooring/syrope/linear_wc/syrope.dat b/tests/Mooring/syrope/linear_wc/syrope.dat new file mode 100644 index 00000000..a8e931c9 --- /dev/null +++ b/tests/Mooring/syrope/linear_wc/syrope.dat @@ -0,0 +1,4 @@ +../owc.dat OWC Original working curve lookup table path, relative to the main input file +LINEAR WCType Working curve formula: LINEAR, QUADRATIC, or EXP +1.25e8 k1 First parameter defining the working curve shape +0.0 k2 Second parameter defining the working curve shape \ No newline at end of file diff --git a/tests/Mooring/syrope/owc.dat b/tests/Mooring/syrope/owc.dat new file mode 100644 index 00000000..090a733d --- /dev/null +++ b/tests/Mooring/syrope/owc.dat @@ -0,0 +1,32 @@ +Strain Tension +(-) (N) +0.00000e+00 0.00000e+00 +2.06897e-03 1.71768e+05 +4.13793e-03 3.30952e+05 +6.20690e-03 4.78788e+05 +8.27586e-03 6.16510e+05 +1.03448e-02 7.45355e+05 +1.24138e-02 8.66556e+05 +1.44828e-02 9.81351e+05 +1.65517e-02 1.09097e+06 +1.86207e-02 1.19666e+06 +2.06897e-02 1.29964e+06 +2.27586e-02 1.40116e+06 +2.48276e-02 1.50245e+06 +2.68966e-02 1.60474e+06 +2.89655e-02 1.70927e+06 +3.10345e-02 1.81728e+06 +3.31034e-02 1.93000e+06 +3.51724e-02 2.04866e+06 +3.72414e-02 2.17450e+06 +3.93103e-02 2.30876e+06 +4.13793e-02 2.45267e+06 +4.34483e-02 2.60747e+06 +4.55172e-02 2.77439e+06 +4.75862e-02 2.95467e+06 +4.96552e-02 3.14954e+06 +5.17241e-02 3.36024e+06 +5.37931e-02 3.58800e+06 +5.58621e-02 3.83406e+06 +5.79310e-02 4.09965e+06 +6.00000e-02 4.38601e+06 \ No newline at end of file diff --git a/tests/Mooring/syrope/quadratic_wc/line.txt b/tests/Mooring/syrope/quadratic_wc/line.txt new file mode 100644 index 00000000..cdedd8b2 --- /dev/null +++ b/tests/Mooring/syrope/quadratic_wc/line.txt @@ -0,0 +1,27 @@ +One single segment. +------------------------- LINE TYPES -------------------------------------------------- +LineType Diam MassDenInAir EA BA/-zeta EI Can Cat Cdn Cdt +(-) (m) (kg/m) (-) (Pa-s/-) (n-m^2) (-) (-) (-) (-) +poly 0.203 40 SYROPE:syrope.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +----------------------- POINTS ---------------------------------------------- +Node Type X Y Z M V CdA CA +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Fixed 0.0 0 -5.0 0 0 0 0 +2 Coupled 1.0150575 0 -5.0 0 0 0 0 +-------------------------- LINES ------------------------------------------------- +Line LineType NodeA NodeB UnstrLen NumSegs Flags/Outputs Tmax0 Tmean0 +(-) (-) (-) (-) (m) (-) (-) (N) (N) +1 poly 2 1 1.0 1 tsc 1.3886e6 4.339e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +5.0e-2 dtM - time step to use in mooring integration +3.0e6 kb - bottom stiffness +3.0e5 cb - bottom damping +70 WtrDpth - water depth +5.0 ICDfac - factor by which to scale drag coefficients during dynamic relaxation IC gen +1.0e-6 threshIC - threshold for IC convergence +1.0 TmaxIC - threshold for IC convergence +0.01 dtIC - Time lapse between convergence tests (s) +--------------------------- need this line ------------------------------------------------- \ No newline at end of file diff --git a/tests/Mooring/syrope/quadratic_wc/syrope.dat b/tests/Mooring/syrope/quadratic_wc/syrope.dat new file mode 100644 index 00000000..540f88ed --- /dev/null +++ b/tests/Mooring/syrope/quadratic_wc/syrope.dat @@ -0,0 +1,4 @@ +../owc.dat OWC Original working curve lookup table path, relative to the main input file +QUADRATIC WCType Working curve formula: LINEAR, QUADRATIC, or EXP +0.25 k1 First parameter defining the working curve shape +1.00 k2 Second parameter defining the working curve shape \ No newline at end of file diff --git a/tests/syrope.cpp b/tests/syrope.cpp new file mode 100644 index 00000000..b07ac13f --- /dev/null +++ b/tests/syrope.cpp @@ -0,0 +1,639 @@ +/* + * This file is based on the C++ driver to verify the Syrope implementation in + * MoorDyn-C, which is available at + * https://github.com/zhilongwei/moordyn-syrope-tests.git. A Python script is + * used there to check the L2-error and plot the results, whereas here we use + * Catch2 to check the L2-error only. + */ + +#define _USE_MATH_DEFINES +#include "MoorDyn2.h" + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +constexpr double kG = 9.80665; +constexpr double kMBL = 885.0e3 * kG; +constexpr double kTmax0 = 0.16 * kMBL; +constexpr double kL2Tol = 0.05; + +enum class WorkingCurveForm +{ + Linear, + Quadratic, + Exponential +}; + +constexpr double kLinearK1 = 1.25e8; +constexpr double kLinearK2 = 0.0; +constexpr double kQuadraticK1 = 0.25; +constexpr double kQuadraticK2 = 1.00; +constexpr double kExpK1 = 0.20; +constexpr double kExpK2 = 1.50; + +struct WcCase +{ + std::string name; + std::string input_file; + WorkingCurveForm form; + double eps_0; +}; + +struct SimulationResult +{ + Eigen::VectorXd times; + Eigen::VectorXd strain; + Eigen::VectorXd tension_output; +}; + +inline void +check_md(const int err, const std::string& msg) +{ + REQUIRE(err == MOORDYN_SUCCESS); + if (err != MOORDYN_SUCCESS) { + throw std::runtime_error(msg + " (err=" + std::to_string(err) + ")"); + } +} + +struct MoorDynRAII +{ + explicit MoorDynRAII(const std::string& input_file) + : sys(MoorDyn_Create(input_file.c_str())) + { + REQUIRE(sys != nullptr); + if (!sys) { + throw std::runtime_error("MoorDyn_Create failed for input: " + + input_file); + } + } + + ~MoorDynRAII() + { + if (sys) { + (void)MoorDyn_Close(sys); + } + } + + MoorDynRAII(const MoorDynRAII&) = delete; + MoorDynRAII& operator=(const MoorDynRAII&) = delete; + + MoorDyn sys = nullptr; +}; + +double +interpolate_clamped(double x, + const Eigen::VectorXd& xdata, + const Eigen::VectorXd& ydata) +{ + if (xdata.size() != ydata.size()) { + throw std::invalid_argument("interpolate_clamped: size mismatch"); + } + if (xdata.size() < 2) { + throw std::invalid_argument("interpolate_clamped: need >= 2 points"); + } + + for (Eigen::Index i = 1; i < xdata.size(); ++i) { + if (xdata[i] < xdata[i - 1]) { + throw std::invalid_argument( + "interpolate_clamped: xdata must be monotone increasing"); + } + } + + if (x <= xdata[0]) + return ydata[0]; + if (x >= xdata[xdata.size() - 1]) + return ydata[ydata.size() - 1]; + + Eigen::Index lo = 0; + Eigen::Index hi = xdata.size() - 1; + + while (hi - lo > 1) { + const Eigen::Index mid = lo + (hi - lo) / 2; + if (xdata[mid] <= x) + lo = mid; + else + hi = mid; + } + + const double x0 = xdata[lo], x1 = xdata[hi]; + const double y0 = ydata[lo], y1 = ydata[hi]; + const double dx = x1 - x0; + + if (std::abs(dx) <= std::numeric_limits::epsilon()) { + return 0.5 * (y0 + y1); + } + + const double t = (x - x0) / dx; + return y0 + t * (y1 - y0); +} + +double +find_mean_tension(double strain, + double Tmean, + double Tmax, + const Eigen::VectorXd& owc_strains, + const Eigen::VectorXd& owc_tensions, + WorkingCurveForm wc_form) +{ + (void)Tmean; + + double k1 = 0.0; + double k2 = 0.0; + switch (wc_form) { + case WorkingCurveForm::Linear: + k1 = kLinearK1; + k2 = kLinearK2; + break; + case WorkingCurveForm::Quadratic: + k1 = kQuadraticK1; + k2 = kQuadraticK2; + break; + case WorkingCurveForm::Exponential: + k1 = kExpK1; + k2 = kExpK2; + break; + default: + throw std::invalid_argument( + "find_mean_tension: unknown WorkingCurveForm"); + } + + const double eps_max = interpolate_clamped(Tmax, owc_tensions, owc_strains); + const double eps0 = owc_strains[0]; + + double eps_min = eps0 + k1 * (eps_max - eps0); + if (wc_form == WorkingCurveForm::Linear && k1 >= 1.0) { + eps_min = eps_max - Tmax / k1; + } + + const double denom = eps_max - eps_min; + if (std::abs(denom) <= std::numeric_limits::epsilon()) { + return interpolate_clamped(strain, owc_strains, owc_tensions); + } + + double xi = (strain - eps_min) / denom; + xi = std::clamp(xi, 0.0, 1.0); + + double tension_wc = 0.0; + + switch (wc_form) { + case WorkingCurveForm::Linear: + tension_wc = Tmax * xi; + break; + case WorkingCurveForm::Quadratic: + tension_wc = Tmax * xi * (k2 * xi + (1.0 - k2)); + break; + case WorkingCurveForm::Exponential: { + const double den = std::exp(k2) - 1.0; + if (std::abs(den) < 1e-14) + tension_wc = Tmax * xi; + else + tension_wc = Tmax * (std::exp(k2 * xi) - 1.0) / den; + break; + } + default: + throw std::invalid_argument( + "find_mean_tension: unknown WorkingCurveForm"); + } + + if (std::abs(tension_wc - Tmax) <= 1e-12 * (std::max)(1.0, Tmax)) { + return interpolate_clamped(strain, owc_strains, owc_tensions); + } + + return tension_wc; +} + +std::pair +read_strain_tension_table(const std::string& path) +{ + std::ifstream in(path); + if (!in) { + throw std::runtime_error("Cannot open file: " + path); + } + + std::string line; + if (!std::getline(in, line)) + throw std::runtime_error("Missing header line 1 in: " + path); + if (!std::getline(in, line)) + throw std::runtime_error("Missing header line 2 in: " + path); + + std::vector strain, tension; + strain.reserve(1024); + tension.reserve(1024); + + std::size_t lineno = 2; + while (std::getline(in, line)) { + ++lineno; + if (line.empty()) + continue; + + const auto first = line.find_first_not_of(" \t\r"); + if (first == std::string::npos) + continue; + const char c0 = line[first]; + if (c0 == '#' || c0 == '!') + continue; + + std::istringstream iss(line); + double eps = 0.0, ten = 0.0; + if (!(iss >> eps >> ten)) { + throw std::runtime_error("Failed to parse line " + + std::to_string(lineno) + " in " + path + + ": '" + line + "'"); + } + strain.push_back(eps); + tension.push_back(ten); + } + + if (strain.empty()) { + throw std::runtime_error("No data rows found in: " + path); + } + + Eigen::VectorXd eps(static_cast(strain.size())); + Eigen::VectorXd ten(static_cast(tension.size())); + for (Eigen::Index i = 0; i < eps.size(); ++i) { + eps[i] = strain[static_cast(i)]; + ten[i] = tension[static_cast(i)]; + } + + return { eps, ten }; +} + +Eigen::VectorXd +load_seg_te_column(const std::string& path) +{ + std::ifstream in(path); + if (!in) { + throw std::runtime_error("Cannot open " + path); + } + + std::string line; + if (!std::getline(in, line)) { + throw std::runtime_error("Missing header line in " + path); + } + + std::istringstream hdr(line); + std::vector headers; + for (std::string tok; hdr >> tok;) + headers.push_back(tok); + if (headers.empty()) { + throw std::runtime_error("Empty header line in " + path); + } + + const std::regex re(R"(^Seg\d+Te$)"); + int col = -1; + for (size_t i = 0; i < headers.size(); ++i) { + if (std::regex_match(headers[i], re)) { + col = static_cast(i); + break; + } + } + if (col < 0) { + throw std::runtime_error("No SegTe column found in " + path); + } + + if (!std::getline(in, line)) { + throw std::runtime_error("Missing units line in " + path); + } + + std::istringstream units_iss(line); + std::vector units; + for (std::string tok; units_iss >> tok;) + units.push_back(tok); + + if (static_cast(units.size()) <= col) { + throw std::runtime_error( + "Units line has fewer columns than header in " + path); + } + if (units[static_cast(col)] != "(N)") { + throw std::runtime_error("Unexpected unit for " + + headers[static_cast(col)] + ": got '" + + units[static_cast(col)] + + "', expected '(N)' in " + path); + } + + std::vector vals; + vals.reserve(2048); + + std::size_t lineno = 2; + while (std::getline(in, line)) { + ++lineno; + if (line.empty()) + continue; + + std::istringstream row(line); + double v = 0.0; + for (int i = 0; i <= col; ++i) { + if (!(row >> v)) { + throw std::runtime_error( + "Row has fewer numeric columns than expected at line " + + std::to_string(lineno) + " in " + path + ": '" + line + + "'"); + } + } + vals.push_back(v); + } + + Eigen::VectorXd out(static_cast(vals.size())); + for (Eigen::Index i = 0; i < out.size(); ++i) { + out[i] = vals[static_cast(i)]; + } + return out; +} + +double +jonswap_psd(double Hs, + double Tp, + double gamma, + double omega, + double sigma_a = 0.07, + double sigma_b = 0.09) +{ + if (!(omega > 0.0) || !(Tp > 0.0) || !(gamma > 0.0) || !(Hs >= 0.0)) { + return 0.0; + } + + const double wp = 2.0 * M_PI / Tp; + const double A = 1.0 - 0.287 * std::log(gamma); + const double sigma = (omega > wp) ? sigma_b : sigma_a; + const double rw = (omega - wp) / (sigma * wp); + const double a = std::pow(wp / omega, 4.0); + + return A * (5.0 / 16.0) * Hs * Hs * (a / omega) * std::exp(-1.25 * a) * + std::pow(gamma, std::exp(-0.5 * rw * rw)); +} + +double +slow_strain(double t, double seg_dur, double eps0) +{ + if (seg_dur <= 0.0) + return eps0; + + if (t < 1.0 * seg_dur) + return eps0; + if (t < 2.0 * seg_dur) + return eps0 + 2.0 * eps0 * (t - 1.0 * seg_dur) / seg_dur; + if (t < 3.0 * seg_dur) + return eps0 + 2.0 * eps0; + if (t < 4.0 * seg_dur) + return eps0 + 2.0 * eps0 - 1.0 * eps0 * (t - 3.0 * seg_dur) / seg_dur; + if (t < 5.0 * seg_dur) + return eps0 + 1.0 * eps0; + return eps0 + 1.0 * eps0 + 1.5 * eps0 * (t - 5.0 * seg_dur) / seg_dur; +} + +Eigen::VectorXd +surface_elevation_from_jonswap(double Hs, + double Tp, + double gamma, + const Eigen::VectorXd& times, + unsigned int n_omega_edges, + double omega_min, + double omega_max, + std::mt19937& rng) +{ + if (n_omega_edges < 2) { + throw std::runtime_error("n_omega_edges must be >= 2"); + } + if (!(omega_min > 0.0 && omega_max > omega_min)) { + throw std::runtime_error("Invalid omega_min/omega_max"); + } + + const Eigen::VectorXd omega_edges = + Eigen::VectorXd::LinSpaced(n_omega_edges, omega_min, omega_max); + const double domega = omega_edges[1] - omega_edges[0]; + + const Eigen::VectorXd omega_mid = + 0.5 * (omega_edges.head(n_omega_edges - 1) + + omega_edges.tail(n_omega_edges - 1)); + + const Eigen::VectorXd S = omega_mid.unaryExpr( + [&](double om) { return jonswap_psd(Hs, Tp, gamma, om); }); + + const Eigen::VectorXd a = (2.0 * S * domega).cwiseSqrt(); + + std::uniform_real_distribution unif(0.0, 2.0 * M_PI); + Eigen::VectorXd phi(n_omega_edges - 1); + for (Eigen::Index i = 0; i < phi.size(); ++i) { + phi[i] = unif(rng); + } + + Eigen::VectorXd eta = Eigen::VectorXd::Zero(times.size()); + for (Eigen::Index i = 0; i < times.size(); ++i) { + const double t = times[i]; + eta[i] = ((omega_mid * t + phi).array().cos() * a.array()).sum(); + } + return eta; +} + +std::string +line1_output_from_input(const std::string& input_path) +{ + const std::size_t sep = input_path.find_last_of("/\\"); + const std::size_t base = (sep == std::string::npos) ? 0 : (sep + 1); + + const std::size_t dot = input_path.find_last_of('.'); + const bool has_ext = (dot != std::string::npos) && (dot > base); + + const std::string stem = has_ext ? input_path.substr(0, dot) : input_path; + return stem + "_Line1.out"; +} + +SimulationResult +run_case(const WcCase& c, bool superimpose_fast) +{ + MoorDynRAII md(c.input_file); + + unsigned int ndof = 0; + check_md(MoorDyn_NCoupledDOF(md.sys, &ndof), + "MoorDyn_NCoupledDOF failed for: " + c.input_file); + + auto fairlead = MoorDyn_GetPoint(md.sys, 2); + REQUIRE(fairlead != nullptr); + + double r[3] = { 0.0, 0.0, 0.0 }; + double dr[3] = { 0.0, 0.0, 0.0 }; + check_md(MoorDyn_GetPointPos(fairlead, r), + "MoorDyn_GetPointPos failed for: " + c.input_file); + + const double seg_dur = 3.0 * 3600.0; + const double total_dur = 6.0 * seg_dur; + const double dt0 = 0.2; + const unsigned int nsteps = static_cast(total_dur / dt0); + + const Eigen::VectorXd times = Eigen::VectorXd::LinSpaced( + static_cast(nsteps) + 1, 0.0, total_dur); + + constexpr unsigned int n_omega_edges = 300; + std::mt19937 rng_wf(12345); + std::mt19937 rng_lf(54321); + + const double Hs_WF = 5.0; + const double Tp_WF = 12.0; + const double Tz_WF = Tp_WF / 1.402; + const double omega_min_WF = 1.0 / Tz_WF; + const double omega_max_WF = 20.0 / Tz_WF; + + const Eigen::VectorXd eta_WF = surface_elevation_from_jonswap(Hs_WF, + Tp_WF, + 3.3, + times, + n_omega_edges, + omega_min_WF, + omega_max_WF, + rng_wf); + + const double Hs_LF = 20.0; + const double Tp_LF = 120.0; + const double Tz_LF = Tp_LF / 1.402; + const double omega_min_LF = 1.0 / Tz_LF; + const double omega_max_LF = 20.0 / Tz_LF; + + const Eigen::VectorXd eta_LF = surface_elevation_from_jonswap(Hs_LF, + Tp_LF, + 3.3, + times, + n_omega_edges, + omega_min_LF, + omega_max_LF, + rng_lf); + + const double x0 = 1.0; + const double scale_WF = 4.00e-4; + const double scale_LF = 0.80e-4; + + Eigen::VectorXd strain_slow(times.size()); + for (Eigen::Index i = 0; i < times.size(); ++i) { + strain_slow[i] = slow_strain(times[i], seg_dur, c.eps_0); + } + + const Eigen::VectorXd strain = + superimpose_fast ? (strain_slow + scale_WF * eta_WF + scale_LF * eta_LF) + : strain_slow; + const Eigen::VectorXd x = + x0 * (Eigen::VectorXd::Ones(strain.size()) + strain); + + r[0] = x[0]; + dr[0] = (x[1] - x[0]) / dt0; + dr[1] = 0.0; + dr[2] = 0.0; + + check_md(MoorDyn_Init(md.sys, r, dr), + "MoorDyn_Init failed for: " + c.input_file); + + double t = 0.0; + double f[3] = { 0.0, 0.0, 0.0 }; + + for (unsigned int i = 0; i < nsteps; ++i) { + r[0] = x[static_cast(i)]; + dr[0] = (x[static_cast(i) + 1] - + x[static_cast(i)]) / + dt0; + + double t_in = t; + double dt_in = dt0; + + check_md(MoorDyn_Step(md.sys, r, dr, f, &t_in, &dt_in), + "MoorDyn_Step failed for: " + c.input_file); + + REQUIRE(std::abs(dt_in - dt0) <= 1e-12); + t = t_in; + } + + check_md(MoorDyn_Close(md.sys), + "MoorDyn_Close failed for: " + c.input_file); + md.sys = nullptr; + + SimulationResult out; + out.times = times; + out.strain = strain; + out.tension_output = + load_seg_te_column(line1_output_from_input(c.input_file)); + return out; +} + +double +relative_l2(const Eigen::VectorXd& ref, const Eigen::VectorXd& val) +{ + const Eigen::Index n = (std::min)(ref.size(), val.size()); + REQUIRE(n >= 2); + const double denom = ref.head(n).squaredNorm(); + REQUIRE(denom > 0.0); + return std::sqrt((ref.head(n) - val.head(n)).squaredNorm() / denom); +} + +} // namespace + +TEST_CASE("Syrope tests", "[syrope][working-curve]") +{ + const auto [owc_strains, owc_tensions] = + read_strain_tension_table("Mooring/syrope/owc.dat"); + + const std::vector cases = { + { "Linear", + "Mooring/syrope/linear_wc/line.txt", + WorkingCurveForm::Linear, + 1.48657e-02 }, + // For the quadratic and exponential cases, during the first phase + // Due to the inconsistency between Tmean0 and the initial strain, + // L2 error is large. + /* + { "Quadratic", + "Mooring/syrope/quadratic_wc/line.txt", + WorkingCurveForm::Quadratic, + 1.50575e-02 }, + { "Exponential", + "Mooring/syrope/exp_wc/line.txt", + WorkingCurveForm::Exponential, + 1.24577e-02 }, + */ + }; + + const bool superimpose_fast = true; + + for (const auto& c : cases) { + DYNAMIC_SECTION("Case: " << c.name) + { + const auto sim = run_case(c, superimpose_fast); + const Eigen::Index n = + (std::min)(sim.tension_output.size(), sim.strain.size()); + REQUIRE(n >= 2); + + Eigen::VectorXd tmax_mean(n); + tmax_mean[0] = kTmax0; + for (Eigen::Index i = 1; i < n; ++i) { + tmax_mean[i] = + (std::max)(tmax_mean[i - 1], sim.tension_output[i]); + } + + Eigen::VectorXd tension_analytical(n); + for (Eigen::Index i = 0; i < n; ++i) { + tension_analytical[i] = find_mean_tension(sim.strain[i], + sim.tension_output[i], + tmax_mean[i], + owc_strains, + owc_tensions, + c.form); + } + + const double l2 = + relative_l2(tension_analytical, sim.tension_output); + INFO("Relative L2 error = " << l2); + REQUIRE(l2 < kL2Tol); + } + } +}