Nektar++
Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | Friends | List of all members
Nektar::ForcingMovingBody Class Reference

#include <ForcingMovingBody.h>

Inheritance diagram for Nektar::ForcingMovingBody:
[legend]

Static Public Member Functions

static SolverUtils::ForcingSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
 Creates an instance of this class. More...
 
- Static Public Member Functions inherited from Nektar::SolverUtils::Forcing
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtrLoad (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
 

Static Public Attributes

static std::string className
 Name of the class. More...
 

Protected Member Functions

virtual void v_InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
 
virtual void v_Apply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time)
 
- Protected Member Functions inherited from Nektar::SolverUtils::Forcing
SOLVER_UTILS_EXPORT Forcing (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation)
 Constructor. More...
 
virtual SOLVER_UTILS_EXPORT void v_PreApply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time)
 
virtual SOLVER_UTILS_EXPORT void v_ApplyCoeff (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time)
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, std::string pName, bool pCache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void EvaluateTimeFunction (LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
 
SOLVER_UTILS_EXPORT void EvaluateTimeFunction (const NekDouble pTime, const LibUtilities::EquationSharedPtr &pEqn, Array< OneD, NekDouble > &pArray)
 

Protected Attributes

GlobalMapping::MappingSharedPtr m_mapping
 
- Protected Attributes inherited from Nektar::SolverUtils::Forcing
LibUtilities::SessionReaderSharedPtr m_session
 Session reader. More...
 
const std::weak_ptr< EquationSystemm_equ
 Weak pointer to equation system using this forcing. More...
 
Array< OneD, Array< OneD, NekDouble > > m_Forcing
 Evaluated forcing function. More...
 
int m_NumVariable
 Number of variables. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtrm_sessionFunctions
 Map of known SessionFunctions. More...
 

Private Member Functions

 ForcingMovingBody (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation)
 
void CheckIsFromFile (const TiXmlElement *pForce)
 
void InitialiseCableModel (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 
void InitialiseFilter (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
 
void Newmark_betaSolver (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
 
void EvaluateStructDynModel (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Hydroforces, NekDouble time)
 
void SetDynEqCoeffMatrix (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 
void RollOver (Array< OneD, Array< OneD, NekDouble >> &input)
 

Private Attributes

int m_movingBodyCalls
 number of times the movbody have been called More...
 
int m_np
 number of planes per processors More...
 
int m_vdim
 vibration dimension More...
 
NekDouble m_structrho
 mass of the cable per unit length More...
 
NekDouble m_structdamp
 damping ratio of the cable More...
 
NekDouble m_lhom
 length ratio of the cable More...
 
NekDouble m_kinvis
 fluid viscous More...
 
NekDouble m_timestep
 time step More...
 
LibUtilities::NektarFFTSharedPtr m_FFT
 
FilterMovingBodySharedPtr m_MovBodyfilter
 
Array< OneD, NekDoublem_Aeroforces
 storage for the cable's force(x,y) variables More...
 
Array< OneD, Array< OneD, NekDouble > > m_MotionVars
 storage for the cable's motion(x,y) variables More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
 fictitious velocity storage More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
 fictitious acceleration storage More...
 
Array< OneD, DNekMatSharedPtrm_CoeffMat_A
 matrices in Newmart-beta method More...
 
Array< OneD, DNekMatSharedPtrm_CoeffMat_B
 matrices in Newmart-beta method More...
 
Array< OneD, std::string > m_funcName
 [0] is displacements, [1] is velocities, [2] is accelerations More...
 
Array< OneD, std::string > m_motion
 motion direction: [0] is 'x' and [1] is 'y' More...
 
Array< OneD, bool > m_IsFromFile
 do determine if the the body motion come from an extern file More...
 
Array< OneD, Array< OneD, NekDouble > > m_zta
 Store the derivatives of motion variables in x-direction. More...
 
Array< OneD, Array< OneD, NekDouble > > m_eta
 Store the derivatives of motion variables in y-direction. More...
 
unsigned int m_outputFrequency
 
Array< OneD, std::ofstream > m_outputStream
 
std::string m_outputFile_fce
 
std::string m_outputFile_mot
 

Friends

class MemoryManager< ForcingMovingBody >
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Forcing
virtual SOLVER_UTILS_EXPORT ~Forcing ()
 
SOLVER_UTILS_EXPORT void InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
 Initialise the forcing object. More...
 
SOLVER_UTILS_EXPORT void Apply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time)
 Apply the forcing. More...
 
SOLVER_UTILS_EXPORT void PreApply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time)
 Change the advection velocity before applying the forcing. For example, subtracting the frame velocity from the advection velocity in the MovingRefercenceFrame. More...
 
SOLVER_UTILS_EXPORT void ApplyCoeff (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time)
 Apply the forcing. More...
 
SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & GetForces ()
 
SOLVER_UTILS_EXPORT Array< OneD, Array< OneD, NekDouble > > & UpdateForces ()
 

Detailed Description

Definition at line 49 of file ForcingMovingBody.h.

Constructor & Destructor Documentation

◆ ForcingMovingBody()

Nektar::ForcingMovingBody::ForcingMovingBody ( const LibUtilities::SessionReaderSharedPtr pSession,
const std::weak_ptr< SolverUtils::EquationSystem > &  pEquation 
)
private

Definition at line 48 of file ForcingMovingBody.cpp.

51  : Forcing(pSession, pEquation)
52 {
53 }
SOLVER_UTILS_EXPORT Forcing(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation)
Constructor.
Definition: Forcing.cpp:50

Member Function Documentation

◆ CheckIsFromFile()

void Nektar::ForcingMovingBody::CheckIsFromFile ( const TiXmlElement *  pForce)
private

Definition at line 1063 of file ForcingMovingBody.cpp.

1064 {
1065 
1066  m_funcName = Array<OneD, std::string>(3);
1067  m_motion = Array<OneD, std::string>(2);
1068  m_motion[0] = "x";
1069  m_motion[1] = "y";
1070 
1071  m_IsFromFile = Array<OneD, bool>(6);
1072  // Loading the x-dispalcement (m_zta) and the y-displacement (m_eta)
1073  // Those two variables are bith functions of z and t and the may come
1074  // from an equation (forced vibration) or from another solver which, given
1075  // the aerodynamic forces at the previous step, calculates the
1076  // displacements.
1077 
1078  // Get the body displacement: m_zta and m_eta
1079  const TiXmlElement *funcNameElmt_D =
1080  pForce->FirstChildElement("DISPLACEMENTS");
1081  ASSERTL0(funcNameElmt_D,
1082  "MOVINGBODYFORCE tag has to specify a function name which "
1083  "prescribes the body displacement as d(z,t).");
1084 
1085  m_funcName[0] = funcNameElmt_D->GetText();
1086  ASSERTL0(m_session->DefinesFunction(m_funcName[0]),
1087  "Function '" + m_funcName[0] + "' not defined.");
1088 
1089  // Get the body velocity of movement: d(m_zta)/dt and d(m_eta)/dt
1090  const TiXmlElement *funcNameElmt_V =
1091  pForce->FirstChildElement("VELOCITIES");
1092  ASSERTL0(funcNameElmt_D,
1093  "MOVINGBODYFORCE tag has to specify a function name which "
1094  "prescribes the body velocity of movement as v(z,t).");
1095 
1096  m_funcName[1] = funcNameElmt_V->GetText();
1097  ASSERTL0(m_session->DefinesFunction(m_funcName[1]),
1098  "Function '" + m_funcName[1] + "' not defined.");
1099 
1100  // Get the body acceleration: dd(m_zta)/ddt and dd(m_eta)/ddt
1101  const TiXmlElement *funcNameElmt_A =
1102  pForce->FirstChildElement("ACCELERATIONS");
1103  ASSERTL0(funcNameElmt_A,
1104  "MOVINGBODYFORCE tag has to specify a function name which "
1105  "prescribes the body acceleration as a(z,t).");
1106 
1107  m_funcName[2] = funcNameElmt_A->GetText();
1108  ASSERTL0(m_session->DefinesFunction(m_funcName[2]),
1109  "Function '" + m_funcName[2] + "' not defined.");
1110 
1112 
1113  // Check Displacement x
1114  vType = m_session->GetFunctionType(m_funcName[0], m_motion[0]);
1116  {
1117  m_IsFromFile[0] = false;
1118  }
1119  else if (vType == LibUtilities::eFunctionTypeFile)
1120  {
1121  m_IsFromFile[0] = true;
1122  }
1123  else
1124  {
1125  ASSERTL0(false, "The displacements in x must be specified via an "
1126  "equation <E> or a file <F>");
1127  }
1128 
1129  // Check Displacement y
1130  vType = m_session->GetFunctionType(m_funcName[0], m_motion[1]);
1132  {
1133  m_IsFromFile[1] = false;
1134  }
1135  else if (vType == LibUtilities::eFunctionTypeFile)
1136  {
1137  m_IsFromFile[1] = true;
1138  }
1139  else
1140  {
1141  ASSERTL0(false, "The displacements in y must be specified via an "
1142  "equation <E> or a file <F>");
1143  }
1144 
1145  // Check Velocity x
1146  vType = m_session->GetFunctionType(m_funcName[1], m_motion[0]);
1148  {
1149  m_IsFromFile[2] = false;
1150  }
1151  else if (vType == LibUtilities::eFunctionTypeFile)
1152  {
1153  m_IsFromFile[2] = true;
1154  }
1155  else
1156  {
1157  ASSERTL0(false, "The velocities in x must be specified via an equation "
1158  "<E> or a file <F>");
1159  }
1160 
1161  // Check Velocity y
1162  vType = m_session->GetFunctionType(m_funcName[1], m_motion[1]);
1164  {
1165  m_IsFromFile[3] = false;
1166  }
1167  else if (vType == LibUtilities::eFunctionTypeFile)
1168  {
1169  m_IsFromFile[3] = true;
1170  }
1171  else
1172  {
1173  ASSERTL0(false, "The velocities in y must be specified via an equation "
1174  "<E> or a file <F>");
1175  }
1176 
1177  // Check Acceleration x
1178  vType = m_session->GetFunctionType(m_funcName[2], m_motion[0]);
1180  {
1181  m_IsFromFile[4] = false;
1182  }
1183  else if (vType == LibUtilities::eFunctionTypeFile)
1184  {
1185  m_IsFromFile[4] = true;
1186  }
1187  else
1188  {
1189  ASSERTL0(false, "The accelerations in x must be specified via an "
1190  "equation <E> or a file <F>");
1191  }
1192 
1193  // Check Acceleration y
1194  vType = m_session->GetFunctionType(m_funcName[2], m_motion[1]);
1196  {
1197  m_IsFromFile[5] = false;
1198  }
1199  else if (vType == LibUtilities::eFunctionTypeFile)
1200  {
1201  m_IsFromFile[5] = true;
1202  }
1203  else
1204  {
1205  ASSERTL0(false, "The accelerations in y must be specified via an "
1206  "equation <E> or a file <F>");
1207  }
1208 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:116

References ASSERTL0, Nektar::LibUtilities::eFunctionTypeExpression, Nektar::LibUtilities::eFunctionTypeFile, m_funcName, m_IsFromFile, m_motion, and Nektar::SolverUtils::Forcing::m_session.

Referenced by v_InitObject().

◆ create()

static SolverUtils::ForcingSharedPtr Nektar::ForcingMovingBody::create ( const LibUtilities::SessionReaderSharedPtr pSession,
const std::weak_ptr< SolverUtils::EquationSystem > &  pEquation,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
const unsigned int &  pNumForcingFields,
const TiXmlElement *  pForce 
)
inlinestatic

Creates an instance of this class.

Definition at line 55 of file ForcingMovingBody.h.

60  {
63  pEquation);
64  p->InitObject(pFields, pNumForcingFields, pForce);
65  return p;
66  }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
SOLVER_UTILS_EXPORT typedef std::shared_ptr< Forcing > ForcingSharedPtr
A shared pointer to an EquationSystem object.
Definition: Forcing.h:52

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::SolverUtils::ForcingSharedPtr, and CellMLToNektar.cellml_metadata::p.

◆ EvaluateStructDynModel()

void Nektar::ForcingMovingBody::EvaluateStructDynModel ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
Array< OneD, NekDouble > &  Hydroforces,
NekDouble  time 
)
private

Definition at line 225 of file ForcingMovingBody.cpp.

228 {
229  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
230  int colrank = vcomm->GetColumnComm()->GetRank();
231  int nproc = vcomm->GetColumnComm()->GetSize();
232 
233  bool homostrip;
234  m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
235 
236  // number of structral modes and number of strips
237  int npts, nstrips;
238  if (!homostrip) // full resolutions
239  {
240  npts = m_session->GetParameter("HomModesZ");
241  }
242  else
243  {
244  m_session->LoadParameter("HomStructModesZ", npts);
245  m_session->LoadParameter("Strip_Z", nstrips);
246 
247  // ASSERTL0(nstrips < = npts,
248  // "the number of struct. modes must be larger than that of the
249  // strips.");
250  }
251 
252  // the hydrodynamic forces
253  Array<OneD, Array<OneD, NekDouble>> fces(2);
254  // forces in x-direction
255  fces[0] = Array<OneD, NekDouble>(npts, 0.0);
256  // forces in y-direction
257  fces[1] = Array<OneD, NekDouble>(npts, 0.0);
258 
259  // fill the force vectors
260  if (colrank == 0)
261  {
262  if (!homostrip) // full resolutions
263  {
264  Vmath::Vcopy(m_np, Hydroforces, 1, fces[0], 1);
265  Vmath::Vcopy(m_np, Hydroforces + m_np, 1, fces[1], 1);
266  }
267  else // strip modelling
268  {
269  fces[0][0] = Hydroforces[0];
270  fces[1][0] = Hydroforces[m_np];
271  }
272 
273  if (!homostrip) // full resolutions
274  {
275  Array<OneD, NekDouble> tmp(2 * m_np);
276  for (int i = 1; i < nproc; ++i)
277  {
278  vcomm->GetColumnComm()->Recv(i, tmp);
279  for (int n = 0; n < m_np; n++)
280  {
281  for (int j = 0; j < 2; ++j)
282  {
283  fces[j][i * m_np + n] = tmp[j * m_np + n];
284  }
285  }
286  }
287  }
288  else // strip modelling
289  // if the body is submerged partly, the fces are filled partly
290  // by the flow induced forces
291  {
292  Array<OneD, NekDouble> tmp(2);
293  for (int i = 1; i < nstrips; ++i)
294  {
295  vcomm->GetColumnComm()->Recv(i, tmp);
296 
297  for (int j = 0; j < 2; ++j)
298  {
299  fces[j][i] = tmp[j];
300  }
301  }
302  }
303  }
304  else
305  {
306  if (!homostrip) // full resolutions
307  {
308  vcomm->GetColumnComm()->Send(0, Hydroforces);
309  }
310  else // strip modelling
311  {
312  for (int i = 1; i < nstrips; ++i)
313  {
314  if (colrank == i)
315  {
316  Array<OneD, NekDouble> tmp(2);
317  tmp[0] = Hydroforces[0];
318  tmp[1] = Hydroforces[m_np];
319  vcomm->GetColumnComm()->Send(0, tmp);
320  }
321  }
322  }
323  }
324 
325  if (colrank == 0)
326  {
327  // Fictitious mass method used to stablize the explicit coupling at
328  // relatively lower mass ratio
329  bool fictmass;
330  m_session->MatchSolverInfo("FictitiousMassMethod", "True", fictmass,
331  false);
332  if (fictmass)
333  {
334  NekDouble fictrho, fictdamp;
335  m_session->LoadParameter("FictMass", fictrho);
336  m_session->LoadParameter("FictDamp", fictdamp);
337 
338  static NekDouble Betaq_Coeffs[2][2] = {{1.0, 0.0}, {2.0, -1.0}};
339 
340  // only consider second order approximation for fictitious variables
341  int intOrder = 2;
342  int nint = min(m_movingBodyCalls + 1, intOrder);
343  int nlevels = m_fV[0].size();
344 
345  for (int i = 0; i < m_motion.size(); ++i)
346  {
347  RollOver(m_fV[i]);
348  RollOver(m_fA[i]);
349 
350  int Voffset = npts;
351  int Aoffset = 2 * npts;
352 
353  Vmath::Vcopy(npts, m_MotionVars[i] + Voffset, 1, m_fV[i][0], 1);
354  Vmath::Vcopy(npts, m_MotionVars[i] + Aoffset, 1, m_fA[i][0], 1);
355 
356  // Extrapolate to n+1
357  Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
358  m_fV[i][nint - 1], 1, m_fV[i][nlevels - 1], 1);
359  Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
360  m_fA[i][nint - 1], 1, m_fA[i][nlevels - 1], 1);
361 
362  for (int n = 0; n < nint - 1; ++n)
363  {
364  Vmath::Svtvp(npts, Betaq_Coeffs[nint - 1][n], m_fV[i][n], 1,
365  m_fV[i][nlevels - 1], 1, m_fV[i][nlevels - 1],
366  1);
367  Vmath::Svtvp(npts, Betaq_Coeffs[nint - 1][n], m_fA[i][n], 1,
368  m_fA[i][nlevels - 1], 1, m_fA[i][nlevels - 1],
369  1);
370  }
371 
372  // Add the fictitious forces on the RHS of the equation
373  Vmath::Svtvp(npts, fictdamp, m_fV[i][nlevels - 1], 1, fces[i],
374  1, fces[i], 1);
375  Vmath::Svtvp(npts, fictrho, m_fA[i][nlevels - 1], 1, fces[i], 1,
376  fces[i], 1);
377  }
378  }
379  }
380  // structural solver is implemented on the root process
381  if (colrank == 0)
382  {
383  // Tensioned cable model is evaluated in wave space
384  for (int n = 0, cn = 1; n < m_vdim; n++, cn--)
385  {
386  Newmark_betaSolver(pFields, fces[cn], m_MotionVars[cn]);
387  }
388  }
389 
390  Array<OneD, NekDouble> Motvars(2 * 2 * m_np);
391  // send physical coeffients to all planes of each processor
392  if (!homostrip) // full resolutions
393  {
394  Array<OneD, NekDouble> tmp(2 * 2 * m_np);
395 
396  if (colrank != 0)
397  {
398  vcomm->GetColumnComm()->Recv(0, tmp);
399  Vmath::Vcopy(2 * 2 * m_np, tmp, 1, Motvars, 1);
400  }
401  else
402  {
403  for (int i = 1; i < nproc; ++i)
404  {
405  for (int j = 0; j < 2; j++) // moving dimensions
406  {
407  for (int k = 0; k < 2; k++) // disp. and vel.
408  {
409  for (int n = 0; n < m_np; n++)
410  {
411  tmp[j * 2 * m_np + k * m_np + n] =
412  m_MotionVars[j][k * npts + i * m_np + n];
413  }
414  }
415  }
416  vcomm->GetColumnComm()->Send(i, tmp);
417  }
418 
419  for (int j = 0; j < 2; j++)
420  {
421  for (int k = 0; k < 2; k++)
422  {
423  for (int n = 0; n < m_np; n++)
424  {
425  tmp[j * 2 * m_np + k * m_np + n] =
426  m_MotionVars[j][k * npts + n];
427  }
428  Vmath::Vcopy(2 * 2 * m_np, tmp, 1, Motvars, 1);
429  }
430  }
431  }
432  }
433  else // strip modelling
434  {
435  Array<OneD, NekDouble> tmp(2 * 2);
436 
437  if (colrank != 0)
438  {
439  for (int j = 1; j < nproc / nstrips; j++)
440  {
441  if (colrank == j * nstrips)
442  {
443  vcomm->GetColumnComm()->Recv(0, tmp);
444 
445  for (int plane = 0; plane < m_np; plane++)
446  {
447  for (int var = 0; var < 2; var++)
448  {
449  for (int k = 0; k < 2; k++)
450  {
451  Motvars[var * 2 * m_np + k * m_np + plane] =
452  tmp[var * 2 + k];
453  }
454  }
455  }
456  }
457  }
458 
459  for (int i = 1; i < nstrips; i++)
460  {
461  for (int j = 0; j < nproc / nstrips; j++)
462  {
463  if (colrank == i + j * nstrips)
464  {
465  vcomm->GetColumnComm()->Recv(0, tmp);
466 
467  for (int plane = 0; plane < m_np; plane++)
468  {
469  for (int var = 0; var < 2; var++)
470  {
471  for (int k = 0; k < 2; k++)
472  {
473  Motvars[var * 2 * m_np + k * m_np + plane] =
474  tmp[var * 2 + k];
475  }
476  }
477  }
478  }
479  }
480  }
481  }
482  else
483  {
484  for (int j = 0; j < 2; ++j)
485  {
486  for (int k = 0; k < 2; ++k)
487  {
488  tmp[j * 2 + k] = m_MotionVars[j][k * npts];
489  }
490  }
491 
492  for (int j = 1; j < nproc / nstrips; j++)
493  {
494  vcomm->GetColumnComm()->Send(j * nstrips, tmp);
495  }
496 
497  for (int plane = 0; plane < m_np; plane++)
498  {
499  for (int j = 0; j < 2; j++)
500  {
501  for (int k = 0; k < 2; ++k)
502  {
503  Motvars[j * 2 * m_np + k * m_np + plane] =
504  m_MotionVars[j][k * npts];
505  }
506  }
507  }
508 
509  for (int i = 1; i < nstrips; ++i)
510  {
511  for (int j = 0; j < 2; ++j)
512  {
513  for (int k = 0; k < 2; ++k)
514  {
515  tmp[j * 2 + k] = m_MotionVars[j][i + k * npts];
516  }
517  }
518 
519  for (int j = 0; j < nproc / nstrips; j++)
520  {
521  vcomm->GetColumnComm()->Send(i + j * nstrips, tmp);
522  }
523  }
524  }
525  }
526 
527  // Set the m_forcing term based on the motion of the cable
528  for (int var = 0; var < 2; var++)
529  {
530  for (int plane = 0; plane < m_np; plane++)
531  {
532  int n = pFields[0]->GetPlane(plane)->GetTotPoints();
533 
534  Array<OneD, NekDouble> tmp;
535 
536  int offset = plane * n;
537  int xoffset = var * m_np + plane;
538  int yoffset = 2 * m_np + xoffset;
539 
540  Vmath::Fill(n, Motvars[xoffset], tmp = m_zta[var] + offset, 1);
541  Vmath::Fill(n, Motvars[yoffset], tmp = m_eta[var] + offset, 1);
542  }
543  }
544 }
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
fictitious velocity storage
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
fictitious acceleration storage
int m_np
number of planes per processors
int m_vdim
vibration dimension
Array< OneD, Array< OneD, NekDouble > > m_MotionVars
storage for the cable's motion(x,y) variables
int m_movingBodyCalls
number of times the movbody have been called
void RollOver(Array< OneD, Array< OneD, NekDouble >> &input)
void Newmark_betaSolver(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
double NekDouble
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255

References Vmath::Fill(), m_eta, m_fA, m_fV, m_motion, m_MotionVars, m_movingBodyCalls, m_np, Nektar::SolverUtils::Forcing::m_session, m_vdim, m_zta, Newmark_betaSolver(), RollOver(), Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

Referenced by v_Apply().

◆ InitialiseCableModel()

void Nektar::ForcingMovingBody::InitialiseCableModel ( const LibUtilities::SessionReaderSharedPtr pSession,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields 
)
private

Definition at line 677 of file ForcingMovingBody.cpp.

680 {
681  m_movingBodyCalls = 0;
682  m_session->LoadParameter("Kinvis", m_kinvis);
683  m_session->LoadParameter("TimeStep", m_timestep, 0.01);
684 
685  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
686  int colrank = vcomm->GetColumnComm()->GetRank();
687  int nproc = vcomm->GetColumnComm()->GetSize();
688 
689  // number of structral modes
690  int npts;
691  bool homostrip;
692  m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
693 
694  if (!homostrip) // full resolutions
695  {
696  npts = m_session->GetParameter("HomModesZ");
697  }
698  else
699  {
700  m_session->LoadParameter("HomStructModesZ", npts);
701  }
702 
703  m_MotionVars = Array<OneD, Array<OneD, NekDouble>>(2);
704  m_MotionVars[0] = Array<OneD, NekDouble>(3 * npts, 0.0);
705  m_MotionVars[1] = Array<OneD, NekDouble>(3 * npts, 0.0);
706 
707  Array<OneD, unsigned int> ZIDs;
708  ZIDs = pFields[0]->GetZIDs();
709  m_np = ZIDs.size();
710 
711  std::string vibtype = m_session->GetSolverInfo("VibrationType");
712 
713  if (boost::iequals(vibtype, "Constrained"))
714  {
715  m_vdim = 1;
716  }
717  else if (boost::iequals(vibtype, "Free"))
718  {
719  m_vdim = 2;
720  }
721  else if (boost::iequals(vibtype, "Forced"))
722  {
723  m_vdim = 0;
724  return;
725  }
726 
727  if (!homostrip)
728  {
729  m_session->LoadParameter("LZ", m_lhom);
730  int nplanes = m_session->GetParameter("HomModesZ");
732  nplanes);
733  }
734  else
735  {
736  int nstrips;
737  NekDouble DistStrip;
738 
739  m_session->LoadParameter("DistStrip", DistStrip);
740  m_session->LoadParameter("Strip_Z", nstrips);
741  m_lhom = nstrips * DistStrip;
743  nstrips);
744  }
745 
746  // load the structural dynamic parameters from xml file
747  m_session->LoadParameter("StructRho", m_structrho);
748  m_session->LoadParameter("StructDamp", m_structdamp, 0.0);
749 
750  // Identify whether the fictitious mass method is active for explicit
751  // coupling of fluid solver and structural dynamics solver
752  bool fictmass;
753  m_session->MatchSolverInfo("FictitiousMassMethod", "True", fictmass, false);
754  if (fictmass)
755  {
756  NekDouble fictrho, fictdamp;
757  m_session->LoadParameter("FictMass", fictrho);
758  m_session->LoadParameter("FictDamp", fictdamp);
759  m_structrho += fictrho;
760  m_structdamp += fictdamp;
761 
762  // Storage array of Struct Velocity and Acceleration used for
763  // extrapolation of fictitious force
764  m_fV = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
765  m_fA = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
766  for (int i = 0; i < m_motion.size(); ++i)
767  {
768  m_fV[i] = Array<OneD, Array<OneD, NekDouble>>(2);
769  m_fA[i] = Array<OneD, Array<OneD, NekDouble>>(2);
770 
771  for (int n = 0; n < 2; ++n)
772  {
773  m_fV[i][n] = Array<OneD, NekDouble>(npts, 0.0);
774  m_fA[i][n] = Array<OneD, NekDouble>(npts, 0.0);
775  }
776  }
777  }
778 
779  // Setting the coefficient matrices for solving structural dynamic ODEs
780  SetDynEqCoeffMatrix(pFields);
781 
782  // Set initial condition for cable's motion
783  int cnt = 0;
784 
785  for (int j = 0; j < m_funcName.size(); j++)
786  {
787  // loading from the specified files through inputstream
788  if (m_IsFromFile[cnt] && m_IsFromFile[cnt + 1])
789  {
790  std::ifstream inputStream;
791  int nzpoints;
792 
793  if (homostrip)
794  {
795  m_session->LoadParameter("HomStructModesZ", nzpoints);
796  }
797  else
798  {
799  nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
800  }
801 
802  if (vcomm->GetRank() == 0)
803  {
804  std::string filename =
805  m_session->GetFunctionFilename(m_funcName[j], m_motion[0]);
806 
807  // Open intputstream for cable motions
808  inputStream.open(filename.c_str());
809 
810  // Import the head string from the file
811  Array<OneD, std::string> tmp(9);
812  for (int n = 0; n < tmp.size(); n++)
813  {
814  inputStream >> tmp[n];
815  }
816 
817  NekDouble time, z_cds;
818  // Import the motion variables from the file
819  for (int n = 0; n < nzpoints; n++)
820  {
821  inputStream >> setprecision(6) >> time;
822  inputStream >> setprecision(6) >> z_cds;
823  inputStream >> setprecision(8) >> m_MotionVars[0][n];
824  inputStream >> setprecision(8) >>
825  m_MotionVars[0][n + nzpoints];
826  inputStream >> setprecision(8) >>
827  m_MotionVars[0][n + 2 * nzpoints];
828  inputStream >> setprecision(8) >> m_MotionVars[1][n];
829  inputStream >> setprecision(8) >>
830  m_MotionVars[1][n + nzpoints];
831  inputStream >> setprecision(8) >>
832  m_MotionVars[1][n + 2 * nzpoints];
833  }
834  // Close inputstream for cable motions
835  inputStream.close();
836  }
837  cnt = cnt + 2;
838  }
839  else // Evaluate from the functions specified in xml file
840  {
841  if (!homostrip)
842  {
843  int nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
844  Array<OneD, NekDouble> z_coords(nzpoints, 0.0);
845  Array<OneD, const NekDouble> pts =
846  pFields[0]->GetHomogeneousBasis()->GetZ();
847 
848  Vmath::Smul(nzpoints, m_lhom / 2.0, pts, 1, z_coords, 1);
849  Vmath::Sadd(nzpoints, m_lhom / 2.0, z_coords, 1, z_coords, 1);
850 
851  Array<OneD, NekDouble> x0(m_np, 0.0);
852  Array<OneD, NekDouble> x1(m_np, 0.0);
853  Array<OneD, NekDouble> x2(m_np, 0.0);
854  for (int plane = 0; plane < m_np; plane++)
855  {
856  x2[plane] = z_coords[ZIDs[plane]];
857  }
858 
859  Array<OneD, NekDouble> tmp(m_np, 0.0);
860  Array<OneD, NekDouble> tmp0(m_np, 0.0);
861  Array<OneD, NekDouble> tmp1(m_np, 0.0);
862  LibUtilities::EquationSharedPtr ffunc0, ffunc1;
863 
864  ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
865  ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
866  ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
867  ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
868 
869  int offset = j * npts;
870  Vmath::Vcopy(m_np, tmp0, 1, tmp = m_MotionVars[0] + offset, 1);
871  Vmath::Vcopy(m_np, tmp1, 1, tmp = m_MotionVars[1] + offset, 1);
872 
873  if (colrank == 0)
874  {
875  for (int i = 1; i < nproc; ++i)
876  {
877  vcomm->GetColumnComm()->Recv(i, tmp0);
878  vcomm->GetColumnComm()->Recv(i, tmp1);
879  Vmath::Vcopy(m_np, tmp0, 1,
880  tmp = m_MotionVars[0] + offset + i * m_np,
881  1);
882  Vmath::Vcopy(m_np, tmp1, 1,
883  tmp = m_MotionVars[1] + offset + i * m_np,
884  1);
885  }
886  }
887  else
888  {
889  vcomm->GetColumnComm()->Send(0, tmp0);
890  vcomm->GetColumnComm()->Send(0, tmp1);
891  }
892  }
893  else
894  {
895  if (colrank == 0)
896  {
897  int nstrips;
898  m_session->LoadParameter("Strip_Z", nstrips);
899 
900  ASSERTL0(
901  m_session->DefinesSolverInfo("USEFFT"),
902  "Fourier transformation of cable motion is currently "
903  "implemented only for FFTW module.");
904 
905  NekDouble DistStrip;
906  m_session->LoadParameter("DistStrip", DistStrip);
907 
908  Array<OneD, NekDouble> x0(npts, 0.0);
909  Array<OneD, NekDouble> x1(npts, 0.0);
910  Array<OneD, NekDouble> x2(npts, 0.0);
911  Array<OneD, NekDouble> tmp(npts, 0.0);
912  Array<OneD, NekDouble> tmp0(npts, 0.0);
913  Array<OneD, NekDouble> tmp1(npts, 0.0);
914  for (int plane = 0; plane < npts; plane++)
915  {
916  x2[plane] = plane * DistStrip;
917  }
918  LibUtilities::EquationSharedPtr ffunc0, ffunc1;
919  ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
920  ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
921  ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
922  ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
923 
924  int offset = j * npts;
925  Vmath::Vcopy(npts, tmp0, 1, tmp = m_MotionVars[0] + offset,
926  1);
927  Vmath::Vcopy(npts, tmp1, 1, tmp = m_MotionVars[1] + offset,
928  1);
929  }
930  }
931 
932  cnt = cnt + 2;
933  }
934  }
935 }
NekDouble m_lhom
length ratio of the cable
NekDouble m_kinvis
fluid viscous
void SetDynEqCoeffMatrix(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
NekDouble m_structrho
mass of the cable per unit length
NekDouble m_timestep
time step
LibUtilities::NektarFFTSharedPtr m_FFT
NekDouble m_structdamp
damping ratio of the cable
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:130
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::LibUtilities::GetNektarFFTFactory(), m_fA, m_FFT, m_funcName, m_fV, m_IsFromFile, m_kinvis, m_lhom, m_motion, m_MotionVars, m_movingBodyCalls, m_np, Nektar::SolverUtils::Forcing::m_session, m_structdamp, m_structrho, m_timestep, m_vdim, Vmath::Sadd(), SetDynEqCoeffMatrix(), Vmath::Smul(), and Vmath::Vcopy().

Referenced by v_InitObject().

◆ InitialiseFilter()

void Nektar::ForcingMovingBody::InitialiseFilter ( const LibUtilities::SessionReaderSharedPtr pSession,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
const TiXmlElement *  pForce 
)
private

Definition at line 1213 of file ForcingMovingBody.cpp.

1217 {
1218  // Get the outputfile name, output frequency and
1219  // the boundary's ID for the cable's wall
1220  std::string typeStr = pForce->Attribute("TYPE");
1221  std::map<std::string, std::string> vParams;
1222 
1223  const TiXmlElement *param = pForce->FirstChildElement("PARAM");
1224  while (param)
1225  {
1226  ASSERTL0(param->Attribute("NAME"),
1227  "Missing attribute 'NAME' for parameter in filter " + typeStr +
1228  "'.");
1229  std::string nameStr = param->Attribute("NAME");
1230 
1231  ASSERTL0(param->GetText(), "Empty value string for param.");
1232  std::string valueStr = param->GetText();
1233 
1234  vParams[nameStr] = valueStr;
1235 
1236  param = param->NextSiblingElement("PARAM");
1237  }
1238 
1239  // Creat the filter for MovingBody, where we performed the calculation of
1240  // fluid forces and write both the aerodynamic forces and motion variables
1241  // into the output files
1243  pSession, m_equ, vParams);
1244 
1245  // Initialise the object of MovingBody filter
1246  m_MovBodyfilter->Initialise(pFields, 0.0);
1247 }
FilterMovingBodySharedPtr m_MovBodyfilter
const std::weak_ptr< EquationSystem > m_equ
Weak pointer to equation system using this forcing.
Definition: Forcing.h:118

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::SolverUtils::Forcing::m_equ, and m_MovBodyfilter.

