PanelLoadLaw.h

#pragma once
 
#include "LocalThroughFlow.h"
#include "PanelLoadTypes.h"
#include "ScreenKF2012.h"
#include "ScreenMF2022.h"
#include "TwineCrossFlow.h"
 
#include <fhsim/ISimObjectLogger.h>
#include <fhsim/simobject/ISimObjectCreator.h>
#include <cmath>
#include <limits>
#include <sstream>
#include <stdexcept>
#include <string>
#include <variant>
 
namespace hydrodynamics
{
 
namespace panel_load_law_detail
{
template<class T>
PanelLoad<T> ZeroLoad(unsigned flags)
{
   PanelLoad<T> load {};
   for (int k = 0; k < 3; k++)
      load.force[k] = T(0.0);
   load.drag = load.lift = load.theta = load.reynolds = T(0.0);
   load.momentumDeficit = load.throughFlowSpeed = load.localCd = T(0.0);
   load.flags = flags;
   return load;
}
} // namespace panel_load_law_detail
 
using PanelLoadLaw = std::variant<LocalThroughFlow, ScreenKF2012, ScreenMF2022, TwineCrossFlow>;
 
template<class T>
PanelLoad<T> Evaluate(const PanelLoadLaw& law, const PanelGeometry<T>& panel, const Netting<T>& netting,
      const Fluid& fluid, const FlowSample<T>& flow)
{
   const T* u       = flow.relativeVelocity;
   const T  speedSq = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
   if (!(speedSq < std::numeric_limits<double>::infinity()))
      return panel_load_law_detail::ZeroLoad<T>(kPanelFlagFlowNotFinite);
   return std::visit([&](const auto& heldLaw) { return Evaluate(heldLaw, panel, netting, fluid, flow); }, law);
}
 
namespace panel_load_law_detail
{
 
inline bool AcceptsAtPanel(const LocalThroughFlow& /*law*/)
{
   return false;
}
 
inline bool AcceptsAtPanel(const ScreenKF2012& law)
{
   return !law.useInduction;
}
 
inline bool AcceptsAtPanel(const ScreenMF2022& /*law*/)
{
   return ScreenMF2022::acceptsAtPanel;
}
 
inline bool AcceptsAtPanel(const TwineCrossFlow& /*law*/)
{
   return true;
}
 
inline std::string FormatValue(double value)
{
   std::ostringstream text;
   text << value;
   return text.str();
}
 
[[noreturn]] inline void RejectParameter(ISimObjectCreator* creator, const std::string& parameter, const std::string& message)
{
   creator->ReportParameterError(parameter, message);
   throw std::invalid_argument(message);
}
 
inline void ReadHarmonics(ISimObjectCreator* creator, const std::string& model, double& a3, double& b4)
{
   creator->GetDoubleParam("HydroA3", &a3, a3);
   creator->GetDoubleParam("HydroB4", &b4, b4);
 
   if (!std::isfinite(a3)) {
      RejectParameter(creator, "HydroA3", "HydroA3 = " + FormatValue(a3) + " is not finite for HydroModel " + model + ".");
   }
   if (!std::isfinite(b4)) {
      RejectParameter(creator, "HydroB4", "HydroB4 = " + FormatValue(b4) + " is not finite for HydroModel " + model + ".");
   }
   if (a3 < kf2012::kMinA3) {
      const std::string problem = "HydroA3 = " + FormatValue(a3) + " is below " + FormatValue(kf2012::kMinA3)
            + " for HydroModel " + model + ".";
      const std::string reason = " The slope of C_D(theta) / cd against cos(theta) at normal incidence, 1 + 8 a3,"
                                 " is then negative, so C_D is no longer largest at normal incidence (KF2012 p. 223).";
      RejectParameter(creator, "HydroA3", problem + reason);
   }
   if (a3 > kf2012::kMaxA3) {
      const std::string problem = "HydroA3 = " + FormatValue(a3) + " is above " + FormatValue(kf2012::kMaxA3)
            + " for HydroModel " + model + ".";
      const std::string reason = " The drag factor cos(theta)(1 - 4 a3 + 4 a3 cos^2 theta) of KF2012 Eq. 14 then"
                                 " turns negative at grazing angles, theta above acos(sqrt((4 a3 - 1)/(4 a3))),"
                                 " as on trawl panels nearly along the tow or cage side panels at 75 to 90 deg"
                                 " to the current: the panel is pushed upstream and its wake speeds the flow up.";
      RejectParameter(creator, "HydroA3", problem + reason);
   }
   if (std::abs(b4) > kf2012::kMaxAbsB4) {
      const std::string problem = "HydroB4 = " + FormatValue(b4) + " is outside -" + FormatValue(kf2012::kMaxAbsB4)
            + " to " + FormatValue(kf2012::kMaxAbsB4) + " for HydroModel " + model + ".";
      const std::string reason = " The lift factor 2 b2 + 4 b4 cos(2 theta) of KF2012 Eq. 14 (b2 = 1) then changes"
                                 " sign near 90 deg (b4 > 0.5) or near 0 deg (b4 < -0.5), so the lift reverses there.";
      RejectParameter(creator, "HydroB4", problem + reason);
   }
}
 
inline TwineCrossFlow ReadTwineCrossFlow(ISimObjectCreator* creator)
{
   TwineCrossFlow law;
   creator->GetDoubleParam("Ct_nominal", &law.ct, law.ct);
   creator->GetDoubleParam("CnKnots_nominal", &law.cnKnots, law.cnKnots);
   return law;
}
 
} // namespace panel_load_law_detail
 
inline bool Accepts(const PanelLoadLaw& law, FlowReference reference)
{
   if (reference == FlowReference::FreeStream) {
      return true;
   }
   return std::visit([](const auto& heldLaw) { return panel_load_law_detail::AcceptsAtPanel(heldLaw); }, law);
}
 
inline Fluid ReadFluid(ISimObjectCreator* creator)
{
   Fluid fluid {1025.0, 1.19e-6};
   creator->GetDoubleParam("Rho", &fluid.rho, fluid.rho);
   creator->GetDoubleParam("Nu", &fluid.nu, fluid.nu);
   const auto rejectUnlessPositive = [creator](const std::string& parameter, double value) {
      if (!(std::isfinite(value) && value > 0.0))
         panel_load_law_detail::RejectParameter(creator, parameter, parameter + " = " + panel_load_law_detail::FormatValue(value) + " must be finite and positive.");
   };
   rejectUnlessPositive("Rho", fluid.rho);
   rejectUnlessPositive("Nu", fluid.nu);
   return fluid;
}
 
inline PanelLoadLaw ReadPanelLoadLaw(ISimObjectCreator* creator)
{
   std::string model;
   creator->GetStringParam("HydroModel", model, "LocalThroughFlow");
 
   if (model == "LocalThroughFlow") {
      LocalThroughFlow law;
      panel_load_law_detail::ReadHarmonics(creator, model, law.a3, law.b4);
      return law;
   }
   if (model == "ScreenKF2012") {
      ScreenKF2012 law;
      panel_load_law_detail::ReadHarmonics(creator, model, law.a3, law.b4);
      creator->GetBoolParam("HydroInduction", &law.useInduction, law.useInduction);
      return law;
   }
   if (model == "ScreenMF2022") {
      return ScreenMF2022 {};
   }
   if (model == "TwineCrossFlow") {
      return panel_load_law_detail::ReadTwineCrossFlow(creator);
   }
 
   const std::string message = "Unknown HydroModel \"" + model + "\"; expected LocalThroughFlow, ScreenKF2012, ScreenMF2022 or TwineCrossFlow.";
   creator->ReportParameterError("HydroModel", message);
   throw std::invalid_argument(message);
}
 
inline bool WarnIfOutsideFitRange(const PanelLoadLaw& law, double solidity, ISimObjectLogger& logger)
{
   if (!std::holds_alternative<ScreenMF2022>(law) || ScreenMF2022::InValidSolidityRange(solidity)) {
      return false;
   }
   std::string message = "Warning: ScreenMF2022 was fitted for solidities 0.18 to 0.36 (MF2022 Eq. 10); the solidity " + std::to_string(solidity) + " is outside that range, so its coefficients are extrapolated.";
   if (solidity < screen_mf2022::kMinSolidity) {
      message += " Below " + std::to_string(screen_mf2022::kMinSolidity) + " the angle fit of Eq. 10 gives a negative drag coefficient for some flow angles, so the solidity is clamped to " + std::to_string(screen_mf2022::kMinSolidity) + ".";
   } else if (solidity > screen_mf2022::kMaxSolidity) {
      message += " Above " + std::to_string(screen_mf2022::kMaxSolidity) + " the solidity is clamped to " + std::to_string(screen_mf2022::kMaxSolidity) + ", so the coefficients are evaluated there, not extrapolated further (MARE-0147).";
   }
   logger.LogParameterInfo("HydroModel", message);
   return true;
}
 
} // namespace hydrodynamics