F/fsi pitch (#1110)

* Add some operators for SixDOF

* Add missing copyrights

* Rename variables and add unit test

* Add hub motion based displacement computations

* Minor comments and syntax updates

* Add unit test for orientation displacements

include/aero/aero_utils/displacements.h

@@ -20,29 +20,49 @@ namespace aero {
 //! WienerMilenkovic parameter
 struct SixDOF
 {
-  // Kind of dangeraous constructor
+  // Kind of dangerous constructor
   SixDOF(double* vec)
-    : translation_({vec[0], vec[1], vec[2]}),
-      rotation_({vec[3], vec[4], vec[5]})
+    : position_({vec[0], vec[1], vec[2]}),
+      orientation_({vec[3], vec[4], vec[5]})
   {
   }
+
   SixDOF(vs::Vector transDisp, vs::Vector rotDisp)
-    : translation_(transDisp), rotation_(rotDisp)
+    : position_(transDisp), orientation_(rotDisp)
   {
   }
-  vs::Vector translation_;
-  vs::Vector rotation_;
+
+  vs::Vector position_;
+  vs::Vector orientation_;
 };
 
+KOKKOS_FORCEINLINE_FUNCTION
+SixDOF
+operator+(const SixDOF& a, const SixDOF& b)
+{
+  // adding b to a, so pushing b wmp onto a stack
+  return SixDOF(
+    a.position_ + b.position_, wmp::push(b.orientation_, a.orientation_));
+}
+
+KOKKOS_FORCEINLINE_FUNCTION
+SixDOF
+operator-(const SixDOF& a, const SixDOF& b)
+{
+  // subtracting b from a, so popping b from the a stack
+  return SixDOF(
+    a.position_ - b.position_, wmp::pop(b.orientation_, a.orientation_));
+}
+
 KOKKOS_FORCEINLINE_FUNCTION
 SixDOF
 linear_interp_total_displacement(
   const SixDOF start, const SixDOF end, const double interpFactor)
 {
   auto transDisp = wmp::linear_interp_translation(
-    start.translation_, end.translation_, interpFactor);
-  auto rotDisp =
-    wmp::linear_interp_rotation(start.rotation_, end.rotation_, interpFactor);
+    start.position_, end.position_, interpFactor);
+  auto rotDisp = wmp::linear_interp_rotation(
+    start.orientation_, end.orientation_, interpFactor);
   return SixDOF(transDisp, rotDisp);
 }
 
@@ -52,9 +72,9 @@ linear_interp_total_velocity(
   const SixDOF start, const SixDOF end, const double interpFactor)
 {
   auto transDisp = wmp::linear_interp_translation(
-    start.translation_, end.translation_, interpFactor);
+    start.position_, end.position_, interpFactor);
   auto rotDisp = wmp::linear_interp_translation(
-    start.translation_, end.translation_, interpFactor);
+    start.position_, end.position_, interpFactor);
   return SixDOF(transDisp, rotDisp);
 }
 
@@ -65,8 +85,20 @@ vs::Vector
 local_aero_coordinates(
   const vs::Vector inertialPos, const SixDOF aeroRefPosition)
 {
-  const auto shift = inertialPos - aeroRefPosition.translation_;
-  return wmp::rotate(aeroRefPosition.rotation_, shift);
+  const auto shift = inertialPos - aeroRefPosition.position_;
+  return wmp::rotate(aeroRefPosition.orientation_, shift);
+}
+
+// TODO(psakiev) test this
+//! Translate coordinate system for SixDOF variable from inertial to a reference
+//! coordinate system
+KOKKOS_FORCEINLINE_FUNCTION
+SixDOF
+local_aero_transformation(const SixDOF& inertialPos, const SixDOF& refPos)
+{
+  return SixDOF(
+    local_aero_coordinates(inertialPos.position_, refPos),
+    wmp::push(refPos.orientation_, inertialPos.orientation_));
 }
 
 //! Convert one array of 6 deflections (transX, transY, transZ, wmX, wmY,