Referenced by v_InitObject().

◆ Newmark_betaSolver()

void Nektar::ForcingMovingBody::Newmark_betaSolver ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
Array< OneD, NekDouble > &  FcePhysinArray,
Array< OneD, NekDouble > &  MotPhysinArray 
)
private

Definition at line 549 of file ForcingMovingBody.cpp.

552 {
553  std::string supptype = m_session->GetSolverInfo("SupportType");
554 
555  int npts = HydroForces.size();
556 
557  Array<OneD, Array<OneD, NekDouble>> fft_i(4);
558  Array<OneD, Array<OneD, NekDouble>> fft_o(4);
559 
560  for (int i = 0; i < 4; i++)
561  {
562  fft_i[i] = Array<OneD, NekDouble>(npts, 0.0);
563  fft_o[i] = Array<OneD, NekDouble>(npts, 0.0);
564  }
565 
566  Vmath::Vcopy(npts, HydroForces, 1, fft_i[0], 1);
567  for (int i = 0; i < 3; i++)
568  {
569  Vmath::Vcopy(npts, BodyMotions + i * npts, 1, fft_i[i + 1], 1);
570  }
571 
572  // Implement Fourier transformation of the motion variables
573  if (boost::iequals(supptype, "Free-Free"))
574  {
575  for (int j = 0; j < 4; ++j)
576  {
577  m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
578  }
579  }
580  else if (boost::iequals(supptype, "Pinned-Pinned"))
581  {
582  // TODO:
583  int N = fft_i[0].size();
584 
585  for (int var = 0; var < 4; var++)
586  {
587  for (int i = 0; i < N; i++)
588  {
589  fft_o[var][i] = 0;
590 
591  for (int k = 0; k < N; k++)
592  {
593  fft_o[var][i] +=
594  fft_i[var][k] * sin(M_PI / (N) * (k + 1 / 2) * (i + 1));
595  }
596  }
597  }
598  }
599  else
600  {
601  ASSERTL0(false, "Unrecognized support type for cable's motion");
602  }
603 
604  // solve the ODE in the wave space
605  for (int i = 0; i < npts; ++i)
606  {
607  int nrows = 3;
608 
609  Array<OneD, NekDouble> tmp0, tmp1, tmp2;
610  tmp0 = Array<OneD, NekDouble>(3, 0.0);
611  tmp1 = Array<OneD, NekDouble>(3, 0.0);
612  tmp2 = Array<OneD, NekDouble>(3, 0.0);
613 
614  for (int var = 0; var < 3; var++)
615  {
616  tmp0[var] = fft_o[var + 1][i];
617  }
618 
619  tmp2[0] = fft_o[0][i];
620 
621  Blas::Dgemv('N', nrows, nrows, 1.0, &(m_CoeffMat_B[i]->GetPtr())[0],
622  nrows, &tmp0[0], 1, 0.0, &tmp1[0], 1);
623  Blas::Dgemv('N', nrows, nrows, 1.0 / m_structrho,
624  &(m_CoeffMat_A[i]->GetPtr())[0], nrows, &tmp2[0], 1, 1.0,
625  &tmp1[0], 1);
626 
627  for (int var = 0; var < 3; var++)
628  {
629  fft_i[var][i] = tmp1[var];
630  }
631  }
632 
633  // get physical coeffients via Backward fourier transformation of wave
634  // coefficients
635  if (boost::iequals(supptype, "Free-Free"))
636  {
637  for (int var = 0; var < 3; var++)
638  {
639  m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
640  }
641  }
642  else if (boost::iequals(supptype, "Pinned-Pinned"))
643  {
644  // TODO:
645  int N = fft_i[0].size();
646 
647  for (int var = 0; var < 3; var++)
648  {
649  for (int i = 0; i < N; i++)
650  {
651  fft_o[var][i] = 0;
652 
653  for (int k = 0; k < N; k++)
654  {
655  fft_o[var][i] += fft_i[var][k] *
656  sin(M_PI / (N) * (k + 1) * (i + 1 / 2)) *
657  2 / N;
658  }
659  }
660  }
661  }
662  else
663  {
664  ASSERTL0(false, "Unrecognized support type for cable's motion");
665  }
666 
667  for (int i = 0; i < 3; i++)
668  {
669  Array<OneD, NekDouble> tmp(npts, 0.0);
670  Vmath::Vcopy(npts, fft_o[i], 1, tmp = BodyMotions + i * npts, 1);
671  }
672 }
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
matrices in Newmart-beta method
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:246

