From 698bc17ba4a2371f49b68265ec67c168eaed1079 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 12 Jan 2026 12:49:46 +0000 Subject: [PATCH 01/14] Syrope model implementation --- source/Line.cpp | 103 ++++++++++++++++++++++++++++- source/Line.hpp | 63 ++++++++++++++++++ source/Misc.hpp | 3 + source/MoorDyn2.cpp | 156 +++++++++++++++++++++++++++++++++++++++----- 4 files changed, 307 insertions(+), 18 deletions(-) diff --git a/source/Line.cpp b/source/Line.cpp index f14b743c..e8e8f764 100644 --- a/source/Line.cpp +++ b/source/Line.cpp @@ -118,6 +118,41 @@ 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 +193,44 @@ 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); @@ -613,6 +676,15 @@ Line::initialize() } } + // 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(stiffzs_, stiffxs_, 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 +1091,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 @@ -1095,6 +1168,34 @@ Line::getStateDeriv(InstanceStateVarView drdt) // 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]) { + Tmean[i] = Tmax[i]; + setWorkingCurve(Tmax[i]); + } + } } // calculate the curvatures and normal vectors (just if needed) diff --git a/source/Line.hpp b/source/Line.hpp index 9ee1b5f5..9b8744d5 100644 --- a/source/Line.hpp +++ b/source/Line.hpp @@ -90,8 +90,21 @@ 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 +143,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 +237,8 @@ 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,27 @@ 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..9bc220b6 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -1178,6 +1178,9 @@ 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..78ac2a5b 100644 --- a/source/MoorDyn2.cpp +++ b/source/MoorDyn2.cpp @@ -996,6 +996,58 @@ moordyn::MoorDyn::ReadInFile() } } + // If Syrope working curve is defined, read it + if ((i = findStartOfSection(in_txt, { "SYROPE WORKING CURVES"})) != -1) { + LOGDBG << " Reading Syrope working curve formulas:" << endl; + + // parse until the next header or the end of the file + while ((in_txt[i].find("---") == string::npos) && + (i < (int)in_txt.size())) { + vector entries = moordyn::str::split(in_txt[i], ' '); + // Expecting 4 entries: LineType (string), WCType (string), k1 (float), and k2 (float) + if (entries.size() != 4) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "4 fields are required, but just " << entries.size() + << " are provided" << endl; + return MOORDYN_INVALID_INPUT; + } + + string line_prop_type = entries[0]; + string wc_formula = entries[1]; + double k1 = atof(entries[2].c_str()); + double k2 = atof(entries[3].c_str()); + + // Loop all line properties to find the matching one + bool line_type_found = false; + for (auto& line_prop : LinePropList) { + if (line_prop->type == line_prop_type) { + line_type_found = true; + + if (wc_formula == "LINEAR") line_prop->SyropeWCForm = 1; + else if (wc_formula == "QUADRATIC") line_prop->SyropeWCForm = 2; + else if (wc_formula == "EXP") line_prop->SyropeWCForm = 3; + else { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "Invalid working curve formula: " << wc_formula << endl; + return MOORDYN_INVALID_INPUT; + } + line_prop->k1 = k1; + line_prop->k2 = k2; + } + } + + LOGDBG << " Line Type: " << line_prop_type + << ", WC Type: " << wc_formula + << ", k1: " << k1 + << ", k2: " << k2 << endl; + i++; + } + } + if ((i = findStartOfSection(in_txt, { "ROD DICTIONARY", "ROD TYPES" })) != -1) { LOGDBG << " Reading rod types:" << endl; @@ -1332,6 +1384,44 @@ moordyn::MoorDyn::ReadInFile() } } + if ((i = findStartOfSection(in_txt, { "SYROPE IC"})) != -1) { + LOGDBG << " Reading Syrope line initial conditions:" << endl; + + // parse until the next header or the end of the file + while ((in_txt[i].find("---") == string::npos) && + (i < (int)in_txt.size())) { + vector entries = moordyn::str::split(in_txt[i], ' '); + if (entries.size() < 3) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "3 fields are required, but just " << entries.size() + << " are provided" << endl; + return MOORDYN_INVALID_INPUT; + } + + int line_id = atoi(entries[0].c_str()); + real Tmax0 = atof(entries[1].c_str()); + real Tmean0 = atof(entries[2].c_str()); + + if (line_id < 1 || line_id > (int)LineList.size()) { + LOGERR << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "There is no line with ID " << line_id << "!" << endl; + return MOORDYN_INVALID_INPUT; + } + + auto* line = LineList[line_id - 1]; + if (line->getElasticModel() == 4) { + line->setInitialTmax(Tmax0); + line->setInitialTmean(Tmean0); + } + + i++; + } + } + if ((i = findStartOfSection(in_txt, { "FAILURE" })) != -1) { LOGDBG << " Reading failure conditions:" << endl; // parse until the next header or the end of the file @@ -1922,29 +2012,61 @@ 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 owc = 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; - } - err = read_curve(EA_stuff[0].c_str(), + + err = read_curve(owc.c_str(), &(obj->EA), &(obj->nEApoints), obj->stiffXs, obj->stiffYs); - if (err) - return nullptr; + 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; + } vector BA_stuff = moordyn::str::split(entries[4], '|'); unsigned int BA_N = BA_stuff.size(); From 689a22d554187218907d1322dce3f5bc2ebaf5ff Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 12 Jan 2026 13:00:19 +0000 Subject: [PATCH 02/14] Optimization of the exponential working curve formula --- source/Line.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/Line.cpp b/source/Line.cpp index e8e8f764..56818117 100644 --- a/source/Line.cpp +++ b/source/Line.cpp @@ -142,7 +142,7 @@ Line::setWorkingCurve(real Tmax) 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)); + stiffys_[I] = Tmax * (1.0 - exp(k2 * xi)) / (1.0 - exp(k2)); else { LOGERR << "Line " << lineId << ": Invalid Syrope working curve formula " << syropeWCForm << "." << endl; From 0136b33c138009f8d0db686c8a30feefb6337988 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 12 Jan 2026 13:34:37 +0000 Subject: [PATCH 03/14] Fixed wrong initial mean tension and wrong working when Tmean < Tmax --- source/Line.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/Line.cpp b/source/Line.cpp index 56818117..5f9c20e4 100644 --- a/source/Line.cpp +++ b/source/Line.cpp @@ -669,7 +669,7 @@ 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 @@ -681,7 +681,7 @@ Line::initialize() 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(stiffzs_, stiffxs_, Tmean[i]) * l[i]; // stretch instead of strain + dl_1[i] = interp(stiffys_, stiffzs_, Tmean[i]) * l[i]; // stretch instead of strain } } @@ -1192,7 +1192,7 @@ Line::getStateDeriv(InstanceStateVarView drdt) } if (Tmean[i] > Tmax[i]) { - Tmean[i] = Tmax[i]; + Tmax[i] = Tmean[i]; setWorkingCurve(Tmax[i]); } } From 88daee9f2897294923ac1cdde0c045980da338f6 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 12 Jan 2026 14:28:56 +0000 Subject: [PATCH 04/14] Documentation for Syrope extension --- docs/inputs.rst | 59 +++++++++++++++++++++++++++++++++++++++++++++ docs/references.rst | 4 +++ 2 files changed, 63 insertions(+) diff --git a/docs/inputs.rst b/docs/inputs.rst index 4ac84b86..951d1828 100644 --- a/docs/inputs.rst +++ b/docs/inputs.rst @@ -333,6 +333,65 @@ needs to be input as a value. Example inputs are below: Polyester ... ... EA_s|EA_d BA_s|BA_d <-- Constant dynamic stiffness method with static and dynamic damping Polyester ... ... EA_s|EA_Dc|EA_D_Lm BA_s|BA_d <-- Load dependent dynamic stiffness method with static and dynamic damping +MoorDyn also 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 (hysteresis). 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 above. +To indicate Syrope, the **first entry must start with** ``SYROPE:``, followed by the name (or path) +of a tabulated file defining the **original working curve**. The second and third entries are +``EA_Dc`` and ``EA_D_Lm`` (as in the load-dependent dynamic stiffness method). 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:owc.txt|EA_Dc|EA_D_Lm BA_s|BA_d + +The file ``owc.txt`` follows the same format as the tabulated non-linear stiffness files described +above (three header lines, then columns for strain and tension). + +In addition to the line-type inputs above, users must also provide (i) the working-curve +parameterization and (ii) Syrope initial conditions for each Syrope line. + +The working-curve parameters are provided in a dedicated table. The columns are: + +- LineType – a string matching a Line Dictionary entry that uses the Syrope model +- WCType – the working-curve type (options: ``LINEAR``, ``QUADRATIC``, ``EXP``) +- k1 – the primary coefficient of the working curve (sets the strain offset where mean tension is zero) +- k2 – the secondary coefficient of the working curve (sets the curve shape) + +Example: + +.. code-block:: none + + -------------------------- SYROPE WORKING CURVES ------------------------------ + LineType WCType k1 k2 + (-) (-) (-) (-) + poly LINEAR k1 k2 + +Initial conditions are specified per line. The **Line** ID must match a line defined in the +*Lines* section and must reference a line type that uses the Syrope model. + +- Tmax – the highest mean tension experienced by the line prior to the current time (N) +- Tmean – the current mean tension at initialization (N) + +In general, ``Tmax0 >= Tmean0``. + +Example: + +.. code-block:: none + + -------------------------- SYROPE IC ------------------------------------------ + Line Tmax Tmean + (-) (N) (N) + 1 Tmax0 Tmean0 + Rod Types ^^^^^^^^^ diff --git a/docs/references.rst b/docs/references.rst index a6e24059..4e83bd7a 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -107,6 +107,10 @@ Viscoelastic approach for non-linear rope behavior: `Hall, Matthew, Brian Duong, and Ericka Lozon, “Streamlined Loads Analysis of Floating Wind Turbines With Fiber Rope Mooring Lines.” In ASME 2023 5th International Offshore Wind Technical Conference, V001T01A029. Exeter, UK: American Society of Mechanical Engineers, 2023. `_ +Syrope model for polyester ropes: + + `Falkenberg, Erik, Åhjem, Vidar, and Limin Yang. "Best Practice for Analysis of Polyester Rope Mooring Systems." Paper presented at the Offshore Technology Conference, Houston, Texas, USA, May 2017. `_ + Updated MoorDyn-OpenFOAM Coupling: `Haifei Chen, Tanausú Almeida Medina, and Jose Luis Cercos-Pita, "CFD simulation of multiple moored floating structures using OpenFOAM: An open-access mooring restraints From acd8d0bcce191bf26fee163eb0f1511e545424b5 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Tue, 13 Jan 2026 15:45:15 +0000 Subject: [PATCH 05/14] Test for Syrope model with slow loading --- tests/Mooring/syrope/owc.dat | 32 ++++ tests/Mooring/syrope/slow_loading.txt | 35 ++++ tests/syrope.cpp | 246 ++++++++++++++++++++++++++ 3 files changed, 313 insertions(+) create mode 100644 tests/Mooring/syrope/owc.dat create mode 100644 tests/Mooring/syrope/slow_loading.txt create mode 100644 tests/syrope.cpp 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/slow_loading.txt b/tests/Mooring/syrope/slow_loading.txt new file mode 100644 index 00000000..05aca1c2 --- /dev/null +++ b/tests/Mooring/syrope/slow_loading.txt @@ -0,0 +1,35 @@ +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:owc.dat|1.53e8|23.12 2.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +-------------------------- SYROPE WORKING CURVES ------------------------------------------------- +LineType WCType k1 k2 +(-) (-) (-) (-) +poly LINEAR 1.25e8 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.0137912 0 -5.0 0 0 0 0 +-------------------------- LINES ------------------------------------------------- +Line LineType NodeA NodeB UnstrLen NumSegs Flags/Outputs +(-) (-) (-) (-) (m) (-) (-) +1 poly 2 1 1.0 1 ts +-------------------------- SYROPE IC ------------------------------------------------- +Line Tmax Tmean +(-) (N) (N) +1 1.302e6 4.339e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +1.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 +100.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/syrope.cpp b/tests/syrope.cpp new file mode 100644 index 00000000..9fab8f86 --- /dev/null +++ b/tests/syrope.cpp @@ -0,0 +1,246 @@ +/** @file syrope.cpp + * @brief Test case for Syrope model. + * Compares MoorDyn simulation output against analytical expectations from given original working curve + * and working curves. + */ + +#include "MoorDyn2.h" + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +constexpr double G = 9.80665; +constexpr double MBL = 885.0e3 * G; +constexpr double TMAX_0 = 0.15 * MBL; +constexpr double K1 = 1.25e8; + +// Linear interpolation on monotone-increasing xdata, clamped to endpoints. +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"); + + 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) +{ + if (Tmean < Tmax) { + const double eps_max = interpolate_clamped(Tmax, owc_tensions, owc_strains); + const double eps_min = eps_max - Tmax / K1; + return Tmax * (strain - eps_min) / (eps_max - eps_min); + } + return interpolate_clamped(strain, owc_strains, owc_tensions); +} + +bool read_three_whitespace_columns_skip_header(const std::string& path, + int header_lines, + std::vector& c0, + std::vector& c1, + std::vector& c2) +{ + std::ifstream in(path); + if (!in.is_open()) return false; + + std::string line; + for (int i = 0; i < header_lines; ++i) { + if (!std::getline(in, line)) return false; + } + + while (std::getline(in, line)) { + std::istringstream iss(line); + double a = 0.0, b = 0.0, d = 0.0; + if (!(iss >> a >> b >> d)) continue; + c0.push_back(a); + c1.push_back(b); + c2.push_back(d); + } + return true; +} + +bool read_two_whitespace_columns(const std::string& path, + std::vector& c0, + std::vector& c1) +{ + std::ifstream in(path); + if (!in.is_open()) return false; + + std::string line; + while (std::getline(in, line)) { + std::istringstream iss(line); + double a = 0.0, b = 0.0; + if (!(iss >> a >> b)) continue; + c0.push_back(a); + c1.push_back(b); + } + return true; +} + +} // namespace + +TEST_CASE("Syrope: simulation runs and working-curve postprocessing is consistent for slow loading") +{ + // --- Run MoorDyn case --- + MoorDyn system = MoorDyn_Create("Mooring/syrope/slow_loading.txt"); + REQUIRE(system); + + unsigned int ndof = 0; + REQUIRE(MoorDyn_NCoupledDOF(system, &ndof) == MOORDYN_SUCCESS); + REQUIRE(ndof == 3); + + double r[3] = {0.0, 0.0, 0.0}; + double dr[3] = {0.0, 0.0, 0.0}; + + MoorDynPoint anchor = MoorDyn_GetPoint(system, 1); + MoorDynPoint fairlead = MoorDyn_GetPoint(system, 2); + MoorDynLine line = MoorDyn_GetLine(system, 1); + REQUIRE(anchor); + REQUIRE(fairlead); + REQUIRE(line); + + double l0 = 0.0; + REQUIRE(MoorDyn_GetLineUnstretchedLength(line, &l0) == MOORDYN_SUCCESS); + + REQUIRE(MoorDyn_GetPointPos(fairlead, r) == MOORDYN_SUCCESS); + std::fill(dr, dr + 3, 0.0); + REQUIRE(MoorDyn_Init(system, r, dr) == MOORDYN_SUCCESS); + + const double tdur = 3600.0 * 3.0; // 3 hours for one loading segment + const double Tdur = tdur * 6.0; // six loading segments + double dt = 10.0; + const unsigned int nsteps = static_cast(Tdur / dt); + + Eigen::VectorXd times = Eigen::VectorXd::LinSpaced(static_cast(nsteps) + 1, 0.0, Tdur); + Eigen::VectorXd strains = Eigen::VectorXd::Zero(static_cast(nsteps) + 1); + + const double eps_0 = 1.37912e-02; + for (unsigned int i = 0; i <= nsteps; ++i) { + const double tt = times[static_cast(i)]; + double eps = eps_0; + + if (tt < tdur) eps = eps_0; + else if (tt < 2.0 * tdur) eps = eps_0 + 2.0 * eps_0 * (tt - tdur) / tdur; + else if (tt < 3.0 * tdur) eps = eps_0 + 2.0 * eps_0; + else if (tt < 4.0 * tdur) eps = eps_0 + 2.0 * eps_0 - 1.0 * eps_0 * (tt - 3.0 * tdur) / tdur; + else if (tt < 5.0 * tdur) eps = eps_0 + 1.0 * eps_0; + else eps = eps_0 + 1.0 * eps_0 + 1.5 * eps_0 * (tt - 5.0 * tdur) / tdur; + + strains[static_cast(i)] = eps; + } + + double t = 0.0; + double flat[3] = {0.0, 0.0, 0.0}; + + for (unsigned int i = 0; i < nsteps; ++i) { + r[0] = 1.0 + strains[static_cast(i)]; + + double t_in = times[static_cast(i)]; + double dt_in = dt; + + INFO("Step i=" << i << " t=" << t_in << " dt=" << dt_in << " r[0]=" << r[0]); + REQUIRE(MoorDyn_Step(system, r, dr, flat, &t_in, &dt_in) == MOORDYN_SUCCESS); + + dt = dt_in; + } + + REQUIRE(MoorDyn_Close(system) == MOORDYN_SUCCESS); + + // --- Read MoorDyn output --- + std::vector timedata; + std::vector tensiondata; // Mean tension instead of instant tension, + // MoorDyn_GetLineNodeTen() returns instant tension + // (sum of the axial stiffness plus the internal damping) + std::vector straindata; + REQUIRE(read_three_whitespace_columns_skip_header("Mooring/syrope/slow_loading_Line1.out", + /*header_lines=*/2, + timedata, tensiondata, straindata)); + + REQUIRE(timedata.size() >= 2); + REQUIRE(tensiondata.size() == timedata.size()); + REQUIRE(straindata.size() == timedata.size()); + + const Eigen::Index n_csv = static_cast(timedata.size()); + const Eigen::VectorXd csv_time = + Eigen::Map(timedata.data(), n_csv); + const Eigen::VectorXd csv_tension = + Eigen::Map(tensiondata.data(), n_csv); + const Eigen::VectorXd csv_strain = + Eigen::Map(straindata.data(), n_csv); + + const Eigen::Index n_use = std::min(n_csv, strains.size()); + REQUIRE(n_use >= 2); + + // --- Read OWC (strain, tension) --- + std::vector owc_strain_data, owc_tension_data; + REQUIRE(read_two_whitespace_columns("Mooring/syrope/owc.dat", owc_strain_data, owc_tension_data)); + REQUIRE(owc_strain_data.size() >= 2); + REQUIRE(owc_strain_data.size() == owc_tension_data.size()); + + const Eigen::VectorXd owc_strains = + Eigen::Map(owc_strain_data.data(), + static_cast(owc_strain_data.size())); + const Eigen::VectorXd owc_tensions = + Eigen::Map(owc_tension_data.data(), + static_cast(owc_tension_data.size())); + + // --- Post-process: Tmax and analytical tension --- + Eigen::VectorXd Tmax = Eigen::VectorXd::Zero(n_use); + Tmax[0] = TMAX_0; + for (Eigen::Index i = 1; i < n_use; ++i) { + Tmax[i] = std::max(Tmax[i - 1], csv_tension[i]); + } + + Eigen::VectorXd anal_tension = Eigen::VectorXd::Zero(n_use); + for (Eigen::Index i = 0; i < n_use; ++i) { + anal_tension[i] = find_mean_tension(strains[i], csv_tension[i], Tmax[i], owc_strains, owc_tensions); + } + + // Relative L2 error + const Eigen::VectorXd abs_err = (csv_tension.head(n_use) - anal_tension).cwiseAbs(); + const double denom = std::max(anal_tension.squaredNorm(), 1e-30); + const double rel_l2 = std::sqrt(abs_err.squaredNorm() / denom); + + REQUIRE(std::isfinite(rel_l2)); + + // For slow loading, expect <1% relative L2 error + REQUIRE(rel_l2 < 0.01); +} \ No newline at end of file From a89e723fcb566e15cdcf61e2b6ea6bdd8be0a5a3 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Thu, 15 Jan 2026 15:59:05 +0000 Subject: [PATCH 06/14] Three working curve formulas in one test --- tests/syrope.cpp | 653 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 498 insertions(+), 155 deletions(-) diff --git a/tests/syrope.cpp b/tests/syrope.cpp index 9fab8f86..0048452b 100644 --- a/tests/syrope.cpp +++ b/tests/syrope.cpp @@ -1,9 +1,3 @@ -/** @file syrope.cpp - * @brief Test case for Syrope model. - * Compares MoorDyn simulation output against analytical expectations from given original working curve - * and working curves. - */ - #include "MoorDyn2.h" #include @@ -14,33 +8,116 @@ #include #include #include +#include +#include #include #include #include +#include #include namespace { -constexpr double G = 9.80665; -constexpr double MBL = 885.0e3 * G; -constexpr double TMAX_0 = 0.15 * MBL; -constexpr double K1 = 1.25e8; +// ----------------------------- +// Constants +// ----------------------------- +constexpr double kG = 9.80665; +constexpr double kMBL = 885.0e3 * kG; + +// initial conditions +constexpr double kTmax0 = 0.15 * kMBL; +constexpr double kTmean0 = 0.05 * kMBL; // currently unused + +// dynamic coefficients (currently unused) +constexpr double kAlpha = 17.667 * kMBL; +constexpr double kBeta = 0.2313 * 100.0; + +enum class WorkingCurveForm { Linear, Quadratic, Exponential }; + +// working curve coefficients +constexpr double kLinearK1 = 1.25e8; +constexpr double kLinearK2 = 0.0; +constexpr double kQuadraticK1 = 0.25; +constexpr double kQuadraticK2 = 0.80; +constexpr double kExpK1 = 0.20; +constexpr double kExpK2 = 1.10; + +// Tolerance for regression (adjust as appropriate once you have baselines) +constexpr double kL2Tol = 0.05; + +// ----------------------------- +// working curve cases +// ----------------------------- +struct WcCase +{ + std::string name; + std::string input_file; + WorkingCurveForm form; + double eps_0; +}; + +// ----------------------------- +// MoorDyn error handling +// ----------------------------- +inline void check_md(const int err, const std::string& msg) +{ + REQUIRE(err == MOORDYN_SUCCESS); + if (err != MOORDYN_SUCCESS) { + throw std::runtime_error(msg); + } +} + +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; -// Linear interpolation on monotone-increasing xdata, clamped to endpoints. + MoorDyn sys = nullptr; +}; + +// ----------------------------- +// Linear interpolation (clamped) +// xdata must be monotone increasing +// ----------------------------- double interpolate_clamped(double x, const Eigen::VectorXd& xdata, const Eigen::VectorXd& ydata) { - if (xdata.size() != ydata.size()) + if (xdata.size() != ydata.size()) { throw std::invalid_argument("interpolate_clamped: size mismatch"); - if (xdata.size() < 2) + } + 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; @@ -51,196 +128,462 @@ double interpolate_clamped(double x, const double y0 = ydata[lo], y1 = ydata[hi]; const double dx = x1 - x0; - if (std::abs(dx) <= std::numeric_limits::epsilon()) + if (std::abs(dx) <= std::numeric_limits::epsilon()) { return 0.5 * (y0 + y1); + } const double t = (x - x0) / dx; return y0 + t * (y1 - y0); } +// ----------------------------- +// Mean tension from working curve +// ----------------------------- double find_mean_tension(double strain, double Tmean, double Tmax, const Eigen::VectorXd& owc_strains, - const Eigen::VectorXd& owc_tensions) + const Eigen::VectorXd& owc_tensions, + WorkingCurveForm wc_form) { - if (Tmean < Tmax) { - const double eps_max = interpolate_clamped(Tmax, owc_tensions, owc_strains); - const double eps_min = eps_max - Tmax / K1; - return Tmax * (strain - eps_min) / (eps_max - eps_min); + 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"); + } + + // If not unloading/reloading, mean tension follows the original working curve + if (Tmean >= Tmax) { + return interpolate_clamped(strain, owc_strains, owc_tensions); + } + + // Unloading/reloading: build working curve between eps_min and eps_max + 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); + + // Linear WC special rule: if k1 >= 1.0 treat as slope (dimensional stiffness) + 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); + + switch (wc_form) { + case WorkingCurveForm::Linear: + return Tmax * xi; + case WorkingCurveForm::Quadratic: + return Tmax * xi * (k2 * xi + (1.0 - k2)); + case WorkingCurveForm::Exponential: { + const double den = std::exp(k2) - 1.0; + if (std::abs(den) < 1e-14) return Tmax * xi; // limit -> linear + return Tmax * (std::exp(k2 * xi) - 1.0) / den; + } + default: + throw std::invalid_argument("find_mean_tension: unknown WorkingCurveForm"); } - return interpolate_clamped(strain, owc_strains, owc_tensions); } -bool read_three_whitespace_columns_skip_header(const std::string& path, - int header_lines, - std::vector& c0, - std::vector& c1, - std::vector& c2) +// ----------------------------- +// Read OWC table: 2 header lines then two numeric columns +// ----------------------------- +std::pair +read_strain_tension_table(const std::string& path) { std::ifstream in(path); - if (!in.is_open()) return false; + if (!in) { + throw std::runtime_error("Cannot open file: " + path); + } std::string line; - for (int i = 0; i < header_lines; ++i) { - if (!std::getline(in, line)) return false; - } + 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 a = 0.0, b = 0.0, d = 0.0; - if (!(iss >> a >> b >> d)) continue; - c0.push_back(a); - c1.push_back(b); - c2.push_back(d); + 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); } - return true; + + 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}; } -bool read_two_whitespace_columns(const std::string& path, - std::vector& c0, - std::vector& c1) +// ----------------------------- +// Read SegTe column from MoorDyn output file +// Assumes: +// 1) header tokens line +// 2) units tokens line, with "(N)" under SegNTe +// ----------------------------- +Eigen::VectorXd load_seg_te_column(const std::string& path) { std::ifstream in(path); - if (!in.is_open()) return false; + if (!in) { + throw std::runtime_error("Cannot open " + path); + } std::string line; + + // Header 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); + } + + // Find first SegTe + 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); + } + + // Units line + 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); + } + + // Data + std::vector vals; + vals.reserve(2048); + + std::size_t lineno = 2; while (std::getline(in, line)) { - std::istringstream iss(line); - double a = 0.0, b = 0.0; - if (!(iss >> a >> b)) continue; - c0.push_back(a); - c1.push_back(b); + ++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); } - return true; + + Eigen::VectorXd out(static_cast(vals.size())); + for (Eigen::Index i = 0; i < out.size(); ++i) { + out[i] = vals[static_cast(i)]; + } + return out; } -} // namespace +// ----------------------------- +// JONSWAP PSD (as in your original code) +// ----------------------------- +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)); +} + +// ----------------------------- +// Slow strain profile +// ----------------------------- +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; +} + +// ----------------------------- +// Surface elevation from JONSWAP spectrum +// Deterministic random phases via RNG seed +// ----------------------------- +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; +} -TEST_CASE("Syrope: simulation runs and working-curve postprocessing is consistent for slow loading") +// Helper: "_Line1.txt" without +std::string line1_output_from_input(const std::string& input_path) { - // --- Run MoorDyn case --- - MoorDyn system = MoorDyn_Create("Mooring/syrope/slow_loading.txt"); - REQUIRE(system); + 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"; +} + +// ----------------------------- +// Run a single case and return relative L2 error +// ----------------------------- +double run_case(const WcCase& c, + const Eigen::VectorXd& owc_strains, + const Eigen::VectorXd& owc_tensions) +{ + MoorDynRAII md(c.input_file); unsigned int ndof = 0; - REQUIRE(MoorDyn_NCoupledDOF(system, &ndof) == MOORDYN_SUCCESS); - REQUIRE(ndof == 3); + check_md(MoorDyn_NCoupledDOF(md.sys, &ndof), "MoorDyn_NCoupledDOF failed for: " + c.input_file); + INFO("ndof=" << ndof); + + auto fairlead = MoorDyn_GetPoint(md.sys, 2); + auto line = MoorDyn_GetLine(md.sys, 1); + REQUIRE(fairlead != nullptr); + REQUIRE(line != 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); + + // time settings + const double seg_dur = 3.0 * 3600.0; + const double total_dur = 6.0 * seg_dur; + const double dt0 = 0.1; + const unsigned int nsteps = static_cast(total_dur / dt0); + + const Eigen::VectorXd times = + Eigen::VectorXd::LinSpaced(static_cast(nsteps) + 1, 0.0, total_dur); + + // deterministic excitation + constexpr unsigned int n_omega_edges = 300; + std::mt19937 rng_wf(12345); + std::mt19937 rng_lf(54321); + + // WF + 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); + + // LF + 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); + + // prescribed strain and x(t) + 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); + } - MoorDynPoint anchor = MoorDyn_GetPoint(system, 1); - MoorDynPoint fairlead = MoorDyn_GetPoint(system, 2); - MoorDynLine line = MoorDyn_GetLine(system, 1); - REQUIRE(anchor); - REQUIRE(fairlead); - REQUIRE(line); + const Eigen::VectorXd strain = strain_slow + scale_WF * eta_WF + scale_LF * eta_LF; + const Eigen::VectorXd x = x0 * (Eigen::VectorXd::Ones(strain.size()) + strain); - double l0 = 0.0; - REQUIRE(MoorDyn_GetLineUnstretchedLength(line, &l0) == MOORDYN_SUCCESS); + // init moordyn + r[0] = x[0]; + dr[0] = (x[1] - x[0]) / dt0; + dr[1] = 0.0; + dr[2] = 0.0; - REQUIRE(MoorDyn_GetPointPos(fairlead, r) == MOORDYN_SUCCESS); - std::fill(dr, dr + 3, 0.0); - REQUIRE(MoorDyn_Init(system, r, dr) == MOORDYN_SUCCESS); + check_md(MoorDyn_Init(md.sys, r, dr), "MoorDyn_Init failed for: " + c.input_file); - const double tdur = 3600.0 * 3.0; // 3 hours for one loading segment - const double Tdur = tdur * 6.0; // six loading segments - double dt = 10.0; - const unsigned int nsteps = static_cast(Tdur / dt); + // stepping (enforce fixed dt) + double t = 0.0; + double f[3] = {0.0, 0.0, 0.0}; - Eigen::VectorXd times = Eigen::VectorXd::LinSpaced(static_cast(nsteps) + 1, 0.0, Tdur); - Eigen::VectorXd strains = Eigen::VectorXd::Zero(static_cast(nsteps) + 1); + 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; - const double eps_0 = 1.37912e-02; - for (unsigned int i = 0; i <= nsteps; ++i) { - const double tt = times[static_cast(i)]; - double eps = eps_0; + double t_in = t; + double dt_in = dt0; - if (tt < tdur) eps = eps_0; - else if (tt < 2.0 * tdur) eps = eps_0 + 2.0 * eps_0 * (tt - tdur) / tdur; - else if (tt < 3.0 * tdur) eps = eps_0 + 2.0 * eps_0; - else if (tt < 4.0 * tdur) eps = eps_0 + 2.0 * eps_0 - 1.0 * eps_0 * (tt - 3.0 * tdur) / tdur; - else if (tt < 5.0 * tdur) eps = eps_0 + 1.0 * eps_0; - else eps = eps_0 + 1.0 * eps_0 + 1.5 * eps_0 * (tt - 5.0 * tdur) / tdur; + check_md(MoorDyn_Step(md.sys, r, dr, f, &t_in, &dt_in), + "MoorDyn_Step failed for: " + c.input_file); - strains[static_cast(i)] = eps; + REQUIRE(std::abs(dt_in - dt0) <= 1e-12); + t = t_in; } - double t = 0.0; - double flat[3] = {0.0, 0.0, 0.0}; + // close to flush outputs + check_md(MoorDyn_Close(md.sys), "MoorDyn_Close failed for: " + c.input_file); + md.sys = nullptr; - for (unsigned int i = 0; i < nsteps; ++i) { - r[0] = 1.0 + strains[static_cast(i)]; - - double t_in = times[static_cast(i)]; - double dt_in = dt; - - INFO("Step i=" << i << " t=" << t_in << " dt=" << dt_in << " r[0]=" << r[0]); - REQUIRE(MoorDyn_Step(system, r, dr, flat, &t_in, &dt_in) == MOORDYN_SUCCESS); - - dt = dt_in; - } - - REQUIRE(MoorDyn_Close(system) == MOORDYN_SUCCESS); - - // --- Read MoorDyn output --- - std::vector timedata; - std::vector tensiondata; // Mean tension instead of instant tension, - // MoorDyn_GetLineNodeTen() returns instant tension - // (sum of the axial stiffness plus the internal damping) - std::vector straindata; - REQUIRE(read_three_whitespace_columns_skip_header("Mooring/syrope/slow_loading_Line1.out", - /*header_lines=*/2, - timedata, tensiondata, straindata)); - - REQUIRE(timedata.size() >= 2); - REQUIRE(tensiondata.size() == timedata.size()); - REQUIRE(straindata.size() == timedata.size()); - - const Eigen::Index n_csv = static_cast(timedata.size()); - const Eigen::VectorXd csv_time = - Eigen::Map(timedata.data(), n_csv); - const Eigen::VectorXd csv_tension = - Eigen::Map(tensiondata.data(), n_csv); - const Eigen::VectorXd csv_strain = - Eigen::Map(straindata.data(), n_csv); - - const Eigen::Index n_use = std::min(n_csv, strains.size()); - REQUIRE(n_use >= 2); - - // --- Read OWC (strain, tension) --- - std::vector owc_strain_data, owc_tension_data; - REQUIRE(read_two_whitespace_columns("Mooring/syrope/owc.dat", owc_strain_data, owc_tension_data)); - REQUIRE(owc_strain_data.size() >= 2); - REQUIRE(owc_strain_data.size() == owc_tension_data.size()); - - const Eigen::VectorXd owc_strains = - Eigen::Map(owc_strain_data.data(), - static_cast(owc_strain_data.size())); - const Eigen::VectorXd owc_tensions = - Eigen::Map(owc_tension_data.data(), - static_cast(owc_tension_data.size())); - - // --- Post-process: Tmax and analytical tension --- - Eigen::VectorXd Tmax = Eigen::VectorXd::Zero(n_use); - Tmax[0] = TMAX_0; - for (Eigen::Index i = 1; i < n_use; ++i) { - Tmax[i] = std::max(Tmax[i - 1], csv_tension[i]); - } - - Eigen::VectorXd anal_tension = Eigen::VectorXd::Zero(n_use); - for (Eigen::Index i = 0; i < n_use; ++i) { - anal_tension[i] = find_mean_tension(strains[i], csv_tension[i], Tmax[i], owc_strains, owc_tensions); - } - - // Relative L2 error - const Eigen::VectorXd abs_err = (csv_tension.head(n_use) - anal_tension).cwiseAbs(); - const double denom = std::max(anal_tension.squaredNorm(), 1e-30); - const double rel_l2 = std::sqrt(abs_err.squaredNorm() / denom); - - REQUIRE(std::isfinite(rel_l2)); - - // For slow loading, expect <1% relative L2 error - REQUIRE(rel_l2 < 0.01); + // read output mean tension + const std::string out_path = line1_output_from_input(c.input_file); + INFO("output=" << out_path); + + Eigen::VectorXd tension_mean_output; + REQUIRE_NOTHROW(tension_mean_output = load_seg_te_column(out_path)); + + const Eigen::Index n = std::min(tension_mean_output.size(), strain.size()); + REQUIRE(n >= 2); + + // preceding highest mean tension + 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], tension_mean_output[i]); + } + + // analytical mean tension from WC + OWC + Eigen::VectorXd tension_analytical(n); + for (Eigen::Index i = 0; i < n; ++i) { + tension_analytical[i] = + find_mean_tension(strain[i], + tension_mean_output[i], + Tmax_mean[i], + owc_strains, + owc_tensions, + c.form); + } + + const double denom = tension_analytical.squaredNorm(); + REQUIRE(denom > 0.0); + + const double l2_rel = + std::sqrt((tension_analytical - tension_mean_output.head(n)).squaredNorm() / denom); + + return l2_rel; +} + +} // namespace + +TEST_CASE("Syrope working curve regression: mean tension matches analytical WC", "[syrope][working-curve]") +{ + // shared OWC table for all cases + const auto [owc_strains, owc_tensions] = read_strain_tension_table("Mooring/syrope/owc.dat"); + + const std::vector cases = { + {"Linear", "Mooring/syrope/linear_wc_input.txt", WorkingCurveForm::Linear, 1.37912e-02}, + {"Quadratic", "Mooring/syrope/quadratic_wc_input.txt", WorkingCurveForm::Quadratic, 1.41587e-02}, + {"Exponential", "Mooring/syrope/exponential_wc_input.txt", WorkingCurveForm::Exponential, 1.18597e-02}, + }; + + for (const auto& c : cases) { + DYNAMIC_SECTION("Case: " << c.name) { + const double l2_rel = run_case(c, owc_strains, owc_tensions); + INFO("Relative L2 error = " << l2_rel); + REQUIRE(l2_rel < kL2Tol); + } + } } \ No newline at end of file From 8e04b7b256e94cb9fbd59f3eee74293243206a47 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Fri, 16 Jan 2026 08:11:30 +0000 Subject: [PATCH 07/14] Formatting --- source/Line.cpp | 81 ++-- source/Line.hpp | 9 +- source/Misc.hpp | 7 +- source/MoorDyn2.cpp | 102 +++-- tests/syrope.cpp | 1025 +++++++++++++++++++++++-------------------- 5 files changed, 671 insertions(+), 553 deletions(-) diff --git a/source/Line.cpp b/source/Line.cpp index 5f9c20e4..27b424b4 100644 --- a/source/Line.cpp +++ b/source/Line.cpp @@ -123,33 +123,42 @@ 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; + 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 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 + 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; + 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"); + 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]); + stiffzs_[I] = stiffxs_[I] - + 1.0 / vbeta * log(1.0 + vbeta / alphaMBL * stiffys_[I]); } } @@ -210,14 +219,21 @@ Line::setup(int number_in, // 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])); + 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"); + 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"); } } @@ -230,7 +246,6 @@ Line::setup(int number_in, 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); @@ -676,12 +691,14 @@ Line::initialize() } } - // If using the Syrope model, initalize the deltaL_1 to slow stretch based on the mean tension and the active working curve + // 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 + dl_1[i] = interp(stiffys_, stiffzs_, Tmean[i]) * + l[i]; // stretch instead of strain } } @@ -1091,7 +1108,7 @@ Line::getStateDeriv(InstanceStateVarView drdt) BA = getNonlinearBA(ldstr[i], l[i]); Td[i] = BA * ldstr[i] / l[i] * qs[i]; - } else if (ElasticMod == ELASTIC_VISCO_CTE || + } else if (ElasticMod == ELASTIC_VISCO_CTE || ElasticMod == ELASTIC_VISCO_MEAN) { // viscoelastic model from // https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018 @@ -1167,26 +1184,34 @@ 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) { + } 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); + 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 + } 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); + 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]; } diff --git a/source/Line.hpp b/source/Line.hpp index 9b8744d5..4691374c 100644 --- a/source/Line.hpp +++ b/source/Line.hpp @@ -105,7 +105,6 @@ class DECLDIR Line final SYROPE_EXP_WC = 3, } syrope_wc_formula; - /** @brief Get the non-linear stiffness. This is interpolated from a * curve provided in the input file. * @param l_stretched The actual length of the segment @@ -237,7 +236,8 @@ 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) + /// 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; @@ -532,10 +532,7 @@ 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); - } + inline void setInitialTmax(moordyn::real Tmax0) { Tmax.assign(N, Tmax0); } /** @brief Set the initial mean tension values * @param Tmean0 Initial Tmean value diff --git a/source/Misc.hpp b/source/Misc.hpp index 9bc220b6..d93013fd 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -1178,9 +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 + 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 78ac2a5b..c930cdab 100644 --- a/source/MoorDyn2.cpp +++ b/source/MoorDyn2.cpp @@ -997,20 +997,21 @@ moordyn::MoorDyn::ReadInFile() } // If Syrope working curve is defined, read it - if ((i = findStartOfSection(in_txt, { "SYROPE WORKING CURVES"})) != -1) { + if ((i = findStartOfSection(in_txt, { "SYROPE WORKING CURVES" })) != -1) { LOGDBG << " Reading Syrope working curve formulas:" << endl; // parse until the next header or the end of the file while ((in_txt[i].find("---") == string::npos) && (i < (int)in_txt.size())) { vector entries = moordyn::str::split(in_txt[i], ' '); - // Expecting 4 entries: LineType (string), WCType (string), k1 (float), and k2 (float) + // Expecting 4 entries: LineType (string), WCType (string), k1 + // (float), and k2 (float) if (entries.size() != 4) { LOGERR << "Error in " << _filepath << " at line " << i + 1 - << ":" << endl - << "'" << in_txt[i] << "'" << endl - << "4 fields are required, but just " << entries.size() - << " are provided" << endl; + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "4 fields are required, but just " << entries.size() + << " are provided" << endl; return MOORDYN_INVALID_INPUT; } @@ -1024,15 +1025,20 @@ moordyn::MoorDyn::ReadInFile() for (auto& line_prop : LinePropList) { if (line_prop->type == line_prop_type) { line_type_found = true; - - if (wc_formula == "LINEAR") line_prop->SyropeWCForm = 1; - else if (wc_formula == "QUADRATIC") line_prop->SyropeWCForm = 2; - else if (wc_formula == "EXP") line_prop->SyropeWCForm = 3; + + if (wc_formula == "LINEAR") + line_prop->SyropeWCForm = 1; + else if (wc_formula == "QUADRATIC") + line_prop->SyropeWCForm = 2; + else if (wc_formula == "EXP") + line_prop->SyropeWCForm = 3; else { - LOGERR << "Error in " << _filepath << " at line " << i + 1 - << ":" << endl - << "'" << in_txt[i] << "'" << endl - << "Invalid working curve formula: " << wc_formula << endl; + LOGERR + << "Error in " << _filepath << " at line " << i + 1 + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "Invalid working curve formula: " << wc_formula + << endl; return MOORDYN_INVALID_INPUT; } line_prop->k1 = k1; @@ -1041,9 +1047,8 @@ moordyn::MoorDyn::ReadInFile() } LOGDBG << " Line Type: " << line_prop_type - << ", WC Type: " << wc_formula - << ", k1: " << k1 - << ", k2: " << k2 << endl; + << ", WC Type: " << wc_formula << ", k1: " << k1 + << ", k2: " << k2 << endl; i++; } } @@ -1384,7 +1389,7 @@ moordyn::MoorDyn::ReadInFile() } } - if ((i = findStartOfSection(in_txt, { "SYROPE IC"})) != -1) { + if ((i = findStartOfSection(in_txt, { "SYROPE IC" })) != -1) { LOGDBG << " Reading Syrope line initial conditions:" << endl; // parse until the next header or the end of the file @@ -1393,22 +1398,22 @@ moordyn::MoorDyn::ReadInFile() vector entries = moordyn::str::split(in_txt[i], ' '); if (entries.size() < 3) { LOGERR << "Error in " << _filepath << " at line " << i + 1 - << ":" << endl - << "'" << in_txt[i] << "'" << endl - << "3 fields are required, but just " << entries.size() - << " are provided" << endl; + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "3 fields are required, but just " << entries.size() + << " are provided" << endl; return MOORDYN_INVALID_INPUT; } int line_id = atoi(entries[0].c_str()); real Tmax0 = atof(entries[1].c_str()); real Tmean0 = atof(entries[2].c_str()); - + if (line_id < 1 || line_id > (int)LineList.size()) { LOGERR << "Error in " << _filepath << " at line " << i + 1 - << ":" << endl - << "'" << in_txt[i] << "'" << endl - << "There is no line with ID " << line_id << "!" << endl; + << ":" << endl + << "'" << in_txt[i] << "'" << endl + << "There is no line with ID " << line_id << "!" << endl; return MOORDYN_INVALID_INPUT; } @@ -2019,14 +2024,15 @@ moordyn::MoorDyn::readLineProps(string inputText, int lineNum) if (is_syrope) { obj->ElasticMod = 4; // Syrope model - LOGDBG << "Line type '" << obj->type - << "' is using the Syrope model for axial tension-strain behavior." - << endl; - + 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; + LOGERR << "A line type EA entry with Syrope model must have 3 " + "(bar-separated) values." + << endl; return nullptr; } @@ -2035,35 +2041,35 @@ moordyn::MoorDyn::readLineProps(string inputText, int lineNum) obj->vbeta = atof(EA_stuff[2].c_str()); err = read_curve(owc.c_str(), - &(obj->EA), - &(obj->nEApoints), - obj->stiffXs, - obj->stiffYs); + &(obj->EA), + &(obj->nEApoints), + obj->stiffXs, + obj->stiffYs); if (err) return nullptr; - } - else { // Not Syrope, other models + } 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->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 + 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; + 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); + &(obj->EA), + &(obj->nEApoints), + obj->stiffXs, + obj->stiffYs); if (err) return nullptr; } diff --git a/tests/syrope.cpp b/tests/syrope.cpp index 0048452b..4d938da7 100644 --- a/tests/syrope.cpp +++ b/tests/syrope.cpp @@ -21,26 +21,31 @@ namespace { // ----------------------------- // Constants // ----------------------------- -constexpr double kG = 9.80665; +constexpr double kG = 9.80665; constexpr double kMBL = 885.0e3 * kG; // initial conditions -constexpr double kTmax0 = 0.15 * kMBL; +constexpr double kTmax0 = 0.15 * kMBL; constexpr double kTmean0 = 0.05 * kMBL; // currently unused // dynamic coefficients (currently unused) constexpr double kAlpha = 17.667 * kMBL; -constexpr double kBeta = 0.2313 * 100.0; +constexpr double kBeta = 0.2313 * 100.0; -enum class WorkingCurveForm { Linear, Quadratic, Exponential }; +enum class WorkingCurveForm +{ + Linear, + Quadratic, + Exponential +}; // working curve coefficients -constexpr double kLinearK1 = 1.25e8; -constexpr double kLinearK2 = 0.0; +constexpr double kLinearK1 = 1.25e8; +constexpr double kLinearK2 = 0.0; constexpr double kQuadraticK1 = 0.25; constexpr double kQuadraticK2 = 0.80; -constexpr double kExpK1 = 0.20; -constexpr double kExpK2 = 1.10; +constexpr double kExpK1 = 0.20; +constexpr double kExpK2 = 1.10; // Tolerance for regression (adjust as appropriate once you have baselines) constexpr double kL2Tol = 0.05; @@ -50,149 +55,172 @@ constexpr double kL2Tol = 0.05; // ----------------------------- struct WcCase { - std::string name; - std::string input_file; - WorkingCurveForm form; - double eps_0; + std::string name; + std::string input_file; + WorkingCurveForm form; + double eps_0; }; // ----------------------------- // MoorDyn error handling // ----------------------------- -inline void check_md(const int err, const std::string& msg) +inline void +check_md(const int err, const std::string& msg) { - REQUIRE(err == MOORDYN_SUCCESS); - if (err != MOORDYN_SUCCESS) { - throw std::runtime_error(msg); - } + REQUIRE(err == MOORDYN_SUCCESS); + if (err != MOORDYN_SUCCESS) { + throw std::runtime_error(msg); + } } 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; + 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; }; // ----------------------------- // Linear interpolation (clamped) // xdata must be monotone increasing // ----------------------------- -double interpolate_clamped(double x, - const Eigen::VectorXd& xdata, - const Eigen::VectorXd& ydata) +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); + 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); } // ----------------------------- // Mean tension from working curve // ----------------------------- -double find_mean_tension(double strain, - double Tmean, - double Tmax, - const Eigen::VectorXd& owc_strains, - const Eigen::VectorXd& owc_tensions, - WorkingCurveForm wc_form) +double +find_mean_tension(double strain, + double Tmean, + double Tmax, + const Eigen::VectorXd& owc_strains, + const Eigen::VectorXd& owc_tensions, + WorkingCurveForm wc_form) { - 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"); - } - - // If not unloading/reloading, mean tension follows the original working curve - if (Tmean >= Tmax) { - return interpolate_clamped(strain, owc_strains, owc_tensions); - } - - // Unloading/reloading: build working curve between eps_min and eps_max - 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); - - // Linear WC special rule: if k1 >= 1.0 treat as slope (dimensional stiffness) - 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); - - switch (wc_form) { - case WorkingCurveForm::Linear: - return Tmax * xi; - case WorkingCurveForm::Quadratic: - return Tmax * xi * (k2 * xi + (1.0 - k2)); - case WorkingCurveForm::Exponential: { - const double den = std::exp(k2) - 1.0; - if (std::abs(den) < 1e-14) return Tmax * xi; // limit -> linear - return Tmax * (std::exp(k2 * xi) - 1.0) / den; - } - default: - throw std::invalid_argument("find_mean_tension: unknown WorkingCurveForm"); - } + 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"); + } + + // If not unloading/reloading, mean tension follows the original working + // curve + if (Tmean >= Tmax) { + return interpolate_clamped(strain, owc_strains, owc_tensions); + } + + // Unloading/reloading: build working curve between eps_min and eps_max + 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); + + // Linear WC special rule: if k1 >= 1.0 treat as slope (dimensional + // stiffness) + 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); + + switch (wc_form) { + case WorkingCurveForm::Linear: + return Tmax * xi; + case WorkingCurveForm::Quadratic: + return Tmax * xi * (k2 * xi + (1.0 - k2)); + case WorkingCurveForm::Exponential: { + const double den = std::exp(k2) - 1.0; + if (std::abs(den) < 1e-14) + return Tmax * xi; // limit -> linear + return Tmax * (std::exp(k2 * xi) - 1.0) / den; + } + default: + throw std::invalid_argument( + "find_mean_tension: unknown WorkingCurveForm"); + } } // ----------------------------- @@ -201,51 +229,57 @@ double find_mean_tension(double strain, 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}; + 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 }; } // ----------------------------- @@ -254,336 +288,391 @@ read_strain_tension_table(const std::string& path) // 1) header tokens line // 2) units tokens line, with "(N)" under SegNTe // ----------------------------- -Eigen::VectorXd load_seg_te_column(const std::string& path) +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; - - // Header 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); - } - - // Find first SegTe - 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); - } - - // Units line - 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); - } - - // Data - 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; + std::ifstream in(path); + if (!in) { + throw std::runtime_error("Cannot open " + path); + } + + std::string line; + + // Header 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); + } + + // Find first SegTe + 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); + } + + // Units line + 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); + } + + // Data + 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; } // ----------------------------- // JONSWAP PSD (as in your original code) // ----------------------------- -double jonswap_psd(double Hs, double Tp, double gamma, double omega, - double sigma_a = 0.07, double sigma_b = 0.09) +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)); + 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)); } // ----------------------------- // Slow strain profile // ----------------------------- -double slow_strain(double t, double seg_dur, double eps0) +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; + 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; } // ----------------------------- // Surface elevation from JONSWAP spectrum // Deterministic random phases via RNG seed // ----------------------------- -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) +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; + 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; } // Helper: "_Line1.txt" without -std::string line1_output_from_input(const std::string& input_path) +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 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::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"; + const std::string stem = has_ext ? input_path.substr(0, dot) : input_path; + return stem + "_Line1.out"; } // ----------------------------- // Run a single case and return relative L2 error // ----------------------------- -double run_case(const WcCase& c, - const Eigen::VectorXd& owc_strains, - const Eigen::VectorXd& owc_tensions) +double +run_case(const WcCase& c, + const Eigen::VectorXd& owc_strains, + const Eigen::VectorXd& owc_tensions) { - MoorDynRAII md(c.input_file); - - unsigned int ndof = 0; - check_md(MoorDyn_NCoupledDOF(md.sys, &ndof), "MoorDyn_NCoupledDOF failed for: " + c.input_file); - INFO("ndof=" << ndof); - - auto fairlead = MoorDyn_GetPoint(md.sys, 2); - auto line = MoorDyn_GetLine(md.sys, 1); - REQUIRE(fairlead != nullptr); - REQUIRE(line != 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); - - // time settings - const double seg_dur = 3.0 * 3600.0; - const double total_dur = 6.0 * seg_dur; - const double dt0 = 0.1; - const unsigned int nsteps = static_cast(total_dur / dt0); - - const Eigen::VectorXd times = - Eigen::VectorXd::LinSpaced(static_cast(nsteps) + 1, 0.0, total_dur); - - // deterministic excitation - constexpr unsigned int n_omega_edges = 300; - std::mt19937 rng_wf(12345); - std::mt19937 rng_lf(54321); - - // WF - 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); - - // LF - 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); - - // prescribed strain and x(t) - 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 = strain_slow + scale_WF * eta_WF + scale_LF * eta_LF; - const Eigen::VectorXd x = x0 * (Eigen::VectorXd::Ones(strain.size()) + strain); - - // init moordyn - 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); - - // stepping (enforce fixed dt) - 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; - } - - // close to flush outputs - check_md(MoorDyn_Close(md.sys), "MoorDyn_Close failed for: " + c.input_file); - md.sys = nullptr; - - // read output mean tension - const std::string out_path = line1_output_from_input(c.input_file); - INFO("output=" << out_path); - - Eigen::VectorXd tension_mean_output; - REQUIRE_NOTHROW(tension_mean_output = load_seg_te_column(out_path)); - - const Eigen::Index n = std::min(tension_mean_output.size(), strain.size()); - REQUIRE(n >= 2); - - // preceding highest mean tension - 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], tension_mean_output[i]); - } - - // analytical mean tension from WC + OWC - Eigen::VectorXd tension_analytical(n); - for (Eigen::Index i = 0; i < n; ++i) { - tension_analytical[i] = - find_mean_tension(strain[i], - tension_mean_output[i], - Tmax_mean[i], - owc_strains, - owc_tensions, - c.form); - } - - const double denom = tension_analytical.squaredNorm(); - REQUIRE(denom > 0.0); - - const double l2_rel = - std::sqrt((tension_analytical - tension_mean_output.head(n)).squaredNorm() / denom); - - return l2_rel; + MoorDynRAII md(c.input_file); + + unsigned int ndof = 0; + check_md(MoorDyn_NCoupledDOF(md.sys, &ndof), + "MoorDyn_NCoupledDOF failed for: " + c.input_file); + INFO("ndof=" << ndof); + + auto fairlead = MoorDyn_GetPoint(md.sys, 2); + auto line = MoorDyn_GetLine(md.sys, 1); + REQUIRE(fairlead != nullptr); + REQUIRE(line != 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); + + // time settings + const double seg_dur = 3.0 * 3600.0; + const double total_dur = 6.0 * seg_dur; + const double dt0 = 0.1; + const unsigned int nsteps = static_cast(total_dur / dt0); + + const Eigen::VectorXd times = Eigen::VectorXd::LinSpaced( + static_cast(nsteps) + 1, 0.0, total_dur); + + // deterministic excitation + constexpr unsigned int n_omega_edges = 300; + std::mt19937 rng_wf(12345); + std::mt19937 rng_lf(54321); + + // WF + 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); + + // LF + 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); + + // prescribed strain and x(t) + 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 = + strain_slow + scale_WF * eta_WF + scale_LF * eta_LF; + const Eigen::VectorXd x = + x0 * (Eigen::VectorXd::Ones(strain.size()) + strain); + + // init moordyn + 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); + + // stepping (enforce fixed dt) + 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; + } + + // close to flush outputs + check_md(MoorDyn_Close(md.sys), + "MoorDyn_Close failed for: " + c.input_file); + md.sys = nullptr; + + // read output mean tension + const std::string out_path = line1_output_from_input(c.input_file); + INFO("output=" << out_path); + + Eigen::VectorXd tension_mean_output; + REQUIRE_NOTHROW(tension_mean_output = load_seg_te_column(out_path)); + + const Eigen::Index n = + std::min(tension_mean_output.size(), strain.size()); + REQUIRE(n >= 2); + + // preceding highest mean tension + 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], tension_mean_output[i]); + } + + // analytical mean tension from WC + OWC + Eigen::VectorXd tension_analytical(n); + for (Eigen::Index i = 0; i < n; ++i) { + tension_analytical[i] = find_mean_tension(strain[i], + tension_mean_output[i], + Tmax_mean[i], + owc_strains, + owc_tensions, + c.form); + } + + const double denom = tension_analytical.squaredNorm(); + REQUIRE(denom > 0.0); + + const double l2_rel = std::sqrt( + (tension_analytical - tension_mean_output.head(n)).squaredNorm() / + denom); + + return l2_rel; } } // namespace -TEST_CASE("Syrope working curve regression: mean tension matches analytical WC", "[syrope][working-curve]") +TEST_CASE("Syrope working curve regression: mean tension matches analytical WC", + "[syrope][working-curve]") { - // shared OWC table for all cases - const auto [owc_strains, owc_tensions] = read_strain_tension_table("Mooring/syrope/owc.dat"); - - const std::vector cases = { - {"Linear", "Mooring/syrope/linear_wc_input.txt", WorkingCurveForm::Linear, 1.37912e-02}, - {"Quadratic", "Mooring/syrope/quadratic_wc_input.txt", WorkingCurveForm::Quadratic, 1.41587e-02}, - {"Exponential", "Mooring/syrope/exponential_wc_input.txt", WorkingCurveForm::Exponential, 1.18597e-02}, - }; - - for (const auto& c : cases) { - DYNAMIC_SECTION("Case: " << c.name) { - const double l2_rel = run_case(c, owc_strains, owc_tensions); - INFO("Relative L2 error = " << l2_rel); - REQUIRE(l2_rel < kL2Tol); - } - } + // shared OWC table for all cases + const auto [owc_strains, owc_tensions] = + read_strain_tension_table("Mooring/syrope/owc.dat"); + + const std::vector cases = { + { "Linear", + "Mooring/syrope/linear_wc_input.txt", + WorkingCurveForm::Linear, + 1.37912e-02 }, + { "Quadratic", + "Mooring/syrope/quadratic_wc_input.txt", + WorkingCurveForm::Quadratic, + 1.41587e-02 }, + { "Exponential", + "Mooring/syrope/exponential_wc_input.txt", + WorkingCurveForm::Exponential, + 1.18597e-02 }, + }; + + for (const auto& c : cases) { + DYNAMIC_SECTION("Case: " << c.name) + { + const double l2_rel = run_case(c, owc_strains, owc_tensions); + INFO("Relative L2 error = " << l2_rel); + REQUIRE(l2_rel < kL2Tol); + } + } } \ No newline at end of file From e489dbd09a64f9138dd553641ebd0ef1c5ac92d1 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 23 Feb 2026 09:00:41 +0000 Subject: [PATCH 08/14] Renamed the input file for the Syrope test --- .../syrope/{slow_loading.txt => fast_and_slow_loading.txt} | 0 tests/syrope.cpp | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) rename tests/Mooring/syrope/{slow_loading.txt => fast_and_slow_loading.txt} (100%) diff --git a/tests/Mooring/syrope/slow_loading.txt b/tests/Mooring/syrope/fast_and_slow_loading.txt similarity index 100% rename from tests/Mooring/syrope/slow_loading.txt rename to tests/Mooring/syrope/fast_and_slow_loading.txt diff --git a/tests/syrope.cpp b/tests/syrope.cpp index 4d938da7..9852669c 100644 --- a/tests/syrope.cpp +++ b/tests/syrope.cpp @@ -376,7 +376,7 @@ load_seg_te_column(const std::string& path) } // ----------------------------- -// JONSWAP PSD (as in your original code) +// JONSWAP PSD // ----------------------------- double jonswap_psd(double Hs, From a7d7fdb17531a178d1045667c8657fbe84702623 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Tue, 3 Mar 2026 11:54:22 +0000 Subject: [PATCH 09/14] Syrope test passed --- tests/CMakeLists.txt | 1 + tests/Mooring/syrope/exp_wc/line.txt | 35 +++ .../line.txt} | 8 +- tests/Mooring/syrope/quadratic_wc/line.txt | 35 +++ tests/syrope.cpp | 272 +++++++----------- 5 files changed, 186 insertions(+), 165 deletions(-) create mode 100644 tests/Mooring/syrope/exp_wc/line.txt rename tests/Mooring/syrope/{fast_and_slow_loading.txt => linear_wc/line.txt} (87%) create mode 100644 tests/Mooring/syrope/quadratic_wc/line.txt 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..d841c79d --- /dev/null +++ b/tests/Mooring/syrope/exp_wc/line.txt @@ -0,0 +1,35 @@ +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:../owc.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +-------------------------- SYROPE WORKING CURVES ------------------------------------------------- +LineType WCType k1 k2 +(-) (-) (-) (-) +poly EXP 0.2 1.5 +----------------------- 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 +(-) (-) (-) (-) (m) (-) (-) +1 poly 2 1 1.0 1 tsc +-------------------------- SYROPE IC ------------------------------------------------- +Line Tmax Tmean +(-) (N) (N) +1 1.3886e6 4.339e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +1.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 +100.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/fast_and_slow_loading.txt b/tests/Mooring/syrope/linear_wc/line.txt similarity index 87% rename from tests/Mooring/syrope/fast_and_slow_loading.txt rename to tests/Mooring/syrope/linear_wc/line.txt index 05aca1c2..2ecc4e8e 100644 --- a/tests/Mooring/syrope/fast_and_slow_loading.txt +++ b/tests/Mooring/syrope/linear_wc/line.txt @@ -2,7 +2,7 @@ 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:owc.dat|1.53e8|23.12 2.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +poly 0.203 40 SYROPE:../owc.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 -------------------------- SYROPE WORKING CURVES ------------------------------------------------- LineType WCType k1 k2 (-) (-) (-) (-) @@ -11,15 +11,15 @@ poly LINEAR 1.25e8 0.0 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.0137912 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 (-) (-) (-) (-) (m) (-) (-) -1 poly 2 1 1.0 1 ts +1 poly 2 1 1.0 1 tsc -------------------------- SYROPE IC ------------------------------------------------- Line Tmax Tmean (-) (N) (N) -1 1.302e6 4.339e5 +1 1.3886e6 4.3394e5 -------------------------- SOLVER OPTIONS--------------------------------------------------- 2 writeLog - Write a log file 0 g - gravitational acceleration diff --git a/tests/Mooring/syrope/quadratic_wc/line.txt b/tests/Mooring/syrope/quadratic_wc/line.txt new file mode 100644 index 00000000..482cdeba --- /dev/null +++ b/tests/Mooring/syrope/quadratic_wc/line.txt @@ -0,0 +1,35 @@ +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:../owc.dat|1.53e8|23.12 5.0e10|1.0e5 0.0 0.0 0.0 0.0 0.0 +-------------------------- SYROPE WORKING CURVES ------------------------------------------------- +LineType WCType k1 k2 +(-) (-) (-) (-) +poly QUADRATIC 0.25 1.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 +(-) (-) (-) (-) (m) (-) (-) +1 poly 2 1 1.0 1 tsc +-------------------------- SYROPE IC ------------------------------------------------- +Line Tmax Tmean +(-) (N) (N) +1 1.3886e6 4.339e5 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +2 writeLog - Write a log file +0 g - gravitational acceleration +1025 rho - water density +1.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 +100.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/syrope.cpp b/tests/syrope.cpp index 9852669c..13e667d8 100644 --- a/tests/syrope.cpp +++ b/tests/syrope.cpp @@ -1,3 +1,10 @@ +/* + * 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. + */ + #include "MoorDyn2.h" #include @@ -18,19 +25,10 @@ namespace { -// ----------------------------- -// Constants -// ----------------------------- constexpr double kG = 9.80665; constexpr double kMBL = 885.0e3 * kG; - -// initial conditions -constexpr double kTmax0 = 0.15 * kMBL; -constexpr double kTmean0 = 0.05 * kMBL; // currently unused - -// dynamic coefficients (currently unused) -constexpr double kAlpha = 17.667 * kMBL; -constexpr double kBeta = 0.2313 * 100.0; +constexpr double kTmax0 = 0.16 * kMBL; +constexpr double kL2Tol = 0.05; enum class WorkingCurveForm { @@ -39,20 +37,13 @@ enum class WorkingCurveForm Exponential }; -// working curve coefficients constexpr double kLinearK1 = 1.25e8; constexpr double kLinearK2 = 0.0; constexpr double kQuadraticK1 = 0.25; -constexpr double kQuadraticK2 = 0.80; +constexpr double kQuadraticK2 = 1.00; constexpr double kExpK1 = 0.20; -constexpr double kExpK2 = 1.10; - -// Tolerance for regression (adjust as appropriate once you have baselines) -constexpr double kL2Tol = 0.05; +constexpr double kExpK2 = 1.50; -// ----------------------------- -// working curve cases -// ----------------------------- struct WcCase { std::string name; @@ -61,15 +52,19 @@ struct WcCase double eps_0; }; -// ----------------------------- -// MoorDyn error handling -// ----------------------------- +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); + throw std::runtime_error(msg + " (err=" + std::to_string(err) + ")"); } } @@ -98,14 +93,10 @@ struct MoorDynRAII MoorDyn sys = nullptr; }; -// ----------------------------- -// Linear interpolation (clamped) -// xdata must be monotone increasing -// ----------------------------- double interpolate_clamped(double x, - const Eigen::VectorXd& xdata, - const Eigen::VectorXd& ydata) + const Eigen::VectorXd& xdata, + const Eigen::VectorXd& ydata) { if (xdata.size() != ydata.size()) { throw std::invalid_argument("interpolate_clamped: size mismatch"); @@ -149,17 +140,16 @@ interpolate_clamped(double x, return y0 + t * (y1 - y0); } -// ----------------------------- -// Mean tension from working curve -// ----------------------------- double find_mean_tension(double strain, - double Tmean, - double Tmax, - const Eigen::VectorXd& owc_strains, - const Eigen::VectorXd& owc_tensions, - WorkingCurveForm wc_form) + 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) { @@ -180,20 +170,10 @@ find_mean_tension(double strain, "find_mean_tension: unknown WorkingCurveForm"); } - // If not unloading/reloading, mean tension follows the original working - // curve - if (Tmean >= Tmax) { - return interpolate_clamped(strain, owc_strains, owc_tensions); - } - - // Unloading/reloading: build working curve between eps_min and eps_max 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); - - // Linear WC special rule: if k1 >= 1.0 treat as slope (dimensional - // stiffness) if (wc_form == WorkingCurveForm::Linear && k1 >= 1.0) { eps_min = eps_max - Tmax / k1; } @@ -206,26 +186,35 @@ find_mean_tension(double strain, 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: - return Tmax * xi; + tension_wc = Tmax * xi; + break; case WorkingCurveForm::Quadratic: - return Tmax * xi * (k2 * xi + (1.0 - k2)); + 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) - return Tmax * xi; // limit -> linear - return Tmax * (std::exp(k2 * xi) - 1.0) / den; + 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; } -// ----------------------------- -// Read OWC table: 2 header lines then two numeric columns -// ----------------------------- std::pair read_strain_tension_table(const std::string& path) { @@ -282,12 +271,6 @@ read_strain_tension_table(const std::string& path) return { eps, ten }; } -// ----------------------------- -// Read SegTe column from MoorDyn output file -// Assumes: -// 1) header tokens line -// 2) units tokens line, with "(N)" under SegNTe -// ----------------------------- Eigen::VectorXd load_seg_te_column(const std::string& path) { @@ -297,8 +280,6 @@ load_seg_te_column(const std::string& path) } std::string line; - - // Header line if (!std::getline(in, line)) { throw std::runtime_error("Missing header line in " + path); } @@ -311,7 +292,6 @@ load_seg_te_column(const std::string& path) throw std::runtime_error("Empty header line in " + path); } - // Find first SegTe const std::regex re(R"(^Seg\d+Te$)"); int col = -1; for (size_t i = 0; i < headers.size(); ++i) { @@ -324,7 +304,6 @@ load_seg_te_column(const std::string& path) throw std::runtime_error("No SegTe column found in " + path); } - // Units line if (!std::getline(in, line)) { throw std::runtime_error("Missing units line in " + path); } @@ -345,7 +324,6 @@ load_seg_te_column(const std::string& path) "', expected '(N)' in " + path); } - // Data std::vector vals; vals.reserve(2048); @@ -361,8 +339,7 @@ load_seg_te_column(const std::string& path) if (!(row >> v)) { throw std::runtime_error( "Row has fewer numeric columns than expected at line " + - std::to_string(lineno) + " in " + path + ": '" + line + - "'"); + std::to_string(lineno) + " in " + path + ": '" + line + "'"); } } vals.push_back(v); @@ -375,16 +352,13 @@ load_seg_te_column(const std::string& path) return out; } -// ----------------------------- -// JONSWAP PSD -// ----------------------------- double jonswap_psd(double Hs, - double Tp, - double gamma, - double omega, - double sigma_a = 0.07, - double sigma_b = 0.09) + 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; @@ -400,9 +374,6 @@ jonswap_psd(double Hs, std::pow(gamma, std::exp(-0.5 * rw * rw)); } -// ----------------------------- -// Slow strain profile -// ----------------------------- double slow_strain(double t, double seg_dur, double eps0) { @@ -422,19 +393,15 @@ slow_strain(double t, double seg_dur, double eps0) return eps0 + 1.0 * eps0 + 1.5 * eps0 * (t - 5.0 * seg_dur) / seg_dur; } -// ----------------------------- -// Surface elevation from JONSWAP spectrum -// Deterministic random phases via RNG seed -// ----------------------------- 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) + 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"); @@ -470,7 +437,6 @@ surface_elevation_from_jonswap(double Hs, return eta; } -// Helper: "_Line1.txt" without std::string line1_output_from_input(const std::string& input_path) { @@ -484,32 +450,24 @@ line1_output_from_input(const std::string& input_path) return stem + "_Line1.out"; } -// ----------------------------- -// Run a single case and return relative L2 error -// ----------------------------- -double +SimulationResult run_case(const WcCase& c, - const Eigen::VectorXd& owc_strains, - const Eigen::VectorXd& owc_tensions) + 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); - INFO("ndof=" << ndof); auto fairlead = MoorDyn_GetPoint(md.sys, 2); - auto line = MoorDyn_GetLine(md.sys, 1); REQUIRE(fairlead != nullptr); - REQUIRE(line != 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); - // time settings const double seg_dur = 3.0 * 3600.0; const double total_dur = 6.0 * seg_dur; const double dt0 = 0.1; @@ -518,12 +476,10 @@ run_case(const WcCase& c, const Eigen::VectorXd times = Eigen::VectorXd::LinSpaced( static_cast(nsteps) + 1, 0.0, total_dur); - // deterministic excitation constexpr unsigned int n_omega_edges = 300; std::mt19937 rng_wf(12345); std::mt19937 rng_lf(54321); - // WF const double Hs_WF = 5.0; const double Tp_WF = 12.0; const double Tz_WF = Tp_WF / 1.402; @@ -539,7 +495,6 @@ run_case(const WcCase& c, omega_max_WF, rng_wf); - // LF const double Hs_LF = 20.0; const double Tp_LF = 120.0; const double Tz_LF = Tp_LF / 1.402; @@ -555,7 +510,6 @@ run_case(const WcCase& c, omega_max_LF, rng_lf); - // prescribed strain and x(t) const double x0 = 1.0; const double scale_WF = 4.00e-4; const double scale_LF = 0.80e-4; @@ -565,21 +519,20 @@ run_case(const WcCase& c, strain_slow[i] = slow_strain(times[i], seg_dur, c.eps_0); } - const Eigen::VectorXd strain = - strain_slow + scale_WF * eta_WF + scale_LF * eta_LF; + 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); - // init moordyn 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); + check_md(MoorDyn_Init(md.sys, r, dr), "MoorDyn_Init failed for: " + c.input_file); - // stepping (enforce fixed dt) double t = 0.0; double f[3] = { 0.0, 0.0, 0.0 }; @@ -599,80 +552,77 @@ run_case(const WcCase& c, t = t_in; } - // close to flush outputs - check_md(MoorDyn_Close(md.sys), - "MoorDyn_Close failed for: " + c.input_file); + check_md(MoorDyn_Close(md.sys), "MoorDyn_Close failed for: " + c.input_file); md.sys = nullptr; - // read output mean tension - const std::string out_path = line1_output_from_input(c.input_file); - INFO("output=" << out_path); - - Eigen::VectorXd tension_mean_output; - REQUIRE_NOTHROW(tension_mean_output = load_seg_te_column(out_path)); + SimulationResult out; + out.times = times; + out.strain = strain; + out.tension_output = load_seg_te_column(line1_output_from_input(c.input_file)); + return out; +} - const Eigen::Index n = - std::min(tension_mean_output.size(), strain.size()); +double +relative_l2(const Eigen::VectorXd& ref, const Eigen::VectorXd& val) +{ + const Eigen::Index n = std::min(ref.size(), val.size()); REQUIRE(n >= 2); - - // preceding highest mean tension - 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], tension_mean_output[i]); - } - - // analytical mean tension from WC + OWC - Eigen::VectorXd tension_analytical(n); - for (Eigen::Index i = 0; i < n; ++i) { - tension_analytical[i] = find_mean_tension(strain[i], - tension_mean_output[i], - Tmax_mean[i], - owc_strains, - owc_tensions, - c.form); - } - - const double denom = tension_analytical.squaredNorm(); + const double denom = ref.head(n).squaredNorm(); REQUIRE(denom > 0.0); - - const double l2_rel = std::sqrt( - (tension_analytical - tension_mean_output.head(n)).squaredNorm() / - denom); - - return l2_rel; + return std::sqrt((ref.head(n) - val.head(n)).squaredNorm() / denom); } } // namespace -TEST_CASE("Syrope working curve regression: mean tension matches analytical WC", - "[syrope][working-curve]") +TEST_CASE("Syrope tests", "[syrope][working-curve]") { - // shared OWC table for all cases const auto [owc_strains, owc_tensions] = read_strain_tension_table("Mooring/syrope/owc.dat"); const std::vector cases = { { "Linear", - "Mooring/syrope/linear_wc_input.txt", + "Mooring/syrope/linear_wc/line.txt", WorkingCurveForm::Linear, - 1.37912e-02 }, + 1.48657e-02 }, { "Quadratic", - "Mooring/syrope/quadratic_wc_input.txt", + "Mooring/syrope/quadratic_wc/line.txt", WorkingCurveForm::Quadratic, - 1.41587e-02 }, + 1.50575e-02 }, { "Exponential", - "Mooring/syrope/exponential_wc_input.txt", + "Mooring/syrope/exp_wc/line.txt", WorkingCurveForm::Exponential, - 1.18597e-02 }, + 1.24577e-02 }, }; + const bool superimpose_fast = true; + for (const auto& c : cases) { DYNAMIC_SECTION("Case: " << c.name) { - const double l2_rel = run_case(c, owc_strains, owc_tensions); - INFO("Relative L2 error = " << l2_rel); - REQUIRE(l2_rel < kL2Tol); + 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); } } -} \ No newline at end of file +} From 41d281dfee790d6688e56029b6a21f99f9069d17 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Wed, 4 Mar 2026 11:43:00 +0000 Subject: [PATCH 10/14] Added a reference for the Syrope model --- docs/references.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/references.rst b/docs/references.rst index 4e83bd7a..72d8cf77 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -109,7 +109,7 @@ Viscoelastic approach for non-linear rope behavior: Syrope model for polyester ropes: - `Falkenberg, Erik, Åhjem, Vidar, and Limin Yang. "Best Practice for Analysis of Polyester Rope Mooring Systems." Paper presented at the Offshore Technology Conference, Houston, Texas, USA, May 2017. `_ + `Wei, Zhilong; Bingham, Harry B.; Shao, Yanlin (2026). ESOMOOR D5.1: Extended MoorDyn solver and validation report. Technical University of Denmark. `_ Updated MoorDyn-OpenFOAM Coupling: From 22818c80302fe8daecd372ca3544e8d25df9d15e Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Tue, 17 Mar 2026 14:48:44 +0000 Subject: [PATCH 11/14] Removed bold texts in Syrope input description --- docs/inputs.rst | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/inputs.rst b/docs/inputs.rst index 951d1828..b31d4899 100644 --- a/docs/inputs.rst +++ b/docs/inputs.rst @@ -333,17 +333,17 @@ needs to be input as a value. Example inputs are below: Polyester ... ... EA_s|EA_d BA_s|BA_d <-- Constant dynamic stiffness method with static and dynamic damping Polyester ... ... EA_s|EA_Dc|EA_D_Lm BA_s|BA_d <-- Load dependent dynamic stiffness method with static and dynamic damping -MoorDyn also 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 +MoorDyn also 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 (hysteresis). 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 above. -To indicate Syrope, the **first entry must start with** ``SYROPE:``, followed by the name (or path) -of a tabulated file defining the **original working curve**. The second and third entries are +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 above. +To indicate Syrope, the first entry must start with ``SYROPE:``, followed by the name (or path) +of a tabulated file defining the original working curve. The second and third entries are ``EA_Dc`` and ``EA_D_Lm`` (as in the load-dependent dynamic stiffness method). Static and dynamic -damping may be provided in **BA** as ``BA_s|BA_d`` (same convention as above). +damping may be provided in BA as ``BA_s|BA_d`` (same convention as above). Example: @@ -375,7 +375,7 @@ Example: (-) (-) (-) (-) poly LINEAR k1 k2 -Initial conditions are specified per line. The **Line** ID must match a line defined in the +Initial conditions are specified per line. The Line ID must match a line defined in the *Lines* section and must reference a line type that uses the Syrope model. - Tmax – the highest mean tension experienced by the line prior to the current time (N) From c3eb0e71b0253242bc5eb3161a0d04d8e3081599 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Tue, 17 Mar 2026 14:54:59 +0000 Subject: [PATCH 12/14] Keep the Syrope reference in Chicago style as the others --- docs/references.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/references.rst b/docs/references.rst index 72d8cf77..b531c35d 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -109,7 +109,7 @@ Viscoelastic approach for non-linear rope behavior: Syrope model for polyester ropes: - `Wei, Zhilong; Bingham, Harry B.; Shao, Yanlin (2026). ESOMOOR D5.1: Extended MoorDyn solver and validation report. Technical University of Denmark. `_ + `Wei, Zhilong, Harry B. Bingham, and Yanlin Shao. 2026. “ESOMOOR D5.1: Extended Moordyn Solver and Validation Report”. Technical University of Denmark. `_ Updated MoorDyn-OpenFOAM Coupling: From b84e560158fc859ed2a750f1307b29cf575ba90f Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Fri, 20 Mar 2026 11:43:00 +0000 Subject: [PATCH 13/14] Use (std::max) and (std::min) instead of std::max and std::min for Windows sake --- tests/syrope.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/syrope.cpp b/tests/syrope.cpp index 13e667d8..0752bf58 100644 --- a/tests/syrope.cpp +++ b/tests/syrope.cpp @@ -208,7 +208,7 @@ find_mean_tension(double strain, "find_mean_tension: unknown WorkingCurveForm"); } - if (std::abs(tension_wc - Tmax) <= 1e-12 * std::max(1.0, Tmax)) { + if (std::abs(tension_wc - Tmax) <= 1e-12 * (std::max)(1.0, Tmax)) { return interpolate_clamped(strain, owc_strains, owc_tensions); } @@ -565,7 +565,7 @@ run_case(const WcCase& c, double relative_l2(const Eigen::VectorXd& ref, const Eigen::VectorXd& val) { - const Eigen::Index n = std::min(ref.size(), val.size()); + const Eigen::Index n = (std::min)(ref.size(), val.size()); REQUIRE(n >= 2); const double denom = ref.head(n).squaredNorm(); REQUIRE(denom > 0.0); @@ -601,13 +601,13 @@ TEST_CASE("Syrope tests", "[syrope][working-curve]") { const auto sim = run_case(c, superimpose_fast); const Eigen::Index n = - std::min(sim.tension_output.size(), sim.strain.size()); + (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]); + tmax_mean[i] = (std::max)(tmax_mean[i - 1], sim.tension_output[i]); } Eigen::VectorXd tension_analytical(n); From 52eb76ded83a83441de625acf58fb13d15483f0b Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Fri, 20 Mar 2026 11:49:49 +0000 Subject: [PATCH 14/14] #define _USE_MATH_DEFINES before including cmath --- tests/syrope.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/syrope.cpp b/tests/syrope.cpp index 0752bf58..f110efdc 100644 --- a/tests/syrope.cpp +++ b/tests/syrope.cpp @@ -5,6 +5,7 @@ * whereas here we use Catch2 to check the L2-error only. */ +#define _USE_MATH_DEFINES #include "MoorDyn2.h" #include