@@ -78,12 +110,12 @@ compute_translational_displacements(
   const SixDOF deflections, const SixDOF referencePos, const vs::Vector cfdPos)
 {
   const auto localPos = local_aero_coordinates(cfdPos, referencePos);
-  const auto delta = cfdPos - referencePos.translation_;
+  const auto delta = cfdPos - referencePos.position_;
   // deflection roations need to be applied from the aerodynamic local frame of
   // reference
   const vs::Vector rotation =
-    wmp::rotate(deflections.rotation_, localPos, true);
-  return deflections.translation_ + rotation - delta;
+    wmp::rotate(deflections.orientation_, localPos, true);
+  return deflections.position_ + rotation - delta;
 }
 
 // TODO(psakiev) this function is a place holder for when we need to add pitch
@@ -117,8 +149,8 @@ compute_mesh_velocity(
   const vs::Vector cfdPos)
 {
   const auto pointLocal = local_aero_coordinates(cfdPos, referencePos);
-  const auto pointRotate = wmp::rotate(totalDis.rotation_, pointLocal);
-  return totalVel.translation_ + (totalVel.rotation_ ^ pointRotate);
+  const auto pointRotate = wmp::rotate(totalDis.orientation_, pointLocal);
+  return totalVel.position_ + (totalVel.orientation_ ^ pointRotate);
 }
 
 } // namespace aero

include/aero/fsi/FSIUtils.h

@@ -0,0 +1,73 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+
+#ifndef FSIUTILS_H
+#define FSIUTILS_H
+#include "aero/aero_utils/displacements.h"
+#include "vs/vector.h"
+#include <Kokkos_Macros.hpp>
+
+namespace fsi {
+
+//! compute displacements from the net motions of the hub (including rotations)
+KOKKOS_FORCEINLINE_FUNCTION
+vs::Vector
+translation_displacements_from_hub_motion(
+  const aero::SixDOF& hubRef,
+  const aero::SixDOF& hubDisp,
+  const aero::SixDOF& bladeRef)
+{
+  const auto refRelToHub =
+    aero::local_aero_coordinates(bladeRef.position_, hubRef);
+  // apply opposite rotation since it is the displacement i.e. what needs to be
+  // added to get to the current position
+  const auto positionDueToRotation =
+    wmp::rotate(hubDisp.orientation_, refRelToHub, true);
+  const auto referencePosition = bladeRef.position_ - hubRef.position_;
+  return positionDueToRotation - referencePosition + hubDisp.position_;
+}
+
+//! assemble the orientation just from changes at the root (no pitch, includes
+//! hub)
+KOKKOS_FORCEINLINE_FUNCTION vs::Vector
+orientation_displacments_from_hub_motion(
+  const aero::SixDOF& rootRef,
+  const aero::SixDOF& rootDisp,
+  const aero::SixDOF& bladeRef)
+{
+  // subtract reference root orientation from the blade reference orientation,
+  // and then rotate the blade orientation so it is in the root frame of
+  // reference
+  const auto rootRelativeRefOrientation = wmp::rotate(
+    rootRef.orientation_,
+    wmp::pop(rootRef.orientation_, bladeRef.orientation_));
+  // subtract orientation displacements at the root
+  const auto rootRelativeDisplacement =
+    wmp::rotate(rootDisp.orientation_, rootRelativeRefOrientation, true);
+  // apply displacement in the reference frame first and then rotate out of root
+  // reference frame
+  return wmp::push(rootDisp.orientation_, rootRelativeDisplacement);
+}
+
+KOKKOS_FORCEINLINE_FUNCTION
+aero::SixDOF
+displacements_from_hub_motion(
+  const aero::SixDOF& hubRef,
+  const aero::SixDOF& hubDisp,
+  const aero::SixDOF& rootRef,
+  const aero::SixDOF& rootDisp,
+  const aero::SixDOF& bladeRef)
+{
+  return aero::SixDOF(
+    translation_displacements_from_hub_motion(hubRef, hubDisp, bladeRef),
+    orientation_displacments_from_hub_motion(rootRef, rootDisp, bladeRef));
+}
+
+} // namespace fsi
+#endif

include/aero/fsi/FSIturbine.h

@@ -1,3 +1,12 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+
 #ifndef FSITURBINE_H
 #define FSITURBINE_H
 

src/aero/fsi/FSIturbine.C

@@ -1,4 +1,14 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+
 #include "aero/fsi/FSIturbine.h"
+#include "aero/fsi/FSIUtils.h"
 #include "utils/ComputeVectorDivergence.h"
 #include <NaluEnv.h>
 
