lsAdvect.hpp

Go to the documentation of this file.

#pragma once
 
#include <lsPreCompileMacros.hpp>
 
#include <functional>
#include <limits>
#include <vector>
 
#include <hrleSparseIterator.hpp>
#include <hrleSparseStarIterator.hpp>
 
#include <vcLogger.hpp>
#include <vcSmartPointer.hpp>
 
#include <lsBooleanOperation.hpp>
#include <lsDomain.hpp>
#include <lsMarkVoidPoints.hpp>
#include <lsReduce.hpp>
 
// Spatial discretization schemes
#include <lsEngquistOsher.hpp>
#include <lsLaxFriedrichs.hpp>
#include <lsLocalLaxFriedrichs.hpp>
#include <lsLocalLaxFriedrichsAnalytical.hpp>
#include <lsLocalLocalLaxFriedrichs.hpp>
#include <lsStencilLocalLaxFriedrichsScalar.hpp>
#include <lsWENO.hpp>
 
// Include implementation of time integration schemes
#include <lsAdvectIntegrationSchemes.hpp>
 
// Velocity accessor
#include <lsVelocityField.hpp>
 
// #define DEBUG_LS_ADVECT_HPP
#ifdef DEBUG_LS_ADVECT_HPP
#include <lsToMesh.hpp>
#include <lsVTKWriter.hpp>
#endif
 