References ASSERTL0, Blas::Dgemv(), m_CoeffMat_A, m_CoeffMat_B, m_FFT, Nektar::SolverUtils::Forcing::m_session, m_structrho, and Vmath::Vcopy().

Referenced by EvaluateStructDynModel().

◆ RollOver()

void Nektar::ForcingMovingBody::RollOver ( Array< OneD, Array< OneD, NekDouble >> &  input)
private

Function to roll time-level storages to the next step layout. The stored data associated with the oldest time-level (not required anymore) are moved to the top, where they will be overwritten as the solution process progresses.

Definition at line 1045 of file ForcingMovingBody.cpp.

1046 {
1047  int nlevels = input.size();
1048 
1049  Array<OneD, NekDouble> tmp;
1050  tmp = input[nlevels - 1];
1051 
1052  for (int n = nlevels - 1; n > 0; --n)
1053  {
1054  input[n] = input[n - 1];
1055  }
1056 
1057  input[0] = tmp;
1058 }

Referenced by EvaluateStructDynModel().

◆ SetDynEqCoeffMatrix()

void Nektar::ForcingMovingBody::SetDynEqCoeffMatrix ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields)
private

Definition at line 940 of file ForcingMovingBody.cpp.

942 {
943  int nplanes;
944 
945  bool homostrip;
946  m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
947 
948  if (!homostrip)
949  {
950  nplanes = m_session->GetParameter("HomModesZ");
951  }
952  else
953  {
954  m_session->LoadParameter("Strip_Z", nplanes);
955  }
956 
957  m_CoeffMat_A = Array<OneD, DNekMatSharedPtr>(nplanes);
958  m_CoeffMat_B = Array<OneD, DNekMatSharedPtr>(nplanes);
959 
960  NekDouble tmp1, tmp2, tmp3;
961  NekDouble tmp4, tmp5, tmp6, tmp7;
962 
963  // load the structural dynamic parameters from xml file
964  NekDouble cabletension;
965  NekDouble bendingstiff;
966  NekDouble structstiff;
967  m_session->LoadParameter("StructStiff", structstiff, 0.0);
968  m_session->LoadParameter("CableTension", cabletension, 0.0);
969  m_session->LoadParameter("BendingStiff", bendingstiff, 0.0);
970 
971  tmp1 = m_timestep * m_timestep;
972  tmp2 = structstiff / m_structrho;
973  tmp3 = m_structdamp / m_structrho;
974  tmp4 = cabletension / m_structrho;
975  tmp5 = bendingstiff / m_structrho;
976 
977  // solve the ODE in the wave space for cable motion to obtain disp, vel and
978  // accel
979 
980  std::string supptype = m_session->GetSolverInfo("SupportType");
981 
982  for (int plane = 0; plane < nplanes; plane++)
983  {
984  int nel = 3;
985  m_CoeffMat_A[plane] =
987  m_CoeffMat_B[plane] =
989 
990  // Initialised to avoid compiler warnings.
991  unsigned int K = 0;
992  NekDouble beta = 0.0;
993 
994  if (boost::iequals(supptype, "Free-Free"))
995  {
996  K = plane / 2;
997  beta = 2.0 * M_PI / m_lhom;
998  }
999  else if (boost::iequals(supptype, "Pinned-Pinned"))
1000  {
1001  K = plane + 1;
1002  beta = M_PI / m_lhom;
1003  }
1004  else
1005  {
1006  ASSERTL0(false, "Unrecognized support type for cable's motion");
1007  }
1008 
1009  tmp6 = beta * K;
1010  tmp6 = tmp6 * tmp6;
1011  tmp7 = tmp6 * tmp6;
1012  tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
1013 
1014  (*m_CoeffMat_A[plane])(0, 0) = tmp7;
1015  (*m_CoeffMat_A[plane])(0, 1) = tmp3;
1016  (*m_CoeffMat_A[plane])(0, 2) = 1.0;
1017  (*m_CoeffMat_A[plane])(1, 0) = 0.0;
1018  (*m_CoeffMat_A[plane])(1, 1) = 1.0;
1019  (*m_CoeffMat_A[plane])(1, 2) = -m_timestep / 2.0;
1020  (*m_CoeffMat_A[plane])(2, 0) = 1.0;
1021  (*m_CoeffMat_A[plane])(2, 1) = 0.0;
1022  (*m_CoeffMat_A[plane])(2, 2) = -tmp1 / 4.0;
1023  (*m_CoeffMat_B[plane])(0, 0) = 0.0;
1024  (*m_CoeffMat_B[plane])(0, 1) = 0.0;
1025  (*m_CoeffMat_B[plane])(0, 2) = 0.0;
1026  (*m_CoeffMat_B[plane])(1, 0) = 0.0;
1027  (*m_CoeffMat_B[plane])(1, 1) = 1.0;
1028  (*m_CoeffMat_B[plane])(1, 2) = m_timestep / 2.0;
1029  (*m_CoeffMat_B[plane])(2, 0) = 1.0;
1030  (*m_CoeffMat_B[plane])(2, 1) = m_timestep;
1031  (*m_CoeffMat_B[plane])(2, 2) = tmp1 / 4.0;
1032 
1033  m_CoeffMat_A[plane]->Invert();
1034  (*m_CoeffMat_B[plane]) =
1035  (*m_CoeffMat_A[plane]) * (*m_CoeffMat_B[plane]);
1036  }
1037 }
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::beta, m_CoeffMat_A, m_CoeffMat_B, m_lhom, Nektar::SolverUtils::Forcing::m_session, m_structdamp, m_structrho, and m_timestep.