@@ -7,7 +17,6 @@
 #include "stk_mesh/base/Field.hpp"
 #include "master_element/MasterElement.h"
 #include "master_element/MasterElementFactory.h"
-#include "aero/aero_utils/displacements.h"
 
 #include "netcdf.h"
 
@@ -1968,8 +1977,8 @@ fsiTurbine::mapDisplacements()
   }
 
   // Now the nacelle
-  stk::mesh::Selector nacsel(stk::mesh::selectUnion(nacelleParts_));
-  const auto& nacbkts = bulk_->get_buckets(stk::topology::NODE_RANK, nacsel);
+  stk::mesh::Selector nacelle(stk::mesh::selectUnion(nacelleParts_));
+  const auto& nacbkts = bulk_->get_buckets(stk::topology::NODE_RANK, nacelle);
   for (auto b : nacbkts) {
     for (size_t in = 0; in < b->size(); in++) {
       auto node = (*b)[in];

src/aero/fsi/OpenfastFSI.C

@@ -1,3 +1,12 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+
 #include "aero/fsi/OpenfastFSI.h"
 #include "aero/fsi/FSIturbine.h"
 #include <NaluParsing.h>

unit_tests/aero/CMakeLists.txt

@@ -1,5 +1,6 @@
 target_sources(${utest_ex_name} PRIVATE
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestWienerMilenkovic.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestDisplacements.C
+   ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestFSIUtils.C
 )
 

unit_tests/aero/UnitTestDisplacements.C

@@ -12,38 +12,84 @@
 #include <aero/aero_utils/displacements.h>
 
 namespace test_displacements {
+//! Test that two WM Params give the same end location for a point
+void
+test_wiener_milenkovic(
+  vs::Vector goldWmp, vs::Vector testWmp, vs::Vector testPoint, double eps)
+{
+  auto goldPnt = wmp::rotate(goldWmp, testPoint);
+  auto testPnt = wmp::rotate(testWmp, testPoint);
+  EXPECT_NEAR(goldPnt.x(), testPnt.x(), eps)
+    << "Gold WMP: " << goldWmp << " testWmp: " << testWmp;
+  EXPECT_NEAR(goldPnt.y(), testPnt.y(), eps)
+    << "Gold WMP: " << goldWmp << " testWmp: " << testWmp;
+  EXPECT_NEAR(goldPnt.z(), testPnt.z(), eps)
+    << "Gold WMP: " << goldWmp << " testWmp: " << testWmp;
+}
+
 TEST(AeroDisplacements, creation_from_pointer)
 {
   std::vector<double> openfastSurrogate(6, 1.0);
   aero::SixDOF disp(openfastSurrogate.data());
   for (int i = 0; i < 3; ++i) {
-    EXPECT_DOUBLE_EQ(openfastSurrogate[i], disp.translation_[i]);
-    EXPECT_DOUBLE_EQ(openfastSurrogate[i + 3], disp.rotation_[i]);
+    EXPECT_DOUBLE_EQ(openfastSurrogate[i], disp.position_[i]);
+    EXPECT_DOUBLE_EQ(openfastSurrogate[i + 3], disp.orientation_[i]);
   }
 }
 
 TEST(AeroDisplacements, creation_from_vs_vector)
 {
   aero::SixDOF disp(vs::Vector::one(), 2.0 * vs::Vector::one());
   for (int i = 0; i < 3; ++i) {
-    EXPECT_DOUBLE_EQ(1.0, disp.translation_[i]);
-    EXPECT_DOUBLE_EQ(2.0, disp.rotation_[i]);
+    EXPECT_DOUBLE_EQ(1.0, disp.position_[i]);
+    EXPECT_DOUBLE_EQ(2.0, disp.orientation_[i]);
   }
 }
 
-//! Test that two WM Params give the same end location for a point
-void
-test_wiener_milenkovic(
-  vs::Vector goldWmp, vs::Vector testWmp, vs::Vector testPoint, double eps)
+TEST(AeroDisplacements, add_six_dof_together)
 {
-  auto goldPnt = wmp::rotate(goldWmp, testPoint);
-  auto testPnt = wmp::rotate(testWmp, testPoint);
-  EXPECT_NEAR(goldPnt.x(), testPnt.x(), eps)
-    << "Gold WMP: " << goldWmp << " testWmp: " << testWmp;
-  EXPECT_NEAR(goldPnt.y(), testPnt.y(), eps)
-    << "Gold WMP: " << goldWmp << " testWmp: " << testWmp;
-  EXPECT_NEAR(goldPnt.z(), testPnt.z(), eps)
-    << "Gold WMP: " << goldWmp << " testWmp: " << testWmp;
+  const auto dispA = vs::Vector::ihat();
+  const auto dispB = vs::Vector::khat();
+
+  const auto orientA =
+    wmp::create_wm_param(vs::Vector::jhat(), utils::radians(15.0));
+  const auto orientB =
+    wmp::create_wm_param(vs::Vector::ihat(), utils::radians(15.0));
+
+  const aero::SixDOF a(dispA, orientA);
+  const aero::SixDOF b(dispB, orientB);
+
+  const auto c = a + b;
+
+  const auto goldTrans = dispA + dispB;
+  // adding b to a, so pushing b wmp onto a stack
+  const auto goldRot = wmp::push(orientB, orientA);
+  for (int i = 0; i < 3; ++i) {
+    EXPECT_DOUBLE_EQ(goldTrans[i], c.position_[i]) << i;
+  }
+  test_wiener_milenkovic(goldRot, c.orientation_, vs::Vector::one(), 1e-12);
+}
+
+TEST(AeroDisplacements, subtract_six_dof)
+{
+  const auto dispA = vs::Vector::ihat();
+  const auto dispB = vs::Vector::khat();
+
+  const auto orientA =
+    wmp::create_wm_param(vs::Vector::jhat(), utils::radians(15.0));
+  const auto orientB =
+    wmp::create_wm_param(vs::Vector::ihat(), utils::radians(15.0));
+
+  const aero::SixDOF a(dispA, orientA);
+  const aero::SixDOF b(dispB, orientB);
+
+  const auto c = a - b;
+  const auto goldTrans = dispA - dispB;
+  const auto goldRot = wmp::pop(orientB, orientA);
+  for (int i = 0; i < 3; ++i) {
+    EXPECT_DOUBLE_EQ(goldTrans[i], c.position_[i]) << i;
+  }
+  test_wiener_milenkovic(goldRot, c.orientation_, vs::Vector::one(), 1e-12);
 }
 
 // Test displacements interpolation along the (1,0,0) axis, and for rotations
@@ -62,7 +108,7 @@ TEST(AeroDisplacements, linear_interp_total_displacements)
     aero::linear_interp_total_displacement(start, end, interpFactor);
 
   for (int i = 0; i < 3; ++i) {
-    EXPECT_DOUBLE_EQ(interpFactor, interpDisp.translation_[i])
+    EXPECT_DOUBLE_EQ(interpFactor, interpDisp.position_[i])
       << "Failed i: " << i;
   }
 
@@ -71,7 +117,7 @@ TEST(AeroDisplacements, linear_interp_total_displacements)
   const auto axisOffset = (vs::Vector::one() - axis) * 1e-3;
   const vs::Vector testPoint = axis * interpFactor + axisOffset;
   auto wmpGold = wmp::create_wm_param(axis, goldAngle);
-  test_wiener_milenkovic(wmpGold, interpDisp.rotation_, testPoint, 1e-10);
+  test_wiener_milenkovic(wmpGold, interpDisp.orientation_, testPoint, 1e-10);
 }
 
 TEST(
@@ -92,7 +138,7 @@ TEST(
   // CFD Pos is using "coordinates" field so it is likely fixed in time. Need to
   // confirm with Ganesh for now we will treat this as an offset from the
   // referencePos in openfast
-  const vs::Vector cfdPos = referencePos.translation_ + vs::Vector::ihat();
+  const vs::Vector cfdPos = referencePos.position_ + vs::Vector::ihat();
 
   const aero::SixDOF deflections(vs::Vector::one() * delta, vs::Vector::zero());
 
@@ -154,7 +200,7 @@ TEST(
   // CFD Pos is using "coordinates" field so it is likely fixed in time. Need to
   // confirm with Ganesh for now we will treat this as an offset from the
   // referencePos in openfast
-  const vs::Vector cfdPos = referencePos.translation_ + vs::Vector::ihat();
+  const vs::Vector cfdPos = referencePos.position_ + vs::Vector::ihat();
 
   const aero::SixDOF deflections(
     vs::Vector::zero(), wmp::create_wm_param(vs::Vector::jhat(), angleDef));