MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
lsfv.cpp
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#include "lsfv.h"
8
9#include <algorithm>
10#include <stack>
11#include <vector>
12#include "FV/fvcartesiansolverxd.h"
13#include "LS/lscartesiansolver.h"
14#include "MEMORY/alloc.h"
15#include "UTIL/functions.h"
16#include "UTIL/kdtree.h"
17#include "globals.h"
18
19#include "globalvariables.h"
20
21#include "coupling.h"
22
23using namespace std;
24
25template <MInt nDim, class SysEqn>
26CouplingLsFv<nDim, SysEqn>::CouplingLsFv(const MInt couplingId, LsSolver* ls, FvCartesianSolver* fv)
27 : Coupling(couplingId), CouplingLS<nDim>(couplingId, ls), CouplingFv<nDim, SysEqn>(couplingId, fv) {
28 TRACE();
29
30 initData();
31 readProperties();
32 checkProperties();
33
34 for(MInt dir = 0; dir < nDim; dir++) {
35 lsSolver().a_meanCoord(dir) = a_meanCoord(dir);
36 }
37}
38
42template <MInt nDim, class SysEqn>
43void CouplingLsFv<nDim, SysEqn>::init() {
44 mAlloc(fvSolver().m_levelSetValues, lsSolver().m_noSets, "fvSolver().m_levelSetValues", AT_);
45 mAlloc(fvSolver().m_curvatureG, lsSolver().m_noSets, "fvSolver().m_curvatureG", AT_);
46
47 // set the time in the lsSolver:
48 lsSolver().m_time = -99;
49 lsSolver().m_time = a_time();
50 // transferGapCellProperty();
51 // computeGCellTimeStep();
52 return;
53}
54
58template <MInt nDim, class SysEqn>
59void CouplingLsFv<nDim, SysEqn>::finalizeCouplerInit() {
60 for(MInt s = 0; s < lsSolver().m_noSets; s++) {
61 fvSolver().m_levelSetValues[s].resize(fvSolver().a_noCells());
62 fvSolver().m_curvatureG[s].resize(fvSolver().a_noCells());
63 }
64
65 transferGapCellProperty();
66 transferLevelSetValues();
67 computeGCellTimeStep();
68}
69
73template <MInt nDim, class SysEqn>
74void CouplingLsFv<nDim, SysEqn>::preCouple(MInt /*recipeStep*/) {
75 if(!lsSolver().m_combustion && !lsSolver().m_levelSetMb) {
76 if(lsSolver().m_semiLagrange) {
77 returnStep_semiLagrange();
78 }
79 }
80}
81
85template <MInt nDim, class SysEqn>
86void CouplingLsFv<nDim, SysEqn>::postCouple(MInt /*recipeStep*/) {
87 lsSolver().m_time = a_time();
88
89 transferLevelSetValues();
90 transferGapCellProperty();
91 testCoupling();
92
93 if(!lsSolver().m_combustion && !lsSolver().m_levelSetMb && !lsSolver().m_semiLagrange && !lsSolver().m_LSSolver) {
94 returnStep();
95 }
96}
97
101template <MInt nDim, class SysEqn>
102void CouplingLsFv<nDim, SysEqn>::postAdaptation() {
103 TRACE();
104
105 for(MInt s = 0; s < lsSolver().m_noSets; s++) {
106 fvSolver().m_levelSetValues[s].clear();
107 fvSolver().m_levelSetValues[s].resize(fvSolver().a_noCells());
108 fvSolver().m_curvatureG[s].clear();
109 fvSolver().m_curvatureG[s].resize(fvSolver().a_noCells());
110 }
111
112 transferLevelSetValues();
113}
114
115
119template <MInt nDim, class SysEqn>
120void CouplingLsFv<nDim, SysEqn>::initData() {
121 TRACE();
122
123 // solver-specific data:
124 m_fvSolverId = fvSolver().m_solverId;
125 m_lsSolverId = lsSolver().m_solverId;
126
127 if(lsSolver().m_closeGaps) {
128 mAlloc(m_hadGapCells, lsSolver().m_noGapRegions, "m_hadGapCells", 0, AT_);
129 mAlloc(m_hasGapCells, lsSolver().m_noGapRegions, "m_hasGapCells", 0, AT_);
130 }
131}
132
136template <MInt nDim, class SysEqn>
137void CouplingLsFv<nDim, SysEqn>::checkProperties() {
138 TRACE();
139
140 ASSERT(lsSolver().a_maxGCellLevel() == fvSolver().maxRefinementLevel(), "");
141}
142
146template <MInt nDim, class SysEqn>
147void CouplingLsFv<nDim, SysEqn>::readProperties() {
148 TRACE();
149
150 // claudia: ugly hardcoded patch for fast checkup of level-set boundary motion - sorry...
162 m_maxVelocity = 14.5;
163 if(Context::propertyExists("maxLevelsetVelocityForTimestep", m_lsSolverId)) {
164 m_maxVelocity =
165 Context::getSolverProperty<MFloat>("maxLevelsetVelocityForTimestep", m_lsSolverId, AT_, &m_maxVelocity);
166 }
167
168 m_cfl = -1;
169 if(Context::propertyExists("cfl", m_lsSolverId)) {
170 m_cfl = Context::getSolverProperty<MFloat>("cfl", m_lsSolverId, AT_);
171 }
172
173 m_G0regionId = -1;
174 if(Context::propertyExists("G0regionId", m_fvSolverId)) {
175 m_G0regionId = Context::getSolverProperty<MInt>("G0regionId", m_fvSolverId, AT_, &m_G0regionId);
176 }
177
178 m_initialCrankAngle = F0;
179 if(Context::propertyExists("initialCrankAngle", m_fvSolverId)) {
180 m_initialCrankAngle =
181 Context::getSolverProperty<MFloat>("initialCrankAngle", m_fvSolverId, AT_, &m_initialCrankAngle);
182 }
183
184 m_timeStepMethod = 8;
185 if(Context::propertyExists("timeStepMethod", m_lsSolverId)) {
186 m_timeStepMethod = Context::getSolverProperty<MInt>("timeStepMethod", m_lsSolverId, AT_);
187 }
188
189 m_solverMethod = Context::getSolverProperty<MString>("solverMethod", m_lsSolverId, AT_);
190
201 m_bandWidthRef = 4;
202 m_bandWidthRefMax = 5;
203 if(Context::propertyExists("bandWidthRef", m_fvSolverId)) {
204 m_bandWidthRef = Context::getSolverProperty<MInt>("bandWidthRef", m_fvSolverId, AT_);
205 m_bandWidthRefMax = Context::getSolverProperty<MInt>("bandWidthRefMax", m_fvSolverId, AT_);
206 }
207}
208
209
213template <MInt nDim, class SysEqn>
214void CouplingLsFv<nDim, SysEqn>::computeGCellTimeStep() {
215 TRACE();
216
217 if(m_timeStepMethod < 6 || m_timeStepMethod > 8) {
218 if(string2enum(m_solverMethod) == MAIA_SEMI_LAGRANGE_LEVELSET
219 || string2enum(m_solverMethod) == MAIA_RUNGE_KUTTA_LEVELSET) {
220 mTerm(1, AT_, "Computation with pure LVS not possible with timeStepMethod other than 6 ,7 or 8! Please check!");
221 }
222 }
223
224 ASSERT(m_cfl > 0, "Couldn't read cfl property for the ls-solver timestepping!");
225
226 if(m_timeStepMethod == 8) {
227 lsSolver().m_timeStep = m_cfl * lsSolver().m_gCellDistance;
228
229 } else if(m_timeStepMethod == 6) {
230 lsSolver().m_timeStep = m_cfl * lsSolver().m_gCellDistance;
231 fvSolver().forceTimeStep(lsSolver().m_timeStep / a_timeRef());
232
233 } else if(m_timeStepMethod == 7) {
234 lsSolver().m_timeStep = m_cfl * lsSolver().m_gCellDistance / m_maxVelocity;
235 fvSolver().forceTimeStep(lsSolver().m_timeStep / a_timeRef());
236
237 } else {
238 lsSolver().m_timeStep = lsTimeStep();
239 }
240}
241
245template <MInt nDim, class SysEqn>
246MBool CouplingLsFv<nDim, SysEqn>::returnStep() {
247 TRACE();
248
249 fvSolver().m_time += lsTimeStep();
250 fvSolver().m_physicalTime += lsTimeStep() * a_timeRef();
251
252 return true;
253}
254
258template <MInt nDim, class SysEqn>
259void CouplingLsFv<nDim, SysEqn>::returnStep_semiLagrange() {
260 TRACE();
261
262 fvSolver().m_time += lsTimeStep() * a_timeRef();
263 fvSolver().m_physicalTime += lsTimeStep() * a_timeRef();
264}
265
266
267template <MInt nDim, class SysEqn>
268MFloat CouplingLsFv<nDim, SysEqn>::interpolateLevelSet(MInt cellId, MFloat* point, MInt set) {
269 MInt interpolationCells[8] = {0, 0, 0, 0, 0, 0, 0, 0};
270 MInt position = 0;
271
272 position = lsSolver().setUpLevelSetInterpolationStencil(cellId, interpolationCells, point);
273
274 // Interpolate level set
275 if(position > -1) {
276 return lsSolver().interpolateLevelSet(interpolationCells, point, set);
277 } else {
278 return lsSolver().a_levelSetFunctionG(cellId, set);
279 }
280}
281
282
283template <MInt nDim, class SysEqn>
284MInt CouplingLsFv<nDim, SysEqn>::noLevelSetFieldData() {
285 MInt noLevelSetFieldData = 0;
286 ASSERT(fvSolver().m_levelSet, "");
287
288 if(lsSolver().m_writeOutAllLevelSetFunctions) {
289 noLevelSetFieldData += noLevelSetFieldData + 2 * a_noSets();
290 } else {
291 noLevelSetFieldData += 1;
292 }
293 if(lsSolver().m_writeOutAllCurvatures) {
294 noLevelSetFieldData += a_noSets();
295 } else {
296 noLevelSetFieldData += 1;
297 }
298
299 return noLevelSetFieldData;
300}
301
302
303template <MInt nDim, class SysEqn>
304void CouplingLsFv<nDim, SysEqn>::setLsInList(MIntScratchSpace&) {
305 TRACE();
306
307 cerr0 << "Setting Sensors for Fv-Solver adaptation!" << endl;
308
309 MInt bandWidthRef = m_bandWidthRef;
310 MInt bandWidthRefmax = m_bandWidthRefMax;
311
312 MIntScratchSpace sendBufferSize(fvSolver().grid().noNeighborDomains(), AT_, "sendBufferSize");
313 MIntScratchSpace receiveBufferSize(fvSolver().grid().noNeighborDomains(), AT_, "receiveBufferSize");
314 MInt listCount = 0;
315 MIntScratchSpace inList(a_noFvGridCells(), AT_, "inList");
316 for(MInt c = 0; c < a_noFvGridCells(); c++)
317 inList[c] = 0;
318
319 // 2) add the g0-Cells and all parents to the list:
320 MInt endSet = a_noSets();
321 if(lsSolver().m_buildCollectedLevelSetFunction) {
322 endSet = 1;
323 }
324
325 for(MInt set = 0; set < endSet; set++) {
326 for(MInt id = 0; id < lsSolver().a_noG0Cells(set); id++) {
327 MInt fvCellId = ls2fvIdParent(a_G0CellId(id, set));
328 if(fvCellId < 0) continue;
329 if(fvSolver().a_isHalo(fvCellId)) continue;
330 inList[fvCellId] = 1;
331 listCount++;
332 MInt parentId = fvSolver().c_parentId(fvCellId);
333 while(parentId > -1) {
334 if(parentId < a_noFvGridCells()) inList[parentId] = 1;
335 parentId = fvSolver().c_parentId(parentId);
336 }
337 }
338 }
339
340 // Exchange the listCount on all Domains
341 MPI_Allreduce(&listCount, &listCount, 1, MPI_INT, MPI_SUM, fvSolver().mpiComm(), AT_, "listCount", "listCount");
342
343 if(listCount == 0) {
344 if(fvSolver().domainId() == 0) cerr << "No G0-Cells found!" << endl;
345 }
346
347 fvSolver().exchangeData(&inList[0], 1);
348
349 for(MInt level = fvSolver().minLevel(); level < fvSolver().maxRefinementLevel(); level++) {
350 if(level == fvSolver().maxRefinementLevel() - 1) bandWidthRef = bandWidthRefmax;
351
352 // 4) Loop over the number of revinement-grid-cells and add those to the list:
353 for(MInt loopMarker = 1; loopMarker < bandWidthRef; loopMarker++) {
354 for(MInt cellId = 0; cellId < a_noFvGridCells(); cellId++) {
355 if(fvSolver().a_level(cellId) != level) continue;
356 if(inList[cellId] != loopMarker) continue;
357
358 // direct neighbor
359 for(MInt n = 0; n < fvSolver().m_noDirs; n++) {
360 MInt nghbrId = fvSolver().c_neighborId(cellId, n, false);
361 if(nghbrId < 0) continue;
362 if(inList[nghbrId] == 0) inList[nghbrId] = loopMarker + 1;
363
364 // diagonal neighbors
365 if(n == 0 || n == 1) {
366 for(MInt nn = 2; nn < fvSolver().m_noDirs; nn++) {
367 MInt nghbrId2 = fvSolver().c_neighborId(nghbrId, nn, false);
368 if(nghbrId2 < 0) continue;
369 if(inList[nghbrId2] == 0) inList[nghbrId2] = loopMarker + 1;
370 }
371 }
372 if(n == 4 || n == 5) {
373 for(MInt nn = 2; nn < 4; nn++) {
374 MInt nghbrId2 = fvSolver().c_neighborId(nghbrId, nn, false);
375 if(nghbrId2 < 0) continue;
376 if(inList[nghbrId2] == 0) inList[nghbrId2] = loopMarker + 1;
377
378 for(MInt nnn = 0; nnn < 1; nnn++) {
379 MInt nghbrId3 = fvSolver().c_neighborId(nghbrId2, nnn, false);
380 if(nghbrId3 < 0) continue;
381 if(inList[nghbrId3] == 0) inList[nghbrId3] = loopMarker + 1;
382 }
383 }
384 }
385 }
386 }
387
388 fvSolver().exchangeData(&inList[0], 1);
389 }
390 }
391}
392
393
397template <MInt nDim, class SysEqn>
398void CouplingLsFv<nDim, SysEqn>::transferLevelSetValues() {
399 TRACE();
400
401 fvSolver().a_noSets() = a_noSets();
402 fvSolver().a_noLevelSetFieldData() = noLevelSetFieldData();
403
404 for(MInt s = 0; s < lsSolver().m_noSets; s++) {
405 fvSolver().m_levelSetValues[s].resize(fvSolver().a_noCells());
406 }
407
408 for(MInt s = 0; s < lsSolver().m_noSets; s++) {
409 for(MInt fc = 0; fc < a_noFvCells(); fc++) {
410 MInt lsId = fv2lsIdParent(fc);
411 fvSolver().a_levelSetFunction(fc, s) = a_outsideGValue();
412 if(lsId > -1) {
413 fvSolver().a_curvatureG(fc, s) = lsSolver().a_curvatureG(lsId, 0);
414 if(fvSolver().a_isInterface(fc)) {
415 // Interplate levelset for cut-cells
416 MFloat point[nDim];
417 for(MInt d = 0; d < nDim; d++) {
418 point[d] = fvSolver().a_coordinate(fc, d);
419 }
420 fvSolver().a_levelSetFunction(fc, s) = interpolateLevelSet(fv2lsIdParent(fc), point, s);
421 } else {
422 fvSolver().a_levelSetFunction(fc, s) = lsSolver().a_levelSetFunctionG(fv2lsIdParent(fc), s);
423 }
424 }
425 }
426 }
427}
428
432template <MInt nDim, class SysEqn>
433void CouplingLsFv<nDim, SysEqn>::transferGapCellProperty() {
434 TRACE();
435
436 if(!lsSolver().m_closeGaps) return;
437 if(lsSolver().m_noGapRegions <= 0) return;
438 if((lsSolver().m_levelSetMb) && fvSolver().m_constructGField) return;
439
440#ifdef COUPLING_DEBUG_
441 testCoupling();
442#endif
443
444 // 1) reset has/had Gap-Cells
445 for(MInt region = 0; region < lsSolver().m_noGapRegions; region++) {
446 m_hadGapCells[region] = 0;
447 m_hasGapCells[region] = 0;
448 }
449
450 // 2) set a_wasGapCell, m_hadGapCells and initialize a_isGapCell
451 for(MInt fc = 0; fc < a_noFvGridCells(); fc++) {
452 if(fvSolver().a_isGapCell(fc)) {
453 MInt lsId = fv2lsId(fc);
454 ASSERT(lsId > -1, "");
455 // skip G0-regions
456 if(m_G0regionId > -1 && a_potentialGapCellClose(lsId) == m_G0regionId) continue;
457 MInt regionId = a_potentialGapCellClose(lsId) - 2;
458 ASSERT(a_potentialGapCellClose(lsId) > 0 && a_potentialGapCellClose(lsId) <= lsSolver().m_noEmbeddedBodies, "");
459 if(regionId > -1 && regionId < lsSolver().m_noGapRegions && globalTimeStep > 0) {
460 // don't set hadGapCells during initialisation, so that initGapClosure is still called
461 // at the first timeStep!
462 m_hadGapCells[regionId]++;
463 }
464 }
465 if(globalTimeStep > 0) {
466 fvSolver().a_wasGapCell(fc) = fvSolver().a_isGapCell(fc);
467 } else {
468 fvSolver().a_wasGapCell(fc) = false;
469 }
470 fvSolver().a_isGapCell(fc) = false;
471 }
472
473 for(MInt fc = a_noFvGridCells(); fc < a_noFvCells(); fc++) {
474 if(globalTimeStep > 0) {
475 fvSolver().a_wasGapCell(fc) = fvSolver().a_isGapCell(fc);
476 } else {
477 fvSolver().a_wasGapCell(fc) = false;
478 }
479 fvSolver().a_isGapCell(fc) = false;
480 }
481
482 // 4) set a_isGapCell and hasGapCells
483 for(MInt gCellId = 0; gCellId < lsSolver().a_noCells(); gCellId++) {
484 if(a_nearGapG(gCellId)) {
485 MInt fvId = ls2fvId(gCellId);
486 if(fvId < 0) {
487 ASSERT(abs(abs(a_levelSetFunctionG(gCellId, 0)) - a_outsideGValue()) < 0.000001,
488 "ERROR, no fv-Cell found for relevant Gap-Cells!");
489 }
490 if(a_potentialGapCellClose(gCellId) > 0 && a_potentialGapCellClose(gCellId) <= lsSolver().m_noEmbeddedBodies) {
491 fvSolver().a_isGapCell(fvId) = true;
492 MInt regionId = a_potentialGapCellClose(gCellId) - 2;
493 if(regionId > -1 && regionId < lsSolver().m_noGapRegions) {
494 m_hasGapCells[regionId]++;
495 }
496 }
497 }
498 }
499 testGapProperty();
500
501 // set isGapCell for g0regionId:
502 if(m_G0regionId > -1) {
503 MInt noG0regionCells = 0;
504 for(MInt gCellId = 0; gCellId < lsSolver().a_noCells(); gCellId++) {
505 if(lsSolver().a_potentialGapCellClose(gCellId) == m_G0regionId && lsSolver().a_gapWidth(gCellId) > 0) {
506 MInt fvId = ls2fvId(gCellId);
507 if(fvId < 0) {
508 ASSERT(abs(abs(a_levelSetFunctionG(gCellId, 0)) - a_outsideGValue()) < 0.000001,
509 "ERROR, no fv-Cell found for relevant Gap-Cells!");
510 }
511 noG0regionCells++;
512 fvSolver().a_isGapCell(fvId) = true;
513 ASSERT(!lsSolver().a_nearGapG(gCellId), "");
514 }
515 }
516#if defined COUPLING_DEBUG_ || !defined NDEBUG
517 MPI_Allreduce(MPI_IN_PLACE, &noG0regionCells, 1, MPI_INT, MPI_SUM, fvSolver().mpiComm(), AT_, "INPLACE",
518 "noG0regionCells");
519 if(lsSolver().domainId() == 0) {
520 cerr << " No of g0-region Cells " << noG0regionCells << endl;
521 }
522#endif
523 }
524
525
526 // at the beginning of the restart
527 if(lsSolver().m_restart && globalTimeStep == fvSolver().m_restartTimeStep) {
528 for(MInt fc = a_noFvGridCells(); fc < a_noFvCells(); fc++) {
529 fvSolver().a_wasGapCell(fc) = fvSolver().a_isGapCell(fc);
530 }
531 for(MInt region = 0; region < lsSolver().m_noGapRegions; region++) {
532 m_hadGapCells[region] = m_hasGapCells[region];
533 }
534 }
535
536 // 3) exchange hadGapCells and hasGapCells for all regions
537 MPI_Allreduce(MPI_IN_PLACE, &m_hadGapCells[0], lsSolver().m_noGapRegions, MPI_INT, MPI_SUM, fvSolver().mpiComm(), AT_,
538 "INPLACE", "m_hadGapCells");
539
540 MPI_Allreduce(MPI_IN_PLACE, &m_hasGapCells[0], lsSolver().m_noGapRegions, MPI_INT, MPI_SUM, fvSolver().mpiComm(), AT_,
541 "INPLACE", "m_hasGapCells");
542
543 fvSolver().exchangeGapInfo();
544
545#if defined COUPLING_DEBUG_ || !defined NDEBUG
546 for(MInt region = 0; region < lsSolver().m_noGapRegions; region++) {
547 if(lsSolver().domainId() == 0) {
548 cerr << globalTimeStep << " region " << region << " has " << m_hasGapCells[region] << " had "
549 << m_hadGapCells[region] << endl;
550 }
551 }
552#endif
553}
554
555
559template <MInt nDim, class SysEqn>
560void CouplingLsFv<nDim, SysEqn>::testGapProperty() {
561 TRACE();
562
563 if(!lsSolver().m_closeGaps) return;
564 if((lsSolver().m_levelSetMb) && fvSolver().m_constructGField) return;
565
566#ifdef COUPLING_DEBUG_
567 for(MInt fc = 0; fc < a_noFvGridCells(); fc++) {
568 if(fvSolver().a_isGapCell(fc) || fvSolver().a_wasGapCell(fc)) {
569 MInt gc = lsSolver().grid().tree().grid2solver(fvSolver().grid().tree().solver2grid(fc));
570 ASSERT(gc > -1, "");
571 ASSERT(lsSolver().a_potentialGapCellClose(gc) > 0
572 && lsSolver().a_potentialGapCellClose(gc) <= fvSolver().m_noEmbeddedBodies,
573 "");
574 }
575 }
576#endif
577}
578
582template <MInt nDim, class SysEqn>
583void CouplingLsFv<nDim, SysEqn>::testCoupling() {
584 TRACE();
585
586 if(lsSolver().m_constructGField) return;
587
588 for(MInt fc = 0; fc < a_noFvGridCells(); fc++) {
589 ASSERT(fvSolver().grid().tree().solver2grid(fc) >= 0, "");
590 ASSERT(fvSolver().grid().solverFlag(fvSolver().grid().tree().solver2grid(fc), m_fvSolverId), "");
591#ifdef COUPLING_DEBUG_
592 if(lsSolver().grid().solverFlag(fvSolver().grid().tree().solver2grid(fc), m_lsSolverId)) {
593 ASSERT(ls2fvId(fv2lsId(fc)) == fc,
594 to_string(fc) + " " + to_string(fv2lsId(fc)) + " " + to_string(ls2fvId(fv2lsId(fc))));
595 }
596#endif
597 }
598 for(MInt cellId = 0; cellId < a_noLsCells(); cellId++) {
599 ASSERT(lsSolver().grid().tree().solver2grid(cellId) >= 0, "");
600 ASSERT(lsSolver().grid().solverFlag(lsSolver().grid().tree().solver2grid(cellId), m_lsSolverId), "");
601#ifdef COUPLING_DEBUG_
602 if(fvSolver().grid().solverFlag(lsSolver().grid().tree().solver2grid(cellId), m_fvSolverId)) {
603 ASSERT(fv2lsId(ls2fvId(cellId)) == cellId,
604 to_string(cellId) + " " + to_string(ls2fvId(cellId)) + " " + to_string(fv2lsId(ls2fvId(cellId))));
605 }
606#endif
607 }
608
609 if(!lsSolver().m_levelSetMb) return;
610
611
612#ifdef COUPLING_DEBUG_
613 for(MInt fc = 0; fc < a_noFvGridCells(); fc++) {
614 for(MInt dir = 0; dir < nDim; dir++) {
615 if(lsSolver().grid().solverFlag(fvSolver().grid().tree().solver2grid(fc), m_lsSolverId)) {
616 ASSERT(abs(fvSolver().c_coordinate(fc, dir) - a_coordinateG(fv2lsId(fc), dir)) < 0.00000001, "");
617 }
618 }
619 }
620#endif
621}
622
623
631template <MInt nDim, class SysEqn>
632void CouplingLsFv<nDim, SysEqn>::computeBodyProperties(MInt returnMode, MFloat* bodyData, MInt body, MFloat time) {
633 TRACE();
634
635 MFloat elapsedTime = time;
636 MFloat angle = F0;
637 MBool& first = m_static_computeBodyProperties_first;
638 MFloat(&amplitude)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_amplitude;
639 MFloat(&freqFactor)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_freqFactor;
640 MFloat(&initialBodyCenter)[m_maxNoEmbeddedBodies * 3] = m_static_computeBodyProperties_initialBodyCenter;
641 MFloat& Strouhal = m_static_computeBodyProperties_Strouhal;
642 MFloat(&mu)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_mu;
643 MFloat(&mu2)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_mu2;
644 MFloat(&liftStartAngle1)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_liftStartAngle1;
645 MFloat(&liftEndAngle1)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_liftEndAngle1;
646 MFloat(&liftStartAngle2)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_liftStartAngle2;
647 MFloat(&liftEndAngle2)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_liftEndAngle2;
648 MFloat(&circleStartAngle)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_circleStartAngle;
649 MFloat(&normal)[m_maxNoEmbeddedBodies * 3] = m_static_computeBodyProperties_normal;
650 MInt(&bodyToFunction)[m_maxNoEmbeddedBodies] = m_static_computeBodyProperties_bodyToFunction;
651 MFloat& omega = m_static_computeBodyProperties_omega;
652 MFloat& rotAngle = m_static_computeBodyProperties_rotAngle;
653
654 if(first) {
655 MInt noEmbeddedBodies = lsSolver().m_noEmbeddedBodies;
656
657 if(noEmbeddedBodies > m_maxNoEmbeddedBodies) {
658 mTerm(1, AT_, "Error in computeBodyProperties: too many embedded Bodies!");
659 }
660
661 // 1: set default values:
662 Strouhal = 0.2;
663 for(MInt k = 0; k < m_maxNoEmbeddedBodies; k++) {
664 amplitude[k] = 0.1;
665 freqFactor[k] = 1.0;
666 bodyToFunction[k] = 1;
667 for(MInt i = 0; i < nDim; i++) {
668 initialBodyCenter[k * nDim + i] = F0;
669 normal[k * nDim + i] = F0;
670 }
671 normal[k * nDim + 0] = 1.0;
672 liftStartAngle1[k] = F0;
673 liftEndAngle1[k] = PI;
674 liftStartAngle2[k] = 3.0 * PI;
675 liftEndAngle2[k] = 4.0 * PI;
676 circleStartAngle[k] = F0;
677 }
678
679 // 2: read Properties
680
693 for(MInt i = 0; i < noEmbeddedBodies; i++)
694 amplitude[i] = Context::getSolverProperty<MFloat>("amplitudes", m_lsSolverId, AT_, &amplitude[i], i);
695
708 for(MInt i = 0; i < noEmbeddedBodies; i++)
709 freqFactor[i] = Context::getSolverProperty<MFloat>("freqFactors", m_lsSolverId, AT_, &freqFactor[i], i);
710
723 for(MInt i = 0; i < noEmbeddedBodies; i++)
724 bodyToFunction[i] =
725 Context::getSolverProperty<MInt>("bodyMovementFunctions", m_lsSolverId, AT_, &bodyToFunction[i], i);
726
727 Strouhal = Context::getSolverProperty<MFloat>("Strouhal", m_lsSolverId, AT_, &Strouhal);
728
739 for(MInt i = 0; i < noEmbeddedBodies; i++)
740 for(MInt j = 0; j < nDim; j++)
741 initialBodyCenter[i * nDim + j] = Context::getSolverProperty<MFloat>(
742 "initialBodyCenters", m_lsSolverId, AT_, &initialBodyCenter[i * nDim + j], i * nDim + j);
743
744 for(MInt i = 0; i < noEmbeddedBodies; i++)
745 for(MInt j = 0; j < nDim; j++)
746 normal[i * nDim + j] = Context::getSolverProperty<MFloat>("bodyMotionNormals", m_lsSolverId, AT_,
747 &normal[i * nDim + j], i * nDim + j);
748
762 for(MInt i = 0; i < noEmbeddedBodies; i++)
763 liftStartAngle1[i] =
764 Context::getSolverProperty<MFloat>("liftStartAngles1", m_lsSolverId, AT_, &liftStartAngle1[i], i);
765
778 for(MInt i = 0; i < noEmbeddedBodies; i++)
779 liftStartAngle2[i] =
780 Context::getSolverProperty<MFloat>("liftStartAngles2", m_lsSolverId, AT_, &liftStartAngle2[i], i);
781
792 for(MInt i = 0; i < noEmbeddedBodies; i++)
793 liftEndAngle1[i] = Context::getSolverProperty<MFloat>("liftEndAngles1", m_lsSolverId, AT_, &liftEndAngle1[i], i);
794
805 for(MInt i = 0; i < noEmbeddedBodies; i++)
806 liftEndAngle2[i] = Context::getSolverProperty<MFloat>("liftEndAngles2", m_lsSolverId, AT_, &liftEndAngle2[i], i);
807
819 for(MInt i = 0; i < noEmbeddedBodies; i++)
820 circleStartAngle[i] =
821 Context::getSolverProperty<MFloat>("circleStartAngles", m_lsSolverId, AT_, &circleStartAngle[i], i);
822
835 rotAngle = 0.0;
836 rotAngle = Context::getSolverProperty<MFloat>("rotAngle", m_lsSolverId, AT_, &rotAngle);
837 rotAngle *= -PI / 180;
838
839 // 3: compute relevant values:
840 const MFloat freq0 = Strouhal * a_UInfinity();
841 const MFloat freq02 = Strouhal;
842 for(MInt k = 0; k < noEmbeddedBodies; k++) {
843 // when using mu : has a dimension!
844 // when using mu2 : dimensionless!
845 mu[k] = freqFactor[k] * freq0 * F2 * PI;
846 mu2[k] = freqFactor[k] * freq02 * F2 * PI;
847 }
848
849 // if bodyMovementFunction is 6 or 7, adjust start and end angles:
850
851 for(MInt i = 0; i < noEmbeddedBodies; i++) {
852 if(bodyToFunction[i] == 6 || bodyToFunction[i] == 7) {
853 liftStartAngle1[i] = liftStartAngle1[i] * PI;
854 liftEndAngle1[i] = liftEndAngle1[i] * PI - liftStartAngle1[i];
855 }
856 }
857
858 omega = freqFactor[body] * sqrt(a_TInfinity()) * a_Ma() * PI / (2.0 * amplitude[body]);
859
860 first = false;
861 }
862
863
864 //--------------------------------
865
866 switch(bodyToFunction[body]) {
867 case 1: // cosine function
868
869 angle = mu[body] * elapsedTime - liftStartAngle1[body];
870 // cerr << " time: " << elapsedTime << " angle: " << angle << endl;
871 switch(returnMode) {
872 case 1: // return body center
873 if(angle > liftEndAngle1[body]) angle = liftEndAngle1[body];
874 if(angle > 0) {
875 bodyData[0] = -amplitude[body] * cos(angle);
876 bodyData[1] = F0;
877 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
878 } else {
879 bodyData[0] = -amplitude[body];
880 bodyData[1] = F0;
881 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
882 }
883 break;
884 case 2: // return body velocity
885 if((angle > 0 && angle < liftEndAngle1[body])) {
886 bodyData[0] = mu[body] * amplitude[body] * sin(angle);
887 // cerr << " time pos vel body "<< body << " is: " << bodyData[0] << " " << bodyData[1] << endl;
888 bodyData[1] = F0;
889 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
890 } else {
891 bodyData[0] = F0;
892 // cerr<< " time neg vel body " << body << " is: " << bodyData[0] << " " << bodyData[1] << endl;
893 bodyData[1] = F0;
894 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
895 }
896 break;
897 case 3: // return body acceleration
898 if((angle > 0 && angle < liftEndAngle1[body])) {
899 bodyData[0] = mu[body] * mu[body] * amplitude[body] * cos(angle);
900 bodyData[1] = F0;
901 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
902 } else {
903 bodyData[0] = F0;
904 bodyData[1] = F0;
905 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
906 }
907 break;
908 default:
909 bodyData[0] = F0;
910 bodyData[1] = F0;
911 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
912 break;
913 }
914 break;
915
916 case 2: { // valve lift shifted quadratic sine
917
918 // cerr << " body position valve. body: " << body << endl;
919 angle = mu2[body] * elapsedTime;
920
921 switch(returnMode) {
922 case 1: // return body center
923 if((angle > liftStartAngle1[body] && angle <= liftEndAngle1[body])
924 || (angle > liftStartAngle2[body] && angle <= liftEndAngle2[body])) {
925 bodyData[0] = amplitude[body] * POW2(sin(mu2[body] * elapsedTime)) * normal[body * nDim + 0];
926 bodyData[1] = amplitude[body] * POW2(sin(mu2[body] * elapsedTime)) * normal[body * nDim + 1];
927 IF_CONSTEXPR(nDim == 3)
928 bodyData[2] = amplitude[body] * POW2(sin(mu2[body] * elapsedTime)) * normal[body * nDim + 2];
929 } else {
930 bodyData[0] = F0;
931 bodyData[1] = F0;
932 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
933 }
934 break;
935 case 2: // return body velocity
936 if((angle > liftStartAngle1[body] && angle <= liftEndAngle1[body])
937 || (angle > liftStartAngle2[body] && angle <= liftEndAngle2[body])) {
938 bodyData[0] = amplitude[body] * mu2[body] * sin(2 * mu2[body] * elapsedTime) * normal[body * nDim + 0];
939 bodyData[1] = amplitude[body] * mu2[body] * sin(2 * mu2[body] * elapsedTime) * normal[body * nDim + 1];
940 IF_CONSTEXPR(nDim == 3)
941 bodyData[2] = amplitude[body] * mu2[body] * sin(2 * mu2[body] * elapsedTime) * normal[body * nDim + 2];
942 } else {
943 bodyData[0] = F0;
944 bodyData[1] = F0;
945 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
946 }
947 break;
948 case 3: // return body acceleration
949 if((angle > liftStartAngle1[body] && angle <= liftEndAngle1[body])
950 || (angle > liftStartAngle2[body] && angle <= liftEndAngle2[body])) {
951 bodyData[0] = 2 * amplitude[body] * mu2[body] * mu2[body] * cos(2 * mu2[body] * elapsedTime)
952 * normal[body * nDim + 0];
953 bodyData[1] = 2 * amplitude[body] * mu2[body] * mu2[body] * cos(2 * mu2[body] * elapsedTime)
954 * normal[body * nDim + 1];
955 IF_CONSTEXPR(nDim == 3)
956 bodyData[2] = 2 * amplitude[body] * mu2[body] * mu2[body] * cos(2 * mu2[body] * elapsedTime)
957 * normal[body * nDim + 2];
958 } else {
959 bodyData[0] = 2 * amplitude[body] * mu2[body] * mu2[body] * normal[body * nDim + 0];
960 bodyData[1] = 2 * amplitude[body] * mu2[body] * mu2[body] * normal[body * nDim + 1];
961 IF_CONSTEXPR(nDim == 3) bodyData[2] = 2 * amplitude[body] * mu2[body] * mu2[body] * normal[body * nDim + 2];
962 }
963 break;
964 default:
965 bodyData[0] = F0;
966 bodyData[1] = F0;
967 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
968 break;
969 }
970
971 break;
972 }
973
974
975 case 3: // piston movement
976
977 // cerr << " body position piston. body: " << body << endl;
978 angle = mu2[body] * elapsedTime;
979
980
981 switch(returnMode) {
982 case 1: // return body center
983 if(elapsedTime > F0) {
984 bodyData[0] = -amplitude[body] * cos(mu2[body] * elapsedTime) * normal[body * nDim + 0];
985 bodyData[1] = -amplitude[body] * cos(mu2[body] * elapsedTime) * normal[body * nDim + 1];
986 IF_CONSTEXPR(nDim == 3)
987 bodyData[2] = -amplitude[body] * cos(mu2[body] * elapsedTime) * normal[body * nDim + 2];
988 } else {
989 bodyData[0] = -amplitude[body] * normal[body * nDim + 0];
990 bodyData[1] = -amplitude[body] * normal[body * nDim + 1];
991 IF_CONSTEXPR(nDim == 3) bodyData[2] = -amplitude[body] * normal[body * nDim + 2];
992 }
993 break;
994 case 2: // return body velocity
995 if(elapsedTime > F0) {
996 bodyData[0] = mu2[body] * amplitude[body] * sin(mu2[body] * elapsedTime) * normal[body * nDim + 0];
997 bodyData[1] = mu2[body] * amplitude[body] * sin(mu2[body] * elapsedTime) * normal[body * nDim + 1];
998 IF_CONSTEXPR(nDim == 3)
999 bodyData[2] = mu2[body] * amplitude[body] * sin(mu2[body] * elapsedTime) * normal[body * nDim + 2];
1000 } else {
1001 bodyData[0] = F0;
1002 bodyData[1] = F0;
1003 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1004 }
1005 break;
1006 case 3: // return body acceleration
1007 if(elapsedTime > F0) {
1008 bodyData[0] =
1009 mu2[body] * mu2[body] * amplitude[body] * cos(mu2[body] * elapsedTime) * normal[body * nDim + 0];
1010 bodyData[1] =
1011 mu2[body] * mu2[body] * amplitude[body] * cos(mu2[body] * elapsedTime) * normal[body * nDim + 1];
1012 IF_CONSTEXPR(nDim == 3)
1013 bodyData[2] =
1014 mu2[body] * mu2[body] * amplitude[body] * cos(mu2[body] * elapsedTime) * normal[body * nDim + 2];
1015 } else {
1016 bodyData[0] = F0;
1017 bodyData[1] = F0;
1018 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1019 }
1020 break;
1021 default:
1022 bodyData[0] = F0;
1023 bodyData[1] = F0;
1024 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1025 break;
1026 }
1027
1028 break;
1029
1030
1031 case 4: // cosine function with normal
1032
1033 angle = mu2[body] * elapsedTime;
1034
1035 switch(returnMode) {
1036 case 1: // return body center
1037 if((angle > liftStartAngle1[body] && angle <= liftEndAngle1[body])
1038 || (angle > liftStartAngle2[body] && angle <= liftEndAngle2[body])) {
1039 bodyData[0] = amplitude[body] * cos(angle) * normal[body * nDim + 0];
1040 bodyData[1] = amplitude[body] * cos(angle) * normal[body * nDim + 1];
1041 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * cos(angle) * normal[body * nDim + 2];
1042 } else {
1043 bodyData[0] = amplitude[body] * normal[body * nDim + 0];
1044 bodyData[1] = amplitude[body] * normal[body * nDim + 1];
1045 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * normal[body * nDim + 2];
1046 }
1047 break;
1048 case 2: // return body velocity
1049 if((angle > liftStartAngle1[body] && angle <= liftEndAngle1[body])
1050 || (angle > liftStartAngle2[body] && angle <= liftEndAngle2[body])) {
1051 bodyData[0] = -mu2[body] * amplitude[body] * sin(angle) * normal[body * nDim + 0];
1052 bodyData[1] = -mu2[body] * amplitude[body] * sin(angle) * normal[body * nDim + 1];
1053 IF_CONSTEXPR(nDim == 3) bodyData[2] = -mu2[body] * amplitude[body] * sin(angle) * normal[body * nDim + 2];
1054 } else {
1055 bodyData[0] = F0;
1056 bodyData[1] = F0;
1057 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1058 }
1059 break;
1060 case 3: // return body acceleration
1061 if((angle > liftStartAngle1[body] && angle <= liftEndAngle1[body])
1062 || (angle > liftStartAngle2[body] && angle <= liftEndAngle2[body])) {
1063 bodyData[0] = -mu2[body] * mu2[body] * amplitude[body] * cos(angle) * normal[body * nDim + 0];
1064 bodyData[1] = -mu2[body] * mu2[body] * amplitude[body] * cos(angle) * normal[body * nDim + 1];
1065 IF_CONSTEXPR(nDim == 3)
1066 bodyData[2] = -mu2[body] * mu2[body] * amplitude[body] * cos(angle) * normal[body * nDim + 2];
1067 } else {
1068 bodyData[0] = F0;
1069 bodyData[1] = F0;
1070 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1071 }
1072 break;
1073 default:
1074 bodyData[0] = F0;
1075 bodyData[1] = F0;
1076 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1077 break;
1078 }
1079 break;
1080
1081
1082 case 5:
1083 // circular motion:
1084 // initialBodyCenter: center of rotation
1085 // amplitude : radius
1086 // only 2d rotation implemented until now!
1087 angle = mu2[body] * elapsedTime + circleStartAngle[body];
1088
1089 switch(returnMode) {
1090 case 1: // return body center
1091 if(elapsedTime > F0) {
1092 bodyData[0] = amplitude[body] * cos(angle);
1093 bodyData[1] = amplitude[body] * sin(angle);
1094 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1095 } else {
1096 bodyData[0] = amplitude[body] * cos(circleStartAngle[body]);
1097 bodyData[1] = amplitude[body] * sin(circleStartAngle[body]);
1098 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1099 }
1100 break;
1101 case 2: // return body velocity
1102 if(elapsedTime > F0) {
1103 bodyData[0] = -amplitude[body] * sin(angle);
1104 bodyData[1] = amplitude[body] * cos(angle);
1105 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1106 } else {
1107 bodyData[0] = F0;
1108 bodyData[1] = F0;
1109 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1110 }
1111 break;
1112 case 3: // return body acceleration
1113 if(elapsedTime > F0) {
1114 bodyData[0] = -amplitude[body] * cos(angle);
1115 bodyData[1] = -amplitude[body] * sin(angle);
1116 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1117 } else {
1118 bodyData[0] = F0;
1119 bodyData[1] = F0;
1120 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1121 }
1122 break;
1123 default:
1124 bodyData[0] = F0;
1125 bodyData[1] = F0;
1126 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1127 break;
1128 }
1129 break;
1130
1131 case 6: // simplified piston motion
1132 {
1133 angle = omega * elapsedTime - liftStartAngle1[body];
1134 if(angle > liftEndAngle1[body]) angle = F0;
1135 // if( domainId() == 0 ) cerr << "angle for piston:" <<angle<< " elapsedTime: " << elapsedTime << endl;
1136
1137 switch(returnMode) {
1138 case 1: // return body center
1139 if(elapsedTime > F0) {
1140 bodyData[0] = -amplitude[body] * cos(angle) * normal[body * nDim + 0];
1141 bodyData[1] = -amplitude[body] * cos(angle) * normal[body * nDim + 1];
1142 IF_CONSTEXPR(nDim == 3) bodyData[2] = -amplitude[body] * cos(angle) * normal[body * nDim + 2];
1143 } else {
1144 bodyData[0] = -amplitude[body] * normal[body * nDim + 0];
1145 bodyData[1] = -amplitude[body] * normal[body * nDim + 1];
1146 IF_CONSTEXPR(nDim == 3) bodyData[2] = -amplitude[body] * normal[body * nDim + 2];
1147 }
1148 break;
1149 case 2: // return body velocity
1150 if(elapsedTime > F0) {
1151 bodyData[0] = omega * amplitude[body] * sin(angle) * normal[body * nDim + 0];
1152 bodyData[1] = omega * amplitude[body] * sin(angle) * normal[body * nDim + 1];
1153 IF_CONSTEXPR(nDim == 3) bodyData[2] = omega * amplitude[body] * sin(angle) * normal[body * nDim + 2];
1154 } else {
1155 bodyData[0] = F0;
1156 bodyData[1] = F0;
1157 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1158 }
1159 break;
1160 case 3: // return body acceleration
1161 if(elapsedTime > F0) {
1162 bodyData[0] = omega * omega * amplitude[body] * cos(angle) * normal[body * nDim + 0];
1163 bodyData[1] = omega * omega * amplitude[body] * cos(angle) * normal[body * nDim + 1];
1164 IF_CONSTEXPR(nDim == 3)
1165 bodyData[2] = omega * omega * amplitude[body] * cos(angle) * normal[body * nDim + 2];
1166 } else {
1167 bodyData[0] = F0;
1168 bodyData[1] = F0;
1169 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1170 }
1171 break;
1172 default:
1173 bodyData[0] = F0;
1174 bodyData[1] = F0;
1175 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1176 break;
1177 }
1178
1179
1180 break;
1181 }
1182
1183 case 7: // valve lift shifted quadratic sine matching piston motion 6
1184
1185 // cerr << " body position valve. body: " << body << endl;
1186 angle = omega * elapsedTime - liftStartAngle1[body];
1187 if(angle > liftEndAngle1[body]) angle = F0;
1188
1189 switch(returnMode) {
1190 case 1: // return body center
1191 if(angle > F0) {
1192 bodyData[0] = amplitude[body] * POW2(sin(angle)) * normal[body * nDim + 0];
1193 bodyData[1] = amplitude[body] * POW2(sin(angle)) * normal[body * nDim + 1];
1194 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * POW2(sin(angle)) * normal[body * nDim + 2];
1195 } else {
1196 bodyData[0] = 0;
1197 bodyData[1] = 0;
1198 IF_CONSTEXPR(nDim == 3) bodyData[2] = 0;
1199 }
1200 break;
1201 case 2: // return body velocity
1202 if(angle > F0) {
1203 bodyData[0] = amplitude[body] * omega * sin(2 * angle) * normal[body * nDim + 0];
1204 bodyData[1] = amplitude[body] * omega * sin(2 * angle) * normal[body * nDim + 1];
1205 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * omega * sin(2 * angle) * normal[body * nDim + 2];
1206 } else {
1207 bodyData[0] = F0;
1208 bodyData[1] = F0;
1209 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1210 }
1211 break;
1212 case 3: // return body acceleration
1213 if(angle > F0) {
1214 bodyData[0] = 2 * amplitude[body] * omega * omega * cos(2 * angle) * normal[body * nDim + 0];
1215 bodyData[1] = 2 * amplitude[body] * omega * omega * cos(2 * angle) * normal[body * nDim + 1];
1216 IF_CONSTEXPR(nDim == 3)
1217 bodyData[2] = 2 * amplitude[body] * omega * omega * cos(2 * angle) * normal[body * nDim + 2];
1218 } else {
1219 bodyData[0] = 2 * amplitude[body] * omega * omega * normal[body * nDim + 0];
1220 bodyData[1] = 2 * amplitude[body] * omega * omega * normal[body * nDim + 1];
1221 IF_CONSTEXPR(nDim == 3) bodyData[2] = 2 * amplitude[body] * omega * omega * normal[body * nDim + 2];
1222 }
1223 break;
1224 default:
1225 bodyData[0] = F0;
1226 bodyData[1] = F0;
1227 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1228 break;
1229 }
1230
1231 break;
1232
1233 case 8:
1234 // translational movement in normal-direction with constant Ma-number (Ma_trans)
1235 // Ma_trans is directly given by the amplitude-factor for each body!
1236 // The body velocity is then based on the free-stream speed of sound.
1237 angle = sqrt(a_TInfinity());
1238
1239 if(std::isnan(angle)) {
1240 cerr << "ERROR in the initialisation of the ls-Solver coupling class!" << endl;
1241 angle = 0;
1242 }
1243
1244 switch(returnMode) {
1245 case 1: // return body center
1246 bodyData[0] = amplitude[body] * angle * elapsedTime * normal[body * nDim + 0];
1247 bodyData[1] = amplitude[body] * angle * elapsedTime * normal[body * nDim + 1];
1248 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * angle * elapsedTime * normal[body * nDim + 2];
1249 break;
1250 case 2: // return body velocity
1251 bodyData[0] = amplitude[body] * angle * normal[body * nDim + 0];
1252 bodyData[1] = amplitude[body] * angle * normal[body * nDim + 1];
1253 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * angle * normal[body * nDim + 2];
1254 break;
1255 case 3: // return body acceleration
1256 bodyData[0] = F0;
1257 bodyData[1] = F0;
1258 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1259 break;
1260
1261 default:
1262 bodyData[0] = F0;
1263 bodyData[1] = F0;
1264 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1265 break;
1266 }
1267
1268 break;
1269
1270 case 9: // sine function with normel (=> periodic motion around the initialBodyCenter!)
1271
1272 angle = mu2[body] * elapsedTime;
1273
1274 switch(returnMode) {
1275 case 1: // return body center
1276 if(angle >= liftStartAngle1[body] && angle <= liftEndAngle1[body]) {
1277 bodyData[0] = amplitude[body] * sin(angle) * normal[body * nDim + 0];
1278 bodyData[1] = amplitude[body] * sin(angle) * normal[body * nDim + 1];
1279 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * sin(angle) * normal[body * nDim + 2];
1280 } else if(angle < liftStartAngle1[body]) {
1281 bodyData[0] = 0;
1282 bodyData[1] = 0;
1283 IF_CONSTEXPR(nDim == 3) bodyData[2] = 0;
1284 } else if(angle > liftEndAngle1[body]) {
1285 bodyData[0] = amplitude[body] * sin(liftEndAngle1[body]) * normal[body * nDim + 0];
1286 bodyData[1] = amplitude[body] * sin(liftEndAngle1[body]) * normal[body * nDim + 1];
1287 IF_CONSTEXPR(nDim == 3) bodyData[2] = amplitude[body] * sin(liftEndAngle1[body]) * normal[body * nDim + 2];
1288 }
1289 break;
1290 case 2: // return body velocity
1291 if(angle > liftStartAngle1[body] && angle <= liftEndAngle1[body]) {
1292 bodyData[0] = mu2[body] * amplitude[body] * cos(angle) * normal[body * nDim + 0];
1293 bodyData[1] = mu2[body] * amplitude[body] * cos(angle) * normal[body * nDim + 1];
1294 IF_CONSTEXPR(nDim == 3) bodyData[2] = mu2[body] * amplitude[body] * cos(angle) * normal[body * nDim + 2];
1295 } else {
1296 bodyData[0] = F0;
1297 bodyData[1] = F0;
1298 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1299 }
1300 break;
1301 case 3: // return body acceleration
1302 if(angle > liftStartAngle1[body] && angle <= liftEndAngle1[body]) {
1303 bodyData[0] = -mu2[body] * mu2[body] * amplitude[body] * sin(angle) * normal[body * nDim + 0];
1304 bodyData[1] = -mu2[body] * mu2[body] * amplitude[body] * sin(angle) * normal[body * nDim + 1];
1305 IF_CONSTEXPR(nDim == 3)
1306 bodyData[2] = -mu2[body] * mu2[body] * amplitude[body] * sin(angle) * normal[body * nDim + 2];
1307 } else {
1308 bodyData[0] = F0;
1309 bodyData[1] = F0;
1310 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1311 }
1312 break;
1313 default:
1314 bodyData[0] = F0;
1315 bodyData[1] = F0;
1316 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1317 break;
1318 }
1319 break;
1320
1321 case 10: {
1322 // piston motion equation:
1323 // reference coordincate system is the cylindeer head!
1324 // l : rod length
1325 // r : crank radius
1326 // TDC : distance at TDC from the cylinder head
1327 // normal : motion normal
1328 // phi : crank angle in radian
1329
1330 const MFloat l = liftStartAngle1[body];
1331 const MFloat TDC = liftEndAngle1[body];
1332 const MFloat r = amplitude[body];
1333
1334 MFloat phi = mu2[body] * elapsedTime;
1335
1336 // consider initial crank-angle:
1337 phi = phi + m_initialCrankAngle * PI / 180;
1338
1339 switch(returnMode) {
1340 case 1: // return body center
1341 bodyData[0] =
1342 normal[body * nDim + 0] * (l + r + TDC - (r * cos(phi) + sqrt(POW2(l) - POW2(r) * POW2(sin(phi)))));
1343 bodyData[1] =
1344 normal[body * nDim + 1] * (l + r + TDC - (r * cos(phi) + sqrt(POW2(l) - POW2(r) * POW2(sin(phi)))));
1345 IF_CONSTEXPR(nDim == 3)
1346 bodyData[2] =
1347 normal[body * nDim + 2] * (l + r + TDC - (r * cos(phi) + sqrt(POW2(l) - POW2(r) * POW2(sin(phi)))));
1348
1349 if(lsSolver().domainId() == 0 && fabs(elapsedTime - a_physicalTime()) < 0.0000000001) {
1350 const MFloat cad = crankAngle(a_physicalTime());
1351 cerr << "Crank-angle-degree : " << cad << endl;
1352 if((cad / 45) > 0.989 && (cad / 45) < 1.01) {
1353 cerr << "Physical-Time : " << a_physicalTime() << endl;
1354 cerr << " " << a_physicalTime() * 0.002898783653689 << endl;
1355 }
1356 cerr << "Piston-position : " << (l + r + TDC - (r * cos(phi) + sqrt(POW2(l) - POW2(r) * POW2(sin(phi)))))
1357 << " in m " << bodyData[1] * 0.075 << endl;
1358 cerr << "Piston-velocity : "
1359 << mu2[body]
1360 * (r * sin(phi) + (POW2(r) * sin(phi) * cos(phi) / (sqrt(POW2(l) - POW2(r) * POW2(sin(phi))))))
1361 << endl;
1362 }
1363
1364 break;
1365 case 2: // return body velocity
1366 bodyData[0] = normal[body * nDim + 0] * mu2[body]
1367 * (r * sin(phi) + (POW2(r) * sin(phi) * cos(phi) / (sqrt(POW2(l) - POW2(r) * POW2(sin(phi))))));
1368
1369 bodyData[1] = normal[body * nDim + 1] * mu2[body]
1370 * (r * sin(phi) + (POW2(r) * sin(phi) * cos(phi) / (sqrt(POW2(l) - POW2(r) * POW2(sin(phi))))));
1371 IF_CONSTEXPR(nDim == 3)
1372 bodyData[2] = normal[body * nDim + 2] * mu2[body]
1373 * (r * sin(phi) + (POW2(r) * sin(phi) * cos(phi) / (sqrt(POW2(l) - POW2(r) * POW2(sin(phi))))));
1374
1375 break;
1376 case 3: // return body acceleration
1377 bodyData[0] =
1378 normal[body * nDim + 0] * mu2[body] * mu2[body]
1379 * (r * cos(phi)
1380 - (POW2(r) * (POW2(cos(phi)) - POW2(sin(phi))) / (sqrt(POW2(l) - POW2(r) * POW2(sin(phi)))))
1381 + (POW4(r) * POW2(sin(phi)) * POW2(cos(phi))) / POW3((sqrt(POW2(l) - POW2(r) * POW2(sin(phi))))));
1382 bodyData[1] =
1383 normal[body * nDim + 1] * mu2[body] * mu2[body]
1384 * (r * cos(phi)
1385 - (POW2(r) * (POW2(cos(phi)) - POW2(sin(phi))) / (sqrt(POW2(l) - POW2(r) * POW2(sin(phi)))))
1386 + (POW4(r) * POW2(sin(phi)) * POW2(cos(phi))) / (POW3(sqrt(POW2(l) - POW2(r) * POW2(sin(phi))))));
1387 IF_CONSTEXPR(nDim == 3)
1388 bodyData[2] =
1389 normal[body * nDim + 2] * mu2[body] * mu2[body]
1390 * (r * cos(phi)
1391 - (POW2(r) * (POW2(cos(phi)) - POW2(sin(phi))) / (sqrt(POW2(l) - POW2(r) * POW2(sin(phi)))))
1392 + (POW4(r) * POW2(sin(phi)) * POW2(cos(phi))) / (POW3(sqrt(POW2(l) - POW2(r) * POW2(sin(phi))))));
1393 break;
1394 default:
1395 bodyData[0] = F0;
1396 bodyData[1] = F0;
1397 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1398 }
1399 break;
1400 }
1401 case 11: // valve motion for single-valve engine:
1402 {
1403 // reference coordincate system is the cylindeer head!
1404 // L : maximum valve lift (2*amplitude[body])
1405 // normal : motion normal
1406 // phi : crank angle in radian
1407 // cad : crank angle degree
1408 // maxNoCycles : maximum number of considered cycles
1409
1410
1411 const MFloat cad = crankAngle(elapsedTime);
1412
1413 MFloat IO = liftStartAngle1[body];
1414 MFloat IC = liftEndAngle1[body];
1415 MFloat EO = liftStartAngle2[body];
1416 MFloat EC = liftEndAngle2[body];
1417
1418 // possible values are:
1419 // IO: -3 // CAD BTDC 3
1420 // IC: 180+47 // CAD ABDC 47
1421 // EO: 3*180-47; // CAD BBDC 47
1422 // EC: 720+3; // CAD ATDC 3
1423
1424 switch(returnMode) {
1425 case 1: // return body center
1426 if(cad < IO) {
1427 bodyData[0] = F0;
1428 bodyData[1] = F0;
1429 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1430 } else if(cad >= 0 && cad < IC) {
1431 MFloat phi_i = (cad - IO) * PI / 180;
1432 MFloat delta_i = (IC - IO) * PI / 180;
1433 MFloat freq_i = 1 / (delta_i / 2 / PI);
1434 MFloat phase = -PI / 2;
1435
1436 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * (1 - sin(freq_i * phi_i - phase));
1437 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * (1 - sin(freq_i * phi_i - phase));
1438 IF_CONSTEXPR(nDim == 3)
1439 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * (1 - sin(freq_i * phi_i - phase));
1440
1441 } else if(cad > IC && cad < EO) {
1442 bodyData[0] = F0;
1443 bodyData[1] = F0;
1444 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1445 } else if(cad >= EO && cad < 720) {
1446 MFloat phi_e = (cad - EO) * PI / 180;
1447 MFloat delta_e = (EC - EO) * PI / 180;
1448 MFloat freq_e = 1 / (delta_e / 2 / PI);
1449 MFloat phase = -PI / 2;
1450
1451 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * (1 - sin(freq_e * phi_e - phase));
1452 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * (1 - sin(freq_e * phi_e - phase));
1453 IF_CONSTEXPR(nDim == 3)
1454 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * (1 - sin(freq_e * phi_e - phase));
1455 } else {
1456 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1457 }
1458
1459 if(lsSolver().domainId() == 0 && fabs(elapsedTime - a_physicalTime()) < 0.0000000001) {
1460 cerr << "Crank-angle-degree : " << cad << endl;
1461 cerr << "Valve-position : " << bodyData[0] << " in mm" << bodyData[0] * 75 << endl;
1462 }
1463
1464 break;
1465 case 2: // return body velocity
1466 if(cad < IO) {
1467 bodyData[0] = F0;
1468 bodyData[1] = F0;
1469 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1470 } else if(cad >= 0 && cad < IC) {
1471 MFloat phi_i = (cad - IO) * PI / 180;
1472 MFloat delta_i = (IC - IO) * PI / 180;
1473 MFloat freq_i = 1 / (delta_i / 2 / PI);
1474 MFloat phase = -PI / 2;
1475
1476 bodyData[0] = -normal[body * nDim + 0] * amplitude[body] * mu2[body] * freq_i * cos(freq_i * phi_i - phase);
1477 bodyData[1] = -normal[body * nDim + 1] * amplitude[body] * mu2[body] * freq_i * cos(freq_i * phi_i - phase);
1478 IF_CONSTEXPR(nDim == 3)
1479 bodyData[2] = -normal[body * nDim + 2] * amplitude[body] * mu2[body] * freq_i * cos(freq_i * phi_i - phase);
1480
1481 } else if(cad > IC && cad < EO) {
1482 bodyData[0] = F0;
1483 bodyData[1] = F0;
1484 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1485 } else if(cad >= EO && cad < 720) {
1486 MFloat phi_e = (cad - EO) * PI / 180;
1487 MFloat delta_e = (EC - EO) * PI / 180;
1488 MFloat freq_e = 1 / (delta_e / 2 / PI);
1489 MFloat phase = -PI / 2;
1490
1491 bodyData[0] = -normal[body * nDim + 0] * amplitude[body] * mu2[body] * freq_e * cos(freq_e * phi_e - phase);
1492 bodyData[1] = -normal[body * nDim + 1] * amplitude[body] * mu2[body] * freq_e * cos(freq_e * phi_e - phase);
1493 IF_CONSTEXPR(nDim == 3)
1494 bodyData[2] = -normal[body * nDim + 2] * amplitude[body] * mu2[body] * freq_e * cos(freq_e * phi_e - phase);
1495 } else {
1496 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1497 }
1498
1499 if(lsSolver().domainId() == 0 && fabs(elapsedTime - a_physicalTime()) < 0.0000000001) {
1500 cerr << "Valve-velocity : " << bodyData[0] << endl;
1501 }
1502
1503 break;
1504 case 3: // return body acceleration
1505 if(cad < IO) {
1506 bodyData[0] = F0;
1507 bodyData[1] = F0;
1508 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1509 } else if(cad >= 0 && cad < IC) {
1510 MFloat phi_i = (cad - IO) * PI / 180;
1511 MFloat delta_i = (IC - IO) * PI / 180;
1512 MFloat freq_i = 1 / (delta_i / 2 / PI);
1513 MFloat phase = -PI / 2;
1514
1515 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq_i)
1516 * sin(freq_i * phi_i - phase);
1517 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq_i)
1518 * sin(freq_i * phi_i - phase);
1519 IF_CONSTEXPR(nDim == 3)
1520 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq_i)
1521 * sin(freq_i * phi_i - phase);
1522
1523 } else if(cad > IC && cad < EO) {
1524 bodyData[0] = F0;
1525 bodyData[1] = F0;
1526 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1527 } else if(cad >= EO && cad < 720) {
1528 MFloat phi_e = (cad - EO) * PI / 180;
1529 MFloat delta_e = (EC - EO) * PI / 180;
1530 MFloat freq_e = 1 / (delta_e / 2 / PI);
1531 MFloat phase = -PI / 2;
1532
1533 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq_e)
1534 * sin(freq_e * phi_e - phase);
1535 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq_e)
1536 * sin(freq_e * phi_e - phase);
1537 IF_CONSTEXPR(nDim == 3)
1538 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq_e)
1539 * sin(freq_e * phi_e - phase);
1540 } else {
1541 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1542 }
1543 break;
1544 default:
1545 bodyData[0] = F0;
1546 bodyData[1] = F0;
1547 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1548 }
1549 break;
1550 }
1551 case 12: // valve motion for inlet valve:
1552 {
1553 // reference coordincate system is the cylindeer head!
1554 // amplitude : dimensionless maximum valve lift
1555 // normal : normalised motion normal
1556 // phi : crank angle in radian
1557 // cad : crank angle degree
1558 // maxNoCycles : maximum number of considered cycles
1559
1560 const MFloat cad = crankAngle(elapsedTime);
1561
1562 MFloat scaleToMeter = 0.075;
1563
1564 // valve-timing in crank-angle-degree
1565
1566 MFloat IO = liftStartAngle1[body];
1567 MFloat IC = liftEndAngle1[body];
1568
1569 // possible values are:
1570 // IO: -24 // 24 CAD BTDC
1571 // IC: 232 // CAD ATDC
1572
1573 MFloat phi = 0;
1574 MFloat delta = -1;
1575 MFloat freq = 0;
1576 MFloat phase = -PI / 2;
1577
1578 switch(returnMode) {
1579 case 1: // return body center
1580 if(cad < IO) {
1581 bodyData[0] = F0;
1582 bodyData[1] = F0;
1583 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1584 } else if(cad >= IO && cad <= IC) {
1585 phi = (cad - IO) * PI / 180;
1586 delta = (IC - IO) * PI / 180;
1587 freq = 1 / (delta / 2 / PI);
1588 phase = -PI / 2;
1589
1590 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * (1 - sin(freq * phi - phase));
1591 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * (1 - sin(freq * phi - phase));
1592 IF_CONSTEXPR(nDim == 3)
1593 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * (1 - sin(freq * phi - phase));
1594
1595 } else if(cad >= 720 + IO) {
1596 phi = (cad - (720 + IO)) * PI / 180;
1597 delta = (IC - IO) * PI / 180;
1598 freq = 1 / (delta / 2 / PI);
1599 phase = -PI / 2;
1600
1601 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * (1 - sin(freq * phi - phase));
1602 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * (1 - sin(freq * phi - phase));
1603 IF_CONSTEXPR(nDim == 3)
1604 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * (1 - sin(freq * phi - phase));
1605
1606 } else if(cad >= IC && cad < 720 + IO) {
1607 bodyData[0] = F0;
1608 bodyData[1] = F0;
1609 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1610 } else {
1611 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1612 }
1613
1614 if(lsSolver().domainId() == 0 && fabs(elapsedTime - a_physicalTime()) < 0.0000000001) {
1615 cerr << "Crank-angle-degree : " << cad << endl;
1616 cerr << "dimensionless-Inlet-Valve-position : " << amplitude[body] * (1 - sin(freq * phi - phase)) << endl;
1617 cerr << "dimensional-Inlet-Valve-position : "
1618 << amplitude[body] * (1 - sin(freq * phi - phase)) * scaleToMeter << " in m " << endl;
1619 cerr << "dimensionless-Inlet-Valve-velocity : "
1620 << -amplitude[body] * mu2[body] * freq * cos(freq * phi - phase) << endl;
1621 }
1622
1623 break;
1624 case 2: // return body velocity
1625 if(cad < IO) {
1626 bodyData[0] = F0;
1627 bodyData[1] = F0;
1628 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1629 } else if(cad >= IO && cad <= IC) {
1630 phi = (cad - IO) * PI / 180;
1631 delta = (IC - IO) * PI / 180;
1632 freq = 1 / (delta / 2 / PI);
1633 phase = -PI / 2;
1634
1635 bodyData[0] = -normal[body * nDim + 0] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1636 bodyData[1] = -normal[body * nDim + 1] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1637 IF_CONSTEXPR(nDim == 3)
1638 bodyData[2] = -normal[body * nDim + 2] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1639
1640 } else if(cad >= 720 + IO) {
1641 phi = (cad - (720 + IO)) * PI / 180;
1642 delta = (IC - IO) * PI / 180;
1643 freq = 1 / (delta / 2 / PI);
1644 phase = -PI / 2;
1645
1646 bodyData[0] = -normal[body * nDim + 0] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1647 bodyData[1] = -normal[body * nDim + 1] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1648 IF_CONSTEXPR(nDim == 3)
1649 bodyData[2] = -normal[body * nDim + 2] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1650
1651 } else if(cad >= IC && cad < 720 + IO) {
1652 bodyData[0] = F0;
1653 bodyData[1] = F0;
1654 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1655 } else {
1656 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1657 }
1658
1659 break;
1660 case 3: // return body acceleration
1661 if(cad < IO) {
1662 bodyData[0] = F0;
1663 bodyData[1] = F0;
1664 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1665 } else if(cad >= IO && cad <= IC) {
1666 phi = (cad - IO) * PI / 180;
1667 delta = (IC - IO) * PI / 180;
1668 freq = 1 / (delta / 2 / PI);
1669 phase = -PI / 2;
1670
1671 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1672 * sin(freq * phi - phase);
1673 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1674 * sin(freq * phi - phase);
1675 IF_CONSTEXPR(nDim == 3)
1676 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1677 * sin(freq * phi - phase);
1678
1679 } else if(cad >= 720 + IO) {
1680 phi = (cad - (720 + IO)) * PI / 180;
1681 delta = (IC - IO) * PI / 180;
1682 freq = 1 / (delta / 2 / PI);
1683 phase = -PI / 2;
1684
1685 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1686 * sin(freq * phi - phase);
1687 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1688 * sin(freq * phi - phase);
1689 IF_CONSTEXPR(nDim == 3)
1690 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1691 * sin(freq * phi - phase);
1692
1693 } else if(cad >= IC && cad < 720 + IO) {
1694 bodyData[0] = F0;
1695 bodyData[1] = F0;
1696 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1697 } else {
1698 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1699 }
1700 break;
1701 default:
1702 bodyData[0] = F0;
1703 bodyData[1] = F0;
1704 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1705 }
1706 break;
1707 }
1708 case 13: // valve motion for outlet valve:
1709 {
1710 // reference coordincate system is the cylindeer head!
1711 // amplitude : dimensionless maximum valve lift
1712 // normal : motion normal
1713 // phi : crank angle in radian
1714 // cad : crank angle degree
1715 // maxNoCycles : maximum number of considered cycles
1716 /*
1717 MFloat cad = mu2[body] * elapsedTime;
1718 MInt maxNoCycles = 20;
1719
1720 for (MInt cycle = maxNoCycles; cycle > 0; cycle-- ) {
1721 if(cad >= 4*PI*cycle) cad = cad - 4*PI*cycle ;
1722 }
1723 cad = cad *180/PI;
1724 */
1725 const MFloat cad = crankAngle(elapsedTime);
1726
1727 // valve-timing in crank-angle-degree
1728
1729 MFloat EO = liftStartAngle1[body];
1730 MFloat EC = liftEndAngle1[body];
1731 MFloat scaleToMeter = 0.075;
1732
1733 // possible values are:
1734 // EO: 480; // CAD ATDC
1735 // EC: 16; // CAD ATDC
1736
1737 MFloat phi = 0;
1738 MFloat delta = -1;
1739 MFloat freq = 0;
1740 MFloat phase = -PI / 2;
1741
1742 switch(returnMode) {
1743 case 1: // return body center
1744 if(cad < EO && cad >= EC) {
1745 bodyData[0] = F0;
1746 bodyData[1] = F0;
1747 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1748 } else if(cad >= EO && cad < 720) {
1749 phi = (cad - EO) * PI / 180;
1750 delta = (720 + EC - EO) * PI / 180;
1751 freq = 1 / (delta / 2 / PI);
1752 phase = -PI / 2;
1753
1754 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * (1 - sin(freq * phi - phase));
1755 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * (1 - sin(freq * phi - phase));
1756 IF_CONSTEXPR(nDim == 3)
1757 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * (1 - sin(freq * phi - phase));
1758
1759 } else if(cad < EC) {
1760 phi = -(EC - cad) * PI / 180;
1761 delta = (720 + EC - EO) * PI / 180;
1762 freq = 1 / (delta / 2 / PI);
1763 phase = -PI / 2;
1764
1765 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * (1 - sin(freq * phi - phase));
1766 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * (1 - sin(freq * phi - phase));
1767 IF_CONSTEXPR(nDim == 3)
1768 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * (1 - sin(freq * phi - phase));
1769
1770 } else {
1771 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1772 }
1773
1774 if(lsSolver().domainId() == 0 && fabs(elapsedTime - a_physicalTime()) < 0.0000000001) {
1775 cerr << "Crank-angle-degree : " << cad << endl;
1776 cerr << "dimensionless-Outlet-Valve-position : " << amplitude[body] * (1 - sin(freq * phi - phase)) << endl;
1777 cerr << "dimensional-Outlet-Valve-position : "
1778 << amplitude[body] * (1 - sin(freq * phi - phase)) * scaleToMeter << " in m " << endl;
1779 cerr << "dimensionless-Outlet-Valve-velocity : "
1780 << -amplitude[body] * mu2[body] * freq * cos(freq * phi - phase) << endl;
1781 }
1782
1783 break;
1784 case 2: // return body velocity
1785 if(cad < EO && cad >= EC) {
1786 bodyData[0] = F0;
1787 bodyData[1] = F0;
1788 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1789 } else if(cad >= EO && cad < 720) {
1790 phi = (cad - EO) * PI / 180;
1791 delta = (720 + EC - EO) * PI / 180;
1792 freq = 1 / (delta / 2 / PI);
1793 phase = -PI / 2;
1794
1795 bodyData[0] = -normal[body * nDim + 0] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1796 bodyData[1] = -normal[body * nDim + 1] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1797 IF_CONSTEXPR(nDim == 3)
1798 bodyData[2] = -normal[body * nDim + 2] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1799
1800
1801 } else if(cad < EC) {
1802 phi = -(EC - cad) * PI / 180;
1803 delta = (720 + EC - EO) * PI / 180;
1804 freq = 1 / (delta / 2 / PI);
1805 phase = -PI / 2;
1806
1807 bodyData[0] = -normal[body * nDim + 0] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1808 bodyData[1] = -normal[body * nDim + 1] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1809 IF_CONSTEXPR(nDim == 3)
1810 bodyData[2] = -normal[body * nDim + 2] * amplitude[body] * mu2[body] * freq * cos(freq * phi - phase);
1811 } else {
1812 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1813 }
1814
1815 break;
1816 case 3: // return body acceleration
1817 if(cad < EO && cad >= EC) {
1818 bodyData[0] = F0;
1819 bodyData[1] = F0;
1820 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1821 } else if(cad >= EO && cad < 720) {
1822 phi = (cad - EO) * PI / 180;
1823 delta = (720 + EC - EO) * PI / 180;
1824 freq = 1 / (delta / 2 / PI);
1825 phase = -PI / 2;
1826
1827 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1828 * sin(freq * phi - phase);
1829 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1830 * sin(freq * phi - phase);
1831 IF_CONSTEXPR(nDim == 3)
1832 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1833 * sin(freq * phi - phase);
1834
1835 } else if(cad < EC) {
1836 phi = -(EC - cad) * PI / 180;
1837 delta = (720 + EC - EO) * PI / 180;
1838 freq = 1 / (delta / 2 / PI);
1839 phase = -PI / 2;
1840
1841 bodyData[0] = normal[body * nDim + 0] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1842 * sin(freq * phi - phase);
1843 bodyData[1] = normal[body * nDim + 1] * amplitude[body] * mu2[body] * mu2[body] * POW2(freq)
1844 * sin(freq * phi - phase);
1845 IF_CONSTEXPR(nDim == 3)
1846 bodyData[2] = normal[body * nDim + 2] * amplitude[body] * mu2[body] * mu2[body]
1847 * (1 + POW2(freq) * sin(freq * phi - phase));
1848 } else {
1849 mTerm(1, AT_, "Unexpected crank-angle-degree!");
1850 }
1851 break;
1852 default:
1853 bodyData[0] = F0;
1854 bodyData[1] = F0;
1855 IF_CONSTEXPR(nDim == 3) bodyData[2] = F0;
1856 }
1857 break;
1858 }
1859 default:
1860 mTerm(1, AT_, "function type not implemented. Please check bodyMovementFunctions property!");
1861 }
1862
1863 // add the initialBodyCenter to the bodyDato for the body-positioning:
1864 if(returnMode == 1) {
1865 for(MInt dir = 0; dir < nDim; dir++) {
1866 bodyData[dir] += initialBodyCenter[body * nDim + dir];
1867 }
1868 }
1869
1870 // rotate the final result around the z-axis
1871
1872 if(rotAngle > 0 || rotAngle < 0) {
1873 MFloat tmp0 = bodyData[0] * cos(rotAngle) + bodyData[1] * sin(rotAngle);
1874 MFloat tmp1 = bodyData[1] * cos(rotAngle) - bodyData[0] * sin(rotAngle);
1875 bodyData[0] = tmp0;
1876 bodyData[1] = tmp1;
1877 bodyData[2] = bodyData[2];
1878 }
1879}
1880
1885template <MInt nDim, class SysEqn>
1886MFloat CouplingLsFv<nDim, SysEqn>::crankAngle(MFloat elapsedTime) {
1887 TRACE();
1888
1889 MFloat& Strouhal = m_static_crankAngle_Strouhal;
1890 MBool& first = m_static_crankAngle_first;
1891
1892 if(first) {
1893 Strouhal = Context::getSolverProperty<MFloat>("Strouhal", m_lsSolverId, AT_, &Strouhal);
1894 first = false;
1895 }
1896
1897 const MFloat mu2 = Strouhal * F2 * PI;
1898 MFloat cad = mu2 * elapsedTime;
1899 const MInt maxNoCycles = 20;
1900
1901 for(MInt cycle = maxNoCycles; cycle > 0; cycle--) {
1902 if(cad >= 4 * PI * cycle) cad = cad - 4 * PI * cycle;
1903 }
1904 cad = cad * 180 / PI;
1905
1906 // consider initial crank angle
1907 cad = cad + m_initialCrankAngle;
1908
1909 return cad;
1910}
1911
1912template class CouplingLsFv<2, FvSysEqnNS<2>>;
1913template class CouplingLsFv<3, FvSysEqnNS<3>>;
1914template class CouplingLsFv<2, FvSysEqnRANS<2, RANSModelConstants<RANS_SA_DV>>>;
1915template class CouplingLsFv<3, FvSysEqnRANS<3, RANSModelConstants<RANS_SA_DV>>>;
1916template class CouplingLsFv<2, FvSysEqnRANS<2, RANSModelConstants<RANS_FS>>>;
1917template class CouplingLsFv<3, FvSysEqnRANS<3, RANSModelConstants<RANS_FS>>>;
1918template class CouplingLsFv<2, FvSysEqnRANS<2, RANSModelConstants<RANS_KOMEGA>>>;
1919template class CouplingLsFv<3, FvSysEqnRANS<3, RANSModelConstants<RANS_KOMEGA>>>;
mAlloc
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
Definition: alloc.h:173
Context::propertyExists
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
Definition: context.cpp:494
CouplingLS< nDim_ >::lsSolver
solverType & lsSolver() const
Definition: coupling.h:188
CouplingLsFv::interpolateLevelSet
MFloat interpolateLevelSet(MInt cellId, MFloat *point, MInt set)
Definition: lsfv.cpp:268
CouplingLsFv::finalizeCouplerInit
void finalizeCouplerInit() override
Definition: lsfv.cpp:59
CouplingLsFv::checkProperties
void checkProperties() override
Checks property-data which is read in by both ls-and Fv-Solver.
Definition: lsfv.cpp:137
CouplingLsFv::testGapProperty
void testGapProperty()
transfers the LevelSetValues from the levelset to the moving boundary Part
Definition: lsfv.cpp:560
CouplingLsFv::returnStep
MBool returnStep()
mimics the behaviour of the rungeKuttaStep() methods with respect to increasing time
Definition: lsfv.cpp:246
CouplingLsFv::init
void init() override
Definition: lsfv.cpp:43
CouplingLsFv::computeBodyProperties
void computeBodyProperties(MInt returnMode, MFloat *bodyData, MInt body, MFloat time)
CouplingLsFv::setLsInList
void setLsInList(MIntScratchSpace &)
Definition: lsfv.cpp:304
CouplingLsFv::nDim
static constexpr MInt nDim
Definition: lsfv.h:29
CouplingLsFv::a_meanCoord
MFloat a_meanCoord(MInt dir) const
Definition: lsfv.h:156
CouplingLsFv::returnStep_semiLagrange
void returnStep_semiLagrange()
mimics the behaviour of the rungeKuttaStep() methods with respect to increasing time
Definition: lsfv.cpp:259
CouplingLsFv::computeGCellTimeStep
void computeGCellTimeStep()
computes the gcell time step
Definition: lsfv.cpp:214
CouplingLsFv::transferGapCellProperty
void transferGapCellProperty()
Sets the gapcell-property.
Definition: lsfv.cpp:433
CouplingLsFv::readProperties
void readProperties() override
CouplingLsFv::noLevelSetFieldData
MInt noLevelSetFieldData()
Definition: lsfv.cpp:284
CouplingLsFv::testCoupling
void testCoupling()
transfers the LevelSetValues from the levelset to the moving boundary Part
Definition: lsfv.cpp:583
CouplingLsFv::preCouple
void preCouple(MInt) override
Definition: lsfv.cpp:74
CouplingLsFv::transferLevelSetValues
void transferLevelSetValues()
Sets the Levelset-Values in fvSolver.
Definition: lsfv.cpp:398
CouplingLsFv::postAdaptation
void postAdaptation() override
finalizeAdaptation
Definition: lsfv.cpp:102
CouplingLsFv::postCouple
void postCouple(MInt recipeStep=0) override
Definition: lsfv.cpp:86
CouplingLsFv::initData
void initData()
Initialize coupling-class-specific Data.
Definition: lsfv.cpp:120
CouplingLsFv::crankAngle
MFloat crankAngle(MFloat)
help-function for engine calculations which returns the crank-angle for a given time
Definition: lsfv.cpp:1886
CouplingLsFv::CouplingLsFv
CouplingLsFv(const MInt couplingId, LsSolver *ls, FvCartesianSolver *fv)
Definition: lsfv.cpp:26
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
string2enum
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
MAIA_RUNGE_KUTTA_LEVELSET
@ MAIA_RUNGE_KUTTA_LEVELSET
Definition: enums.h:66
MAIA_SEMI_LAGRANGE_LEVELSET
@ MAIA_SEMI_LAGRANGE_LEVELSET
Definition: enums.h:67
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
POW4
constexpr Real POW4(const Real x)
Definition: functions.h:127
POW3
constexpr Real POW3(const Real x)
Definition: functions.h:123
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119
globalTimeStep
MInt globalTimeStep
Definition: globalvariables.cpp:33
cerr0
std::ostream cerr0
MInt
int32_t MInt
Definition: maiatypes.h:62
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
MPI_Allreduce
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
Definition: mpioverride.cpp:952
maia::filter::slope::cos
T cos(const T a, const T b, const T x)
Cosine slope filter.
Definition: filter.h:125
maia::math::crankAngle
MFloat crankAngle(const MFloat time, const MFloat Strouhal, const MFloat offset, const MInt mode)
help-function for engine calculations which returns the crank-angle for a given time,...
Definition: maiamath.h:506