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

void v_InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
 
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) override
 
- 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 ~Forcing ()=default
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)=0
 
virtual SOLVER_UTILS_EXPORT 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)=0
 
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

size_t m_movingBodyCalls
 number of times the movbody have been called More...
 
size_t m_np
 number of planes per processors More...
 
size_t 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
 

Friends

class MemoryManager< ForcingMovingBody >
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::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 47 of file ForcingMovingBody.cpp.

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

Member Function Documentation

◆ CheckIsFromFile()

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

Definition at line 1065 of file ForcingMovingBody.cpp.

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

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:53

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 226 of file ForcingMovingBody.cpp.

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

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 678 of file ForcingMovingBody.cpp.

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

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 1215 of file ForcingMovingBody.cpp.

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

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 550 of file ForcingMovingBody.cpp.

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

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 1047 of file ForcingMovingBody.cpp.

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

Referenced by EvaluateStructDynModel().

◆ SetDynEqCoeffMatrix()

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

Definition at line 942 of file ForcingMovingBody.cpp.

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

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 
)
overrideprotectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 102 of file ForcingMovingBody.cpp.

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

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 
)
overrideprotectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 54 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 size_t phystot = pFields[0]->GetTotPoints();
95 for (size_t 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:264

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.
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:41

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 129 of file ForcingMovingBody.h.

◆ m_CoeffMat_A

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

matrices in Newmart-beta method

Definition at line 137 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 139 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 149 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 135 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

◆ m_FFT

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

Definition at line 125 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 141 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 133 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 145 of file ForcingMovingBody.h.

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

◆ m_kinvis

NekDouble Nektar::ForcingMovingBody::m_kinvis
private

fluid viscous

Definition at line 122 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel().

◆ m_lhom

NekDouble Nektar::ForcingMovingBody::m_lhom
private

length ratio of the cable

Definition at line 121 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 143 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 131 of file ForcingMovingBody.h.

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

◆ m_MovBodyfilter

FilterMovingBodySharedPtr Nektar::ForcingMovingBody::m_MovBodyfilter
private

Definition at line 127 of file ForcingMovingBody.h.

Referenced by InitialiseFilter(), and v_Apply().

◆ m_movingBodyCalls

size_t Nektar::ForcingMovingBody::m_movingBodyCalls
private

number of times the movbody have been called

Definition at line 115 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

◆ m_np

size_t Nektar::ForcingMovingBody::m_np
private

number of planes per processors

Definition at line 116 of file ForcingMovingBody.h.

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

◆ m_outputFrequency

unsigned int Nektar::ForcingMovingBody::m_outputFrequency
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 120 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 119 of file ForcingMovingBody.h.

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

◆ m_timestep

NekDouble Nektar::ForcingMovingBody::m_timestep
private

time step

Definition at line 123 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

◆ m_vdim

size_t Nektar::ForcingMovingBody::m_vdim
private

vibration dimension

Definition at line 117 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 147 of file ForcingMovingBody.h.

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