namespace viennals {
 
using namespace viennacore;

template <class T, int D> class Advect {
  using ConstSparseIterator =
      viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType>;
  using hrleIndexType = viennahrle::IndexType;
 
  // Allow the time integration struct to access private members
  friend struct lsInternal::AdvectTimeIntegration<T, D>;
 
  std::vector<SmartPointer<Domain<T, D>>> levelSets;
  SmartPointer<VelocityField<T>> velocities = nullptr;
  SpatialSchemeEnum spatialScheme = SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER;
  TemporalSchemeEnum temporalScheme = TemporalSchemeEnum::FORWARD_EULER;
  double timeStepRatio = 0.4999;
  double dissipationAlpha = 1.0;
  bool calculateNormalVectors = true;
  bool ignoreVoids = false;
  double advectionTime = 0.;
  bool performOnlySingleStep = false;
  double advectedTime = 0.;
  unsigned numberOfTimeSteps = 0;
  bool saveAdvectionVelocities = false;
  bool updatePointData = true;
  bool checkDissipation = true;
  double integrationCutoff = 0.5;
  bool adaptiveTimeStepping = false;
  unsigned adaptiveTimeStepSubdivisions = 20;
  static constexpr double wrappingLayerEpsilon = 1e-4;
  std::vector<SmartPointer<Domain<T, D>>> initialLevelSets;
  std::function<bool(SmartPointer<Domain<T, D>>)> velocityUpdateCallback =
      nullptr;
 
  // this vector will hold the maximum time step for each point and the
  // corresponding velocity
  std::vector<std::vector<std::pair<std::pair<T, T>, T>>> storedRates;
  double currentTimeStep = -1.;
 
  VectorType<T, D> findGlobalAlphas() const {
 
    auto &topDomain = levelSets.back()->getDomain();
    auto &grid = levelSets.back()->getGrid();
 
    const T gridDelta = grid.getGridDelta();
    const T deltaPos = gridDelta;
    const T deltaNeg = -gridDelta;
 
    VectorType<T, D> finalAlphas{};
 
#pragma omp parallel for
    for (unsigned p = 0; p < levelSets.back()->getNumberOfSegments(); ++p) {
      VectorType<T, D> localAlphas{};
 
      viennahrle::Index<D> startVector =
          (p == 0) ? grid.getMinGridPoint()
                   : topDomain.getSegmentation()[p - 1];
      viennahrle::Index<D> endVector =
          (p != static_cast<int>(topDomain.getNumberOfSegments() - 1))
              ? topDomain.getSegmentation()[p]
              : grid.incrementIndices(grid.getMaxGridPoint());
 
      // an iterator for each level set
      std::vector<ConstSparseIterator> iterators;
      for (auto const &ls : levelSets) {
        iterators.emplace_back(ls->getDomain());
      }
 
      // neighborIterator for the top level set
      viennahrle::ConstSparseStarIterator<typename Domain<T, D>::DomainType, 1>
          neighborIterator(topDomain);
 
      for (ConstSparseIterator it(topDomain, startVector);
           it.getStartIndices() < endVector; ++it) {
 
        if (!it.isDefined() || std::abs(it.getValue()) > integrationCutoff)
          continue;
 
        const T value = it.getValue();
        const auto indices = it.getStartIndices();
 
        // check if there is any other levelset at the same point:
        // if yes, take the velocity of the lowest levelset
        for (unsigned lowerLevelSetId = 0; lowerLevelSetId < levelSets.size();
             ++lowerLevelSetId) {
          // put iterator to same position as the top levelset
          iterators[lowerLevelSetId].goToIndicesSequential(indices);
 
          // if the lower surface is actually outside, i.e. its LS value
          // is lower or equal
          if (iterators[lowerLevelSetId].getValue() <=
              value + wrappingLayerEpsilon) {
 
            // move neighborIterator to current position
            neighborIterator.goToIndicesSequential(indices);
 
            Vec3D<T> coords{};
            for (unsigned i = 0; i < D; ++i) {
              coords[i] = indices[i] * gridDelta;
            }
 
            Vec3D<T> normal{};
            T normalModulus = 0.;
            for (unsigned i = 0; i < D; ++i) {
              const T phiPos = neighborIterator.getNeighbor(i).getValue();
              const T phiNeg = neighborIterator.getNeighbor(i + D).getValue();
 
              normal[i] = phiPos - phiNeg;
              normalModulus += normal[i] * normal[i];
            }
            normalModulus = 1. / std::sqrt(normalModulus);
            for (unsigned i = 0; i < D; ++i)
              normal[i] *= normalModulus;
 
            T scaVel = velocities->getScalarVelocity(
                coords, lowerLevelSetId, normal,
                neighborIterator.getCenter().getPointId());
            auto vecVel = velocities->getVectorVelocity(
                coords, lowerLevelSetId, normal,
                neighborIterator.getCenter().getPointId());
 
            for (unsigned i = 0; i < D; ++i) {
              T tempAlpha = std::abs((scaVel + vecVel[i]) * normal[i]);
              localAlphas[i] = std::max(localAlphas[i], tempAlpha);
            }
 
            // exit material loop
            break;
          }
        }
      }
 
#pragma omp critical
      {
        for (unsigned i = 0; i < D; ++i) {
          finalAlphas[i] = std::max(finalAlphas[i], localAlphas[i]);
        }
      }
    } // end of parallel section
 
    return finalAlphas;
  }
 
  // Helper function for linear combination:
  // target = wTarget * target + wSource * source
  bool combineLevelSets(T wTarget, T wSource) {
    // Calculate required expansion width based on CFL and RK steps
    int steps = 1;
    if (temporalScheme == TemporalSchemeEnum::RUNGE_KUTTA_2ND_ORDER) {
      steps = 2;
    } else if (temporalScheme == TemporalSchemeEnum::RUNGE_KUTTA_3RD_ORDER) {
      steps = 3;
    }
 
    // Expand both level sets to ensure sufficient overlap
    int expansionWidth = std::ceil(2.0 * steps * timeStepRatio + 1);
    viennals::Expand<T, D>(levelSets.back(), expansionWidth).apply();
    viennals::Expand<T, D>(initialLevelSets.back(), expansionWidth).apply();
 
    bool movedDown = false;
 
    viennals::BooleanOperation<T, D> op(levelSets.back(),
                                        initialLevelSets.back(),
                                        viennals::BooleanOperationEnum::CUSTOM);
    op.setBooleanOperationComparator(
        [wTarget, wSource, &movedDown](const T &a, const T &b) {
          if (a != Domain<T, D>::POS_VALUE && a != Domain<T, D>::NEG_VALUE &&
              b != Domain<T, D>::POS_VALUE && b != Domain<T, D>::NEG_VALUE) {
            T res = wSource * a + wTarget * b;
            if (res > b)
              movedDown = true;
            return std::make_pair(res, true);
          }
          if (a == Domain<T, D>::POS_VALUE || a == Domain<T, D>::NEG_VALUE)
            return std::make_pair(a, false);
          return std::make_pair(b, false);
        });
    op.apply();
 
    rebuildLS();
 
    return movedDown;
  }
 
  void rebuildLS() {
    // TODO: this function uses Manhattan distances for renormalisation,
    // since this is the quickest. For visualisation applications, better
    // renormalisation is needed, so it might be good to implement
    // Euler distance renormalisation as an option
    auto &grid = levelSets.back()->getGrid();
    auto newlsDomain = SmartPointer<Domain<T, D>>::New(grid);
    auto &newDomain = newlsDomain->getDomain();
    auto &domain = levelSets.back()->getDomain();
 
    // Determine cutoff and width based on discretization scheme to avoid
    // immediate re-expansion
    T cutoff = 1.0;
    int finalWidth = 2;
    if (spatialScheme ==
        SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      cutoff = 1.5;
      finalWidth = 3;
    }
 
    newDomain.initialize(domain.getNewSegmentation(),
                         domain.getAllocation() *
                             (2.0 / levelSets.back()->getLevelSetWidth()));
 
    const bool updateData = updatePointData;
    // save how data should be transferred to new level set
    // list of indices into the old pointData vector
    std::vector<std::vector<unsigned>> newDataSourceIds;
    if (updateData)
      newDataSourceIds.resize(newDomain.getNumberOfSegments());
 
#ifdef DEBUG_LS_ADVECT_HPP
    {
      auto mesh = SmartPointer<Mesh<T>>::New();
      ToMesh<T, D>(levelSets.back(), mesh).apply();
      VTKWriter<T>(mesh, "Advect_beforeRebuild.vtk").apply();
    }
#endif
 
#pragma omp parallel for
    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
      auto &domainSegment = newDomain.getDomainSegment(p);
 
      viennahrle::Index<D> startVector =
          (p == 0) ? grid.getMinGridPoint()
                   : newDomain.getSegmentation()[p - 1];
 
      viennahrle::Index<D> endVector =
          (p != newDomain.getNumberOfSegments() - 1)
              ? newDomain.getSegmentation()[p]
              : grid.incrementIndices(grid.getMaxGridPoint());
 
      // reserve a bit more to avoid reallocation
      // would expect number of points to roughly double
      if (updateData)
        newDataSourceIds[p].reserve(2.5 * domainSegment.getNumberOfPoints());
 
      for (viennahrle::ConstSparseStarIterator<
               typename Domain<T, D>::DomainType, 1>
               it(domain, startVector);
           it.getIndices() < endVector; ++it) {
 
        // if the center is an active grid point
        // <1.0 since it could have been change by 0.5 max
        if (std::abs(it.getCenter().getValue()) <= 1.0) {
 
          int k = 0;
          for (; k < 2 * D; k++)
            if (std::signbit(it.getNeighbor(k).getValue() - 1e-7) !=
                std::signbit(it.getCenter().getValue() + 1e-7))
              break;
 
          // if there is at least one neighbor of opposite sign
          if (k != 2 * D) {
            if (it.getCenter().getDefinedValue() > 0.5) {
              int j = 0;
              for (; j < 2 * D; j++) {
                if (std::abs(it.getNeighbor(j).getValue()) <= 1.0)
                  if (it.getNeighbor(j).getDefinedValue() < -0.5)
                    break;
              }
              if (j == 2 * D) {
                domainSegment.insertNextDefinedPoint(
                    it.getIndices(), it.getCenter().getDefinedValue());
                if (updateData)
                  newDataSourceIds[p].push_back(it.getCenter().getPointId());
                // if there is at least one active grid point, which is < -0.5
              } else {
                domainSegment.insertNextDefinedPoint(it.getIndices(), 0.5);
                if (updateData)
                  newDataSourceIds[p].push_back(it.getNeighbor(j).getPointId());
              }
            } else if (it.getCenter().getDefinedValue() < -0.5) {
              int j = 0;
              for (; j < 2 * D; j++) {
                if (std::abs(it.getNeighbor(j).getValue()) <= 1.0)
                  if (it.getNeighbor(j).getDefinedValue() > 0.5)
                    break;
              }
 
              if (j == 2 * D) {
                domainSegment.insertNextDefinedPoint(
                    it.getIndices(), it.getCenter().getDefinedValue());
                if (updateData)
                  newDataSourceIds[p].push_back(it.getCenter().getPointId());
                // if there is at least one active grid point, which is > 0.5
              } else {
                domainSegment.insertNextDefinedPoint(it.getIndices(), -0.5);
                if (updateData)
                  newDataSourceIds[p].push_back(it.getNeighbor(j).getPointId());
              }
            } else {
              domainSegment.insertNextDefinedPoint(
                  it.getIndices(), it.getCenter().getDefinedValue());
              if (updateData)
                newDataSourceIds[p].push_back(it.getCenter().getPointId());
            }
          } else {
            domainSegment.insertNextUndefinedPoint(
                it.getIndices(), (it.getCenter().getDefinedValue() < 0)
                                     ? Domain<T, D>::NEG_VALUE
                                     : Domain<T, D>::POS_VALUE);
          }
 
        } else { // if the center is not an active grid point
          if (it.getCenter().getValue() >= 0) {
            int usedNeighbor = -1;
            T distance = Domain<T, D>::POS_VALUE;
            for (int i = 0; i < 2 * D; i++) {
              T value = it.getNeighbor(i).getValue();
              if (std::abs(value) <= 1.0 && (value < 0.)) {
                if (distance > value + 1.0) {
                  distance = value + 1.0;
                  usedNeighbor = i;
                }
              }
            }
 
            if (distance <= cutoff) {
              domainSegment.insertNextDefinedPoint(it.getIndices(), distance);
              if (updateData)
                newDataSourceIds[p].push_back(
                    it.getNeighbor(usedNeighbor).getPointId());
            } else {
              domainSegment.insertNextUndefinedPoint(it.getIndices(),
                                                     Domain<T, D>::POS_VALUE);
            }
 
          } else {
            int usedNeighbor = -1;
            T distance = Domain<T, D>::NEG_VALUE;
            for (int i = 0; i < 2 * D; i++) {
              T value = it.getNeighbor(i).getValue();
              if (std::abs(value) <= 1.0 && (value > 0)) {
                if (distance < value - 1.0) {
                  // distance = std::max(distance, value - T(1.0));
                  distance = value - 1.0;
                  usedNeighbor = i;
                }
              }
            }
 
            if (distance >= -cutoff) {
              domainSegment.insertNextDefinedPoint(it.getIndices(), distance);
              if (updateData)
                newDataSourceIds[p].push_back(
                    it.getNeighbor(usedNeighbor).getPointId());
            } else {
              domainSegment.insertNextUndefinedPoint(it.getIndices(),
                                                     Domain<T, D>::NEG_VALUE);
            }
          }
        }
      }
    }
 
    // now copy old data into new level set
    if (updateData) {
      auto &pointData = levelSets.back()->getPointData();
      newlsDomain->getPointData().translateFromMultiData(pointData,
                                                         newDataSourceIds);
    }
 
    newDomain.finalize();
    newDomain.segment();
    levelSets.back()->deepCopy(newlsDomain);
    levelSets.back()->finalize(finalWidth);
  }