Referenced by InitialiseCableModel().

◆ v_Apply()

void Nektar::ForcingMovingBody::v_Apply ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble >> &  inarray,
Array< OneD, Array< OneD, NekDouble >> &  outarray,
const NekDouble time 
)
protectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 102 of file ForcingMovingBody.cpp.

106 {
107  // Update the forces from the calculation of fluid field, which is
108  // implemented in the movingbody filter
109  Array<OneD, NekDouble> Hydroforces(2 * m_np, 0.0);
110  m_MovBodyfilter->UpdateForce(m_session, pFields, Hydroforces, time);
111 
112  // for "free" type (m_vdim = 2), the cable vibrates both in streamwise and
113  // crossflow dimections, for "constrained" type (m_vdim = 1), the cable only
114  // vibrates in crossflow dimection, and for "forced" type (m_vdim = 0), the
115  // calbe vibrates specifically along a given function or file
116  if (m_vdim == 1 || m_vdim == 2)
117  {
118  // For free vibration case, displacements, velocities and acceleartions
119  // are obtained through solving structure dynamic model
120  EvaluateStructDynModel(pFields, Hydroforces, time);
121 
122  // Convert result to format required by mapping
123  int physTot = pFields[0]->GetTotPoints();
124  Array<OneD, Array<OneD, NekDouble>> coords(3);
125  Array<OneD, Array<OneD, NekDouble>> coordsVel(3);
126  for (int i = 0; i < 3; i++)
127  {
128  coords[i] = Array<OneD, NekDouble>(physTot, 0.0);
129  coordsVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
130  }
131  // Get original coordinates
132  pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
133 
134  // Add displacement to coordinates
135  Vmath::Vadd(physTot, coords[0], 1, m_zta[0], 1, coords[0], 1);
136  Vmath::Vadd(physTot, coords[1], 1, m_eta[0], 1, coords[1], 1);
137  // Copy velocities
138  Vmath::Vcopy(physTot, m_zta[1], 1, coordsVel[0], 1);
139  Vmath::Vcopy(physTot, m_eta[1], 1, coordsVel[1], 1);
140 
141  // Update mapping
142  m_mapping->UpdateMapping(time, coords, coordsVel);
143  }
144  else if (m_vdim == 0)
145  {
146  // For forced vibration case, load from specified file or function
147  int cnt = 0;
148  for (int j = 0; j < m_funcName.size(); j++)
149  {
150  if (m_IsFromFile[cnt] && m_IsFromFile[cnt + 1])
151  {
152  ASSERTL0(false, "Motion loading from file needs specific "
153  "implementation: Work in Progress!");
154  }
155  else
156  {
157  GetFunction(pFields, m_session, m_funcName[j], true)
158  ->Evaluate(m_motion[0], m_zta[j], time);
159  GetFunction(pFields, m_session, m_funcName[j], true)
160  ->Evaluate(m_motion[1], m_eta[j], time);
161  cnt = cnt + 2;
162  }
163  }
164 
165  // Update mapping
166  m_mapping->UpdateMapping(time);
167 
168  // Convert result from mapping
169  int physTot = pFields[0]->GetTotPoints();
170  Array<OneD, Array<OneD, NekDouble>> coords(3);
171  Array<OneD, Array<OneD, NekDouble>> coordsVel(3);
172  for (int i = 0; i < 3; i++)
173  {
174  coords[i] = Array<OneD, NekDouble>(physTot, 0.0);
175  coordsVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
176  }
177  // Get original coordinates
178  pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
179 
180  // Get Coordinates and coord velocity from mapping
181  m_mapping->GetCartesianCoordinates(m_zta[0], m_eta[0], coords[2]);
182  m_mapping->GetCoordVelocity(coordsVel);
183 
184  // Calculate displacement
185  Vmath::Vsub(physTot, m_zta[0], 1, coords[0], 1, m_zta[0], 1);
186  Vmath::Vsub(physTot, m_eta[0], 1, coords[1], 1, m_eta[0], 1);
187 
188  Vmath::Vcopy(physTot, coordsVel[0], 1, m_zta[1], 1);
189  Vmath::Vcopy(physTot, coordsVel[1], 1, m_eta[1], 1);
190 
191  for (int var = 0; var < 3; var++)
192  {
193  for (int plane = 0; plane < m_np; plane++)
194  {
195  int n = pFields[0]->GetPlane(plane)->GetTotPoints();
196  int offset = plane * n;
197  int Offset = var * m_np + plane;
198 
199  m_MotionVars[0][Offset] = m_zta[var][offset];
200  m_MotionVars[1][Offset] = m_eta[var][offset];
201  }
202  }
203  }
204  else
205  {
206  ASSERTL0(false, "Unrecogized vibration type for cable's dynamic model");
207  }
208 
209  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
210  int colrank = vcomm->GetColumnComm()->GetRank();
211  // Pass the variables of the cable's motion to the movingbody filter
212  if (colrank == 0)
213  {
214  int n = m_MotionVars[0].size();
215  Array<OneD, NekDouble> tmpArray(2 * n), tmp(n);
216  Vmath::Vcopy(n, m_MotionVars[0], 1, tmpArray, 1);
217  Vmath::Vcopy(n, m_MotionVars[1], 1, tmp = tmpArray + n, 1);
218  m_MovBodyfilter->UpdateMotion(m_session, pFields, tmpArray, time);
219  }
220 }
GlobalMapping::MappingSharedPtr m_mapping
void EvaluateStructDynModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Hydroforces, NekDouble time)
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, std::string pName, bool pCache=false)
Get a SessionFunction by name.
Definition: Forcing.cpp:190
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419

