Skip to content

Commit

Permalink
Wave Body: improve overrides - 8 (#71)
Browse files Browse the repository at this point in the history 
- Remove hardcoded frequency index used to lookup hydro coeffs from hdf5 data.
- Coefficients from hdf5 are linearly interpolated and scaled given the frequency in the simulation params.

Signed-off-by: Rhys Mainwaring <[email protected]>

Signed-off-by: Rhys Mainwaring <[email protected]>
  • Loading branch information
@srmainwaring
srmainwaring authored Oct 9, 2022
1 parent 1c66141 commit f9ab145
Show file tree
Hide file tree
Showing 4 changed files with 363 additions and 46 deletions.
4 changes: 2 additions & 2 deletions gz-waves-models/worlds/ellipsoid_buoy.sdf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
<world name="waves">
<physics name="1ms" type="ignore">
<max_step_size>0.001</max_step_size>
<real_time_factor>-1</real_time_factor>
<real_time_factor>1</real_time_factor>
</physics>

<plugin filename="gz-sim-physics-system"
Expand Down Expand Up @@ -110,7 +110,7 @@
<waves>
<regular>
<period>6</period>
<height>2</height>
<height>4</height>
<direction>0</direction>
<phase>0</phase>
</regular>
Expand Down
211 changes: 167 additions & 44 deletions gz-waves/src/systems/linear_wave_body/LinearWaveBody.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@
#include <vector>
#include <string>

#include <mlinterp>

using namespace gz;
using namespace sim;
using namespace systems;
Expand Down Expand Up @@ -143,7 +145,7 @@ class gz::sim::systems::LinearWaveBodyPrivate
///
/// 2. Improve pose notation to emphasise frame
/// (e.g. Peter Corke's notation).
/// Done. Use Kane/monogram (Drake) notation.
/// Use Kane/monogram (Drake) notation. Done
///
/// 3. Set waterplane pose in parameters. Done.
///
Expand Down Expand Up @@ -172,6 +174,11 @@ class gz::sim::systems::LinearWaveBodyPrivate
/// map[forceName] -> { enable, debug, publish, topic, ... }
/// and perhaps a bitmasks for switching features on/off
///
/// 10. Interpolate hdf5 date for constant coefficient case. Done.
///
/// 11. Optimise interpolation to eliminate unnecessary copies.
///
///

/// \brief WEC-Sim BEMIO hydro data structure (read from HDF5 file)
///
Expand Down Expand Up @@ -594,6 +601,78 @@ LinearWaveBody::~LinearWaveBody()
{
}

namespace
{
/// \todo Optimise to reduce / eliminate copies.
/// The inefficiency is because we read data
/// in Eigen Vectors and Matrices and represent
/// the coefficients as a vector of matrices rathen
/// a matrix of time series suitable for interpolation.

/// \brief Interpolation helpers
///
/// \param _Xd data vector
/// \param _Xi interpolation point
/// \param _Yd data vector of matrix
/// \param _Yi interpolation matrix
template <typename Matrix>
void Interp(
const Eigen::VectorXd &_Xd,
double _Xi,
const std::vector<Matrix> &_Yd,
Matrix *_Yi)
{
// check there is data to interpolate
if (_Xd.size() == 0)
return;

// get dimensions
size_t nf = _Xd.size();
size_t n1 = _Yd[0].rows();
size_t n2 = _Yd[0].cols();

// convert to std::vector
std::vector<double> Xdd;
for (size_t k=0; k<nf; ++k)
{
Xdd.push_back(_Xd[k]);
}

// reshape to (n1 * n2) x 1 x Nf vectors for interpolation
std::vector<std::vector<double>> Ydd(n1*n2);
for (size_t i=0; i<n1; ++i)
{
for (size_t j=0; j<n2; ++j)
{
size_t idx = j*n1 + i;
for (size_t k=0; k<nf; ++k)
{
auto& m = _Yd[k];
double val = m(i, j);
Ydd[idx].push_back(val);
}
}
}

// interpolate
size_t ni = 1;
double yi[] = { 0.0 };
double xi[] = { _Xi };
for (size_t i=0; i<n1; ++i)
{
for (size_t j=0; j<n2; ++j)
{
size_t idx = j*n1 + i;
size_t nd[] = { nf };
double* xd = Xdd.data();
double* yd = Ydd[idx].data();
mlinterp::interp(nd, ni, yd, yi, xd, xi);
(*_Yi)(i, j) = yi[0];
}
}
};
}

/////////////////////////////////////////////////
/// \brief Configure the model
///
Expand Down Expand Up @@ -784,11 +863,6 @@ void LinearWaveBody::Configure(const Entity &_entity,
// 5a. interpolate hydro coeffs from hdf5 data.
if (_sdf->HasElement("hdf5_file"))
{
/// \todo look up omega and interpolate
/// index of nearest value of omega in hydro data
size_t iw = 7;
double w = this->dataPtr->hydroData.w[iw];

/// \note we allow g and rho to override hdf5 values
double g = this->dataPtr->simParams.g;
double rho = this->dataPtr->simParams.rho;
Expand All @@ -797,28 +871,69 @@ void LinearWaveBody::Configure(const Entity &_entity,
this->dataPtr->hydroForceCoeffs.K_hs =
this->dataPtr->hydroData.K_hs * rho * g;

// radiation
this->dataPtr->hydroForceCoeffs.A =
this->dataPtr->hydroData.A[iw] * rho;

this->dataPtr->hydroForceCoeffs.B =
this->dataPtr->hydroData.B[iw] * rho * w;

// excitation
this->dataPtr->hydroForceCoeffs.ex_re =
this->dataPtr->hydroData.ex_re[iw] * rho * w;
this->dataPtr->hydroForceCoeffs.ex_im =
this->dataPtr->hydroData.ex_im[iw] * rho * w;

this->dataPtr->hydroForceCoeffs.fk_re =
this->dataPtr->hydroData.fk_re[iw] * rho * w;
this->dataPtr->hydroForceCoeffs.fk_im =
this->dataPtr->hydroData.fk_im[iw] * rho * w;

this->dataPtr->hydroForceCoeffs.sc_re =
this->dataPtr->hydroData.sc_re[iw] * rho * w;
this->dataPtr->hydroForceCoeffs.sc_im =
this->dataPtr->hydroData.sc_im[iw] * rho * w;
// interpolate and scale.
double w = this->dataPtr->simParams.w;

// radiation added mass [Ndof, Ndof, Nf]
Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.A, // Yd
&this->dataPtr->hydroForceCoeffs.A // Yi
);
this->dataPtr->hydroForceCoeffs.A *= rho;

// radiation damping [Ndof, Ndof, Nf]
Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.B, // Yd
&this->dataPtr->hydroForceCoeffs.B // Yi
);
this->dataPtr->hydroForceCoeffs.B *= rho * w;

// excitation combined [Ndof, Nh, Nf]
Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.ex_re, // Yd
&this->dataPtr->hydroForceCoeffs.ex_re // Yi
);
this->dataPtr->hydroForceCoeffs.ex_re *= rho * g;

Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.ex_im, // Yd
&this->dataPtr->hydroForceCoeffs.ex_im // Yi
);
this->dataPtr->hydroForceCoeffs.ex_im *= rho * g;