  template <class DiscretizationSchemeType>
  double integrateTime(DiscretizationSchemeType spatialScheme,
                       double maxTimeStep) {
 
    auto &topDomain = levelSets.back()->getDomain();
    auto &grid = levelSets.back()->getGrid();
 
    typename PointData<T>::ScalarDataType *voidMarkerPointer;
    if (ignoreVoids) {
      MarkVoidPoints<T, D>(levelSets.back()).apply();
      auto &pointData = levelSets.back()->getPointData();
      voidMarkerPointer =
          pointData.getScalarData(MarkVoidPoints<T, D>::voidPointLabel, true);
      if (voidMarkerPointer == nullptr) {
        VIENNACORE_LOG_WARNING("Advect: Cannot find void point markers. Not "
                               "ignoring void points.");
        ignoreVoids = false;
      }
    }
    const bool ignoreVoidPoints = ignoreVoids;
    const bool useAdaptiveTimeStepping = adaptiveTimeStepping;
    const auto adaptiveFactor = 1.0 / adaptiveTimeStepSubdivisions;
 
    if (!storedRates.empty()) {
      VIENNACORE_LOG_WARNING("Advect: Overwriting previously stored rates.");
    }
 
    storedRates.resize(topDomain.getNumberOfSegments());
 
#pragma omp parallel for
    for (unsigned p = 0; p < topDomain.getNumberOfSegments(); ++p) {
 
      viennahrle::Index<D> startVector =
          (p == 0) ? grid.getMinGridPoint()
                   : topDomain.getSegmentation()[p - 1];
 
      viennahrle::Index<D> endVector =
          (p != topDomain.getNumberOfSegments() - 1)
              ? topDomain.getSegmentation()[p]
              : grid.incrementIndices(grid.getMaxGridPoint());
 
      double tempMaxTimeStep = maxTimeStep;
      // store the rates and value of underneath LS for this segment
      auto &rates = storedRates[p];
      rates.reserve(
          topDomain.getNumberOfPoints() /
              static_cast<double>((levelSets.back())->getNumberOfSegments()) +
          10);
 
      // an iterator for each level set
      std::vector<ConstSparseIterator> iterators;
      for (auto const &ls : levelSets) {
        iterators.emplace_back(ls->getDomain());
      }
 
      DiscretizationSchemeType scheme(spatialScheme);
 
      for (ConstSparseIterator it(topDomain, startVector);
           it.getStartIndices() < endVector; it.next()) {
 
        if (!it.isDefined() || std::abs(it.getValue()) > integrationCutoff)
          continue;
 
        T value = it.getValue();
        double maxStepTime = 0;
        double cfl = timeStepRatio;
 
        for (int currentLevelSetId = levelSets.size() - 1;
             currentLevelSetId >= 0; --currentLevelSetId) {
 
          std::pair<T, T> gradNDissipation;
 
          if (!(ignoreVoidPoints && (*voidMarkerPointer)[it.getPointId()])) {
            // check if there is any other levelset at the same point:
            // if yes, take the velocity of the lowest levelset
            for (unsigned lowerLevelSetId = 0;
                 lowerLevelSetId < levelSets.size(); ++lowerLevelSetId) {
              // put iterator to same position as the top levelset
              iterators[lowerLevelSetId].goToIndicesSequential(
                  it.getStartIndices());
 
              // if the lower surface is actually outside, i.e. its LS value
              // is lower or equal
              if (iterators[lowerLevelSetId].getValue() <=
                  value + wrappingLayerEpsilon) {
                gradNDissipation =
                    scheme(it.getStartIndices(), lowerLevelSetId);
                break;
              }
            }
          }
 
          T velocity = gradNDissipation.first - gradNDissipation.second;
          if (velocity > 0.) {
            // Case 1: Growth / Deposition (Velocity > 0)
            // Limit the time step based on the standard CFL condition.
            maxStepTime += cfl / velocity;
            rates.emplace_back(gradNDissipation,
                               -std::numeric_limits<T>::max());
            break;
          } else if (velocity == 0.) {
            // Case 2: Static (Velocity == 0)
            // No time step limit imposed by this point.
            maxStepTime = std::numeric_limits<T>::max();
            rates.emplace_back(gradNDissipation, std::numeric_limits<T>::max());
            break;
          } else {
            // Case 3: Etching (Velocity < 0)
            // Retrieve the interface location of the underlying material.
            T valueBelow;
            if (currentLevelSetId > 0) {
              iterators[currentLevelSetId - 1].goToIndicesSequential(
                  it.getStartIndices());
              valueBelow = iterators[currentLevelSetId - 1].getValue();
            } else {
              valueBelow = std::numeric_limits<T>::max();
            }
            // Calculate the top material thickness
            T difference = std::abs(valueBelow - value);
 
            if (difference >= cfl) {
              // Sub-case 3a: Standard Advection
              // Far from interface: Use full CFL time step.
              maxStepTime -= cfl / velocity;
              rates.emplace_back(gradNDissipation,
                                 std::numeric_limits<T>::max());
              break;
            } else {
              // Sub-case 3b: Interface Interaction
              // Use adaptiveFactor as threshold.
              if (useAdaptiveTimeStepping &&
                  difference > adaptiveFactor * cfl) {
                // Adaptive Sub-stepping:
                // Approaching boundary: Force small steps to gather
                // flux statistics and prevent numerical overshoot ("Soft
                // Landing").
                maxStepTime -= adaptiveFactor * cfl / velocity;
                rates.emplace_back(gradNDissipation,
                                   std::numeric_limits<T>::max());
                break;
              } else {
                // Terminal Step:
                // Within tolerance: Snap to boundary, consume budget, and
                // switch material.
                cfl -= difference;
                value = valueBelow;
                maxStepTime -= difference / velocity;
                rates.emplace_back(gradNDissipation, valueBelow);
              }
            }
          }
        }
 
        if (maxStepTime < tempMaxTimeStep)
          tempMaxTimeStep = maxStepTime;
      }
 
#pragma omp critical
      {
        // If a Lax Friedrichs scheme is selected the time step is
        // reduced depending on the dissipation coefficients
        // For Engquist Osher scheme this function is empty.
        scheme.reduceTimeStepHamiltonJacobi(
            tempMaxTimeStep, levelSets.back()->getGrid().getGridDelta());
 
        // set global timestep maximum
        if (tempMaxTimeStep < maxTimeStep)
          maxTimeStep = tempMaxTimeStep;
      }
    } // end of parallel section
 
    // maxTimeStep is now the maximum time step possible for all points
    // and rates are stored in a vector
    return maxTimeStep;
  }