References ASSERTL0, EvaluateStructDynModel(), Nektar::SolverUtils::Forcing::GetFunction(), m_eta, m_funcName, m_IsFromFile, m_mapping, m_motion, m_MotionVars, m_MovBodyfilter, m_np, Nektar::SolverUtils::Forcing::m_session, m_vdim, m_zta, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vsub().

◆ v_InitObject()

void Nektar::ForcingMovingBody::v_InitObject ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
const unsigned int &  pNumForcingFields,
const TiXmlElement *  pForce 
)
protectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 55 of file ForcingMovingBody.cpp.

58 {
59  // Just 3D homogenous 1D problems can use this techinque
60  ASSERTL0(pFields[0]->GetExpType() == MultiRegions::e3DH1D,
61  "Moving body implemented just for 3D Homogenous 1D expansions.");
62 
63  // At this point we know in the xml file where those quantities
64  // are declared (equation or file) - via a function name which is now
65  // stored in funcNameD etc. We need now to fill in with this info the
66  // m_zta and m_eta vectors (actuallythey are matrices) Array to control if
67  // the motion is determined by an equation or is from a file.(not Nektar++)
68  // check if we need to load a file or we have an equation
69  CheckIsFromFile(pForce);
70 
71  // Initialise movingbody filter
72  InitialiseFilter(m_session, pFields, pForce);
73 
74  // Initialise the cable model
76 
77  // Load mapping
79  m_mapping->SetTimeDependent(true);
80 
81  if (m_vdim > 0)
82  {
83  m_mapping->SetFromFunction(false);
84  }
85 
86  m_zta = Array<OneD, Array<OneD, NekDouble>>(3);
87  m_eta = Array<OneD, Array<OneD, NekDouble>>(3);
88  // What are this bi-dimensional vectors
89  // ------------------------------------------ m_zta[0] = m_zta | m_eta[0] =
90  // m_eta | m_zta[1] = d(m_zta)/dt |
91  // m_eta[1] = d(m_eta)/dt | m_zta[2] = dd(m_zta)/ddtt |
92  // m_eta[2] = dd(m_eta)/ddtt |
93  //--------------------------------------------------------------------------------
94  int phystot = pFields[0]->GetTotPoints();
95  for (int i = 0; i < m_zta.size(); i++)
96  {
97  m_zta[i] = Array<OneD, NekDouble>(phystot, 0.0);
98  m_eta[i] = Array<OneD, NekDouble>(phystot, 0.0);
99  }
100 }
void CheckIsFromFile(const TiXmlElement *pForce)
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void InitialiseFilter(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:269

References ASSERTL0, CheckIsFromFile(), Nektar::MultiRegions::e3DH1D, InitialiseCableModel(), InitialiseFilter(), Nektar::GlobalMapping::Mapping::Load(), m_eta, m_mapping, Nektar::SolverUtils::Forcing::m_session, m_vdim, and m_zta.

Friends And Related Function Documentation

◆ MemoryManager< ForcingMovingBody >

friend class MemoryManager< ForcingMovingBody >
friend

Definition at line 1 of file ForcingMovingBody.h.

Member Data Documentation

◆ className

std::string Nektar::ForcingMovingBody::className
static
Initial value:
=
"MovingBody", ForcingMovingBody::create, "Moving Body Forcing")
static SolverUtils::ForcingSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
Creates an instance of this class.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44