// excitation Froude-Krylov [Ndof, Nh, Nf]
Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.fk_re, // Yd
&this->dataPtr->hydroForceCoeffs.fk_re // Yi
);
this->dataPtr->hydroForceCoeffs.fk_re *= rho * g;

Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.fk_im, // Yd
&this->dataPtr->hydroForceCoeffs.fk_im // Yi
);
this->dataPtr->hydroForceCoeffs.fk_im *= rho * g;

// excitation scattering [Ndof, Nh, Nf]
Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.sc_re, // Yd
&this->dataPtr->hydroForceCoeffs.sc_re // Yi
);
this->dataPtr->hydroForceCoeffs.sc_re *= rho * g;

Interp(
this->dataPtr->hydroData.w, w, // Xd, Xi
this->dataPtr->hydroData.sc_im, // Yd
&this->dataPtr->hydroForceCoeffs.sc_im // Yi
);
this->dataPtr->hydroForceCoeffs.sc_im *= rho * g;
}

//////////
Expand Down Expand Up @@ -1120,9 +1235,8 @@ void LinearWaveBody::Configure(const Entity &_entity,
}

// display wave parameters
if (this->dataPtr->waveOverrideOn)
{
gzmsg << "Wave Overrides\n"
gzmsg << "Wave Parameters\n"
<< "<wave_period>: "
<< this->dataPtr->waves.period << "\n"
<< "<wave_height>: "
Expand All @@ -1133,10 +1247,22 @@ void LinearWaveBody::Configure(const Entity &_entity,
<< this->dataPtr->waves.phase << "\n";
}

// display simulation parameters
{
gzmsg << "Simulation Parameters\n"
<< "<gravity>: "
<< this->dataPtr->simParams.g << "\n"
<< "<fluid_density>: "
<< this->dataPtr->simParams.rho << "\n"
<< "T: "
<< this->dataPtr->simParams.T << "\n"
<< "w: "
<< this->dataPtr->simParams.w << "\n";
}

// display geometry parameters
if (this->dataPtr->geometryOverrideOn)
{
gzmsg << "Geometry Overrides\n"
gzmsg << "Geometry Parameters\n"
<< "<center_of_waterplane>: "
<< this->dataPtr->geometry.p_BoBwp_B << "\n"
<< "<center_of_buoyancy>: "
Expand All @@ -1146,39 +1272,36 @@ void LinearWaveBody::Configure(const Entity &_entity,
}

// display hydrostatic parameters
if (this->dataPtr->hydrostaticOverrideOn)
{
gzmsg << "Hydrostatic Coefficient Overrides\n"
gzmsg << "Hydrostatic Coefficients\n"
<< "<linear_restoring>:\n"
<< this->dataPtr->hydroForceCoeffs.K_hs << "\n";
}

// display radiation parameters
if (this->dataPtr->radiationOverrideOn)
{
gzmsg << "Radiation Coefficient Overrides\n"
gzmsg << "Radiation Coefficients\n"
<< "<added_mass>:\n"
<< this->dataPtr->hydroForceCoeffs.A << "\n"
<< "<damping>:\n"
<< this->dataPtr->hydroForceCoeffs.B << "\n";
}

// display excitation parameters
if (this->dataPtr->excitationOverrideOn)
{
gzmsg << "Excitation Coefficient Overrides\n"
gzmsg << "Excitation Coefficients\n"
<< "<ex_re>:\n"
<< this->dataPtr->hydroForceCoeffs.ex_re << "\n"
<< this->dataPtr->hydroForceCoeffs.ex_re.transpose() << "\n"
<< "<ex_im>:\n"
<< this->dataPtr->hydroForceCoeffs.ex_im << "\n"
<< this->dataPtr->hydroForceCoeffs.ex_im.transpose() << "\n"
<< "<fk_re>:\n"
<< this->dataPtr->hydroForceCoeffs.fk_re << "\n"
<< this->dataPtr->hydroForceCoeffs.fk_re.transpose() << "\n"
<< "<fk_im>:\n"
<< this->dataPtr->hydroForceCoeffs.fk_im << "\n"
<< this->dataPtr->hydroForceCoeffs.fk_im.transpose() << "\n"
<< "<sc_re>:\n"
<< this->dataPtr->hydroForceCoeffs.sc_re << "\n"
<< this->dataPtr->hydroForceCoeffs.sc_re.transpose() << "\n"
<< "<sc_im>:\n"
<< this->dataPtr->hydroForceCoeffs.sc_im << "\n";
<< this->dataPtr->hydroForceCoeffs.sc_im.transpose() << "\n";
}

/// Radiation added mass - constant coeffients
Expand Down
25 changes: 25 additions & 0 deletions gz-waves/src/systems/linear_wave_body/mlinterp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
/**
* MIT License
*
* Copyright (c) 2017 Parsiad Azimzadeh
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/

#include <mlinterp.hpp>
Loading

0 comments on commit f9ab145

Please sign in to comment.