  void computeRates(double maxTimeStep = std::numeric_limits<double>::max()) {
    prepareLS();
    if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
      auto is = lsInternal::EngquistOsher<T, D, 1>(levelSets.back(), velocities,
                                                   calculateNormalVectors);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_2ND_ORDER) {
      auto is = lsInternal::EngquistOsher<T, D, 2>(levelSets.back(), velocities,
                                                   calculateNormalVectors);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER) {
      auto alphas = findGlobalAlphas();
      auto is = lsInternal::LaxFriedrichs<T, D, 1>(levelSets.back(), velocities,
                                                   dissipationAlpha, alphas,
                                                   calculateNormalVectors);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER) {
      auto alphas = findGlobalAlphas();
      auto is = lsInternal::LaxFriedrichs<T, D, 2>(levelSets.back(), velocities,
                                                   dissipationAlpha, alphas,
                                                   calculateNormalVectors);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER) {
      auto is = lsInternal::LocalLaxFriedrichsAnalytical<T, D, 1>(
          levelSets.back(), velocities);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      auto is = lsInternal::LocalLocalLaxFriedrichs<T, D, 1>(
          levelSets.back(), velocities, dissipationAlpha);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
      auto is = lsInternal::LocalLocalLaxFriedrichs<T, D, 2>(
          levelSets.back(), velocities, dissipationAlpha);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      auto is = lsInternal::LocalLaxFriedrichs<T, D, 1>(
          levelSets.back(), velocities, dissipationAlpha);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
      auto is = lsInternal::LocalLaxFriedrichs<T, D, 2>(
          levelSets.back(), velocities, dissipationAlpha);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme ==
               SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      auto is = lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>(
          levelSets.back(), velocities, dissipationAlpha);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme == SpatialSchemeEnum::WENO_3RD_ORDER) {
      // Instantiate WENO with order 3
      auto is = lsInternal::WENO<T, D, 3>(levelSets.back(), velocities,
                                          dissipationAlpha);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else if (spatialScheme == SpatialSchemeEnum::WENO_5TH_ORDER) {
      // Instantiate WENO with order 5
      auto is = lsInternal::WENO<T, D, 5>(levelSets.back(), velocities,
                                          dissipationAlpha);
      currentTimeStep = integrateTime(is, maxTimeStep);
    } else {
      VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
      currentTimeStep = -1.;
    }
  }
 