Name of the class.

Definition at line 69 of file ForcingMovingBody.h.

◆ m_Aeroforces

Array<OneD, NekDouble> Nektar::ForcingMovingBody::m_Aeroforces
private

storage for the cable's force(x,y) variables

Definition at line 128 of file ForcingMovingBody.h.

◆ m_CoeffMat_A

Array<OneD, DNekMatSharedPtr> Nektar::ForcingMovingBody::m_CoeffMat_A
private

matrices in Newmart-beta method

Definition at line 136 of file ForcingMovingBody.h.

Referenced by Newmark_betaSolver(), and SetDynEqCoeffMatrix().

◆ m_CoeffMat_B

Array<OneD, DNekMatSharedPtr> Nektar::ForcingMovingBody::m_CoeffMat_B
private

matrices in Newmart-beta method

Definition at line 138 of file ForcingMovingBody.h.

Referenced by Newmark_betaSolver(), and SetDynEqCoeffMatrix().

◆ m_eta

Array<OneD, Array<OneD, NekDouble> > Nektar::ForcingMovingBody::m_eta
private

Store the derivatives of motion variables in y-direction.

Definition at line 148 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), v_Apply(), and v_InitObject().

◆ m_fA

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::ForcingMovingBody::m_fA
private

fictitious acceleration storage

