64template <
class T,
int D>
class Advect {
65 using ConstSparseIterator =
66 viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType>;
67 using hrleIndexType = viennahrle::IndexType;
69 std::vector<SmartPointer<Domain<T, D>>> levelSets;
70 SmartPointer<VelocityField<T>> velocities =
nullptr;
73 double timeStepRatio = 0.4999;
74 double dissipationAlpha = 1.0;
75 bool calculateNormalVectors =
true;
76 bool ignoreVoids =
false;
77 double advectionTime = 0.;
78 bool performOnlySingleStep =
false;
79 double advectedTime = 0.;
80 unsigned numberOfTimeSteps = 0;
81 bool saveAdvectionVelocities =
false;
82 bool updatePointData =
true;
83 bool checkDissipation =
true;
84 static constexpr double wrappingLayerEpsilon = 1e-4;
88 std::vector<std::vector<std::pair<std::pair<T, T>,
T>>> storedRates;
89 double currentTimeStep = -1.;
91 VectorType<T, 3> findGlobalAlphas()
const {
93 auto &topDomain = levelSets.back()->getDomain();
94 auto &grid = levelSets.back()->getGrid();
96 const T gridDelta = grid.getGridDelta();
97 const T deltaPos = gridDelta;
98 const T deltaNeg = -gridDelta;
100 VectorType<T, 3> finalAlphas = {0., 0., 0.};
102#pragma omp parallel num_threads((levelSets.back())->getNumberOfSegments())
104 VectorType<T, 3> localAlphas = {0., 0., 0.};
107 p = omp_get_thread_num();
109 viennahrle::Index<D> startVector =
110 (p == 0) ? grid.getMinGridPoint()
111 : topDomain.getSegmentation()[p - 1];
112 viennahrle::Index<D> endVector =
113 (p !=
static_cast<int>(topDomain.getNumberOfSegments() - 1))
114 ? topDomain.getSegmentation()[p]
115 : grid.incrementIndices(grid.getMaxGridPoint());
118 std::vector<ConstSparseIterator> iterators;
119 for (
auto const &
ls : levelSets) {
120 iterators.emplace_back(
ls->getDomain());
124 viennahrle::ConstSparseStarIterator<typename Domain<T, D>::DomainType, 1>
125 neighborIterator(topDomain);
127 for (ConstSparseIterator it(topDomain, startVector);
128 it.getStartIndices() < endVector; ++it) {
130 if (!it.isDefined() || std::abs(it.getValue()) > 0.5)
133 const T value = it.getValue();
134 const auto indices = it.getStartIndices();
138 for (
unsigned lowerLevelSetId = 0; lowerLevelSetId < levelSets.size();
141 iterators[lowerLevelSetId].goToIndicesSequential(indices);
145 if (iterators[lowerLevelSetId].getValue() <=
146 value + wrappingLayerEpsilon) {
149 neighborIterator.goToIndicesSequential(indices);
152 for (
unsigned i = 0; i <
D; ++i) {
153 coords[i] = indices[i] * gridDelta;
156 Vec3D<T> normal = {};
157 T normalModulus = 0.;
158 for (
unsigned i = 0; i <
D; ++i) {
159 const T phiPos = neighborIterator.getNeighbor(i).getValue();
160 const T phiNeg = neighborIterator.getNeighbor(i +
D).getValue();
162 T diffPos = (phiPos - value) / deltaPos;
163 T diffNeg = (phiNeg - value) / deltaNeg;
165 normal[i] = (diffNeg + diffPos) * 0.5;
166 normalModulus += normal[i] * normal[i];
168 normalModulus = std::sqrt(normalModulus);
169 for (
unsigned i = 0; i <
D; ++i)
170 normal[i] /= normalModulus;
172 T scaVel = velocities->getScalarVelocity(
173 coords, lowerLevelSetId, normal,
174 neighborIterator.getCenter().getPointId());
175 auto vecVel = velocities->getVectorVelocity(
176 coords, lowerLevelSetId, normal,
177 neighborIterator.getCenter().getPointId());
179 for (
unsigned i = 0; i <
D; ++i) {
180 T tempAlpha = std::abs((scaVel + vecVel[i]) * normal[i]);
181 localAlphas[i] = std::max(localAlphas[i], tempAlpha);
192 for (
unsigned i = 0; i <
D; ++i) {
193 finalAlphas[i] = std::max(finalAlphas[i], localAlphas[i]);
206 auto &grid = levelSets.back()->getGrid();
207 auto newlsDomain = SmartPointer<Domain<T, D>>::New(grid);
208 auto &newDomain = newlsDomain->getDomain();
209 auto &domain = levelSets.back()->getDomain();
211 newDomain.initialize(domain.getNewSegmentation(),
212 domain.getAllocation() *
213 (2.0 / levelSets.back()->getLevelSetWidth()));
215 const bool updateData = updatePointData;
218 std::vector<std::vector<unsigned>> newDataSourceIds;
220 newDataSourceIds.resize(newDomain.getNumberOfSegments());
222#ifdef DEBUG_LS_ADVECT_HPP
224 auto mesh = SmartPointer<Mesh<T>>::New();
230#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
234 p = omp_get_thread_num();
236 auto &domainSegment = newDomain.getDomainSegment(p);
238 viennahrle::Index<D> startVector =
239 (p == 0) ? grid.getMinGridPoint()
240 : newDomain.getSegmentation()[p - 1];
242 viennahrle::Index<D> endVector =
243 (p !=
static_cast<int>(newDomain.getNumberOfSegments() - 1))
244 ? newDomain.getSegmentation()[p]
245 : grid.incrementIndices(grid.getMaxGridPoint());
250 newDataSourceIds[p].reserve(2.5 * domainSegment.getNumberOfPoints());
253 it(domain, startVector);
254 it.getIndices() < endVector; ++it) {
258 if (std::abs(it.getCenter().getValue()) <= 1.0) {
261 for (; k < 2 *
D; k++)
262 if (std::signbit(it.getNeighbor(k).getValue() - 1e-7) !=
263 std::signbit(it.getCenter().getValue() + 1e-7))
268 if (it.getCenter().getDefinedValue() > 0.5) {
270 for (; j < 2 *
D; j++) {
271 if (std::abs(it.getNeighbor(j).getValue()) <= 1.0)
272 if (it.getNeighbor(j).getDefinedValue() < -0.5)
276 domainSegment.insertNextDefinedPoint(
277 it.getIndices(), it.getCenter().getDefinedValue());
279 newDataSourceIds[p].push_back(it.getCenter().getPointId());
282 domainSegment.insertNextDefinedPoint(it.getIndices(), 0.5);
284 newDataSourceIds[p].push_back(it.getNeighbor(j).getPointId());
286 }
else if (it.getCenter().getDefinedValue() < -0.5) {
288 for (; j < 2 *
D; j++) {
289 if (std::abs(it.getNeighbor(j).getValue()) <= 1.0)
290 if (it.getNeighbor(j).getDefinedValue() > 0.5)
295 domainSegment.insertNextDefinedPoint(
296 it.getIndices(), it.getCenter().getDefinedValue());
298 newDataSourceIds[p].push_back(it.getCenter().getPointId());
301 domainSegment.insertNextDefinedPoint(it.getIndices(), -0.5);
303 newDataSourceIds[p].push_back(it.getNeighbor(j).getPointId());
306 domainSegment.insertNextDefinedPoint(
307 it.getIndices(), it.getCenter().getDefinedValue());
309 newDataSourceIds[p].push_back(it.getCenter().getPointId());
312 domainSegment.insertNextUndefinedPoint(
313 it.getIndices(), (it.getCenter().getDefinedValue() < 0)
319 if (it.getCenter().getValue() >= 0) {
320 int usedNeighbor = -1;
322 for (
int i = 0; i < 2 *
D; i++) {
323 T value = it.getNeighbor(i).getValue();
324 if (std::abs(value) <= 1.0 && (value < 0.)) {
325 if (distance > value + 1.0) {
326 distance = value + 1.0;
332 if (distance <= 1.) {
333 domainSegment.insertNextDefinedPoint(it.getIndices(), distance);
335 newDataSourceIds[p].push_back(
336 it.getNeighbor(usedNeighbor).getPointId());
338 domainSegment.insertNextUndefinedPoint(it.getIndices(),
343 int usedNeighbor = -1;
345 for (
int i = 0; i < 2 *
D; i++) {
346 T value = it.getNeighbor(i).getValue();
347 if (std::abs(value) <= 1.0 && (value > 0)) {
348 if (distance < value - 1.0) {
350 distance = value - 1.0;
356 if (distance >= -1.) {
357 domainSegment.insertNextDefinedPoint(it.getIndices(), distance);
359 newDataSourceIds[p].push_back(
360 it.getNeighbor(usedNeighbor).getPointId());
362 domainSegment.insertNextUndefinedPoint(it.getIndices(),
372 auto &pointData = levelSets.back()->getPointData();
373 newlsDomain->getPointData().translateFromMultiData(pointData,
377 newDomain.finalize();
379 levelSets.back()->deepCopy(newlsDomain);
380 levelSets.back()->finalize(2);
385 double advect(
double maxTimeStep) {
387 if (levelSets.empty()) {
388 Logger::getInstance()
389 .addWarning(
"No level sets passed to Advect. Not advecting.")
391 return std::numeric_limits<double>::max();
393 if (velocities ==
nullptr) {
394 Logger::getInstance()
395 .addWarning(
"No velocity field passed to Advect. Not advecting.")
397 return std::numeric_limits<double>::max();
400 if (currentTimeStep < 0. || storedRates.empty())
412 if (integrationScheme !=
414 for (
unsigned i = 0; i < levelSets.size() - 1; ++i) {
427 template <
class IntegrationSchemeType>
428 double integrateTime(IntegrationSchemeType IntegrationScheme,
429 double maxTimeStep) {
431 auto &topDomain = levelSets.back()->getDomain();
432 auto &grid = levelSets.back()->getGrid();
437 auto &pointData = levelSets.back()->getPointData();
440 if (voidMarkerPointer ==
nullptr) {
441 Logger::getInstance()
442 .addWarning(
"Advect: Cannot find void point markers. Not "
443 "ignoring void points.")
448 const bool ignoreVoidPoints = ignoreVoids;
450 if (!storedRates.empty()) {
451 Logger::getInstance()
452 .addWarning(
"Advect: Overwriting previously stored rates.")
456 storedRates.resize(topDomain.getNumberOfSegments());
458#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
462 p = omp_get_thread_num();
464 viennahrle::Index<D> startVector =
465 (p == 0) ? grid.getMinGridPoint()
466 : topDomain.getSegmentation()[p - 1];
468 viennahrle::Index<D> endVector =
469 (p !=
static_cast<int>(topDomain.getNumberOfSegments() - 1))
470 ? topDomain.getSegmentation()[p]
471 : grid.incrementIndices(grid.getMaxGridPoint());
473 double tempMaxTimeStep = maxTimeStep;
474 auto &tempRates = storedRates[p];
476 topDomain.getNumberOfPoints() /
477 static_cast<double>((levelSets.back())->getNumberOfSegments()) +
481 std::vector<ConstSparseIterator> iterators;
482 for (
auto const &
ls : levelSets) {
483 iterators.emplace_back(
ls->getDomain());
486 IntegrationSchemeType scheme(IntegrationScheme);
488 for (ConstSparseIterator it(topDomain, startVector);
489 it.getStartIndices() < endVector; ++it) {
491 if (!it.isDefined() || std::abs(it.getValue()) > 0.5)
494 T value = it.getValue();
495 double maxStepTime = 0;
496 double cfl = timeStepRatio;
498 for (
int currentLevelSetId = levelSets.size() - 1;
499 currentLevelSetId >= 0; --currentLevelSetId) {
501 std::pair<T, T> gradNDissipation;
503 if (!(ignoreVoidPoints && (*voidMarkerPointer)[it.getPointId()])) {
506 for (
unsigned lowerLevelSetId = 0;
507 lowerLevelSetId < levelSets.size(); ++lowerLevelSetId) {
509 iterators[lowerLevelSetId].goToIndicesSequential(
510 it.getStartIndices());
514 if (iterators[lowerLevelSetId].getValue() <=
515 value + wrappingLayerEpsilon) {
517 scheme(it.getStartIndices(), lowerLevelSetId);
525 if (currentLevelSetId > 0) {
526 iterators[currentLevelSetId - 1].goToIndicesSequential(
527 it.getStartIndices());
528 valueBelow = iterators[currentLevelSetId - 1].getValue();
530 valueBelow = std::numeric_limits<T>::max();
535 T velocity = gradNDissipation.first - gradNDissipation.second;
537 maxStepTime += cfl / velocity;
538 tempRates.push_back(std::make_pair(gradNDissipation,
539 -std::numeric_limits<T>::max()));
542 }
else if (velocity == 0.) {
543 maxStepTime = std::numeric_limits<T>::max();
544 tempRates.push_back(std::make_pair(gradNDissipation,
545 std::numeric_limits<T>::max()));
550 T difference = std::abs(valueBelow - value);
552 if (difference >= cfl) {
553 maxStepTime -= cfl / velocity;
554 tempRates.push_back(std::make_pair(
555 gradNDissipation, std::numeric_limits<T>::max()));
558 maxStepTime -= difference / velocity;
561 tempRates.push_back(std::make_pair(gradNDissipation, valueBelow));
569 if (maxStepTime < tempMaxTimeStep)
570 tempMaxTimeStep = maxStepTime;
578 scheme.reduceTimeStepHamiltonJacobi(
579 tempMaxTimeStep, levelSets.back()->getGrid().getGridDelta());
582 if (tempMaxTimeStep < maxTimeStep)
583 maxTimeStep = tempMaxTimeStep;
595 if (timeStepRatio >= 0.5) {
596 Logger::getInstance()
597 .addWarning(
"Integration time step ratio should be smaller than 0.5. "
598 "Advection might fail!")
602 auto &topDomain = levelSets.back()->getDomain();
603 auto &grid = levelSets.back()->getGrid();
605 assert(currentTimeStep >= 0. &&
"No time step set!");
606 assert(storedRates.size() == topDomain.getNumberOfSegments());
612 const bool saveVelocities = saveAdvectionVelocities;
613 std::vector<std::vector<double>> dissipationVectors(
614 levelSets.back()->getNumberOfSegments());
615 std::vector<std::vector<double>> velocityVectors(
616 levelSets.back()->getNumberOfSegments());
618 const bool checkDiss = checkDissipation;
620#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
624 p = omp_get_thread_num();
626 auto itRS = storedRates[p].cbegin();
627 auto &segment = topDomain.getDomainSegment(p);
628 const unsigned maxId = segment.getNumberOfPoints();
630 if (saveVelocities) {
631 velocityVectors[p].resize(maxId);
632 dissipationVectors[p].resize(maxId);
635 for (
unsigned localId = 0; localId < maxId; ++localId) {
636 T &value = segment.definedValues[localId];
637 double time = currentTimeStep;
642 auto const [gradient, dissipation] = itRS->first;
643 T velocity = gradient - dissipation;
645 if (checkDiss && (gradient < 0 && velocity > 0) ||
646 (gradient > 0 && velocity < 0)) {
650 T rate = time * velocity;
651 while (std::abs(itRS->second - value) < std::abs(rate)) {
652 time -= std::abs((itRS->second - value) / velocity);
653 value = itRS->second;
657 velocity = itRS->first.first - itRS->first.second;
658 if (checkDiss && (itRS->first.first < 0 && velocity > 0) ||
659 (itRS->first.first > 0 && velocity < 0)) {
662 rate = time * velocity;
668 if (saveVelocities) {
669 velocityVectors[p][localId] = rate;
670 dissipationVectors[p][localId] = itRS->first.second;
676 while (std::abs(itRS->second) != std::numeric_limits<T>::max())
684 if (saveVelocities) {
685 auto &pointData = levelSets.back()->getPointData();
690 for (
unsigned i = 0; i < velocityVectors.size(); ++i) {
691 vels.insert(vels.end(),
692 std::make_move_iterator(velocityVectors[i].begin()),
693 std::make_move_iterator(velocityVectors[i].end()));
694 diss.insert(diss.end(),
695 std::make_move_iterator(dissipationVectors[i].begin()),
696 std::make_move_iterator(dissipationVectors[i].end()));
698 pointData.insertReplaceScalarData(std::move(vels),
velocityLabel);
713 levelSets.push_back(passedlsDomain);
718 levelSets.push_back(passedlsDomain);
719 velocities = passedVelocities;
724 : levelSets(passedlsDomains) {
725 velocities = passedVelocities;
731 levelSets.push_back(passedlsDomain);
737 velocities = passedVelocities;
796 integrationScheme = scheme;
816 if (levelSets.empty()) {
817 Logger::getInstance()
818 .addWarning(
"No level sets passed to Advect.")
825 }
else if (integrationScheme ==
828 }
else if (integrationScheme ==
831 }
else if (integrationScheme ==
834 }
else if (integrationScheme ==
839 }
else if (integrationScheme ==
842 }
else if (integrationScheme ==
845 }
else if (integrationScheme ==
848 }
else if (integrationScheme ==
851 }
else if (integrationScheme ==
856 Logger::getInstance()
857 .addWarning(
"Advect: Integration scheme not found.")
864 if (advectionTime == 0.) {
865 advectedTime = advect(std::numeric_limits<double>::max());
866 numberOfTimeSteps = 1;
868 double currentTime = 0.0;
869 numberOfTimeSteps = 0;
870 while (currentTime < advectionTime) {
871 currentTime += advect(advectionTime - currentTime);
873 if (performOnlySingleStep)
876 advectedTime = currentTime;
887 calculateNormalVectors);
888 currentTimeStep = integrateTime(is, maxTimeStep);
889 }
else if (integrationScheme ==
892 calculateNormalVectors);
893 currentTimeStep = integrateTime(is, maxTimeStep);
894 }
else if (integrationScheme ==
896 auto alphas = findGlobalAlphas();
898 dissipationAlpha, alphas,
899 calculateNormalVectors);
900 currentTimeStep = integrateTime(is, maxTimeStep);
901 }
else if (integrationScheme ==
903 auto alphas = findGlobalAlphas();
905 dissipationAlpha, alphas,
906 calculateNormalVectors);
907 currentTimeStep = integrateTime(is, maxTimeStep);
908 }
else if (integrationScheme ==
912 levelSets.back(), velocities);
913 currentTimeStep = integrateTime(is, maxTimeStep);
914 }
else if (integrationScheme ==
917 levelSets.back(), velocities, dissipationAlpha);
918 currentTimeStep = integrateTime(is, maxTimeStep);
919 }
else if (integrationScheme ==
922 levelSets.back(), velocities, dissipationAlpha);
923 currentTimeStep = integrateTime(is, maxTimeStep);
924 }
else if (integrationScheme ==
927 levelSets.back(), velocities, dissipationAlpha);
928 currentTimeStep = integrateTime(is, maxTimeStep);
929 }
else if (integrationScheme ==
932 levelSets.back(), velocities, dissipationAlpha);
933 currentTimeStep = integrateTime(is, maxTimeStep);
934 }
else if (integrationScheme ==
937 levelSets.back(), velocities, dissipationAlpha);
938 currentTimeStep = integrateTime(is, maxTimeStep);
940 Logger::getInstance()
941 .addWarning(
"Advect: Integration scheme not found. Not integrating.")
943 currentTimeStep = -1.;