  // Level Sets below are also considered in order to adjust the advection
  // depth accordingly if there would be a material change.
  void updateLevelSet(double dt) {
    if (timeStepRatio >= 0.5) {
      VIENNACORE_LOG_WARNING(
          "Integration time step ratio should be smaller than 0.5. "
          "Advection might fail!");
    }
 
    auto &topDomain = levelSets.back()->getDomain();
 
    assert(dt >= 0. && "No time step set!");
    assert(storedRates.size() == topDomain.getNumberOfSegments());
 
    // reduce to one layer thickness and apply new values directly to the
    // domain segments --> DO NOT CHANGE SEGMENTATION HERE (true parameter)
    Reduce<T, D>(levelSets.back(), 1, true).apply();
 
    const bool saveVelocities = saveAdvectionVelocities;
    std::vector<std::vector<double>> dissipationVectors(
        levelSets.back()->getNumberOfSegments());
    std::vector<std::vector<double>> velocityVectors(
        levelSets.back()->getNumberOfSegments());
 
    const bool checkDiss = checkDissipation;
 
#pragma omp parallel for
    for (unsigned p = 0; p < topDomain.getNumberOfSegments(); ++p) {
 
      auto itRS = storedRates[p].cbegin();
      auto &segment = topDomain.getDomainSegment(p);
      const unsigned maxId = segment.getNumberOfPoints();
 
      if (saveVelocities) {
        velocityVectors[p].resize(maxId);
        dissipationVectors[p].resize(maxId);
      }
 
      for (unsigned localId = 0; localId < maxId; ++localId) {
        T &value = segment.definedValues[localId];
 
        // Skip points that were not part of computeRates (outer layers)
        if (std::abs(value) > integrationCutoff)
          continue;
 
        double time = dt;
 
        // if there is a change in materials during one time step, deduct
        // the time taken to advect up to the end of the top material and
        // set the LS value to the one below
        auto const [gradient, dissipation] = itRS->first;
        T velocity = gradient - dissipation;
        // check if dissipation is too high and would cause a change in
        // direction of the velocity
        if (checkDiss && ((gradient < 0 && velocity > 0) ||
                          (gradient > 0 && velocity < 0))) {
          velocity = 0;
        }
 
        T rate = time * velocity;
        while (std::abs(itRS->second - value) < std::abs(rate)) {
          time -= std::abs((itRS->second - value) / velocity);
          value = itRS->second;
          ++itRS; // advance the TempStopRates iterator by one
 
          // recalculate velocity and rate
          velocity = itRS->first.first - itRS->first.second;
          if (checkDiss && ((itRS->first.first < 0 && velocity > 0) ||
                            (itRS->first.first > 0 && velocity < 0))) {
            velocity = 0;
          }
          rate = time * velocity;
        }
 
        // now deduct the velocity times the time step we take
        value -= rate;
 
        if (saveVelocities) {
          velocityVectors[p][localId] = rate;
          dissipationVectors[p][localId] = itRS->first.second;
        }
 
        // this is run when two materials are close but the velocity is too slow
        // to actually reach the second material, to get rid of the extra
        // entry in the TempRatesStop
        while (std::abs(itRS->second) != std::numeric_limits<T>::max())
          ++itRS;
 
        // advance the TempStopRates iterator by one
        ++itRS;
      }
    } // end of parallel section
 
    if (saveVelocities) {
      auto &pointData = levelSets.back()->getPointData();
 
      typename PointData<T>::ScalarDataType vels;
      typename PointData<T>::ScalarDataType diss;
 
      for (unsigned i = 0; i < velocityVectors.size(); ++i) {
        vels.insert(vels.end(),
                    std::make_move_iterator(velocityVectors[i].begin()),
                    std::make_move_iterator(velocityVectors[i].end()));
        diss.insert(diss.end(),
                    std::make_move_iterator(dissipationVectors[i].begin()),
                    std::make_move_iterator(dissipationVectors[i].end()));
      }
      pointData.insertReplaceScalarData(std::move(vels), velocityLabel);
      pointData.insertReplaceScalarData(std::move(diss), dissipationLabel);
    }
 
    // clear the stored rates since surface has changed
    storedRates.clear();
  }
 