Definition at line 134 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

◆ m_FFT

LibUtilities::NektarFFTSharedPtr Nektar::ForcingMovingBody::m_FFT
private

Definition at line 124 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and Newmark_betaSolver().

◆ m_funcName

Array<OneD, std::string> Nektar::ForcingMovingBody::m_funcName
private

[0] is displacements, [1] is velocities, [2] is accelerations

Definition at line 140 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), InitialiseCableModel(), and v_Apply().

◆ m_fV

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::ForcingMovingBody::m_fV
private

fictitious velocity storage

Definition at line 132 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

◆ m_IsFromFile

Array<OneD, bool> Nektar::ForcingMovingBody::m_IsFromFile
private

do determine if the the body motion come from an extern file

Definition at line 144 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), InitialiseCableModel(), and v_Apply().

◆ m_kinvis

NekDouble Nektar::ForcingMovingBody::m_kinvis
private

fluid viscous

Definition at line 121 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel().

◆ m_lhom

NekDouble Nektar::ForcingMovingBody::m_lhom
private

length ratio of the cable

Definition at line 120 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

◆ m_mapping

GlobalMapping::MappingSharedPtr Nektar::ForcingMovingBody::m_mapping
protected

Definition at line 73 of file ForcingMovingBody.h.

Referenced by v_Apply(), and v_InitObject().