  void adjustLowerLayers() {
    // Adjust all level sets below the advected one
    if (spatialScheme !=
        viennals::SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      for (unsigned i = 0; i < levelSets.size() - 1; ++i) {
        viennals::BooleanOperation<T, D>(
            levelSets[i], levelSets.back(),
            viennals::BooleanOperationEnum::INTERSECT)
            .apply();
      }
    }
  }

  double advect(double maxTimeStep) {
    switch (temporalScheme) {
    case TemporalSchemeEnum::RUNGE_KUTTA_2ND_ORDER:
      return lsInternal::AdvectTimeIntegration<T, D>::evolveRungeKutta2(
          *this, maxTimeStep);
    case TemporalSchemeEnum::RUNGE_KUTTA_3RD_ORDER:
      return lsInternal::AdvectTimeIntegration<T, D>::evolveRungeKutta3(
          *this, maxTimeStep);
    case TemporalSchemeEnum::FORWARD_EULER:
    default:
      return lsInternal::AdvectTimeIntegration<T, D>::evolveForwardEuler(
          *this, maxTimeStep);
    }
  }
 
public:
  static constexpr char velocityLabel[] = "AdvectionVelocities";
  static constexpr char dissipationLabel[] = "Dissipation";
 
  Advect() = default;
 
  explicit Advect(SmartPointer<Domain<T, D>> passedlsDomain) {
    levelSets.push_back(passedlsDomain);
  }
 
  Advect(SmartPointer<Domain<T, D>> passedlsDomain,
         SmartPointer<VelocityField<T>> passedVelocities) {
    levelSets.push_back(passedlsDomain);
    velocities = passedVelocities;
  }
 
  Advect(std::vector<SmartPointer<Domain<T, D>>> passedlsDomains,
         SmartPointer<VelocityField<T>> passedVelocities)
      : levelSets(passedlsDomains) {
    velocities = passedVelocities;
  }

  void insertNextLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
    levelSets.push_back(passedlsDomain);
  }
 
  void clearLevelSets() { levelSets.clear(); }

  void setVelocityField(SmartPointer<VelocityField<T>> passedVelocities) {
    velocities = passedVelocities;
  }

  void setAdvectionTime(double time) { advectionTime = time; }

  void setSingleStep(bool singleStep) { performOnlySingleStep = singleStep; }

  void setTimeStepRatio(const double &cfl) { timeStepRatio = cfl; }

  void setCalculateNormalVectors(bool cnv) { calculateNormalVectors = cnv; }

  void setIgnoreVoids(bool iV) { ignoreVoids = iV; }

  void setAdaptiveTimeStepping(bool aTS = true, unsigned subdivisions = 20) {
    adaptiveTimeStepping = aTS;
    if (subdivisions < 1) {
      VIENNACORE_LOG_WARNING("Advect: Adaptive time stepping subdivisions must "
                             "be at least 1. Setting to 1.");
      subdivisions = 1;
    }
    adaptiveTimeStepSubdivisions = subdivisions;
  }

  void setSaveAdvectionVelocities(bool sAV) { saveAdvectionVelocities = sAV; }

  double getAdvectedTime() const { return advectedTime; }

  double getCurrentTimeStep() const { return currentTimeStep; }

  unsigned getNumberOfTimeSteps() const { return numberOfTimeSteps; }

  double getTimeStepRatio() const { return timeStepRatio; }

  bool getCalculateNormalVectors() const { return calculateNormalVectors; }

  void setSpatialScheme(SpatialSchemeEnum scheme) { spatialScheme = scheme; }
 
  // Deprecated and will be removed in future versions:
  // use setSpatialScheme instead
  [[deprecated("Use setSpatialScheme instead")]] void
  setIntegrationScheme(IntegrationSchemeEnum scheme) {
    VIENNACORE_LOG_WARNING(
        "Advect::setIntegrationScheme is deprecated and will be removed in "
        "future versions. Use setSpatialScheme instead.");
    spatialScheme = scheme;
  }

  void setTemporalScheme(TemporalSchemeEnum scheme) { temporalScheme = scheme; }

  void setDissipationAlpha(const double &a) { dissipationAlpha = a; }
 
  // Sets the velocity to 0 if the dissipation is too high
  void setCheckDissipation(bool check) { checkDissipation = check; }

  void setUpdatePointData(bool update) { updatePointData = update; }

  void setVelocityUpdateCallback(
      std::function<bool(SmartPointer<Domain<T, D>>)> callback) {
    velocityUpdateCallback = callback;
  }
 
  // Prepare the levelset for advection, based on the provided spatial
  // discretization scheme.
  void prepareLS() {
    // check whether a level set and velocities have been given
    if (levelSets.empty()) {
      VIENNACORE_LOG_ERROR("No level sets passed to Advect.");
      return;
    }
 
    if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
      lsInternal::EngquistOsher<T, D, 1>::prepareLS(levelSets.back());
    } else if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_2ND_ORDER) {
      lsInternal::EngquistOsher<T, D, 2>::prepareLS(levelSets.back());
    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER) {
      lsInternal::LaxFriedrichs<T, D, 1>::prepareLS(levelSets.back());
    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER) {
      lsInternal::LaxFriedrichs<T, D, 2>::prepareLS(levelSets.back());
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER) {
      lsInternal::LocalLaxFriedrichsAnalytical<T, D, 1>::prepareLS(
          levelSets.back());
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      lsInternal::LocalLocalLaxFriedrichs<T, D, 1>::prepareLS(levelSets.back());
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
      lsInternal::LocalLocalLaxFriedrichs<T, D, 2>::prepareLS(levelSets.back());
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      lsInternal::LocalLaxFriedrichs<T, D, 1>::prepareLS(levelSets.back());
    } else if (spatialScheme ==
               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
      lsInternal::LocalLaxFriedrichs<T, D, 2>::prepareLS(levelSets.back());
    } else if (spatialScheme ==
               SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
      lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>::prepareLS(
          levelSets.back());
    } else if (spatialScheme == SpatialSchemeEnum::WENO_3RD_ORDER) {
      lsInternal::WENO<T, D, 3>::prepareLS(levelSets.back());
    } else if (spatialScheme == SpatialSchemeEnum::WENO_5TH_ORDER) {
      lsInternal::WENO<T, D, 5>::prepareLS(levelSets.back());
    } else {
      VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
    }
  }
 
  void apply() {
    // check whether a level set and velocities have been given
    if (levelSets.empty()) {
      VIENNACORE_LOG_ERROR("No level sets passed to Advect. Not advecting.");
      return;
    }
    if (velocities == nullptr) {
      VIENNACORE_LOG_ERROR(
          "No velocity field passed to Advect. Not advecting.");
      return;
    }
 
    if (advectionTime == 0.) {
      advectedTime = advect(std::numeric_limits<double>::max());
      numberOfTimeSteps = 1;
    } else {
      double currentTime = 0.0;
      numberOfTimeSteps = 0;
      while (currentTime < advectionTime) {
        currentTime += advect(advectionTime - currentTime);
        ++numberOfTimeSteps;
        if (performOnlySingleStep)
          break;
      }
      advectedTime = currentTime;
    }
  }
};
 
// add all template specializations for this class
PRECOMPILE_PRECISION_DIMENSION(Advect)
 
} // namespace viennals