◆ m_motion

Array<OneD, std::string> Nektar::ForcingMovingBody::m_motion
private

motion direction: [0] is 'x' and [1] is 'y'

Definition at line 142 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), EvaluateStructDynModel(), InitialiseCableModel(), and v_Apply().

◆ m_MotionVars

Array<OneD, Array<OneD, NekDouble> > Nektar::ForcingMovingBody::m_MotionVars
private

storage for the cable's motion(x,y) variables

Definition at line 130 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), and v_Apply().

◆ m_MovBodyfilter

FilterMovingBodySharedPtr Nektar::ForcingMovingBody::m_MovBodyfilter
private

Definition at line 126 of file ForcingMovingBody.h.

Referenced by InitialiseFilter(), and v_Apply().

◆ m_movingBodyCalls

int Nektar::ForcingMovingBody::m_movingBodyCalls
private

number of times the movbody have been called

Definition at line 114 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

◆ m_np

int Nektar::ForcingMovingBody::m_np
private

number of planes per processors

Definition at line 115 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), and v_Apply().

◆ m_outputFile_fce

std::string Nektar::ForcingMovingBody::m_outputFile_fce
private

Definition at line 152 of file ForcingMovingBody.h.

◆ m_outputFile_mot

std::string Nektar::ForcingMovingBody::m_outputFile_mot
private

Definition at line 153 of file ForcingMovingBody.h.

◆ m_outputFrequency

unsigned int Nektar::ForcingMovingBody::m_outputFrequency
private

Definition at line 150 of file ForcingMovingBody.h.

◆ m_outputStream

Array<OneD, std::ofstream> Nektar::ForcingMovingBody::m_outputStream
private

Definition at line 151 of file ForcingMovingBody.h.

◆ m_structdamp

NekDouble Nektar::ForcingMovingBody::m_structdamp
private

damping ratio of the cable

Definition at line 119 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

◆ m_structrho

NekDouble Nektar::ForcingMovingBody::m_structrho
private

mass of the cable per unit length

Definition at line 118 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), Newmark_betaSolver(), and SetDynEqCoeffMatrix().

◆ m_timestep

NekDouble Nektar::ForcingMovingBody::m_timestep
private

time step

Definition at line 122 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

◆ m_vdim

int Nektar::ForcingMovingBody::m_vdim
private

vibration dimension

Definition at line 116 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), v_Apply(), and v_InitObject().

◆ m_zta

Array<OneD, Array<OneD, NekDouble> > Nektar::ForcingMovingBody::m_zta
private

Store the derivatives of motion variables in x-direction.

Definition at line 146 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), v_Apply(), and v_InitObject().