path: root/extern/mantaflow/preprocessed/plugin/advection.cpp

// DO NOT EDIT !
// This file is generated using the MantaFlow preprocessor (prep generate).

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011-2015 Tobias Pfaff, Nils Thuerey
 *
 * This program is free software, distributed under the terms of the
 * Apache License, Version 2.0
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Plugins for advection
 *
 ******************************************************************************/

#include "vectorbase.h"
#include "grid.h"
#include "kernel.h"
#include <limits>

using namespace std;

namespace Manta {

//! Semi-Lagrange interpolation kernel

template<class T> struct SemiLagrange : public KernelBase {
  SemiLagrange(const FlagGrid &flags,
               const MACGrid &vel,
               Grid<T> &dst,
               const Grid<T> &src,
               Real dt,
               bool isLevelset,
               int orderSpace,
               int orderTrace)
      : KernelBase(&flags, 1),
        flags(flags),
        vel(vel),
        dst(dst),
        src(src),
        dt(dt),
        isLevelset(isLevelset),
        orderSpace(orderSpace),
        orderTrace(orderTrace)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 Grid<T> &dst,
                 const Grid<T> &src,
                 Real dt,
                 bool isLevelset,
                 int orderSpace,
                 int orderTrace)
  {
    if (orderTrace == 1) {
      // traceback position
      Vec3 pos = Vec3(i + 0.5f, j + 0.5f, k + 0.5f) - vel.getCentered(i, j, k) * dt;
      dst(i, j, k) = src.getInterpolatedHi(pos, orderSpace);
    }
    else if (orderTrace == 2) {
      // backtracing using explicit midpoint
      Vec3 p0 = Vec3(i + 0.5f, j + 0.5f, k + 0.5f);
      Vec3 p1 = p0 - vel.getCentered(i, j, k) * dt * 0.5;
      Vec3 p2 = p0 - vel.getInterpolated(p1) * dt;
      dst(i, j, k) = src.getInterpolatedHi(p2, orderSpace);
    }
    else {
      assertMsg(false, "Unknown backtracing order " << orderTrace);
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return vel;
  }
  typedef MACGrid type1;
  inline Grid<T> &getArg2()
  {
    return dst;
  }
  typedef Grid<T> type2;
  inline const Grid<T> &getArg3()
  {
    return src;
  }
  typedef Grid<T> type3;
  inline Real &getArg4()
  {
    return dt;
  }
  typedef Real type4;
  inline bool &getArg5()
  {
    return isLevelset;
  }
  typedef bool type5;
  inline int &getArg6()
  {
    return orderSpace;
  }
  typedef int type6;
  inline int &getArg7()
  {
    return orderTrace;
  }
  typedef int type7;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 1; j < _maxY; j++)
            for (int i = 1; i < _maxX; i++)
              op(i, j, k, flags, vel, dst, src, dt, isLevelset, orderSpace, orderTrace);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, vel, dst, src, dt, isLevelset, orderSpace, orderTrace);
      }
    }
  }
  const FlagGrid &flags;
  const MACGrid &vel;
  Grid<T> &dst;
  const Grid<T> &src;
  Real dt;
  bool isLevelset;
  int orderSpace;
  int orderTrace;
};

//! Semi-Lagrange interpolation kernel for MAC grids

struct SemiLagrangeMAC : public KernelBase {
  SemiLagrangeMAC(const FlagGrid &flags,
                  const MACGrid &vel,
                  MACGrid &dst,
                  const MACGrid &src,
                  Real dt,
                  int orderSpace,
                  int orderTrace)
      : KernelBase(&flags, 1),
        flags(flags),
        vel(vel),
        dst(dst),
        src(src),
        dt(dt),
        orderSpace(orderSpace),
        orderTrace(orderTrace)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 MACGrid &dst,
                 const MACGrid &src,
                 Real dt,
                 int orderSpace,
                 int orderTrace)
  {
    if (orderTrace == 1) {
      // get currect velocity at MAC position
      // no need to shift xpos etc. as lookup field is also shifted
      Vec3 xpos = Vec3(i + 0.5f, j + 0.5f, k + 0.5f) - vel.getAtMACX(i, j, k) * dt;
      Real vx = src.getInterpolatedComponentHi<0>(xpos, orderSpace);
      Vec3 ypos = Vec3(i + 0.5f, j + 0.5f, k + 0.5f) - vel.getAtMACY(i, j, k) * dt;
      Real vy = src.getInterpolatedComponentHi<1>(ypos, orderSpace);
      Vec3 zpos = Vec3(i + 0.5f, j + 0.5f, k + 0.5f) - vel.getAtMACZ(i, j, k) * dt;
      Real vz = src.getInterpolatedComponentHi<2>(zpos, orderSpace);

      dst(i, j, k) = Vec3(vx, vy, vz);
    }
    else if (orderTrace == 2) {
      Vec3 p0 = Vec3(i + 0.5, j + 0.5, k + 0.5);
      Vec3 xp0 = Vec3(i, j + 0.5f, k + 0.5f);
      Vec3 xp1 = xp0 - src.getAtMACX(i, j, k) * dt * 0.5;
      Vec3 xp2 = p0 - src.getInterpolated(xp1) * dt;
      Real vx = src.getInterpolatedComponentHi<0>(xp2, orderSpace);
      Vec3 yp0 = Vec3(i + 0.5f, j, k + 0.5f);
      Vec3 yp1 = yp0 - src.getAtMACY(i, j, k) * dt * 0.5;
      Vec3 yp2 = p0 - src.getInterpolated(yp1) * dt;
      Real vy = src.getInterpolatedComponentHi<1>(yp2, orderSpace);
      Vec3 zp0 = Vec3(i + 0.5f, j + 0.5f, k);
      Vec3 zp1 = zp0 - src.getAtMACZ(i, j, k) * dt * 0.5;
      Vec3 zp2 = p0 - src.getInterpolated(zp1) * dt;
      Real vz = src.getInterpolatedComponentHi<2>(zp2, orderSpace);

      dst(i, j, k) = Vec3(vx, vy, vz);
    }
    else {
      assertMsg(false, "Unknown backtracing order " << orderTrace);
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return vel;
  }
  typedef MACGrid type1;
  inline MACGrid &getArg2()
  {
    return dst;
  }
  typedef MACGrid type2;
  inline const MACGrid &getArg3()
  {
    return src;
  }
  typedef MACGrid type3;
  inline Real &getArg4()
  {
    return dt;
  }
  typedef Real type4;
  inline int &getArg5()
  {
    return orderSpace;
  }
  typedef int type5;
  inline int &getArg6()
  {
    return orderTrace;
  }
  typedef int type6;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 1; j < _maxY; j++)
            for (int i = 1; i < _maxX; i++)
              op(i, j, k, flags, vel, dst, src, dt, orderSpace, orderTrace);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, vel, dst, src, dt, orderSpace, orderTrace);
      }
    }
  }
  const FlagGrid &flags;
  const MACGrid &vel;
  MACGrid &dst;
  const MACGrid &src;
  Real dt;
  int orderSpace;
  int orderTrace;
};

//! Kernel: Correct based on forward and backward SL steps (for both centered & mac grids)

template<class T> struct MacCormackCorrect : public KernelBase {
  MacCormackCorrect(const FlagGrid &flags,
                    Grid<T> &dst,
                    const Grid<T> &old,
                    const Grid<T> &fwd,
                    const Grid<T> &bwd,
                    Real strength,
                    bool isLevelSet,
                    bool isMAC = false)
      : KernelBase(&flags, 0),
        flags(flags),
        dst(dst),
        old(old),
        fwd(fwd),
        bwd(bwd),
        strength(strength),
        isLevelSet(isLevelSet),
        isMAC(isMAC)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const FlagGrid &flags,
                 Grid<T> &dst,
                 const Grid<T> &old,
                 const Grid<T> &fwd,
                 const Grid<T> &bwd,
                 Real strength,
                 bool isLevelSet,
                 bool isMAC = false)
  {
    dst[idx] = fwd[idx];

    if (flags.isFluid(idx)) {
      // only correct inside fluid region; note, strenth of correction can be modified here
      dst[idx] += strength * 0.5 * (old[idx] - bwd[idx]);
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<T> &getArg1()
  {
    return dst;
  }
  typedef Grid<T> type1;
  inline const Grid<T> &getArg2()
  {
    return old;
  }
  typedef Grid<T> type2;
  inline const Grid<T> &getArg3()
  {
    return fwd;
  }
  typedef Grid<T> type3;
  inline const Grid<T> &getArg4()
  {
    return bwd;
  }
  typedef Grid<T> type4;
  inline Real &getArg5()
  {
    return strength;
  }
  typedef Real type5;
  inline bool &getArg6()
  {
    return isLevelSet;
  }
  typedef bool type6;
  inline bool &getArg7()
  {
    return isMAC;
  }
  typedef bool type7;
  void runMessage(){};
  void run()
  {
    const IndexInt _sz = size;
#pragma omp parallel
    {

#pragma omp for
      for (IndexInt i = 0; i < _sz; i++)
        op(i, flags, dst, old, fwd, bwd, strength, isLevelSet, isMAC);
    }
  }
  const FlagGrid &flags;
  Grid<T> &dst;
  const Grid<T> &old;
  const Grid<T> &fwd;
  const Grid<T> &bwd;
  Real strength;
  bool isLevelSet;
  bool isMAC;
};

//! Kernel: Correct based on forward and backward SL steps (for both centered & mac grids)

template<class T> struct MacCormackCorrectMAC : public KernelBase {
  MacCormackCorrectMAC(const FlagGrid &flags,
                       Grid<T> &dst,
                       const Grid<T> &old,
                       const Grid<T> &fwd,
                       const Grid<T> &bwd,
                       Real strength,
                       bool isLevelSet,
                       bool isMAC = false)
      : KernelBase(&flags, 0),
        flags(flags),
        dst(dst),
        old(old),
        fwd(fwd),
        bwd(bwd),
        strength(strength),
        isLevelSet(isLevelSet),
        isMAC(isMAC)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 Grid<T> &dst,
                 const Grid<T> &old,
                 const Grid<T> &fwd,
                 const Grid<T> &bwd,
                 Real strength,
                 bool isLevelSet,
                 bool isMAC = false)
  {
    bool skip[3] = {false, false, false};

    if (!flags.isFluid(i, j, k))
      skip[0] = skip[1] = skip[2] = true;
    if (isMAC) {
      if ((i > 0) && (!flags.isFluid(i - 1, j, k)))
        skip[0] = true;
      if ((j > 0) && (!flags.isFluid(i, j - 1, k)))
        skip[1] = true;
      if ((k > 0) && (!flags.isFluid(i, j, k - 1)))
        skip[2] = true;
    }

    for (int c = 0; c < 3; ++c) {
      if (skip[c]) {
        dst(i, j, k)[c] = fwd(i, j, k)[c];
      }
      else {
        // perform actual correction with given strength
        dst(i, j, k)[c] = fwd(i, j, k)[c] + strength * 0.5 * (old(i, j, k)[c] - bwd(i, j, k)[c]);
      }
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<T> &getArg1()
  {
    return dst;
  }
  typedef Grid<T> type1;
  inline const Grid<T> &getArg2()
  {
    return old;
  }
  typedef Grid<T> type2;
  inline const Grid<T> &getArg3()
  {
    return fwd;
  }
  typedef Grid<T> type3;
  inline const Grid<T> &getArg4()
  {
    return bwd;
  }
  typedef Grid<T> type4;
  inline Real &getArg5()
  {
    return strength;
  }
  typedef Real type5;
  inline bool &getArg6()
  {
    return isLevelSet;
  }
  typedef bool type6;
  inline bool &getArg7()
  {
    return isMAC;
  }
  typedef bool type7;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 0; j < _maxY; j++)
            for (int i = 0; i < _maxX; i++)
              op(i, j, k, flags, dst, old, fwd, bwd, strength, isLevelSet, isMAC);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 0; j < _maxY; j++)
          for (int i = 0; i < _maxX; i++)
            op(i, j, k, flags, dst, old, fwd, bwd, strength, isLevelSet, isMAC);
      }
    }
  }
  const FlagGrid &flags;
  Grid<T> &dst;
  const Grid<T> &old;
  const Grid<T> &fwd;
  const Grid<T> &bwd;
  Real strength;
  bool isLevelSet;
  bool isMAC;
};

// Helper to collect min/max in a template
template<class T> inline void getMinMax(T &minv, T &maxv, const T &val)
{
  if (val < minv)
    minv = val;
  if (val > maxv)
    maxv = val;
}
template<> inline void getMinMax<Vec3>(Vec3 &minv, Vec3 &maxv, const Vec3 &val)
{
  getMinMax(minv.x, maxv.x, val.x);
  getMinMax(minv.y, maxv.y, val.y);
  getMinMax(minv.z, maxv.z, val.z);
}

//! detect out of bounds value
template<class T> inline bool cmpMinMax(T &minv, T &maxv, const T &val)
{
  if (val < minv)
    return true;
  if (val > maxv)
    return true;
  return false;
}
template<> inline bool cmpMinMax<Vec3>(Vec3 &minv, Vec3 &maxv, const Vec3 &val)
{
  return (cmpMinMax(minv.x, maxv.x, val.x) | cmpMinMax(minv.y, maxv.y, val.y) |
          cmpMinMax(minv.z, maxv.z, val.z));
}

#define checkFlag(x, y, z) (flags((x), (y), (z)) & (FlagGrid::TypeFluid | FlagGrid::TypeEmpty))

//! Helper function for clamping non-mac grids (those have specialized per component version below)
//  Note - 2 clamp modes, a sharper one (default, clampMode 1, also uses backward step),
//         and a softer version (clampMode 2) that is recommended in Andy's paper
template<class T>
inline T doClampComponent(const Vec3i &gridSize,
                          const FlagGrid &flags,
                          T dst,
                          const Grid<T> &orig,
                          const T fwd,
                          const Vec3 &pos,
                          const Vec3 &vel,
                          const int clampMode)
{
  T minv(std::numeric_limits<Real>::max()), maxv(-std::numeric_limits<Real>::max());
  bool haveFl = false;

  // forward (and optionally) backward
  Vec3i positions[2];
  int numPos = 1;
  positions[0] = toVec3i(pos - vel);
  if (clampMode == 1) {
    numPos = 2;
    positions[1] = toVec3i(pos + vel);
  }

  for (int l = 0; l < numPos; ++l) {
    Vec3i &currPos = positions[l];

    // clamp lookup to grid
    const int i0 = clamp(currPos.x, 0, gridSize.x - 1);  // note! gridsize already has -1 from call
    const int j0 = clamp(currPos.y, 0, gridSize.y - 1);
    const int k0 = clamp(currPos.z, 0, (orig.is3D() ? (gridSize.z - 1) : 1));
    const int i1 = i0 + 1, j1 = j0 + 1, k1 = (orig.is3D() ? (k0 + 1) : k0);

    // find min/max around source pos
    if (checkFlag(i0, j0, k0)) {
      getMinMax(minv, maxv, orig(i0, j0, k0));
      haveFl = true;
    }
    if (checkFlag(i1, j0, k0)) {
      getMinMax(minv, maxv, orig(i1, j0, k0));
      haveFl = true;
    }
    if (checkFlag(i0, j1, k0)) {
      getMinMax(minv, maxv, orig(i0, j1, k0));
      haveFl = true;
    }
    if (checkFlag(i1, j1, k0)) {
      getMinMax(minv, maxv, orig(i1, j1, k0));
      haveFl = true;
    }

    if (orig.is3D()) {
      if (checkFlag(i0, j0, k1)) {
        getMinMax(minv, maxv, orig(i0, j0, k1));
        haveFl = true;
      }
      if (checkFlag(i1, j0, k1)) {
        getMinMax(minv, maxv, orig(i1, j0, k1));
        haveFl = true;
      }
      if (checkFlag(i0, j1, k1)) {
        getMinMax(minv, maxv, orig(i0, j1, k1));
        haveFl = true;
      }
      if (checkFlag(i1, j1, k1)) {
        getMinMax(minv, maxv, orig(i1, j1, k1));
        haveFl = true;
      }
    }
  }

  if (!haveFl)
    return fwd;
  if (clampMode == 1) {
    dst = clamp(dst, minv, maxv);  // hard clamp
  }
  else {
    if (cmpMinMax(minv, maxv, dst))
      dst = fwd;  // recommended in paper, "softer"
  }
  return dst;
}

//! Helper function for clamping MAC grids, slight differences in flag checks
//  similar to scalar version, just uses single component c of vec3 values
//  for symmetry, reverts to first order near boundaries for clampMode 2
template<int c>
inline Real doClampComponentMAC(const FlagGrid &flags,
                                const Vec3i &gridSize,
                                Real dst,
                                const MACGrid &orig,
                                Real fwd,
                                const Vec3 &pos,
                                const Vec3 &vel,
                                const int clampMode)
{
  Real minv = std::numeric_limits<Real>::max(), maxv = -std::numeric_limits<Real>::max();
  // bool haveFl = false;

  // forward (and optionally) backward
  Vec3i positions[2];
  int numPos = 1;
  positions[0] = toVec3i(pos - vel);
  if (clampMode == 1) {
    numPos = 2;
    positions[1] = toVec3i(pos + vel);
  }

  Vec3i oPos = toVec3i(pos);
  Vec3i nbPos = oPos;
  nbPos[c] -= 1;
  if (clampMode == 2 &&
      (!(checkFlag(oPos.x, oPos.y, oPos.z) && checkFlag(nbPos.x, nbPos.y, nbPos.z))))
    return fwd;  // replaces haveFl check

  for (int l = 0; l < numPos; ++l) {
    Vec3i &currPos = positions[l];

    const int i0 = clamp(currPos.x, 0, gridSize.x - 1);  // note! gridsize already has -1 from call
    const int j0 = clamp(
        currPos.y, 0, gridSize.y - 1);  // but we need a clamp to -2 for the +1 offset below
    const int k0 = clamp(currPos.z, 0, (orig.is3D() ? (gridSize.z - 1) : 0));
    const int i1 = i0 + 1, j1 = j0 + 1, k1 = (orig.is3D() ? (k0 + 1) : k0);

    // find min/max around source pos
    getMinMax(minv, maxv, orig(i0, j0, k0)[c]);
    getMinMax(minv, maxv, orig(i1, j0, k0)[c]);
    getMinMax(minv, maxv, orig(i0, j1, k0)[c]);
    getMinMax(minv, maxv, orig(i1, j1, k0)[c]);

    if (orig.is3D()) {
      getMinMax(minv, maxv, orig(i0, j0, k1)[c]);
      getMinMax(minv, maxv, orig(i1, j0, k1)[c]);
      getMinMax(minv, maxv, orig(i0, j1, k1)[c]);
      getMinMax(minv, maxv, orig(i1, j1, k1)[c]);
    }
  }

  if (clampMode == 1) {
    dst = clamp(dst, minv, maxv);  // hard clamp
  }
  else {
    if (cmpMinMax(minv, maxv, dst))
      dst = fwd;  // recommended in paper, "softer"
  }
  return dst;
}

#undef checkFlag

//! Kernel: Clamp obtained value to min/max in source area, and reset values that point out of grid
//! or into boundaries
//          (note - MAC grids are handled below)

template<class T> struct MacCormackClamp : public KernelBase {
  MacCormackClamp(const FlagGrid &flags,
                  const MACGrid &vel,
                  Grid<T> &dst,
                  const Grid<T> &orig,
                  const Grid<T> &fwd,
                  Real dt,
                  const int clampMode)
      : KernelBase(&flags, 1),
        flags(flags),
        vel(vel),
        dst(dst),
        orig(orig),
        fwd(fwd),
        dt(dt),
        clampMode(clampMode)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 Grid<T> &dst,
                 const Grid<T> &orig,
                 const Grid<T> &fwd,
                 Real dt,
                 const int clampMode)
  {
    T dval = dst(i, j, k);
    Vec3i gridUpper = flags.getSize() - 1;

    dval = doClampComponent<T>(gridUpper,
                               flags,
                               dval,
                               orig,
                               fwd(i, j, k),
                               Vec3(i, j, k),
                               vel.getCentered(i, j, k) * dt,
                               clampMode);

    if (1 && clampMode == 1) {
      // lookup forward/backward , round to closest NB
      Vec3i posFwd = toVec3i(Vec3(i, j, k) + Vec3(0.5, 0.5, 0.5) - vel.getCentered(i, j, k) * dt);
      Vec3i posBwd = toVec3i(Vec3(i, j, k) + Vec3(0.5, 0.5, 0.5) + vel.getCentered(i, j, k) * dt);

      // test if lookups point out of grid or into obstacle (note doClampComponent already checks
      // sides, below is needed for valid flags access)
      if (posFwd.x < 0 || posFwd.y < 0 || posFwd.z < 0 || posBwd.x < 0 || posBwd.y < 0 ||
          posBwd.z < 0 || posFwd.x > gridUpper.x || posFwd.y > gridUpper.y ||
          ((posFwd.z > gridUpper.z) && flags.is3D()) || posBwd.x > gridUpper.x ||
          posBwd.y > gridUpper.y || ((posBwd.z > gridUpper.z) && flags.is3D()) ||
          flags.isObstacle(posFwd) || flags.isObstacle(posBwd)) {
        dval = fwd(i, j, k);
      }
    }
    // clampMode 2 handles flags in doClampComponent call

    dst(i, j, k) = dval;
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return vel;
  }
  typedef MACGrid type1;
  inline Grid<T> &getArg2()
  {
    return dst;
  }
  typedef Grid<T> type2;
  inline const Grid<T> &getArg3()
  {
    return orig;
  }
  typedef Grid<T> type3;
  inline const Grid<T> &getArg4()
  {
    return fwd;
  }
  typedef Grid<T> type4;
  inline Real &getArg5()
  {
    return dt;
  }
  typedef Real type5;
  inline const int &getArg6()
  {
    return clampMode;
  }
  typedef int type6;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 1; j < _maxY; j++)
            for (int i = 1; i < _maxX; i++)
              op(i, j, k, flags, vel, dst, orig, fwd, dt, clampMode);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, vel, dst, orig, fwd, dt, clampMode);
      }
    }
  }
  const FlagGrid &flags;
  const MACGrid &vel;
  Grid<T> &dst;
  const Grid<T> &orig;
  const Grid<T> &fwd;
  Real dt;
  const int clampMode;
};

//! Kernel: same as MacCormackClamp above, but specialized version for MAC grids

struct MacCormackClampMAC : public KernelBase {
  MacCormackClampMAC(const FlagGrid &flags,
                     const MACGrid &vel,
                     MACGrid &dst,
                     const MACGrid &orig,
                     const MACGrid &fwd,
                     Real dt,
                     const int clampMode)
      : KernelBase(&flags, 1),
        flags(flags),
        vel(vel),
        dst(dst),
        orig(orig),
        fwd(fwd),
        dt(dt),
        clampMode(clampMode)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 MACGrid &dst,
                 const MACGrid &orig,
                 const MACGrid &fwd,
                 Real dt,
                 const int clampMode)
  {
    Vec3 pos(i, j, k);
    Vec3 dval = dst(i, j, k);
    Vec3 dfwd = fwd(i, j, k);
    Vec3i gridUpper = flags.getSize() - 1;

    dval.x = doClampComponentMAC<0>(
        flags, gridUpper, dval.x, orig, dfwd.x, pos, vel.getAtMACX(i, j, k) * dt, clampMode);
    dval.y = doClampComponentMAC<1>(
        flags, gridUpper, dval.y, orig, dfwd.y, pos, vel.getAtMACY(i, j, k) * dt, clampMode);
    if (flags.is3D())
      dval.z = doClampComponentMAC<2>(
          flags, gridUpper, dval.z, orig, dfwd.z, pos, vel.getAtMACZ(i, j, k) * dt, clampMode);

    // note - the MAC version currently does not check whether source points were inside an
    // obstacle! (unlike centered version) this would need to be done for each face separately to
    // stay symmetric...

    dst(i, j, k) = dval;
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return vel;
  }
  typedef MACGrid type1;
  inline MACGrid &getArg2()
  {
    return dst;
  }
  typedef MACGrid type2;
  inline const MACGrid &getArg3()
  {
    return orig;
  }
  typedef MACGrid type3;
  inline const MACGrid &getArg4()
  {
    return fwd;
  }
  typedef MACGrid type4;
  inline Real &getArg5()
  {
    return dt;
  }
  typedef Real type5;
  inline const int &getArg6()
  {
    return clampMode;
  }
  typedef int type6;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 1; j < _maxY; j++)
            for (int i = 1; i < _maxX; i++)
              op(i, j, k, flags, vel, dst, orig, fwd, dt, clampMode);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, vel, dst, orig, fwd, dt, clampMode);
      }
    }
  }
  const FlagGrid &flags;
  const MACGrid &vel;
  MACGrid &dst;
  const MACGrid &orig;
  const MACGrid &fwd;
  Real dt;
  const int clampMode;
};

//! template function for performing SL advection
//! (Note boundary width only needed for specialization for MAC grids below)
template<class GridType>
void fnAdvectSemiLagrange(FluidSolver *parent,
                          const FlagGrid &flags,
                          const MACGrid &vel,
                          GridType &orig,
                          int order,
                          Real strength,
                          int orderSpace,
                          int clampMode,
                          int orderTrace)
{
  typedef typename GridType::BASETYPE T;

  Real dt = parent->getDt();
  bool levelset = orig.getType() & GridBase::TypeLevelset;

  // forward step
  GridType *fwd = new GridType(parent, true, false, false);
  SemiLagrange<T>(flags, vel, *fwd, orig, dt, levelset, orderSpace, orderTrace);

  if (order == 1) {
#if OPENMP && OPENMP_OFFLOAD
    orig.copyFrom(*fwd, true, false);
#else
    orig.swap(*fwd);
#endif
  }
  else if (order == 2) {  // MacCormack
    GridType bwd(parent);
    GridType *newGrid = new GridType(parent, true, false, false);

    // bwd <- backwards step
    SemiLagrange<T>(flags, vel, bwd, *fwd, -dt, levelset, orderSpace, orderTrace);

    // newGrid <- compute correction
    MacCormackCorrect<T>(flags, *newGrid, orig, *fwd, bwd, strength, levelset);

    // clamp values
    MacCormackClamp<T>(flags, vel, *newGrid, orig, *fwd, dt, clampMode);

#if OPENMP && OPENMP_OFFLOAD
    orig.copyFrom(*newGrid, true, false);
#else
    orig.swap(*newGrid);
#endif
    if (newGrid)
      delete newGrid;
  }
  if (fwd)
    delete fwd;
}

// outflow functions

//! calculate local propagation velocity for cell (i,j,k)
Vec3 getBulkVel(const FlagGrid &flags, const MACGrid &vel, int i, int j, int k)
{
  Vec3 avg = Vec3(0.);
  int count = 0;
  int size = 1;  // stencil size
  int nmax = (flags.is3D() ? size : 0);
  // average the neighboring fluid / outflow cell's velocity
  for (int n = -nmax; n <= nmax; n++) {
    for (int m = -size; m <= size; m++) {
      for (int l = -size; l <= size; l++) {
        if (flags.isInBounds(Vec3i(i + l, j + m, k + n)) &&
            (flags.isFluid(i + l, j + m, k + n) || flags.isOutflow(i + l, j + m, k + n))) {
          avg += vel(i + l, j + m, k + n);
          count++;
        }
      }
    }
  }
  return count > 0 ? avg / count : avg;
}

//! extrapolate normal velocity components into outflow cell
struct extrapolateVelConvectiveBC : public KernelBase {
  extrapolateVelConvectiveBC(const FlagGrid &flags,
                             const MACGrid &vel,
                             MACGrid &velDst,
                             const MACGrid &velPrev,
                             Real timeStep)
      : KernelBase(&flags, 0),
        flags(flags),
        vel(vel),
        velDst(velDst),
        velPrev(velPrev),
        timeStep(timeStep)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 MACGrid &velDst,
                 const MACGrid &velPrev,
                 Real timeStep)
  {
    if (flags.isOutflow(i, j, k)) {
      const Vec3 bulkVel = getBulkVel(flags, vel, i, j, k);
      const int dim = flags.is3D() ? 3 : 2;
      const Vec3i cur = Vec3i(i, j, k);
      Vec3i low, up, flLow, flUp;
      int cnt = 0;
      // iterate over each velocity component x, y, z
      for (int c = 0; c < dim; c++) {
        low = up = flLow = flUp = cur;
        Real factor = timeStep *
                      max((Real)1.0, abs(bulkVel[c]));  // prevent the extrapolated velocity from
                                                        // exploding when bulk velocity below 1
        low[c] = flLow[c] = cur[c] - 1;
        up[c] = flUp[c] = cur[c] + 1;
        // iterate over bWidth to allow for extrapolation into more distant outflow cells;
        // hard-coded extrapolation distance of two cells
        for (int d = 0; d < 2; d++) {
          bool extrapolateFromLower = flags.isInBounds(flLow) && flags.isFluid(flLow);
          bool extrapolateFromUpper = flags.isInBounds(flUp) && flags.isFluid(flUp);
          if (extrapolateFromLower || extrapolateFromUpper) {
            if (extrapolateFromLower) {
              velDst(i, j, k) += ((vel(i, j, k) - velPrev(i, j, k)) / factor) + vel(low);
              cnt++;
            }
            if (extrapolateFromUpper) {
              // check for cells equally far away from two fluid cells -> average value between
              // both sides
              velDst(i, j, k) += ((vel(i, j, k) - velPrev(i, j, k)) / factor) + vel(up);
              cnt++;
            }
            break;
          }
          flLow[c]--;
          flUp[c]++;
        }
      }
      if (cnt > 0)
        velDst(i, j, k) /= cnt;
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return vel;
  }
  typedef MACGrid type1;
  inline MACGrid &getArg2()
  {
    return velDst;
  }
  typedef MACGrid type2;
  inline const MACGrid &getArg3()
  {
    return velPrev;
  }
  typedef MACGrid type3;
  inline Real &getArg4()
  {
    return timeStep;
  }
  typedef Real type4;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 0; j < _maxY; j++)
            for (int i = 0; i < _maxX; i++)
              op(i, j, k, flags, vel, velDst, velPrev, timeStep);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 0; j < _maxY; j++)
          for (int i = 0; i < _maxX; i++)
            op(i, j, k, flags, vel, velDst, velPrev, timeStep);
      }
    }
  }
  const FlagGrid &flags;
  const MACGrid &vel;
  MACGrid &velDst;
  const MACGrid &velPrev;
  Real timeStep;
};

//! copy extrapolated velocity components
struct copyChangedVels : public KernelBase {
  copyChangedVels(const FlagGrid &flags, const MACGrid &velDst, MACGrid &vel)
      : KernelBase(&flags, 0), flags(flags), velDst(velDst), vel(vel)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, const FlagGrid &flags, const MACGrid &velDst, MACGrid &vel)
  {
    if (flags.isOutflow(i, j, k))
      vel(i, j, k) = velDst(i, j, k);
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return velDst;
  }
  typedef MACGrid type1;
  inline MACGrid &getArg2()
  {
    return vel;
  }
  typedef MACGrid type2;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 0; j < _maxY; j++)
            for (int i = 0; i < _maxX; i++)
              op(i, j, k, flags, velDst, vel);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 0; j < _maxY; j++)
          for (int i = 0; i < _maxX; i++)
            op(i, j, k, flags, velDst, vel);
      }
    }
  }
  const FlagGrid &flags;
  const MACGrid &velDst;
  MACGrid &vel;
};

//! extrapolate normal velocity components into open boundary cells (marked as outflow cells)
void applyOutflowBC(const FlagGrid &flags, MACGrid &vel, const MACGrid &velPrev, double timeStep)
{
  MACGrid velDst(vel.getParent());  // do not overwrite vel while it is read
  extrapolateVelConvectiveBC(flags, vel, velDst, velPrev, max(1.0, timeStep * 4));
  copyChangedVels(flags, velDst, vel);
}

// advection helpers

//! prevent parts of the surface getting "stuck" in obstacle regions
struct knResetPhiInObs : public KernelBase {
  knResetPhiInObs(const FlagGrid &flags, Grid<Real> &sdf)
      : KernelBase(&flags, 0), flags(flags), sdf(sdf)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, const FlagGrid &flags, Grid<Real> &sdf)
  {
    if (flags.isObstacle(i, j, k) && (sdf(i, j, k) < 0.)) {
      sdf(i, j, k) = 0.1;
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<Real> &getArg1()
  {
    return sdf;
  }
  typedef Grid<Real> type1;
  void runMessage(){};
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {

#pragma omp parallel
      {

#pragma omp for
        for (int k = minZ; k < maxZ; k++)
          for (int j = 0; j < _maxY; j++)
            for (int i = 0; i < _maxX; i++)
              op(i, j, k, flags, sdf);
      }
    }
    else {
      const int k = 0;
#pragma omp parallel
      {

#pragma omp for
        for (int j = 0; j < _maxY; j++)
          for (int i = 0; i < _maxX; i++)
            op(i, j, k, flags, sdf);
      }
    }
  }
  const FlagGrid &flags;
  Grid<Real> &sdf;
};
void resetPhiInObs(const FlagGrid &flags, Grid<Real> &sdf)
{
  knResetPhiInObs(flags, sdf);
}
static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "resetPhiInObs", !noTiming);
    PyObject *_retval = nullptr;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      Grid<Real> &sdf = *_args.getPtr<Grid<Real>>("sdf", 1, &_lock);
      _retval = getPyNone();
      resetPhiInObs(flags, sdf);
      _args.check();
    }
    pbFinalizePlugin(parent, "resetPhiInObs", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("resetPhiInObs", e.what());
    return 0;
  }
}
static const Pb::Register _RP_resetPhiInObs("", "resetPhiInObs", _W_0);
extern "C" {
void PbRegister_resetPhiInObs()
{
  KEEP_UNUSED(_RP_resetPhiInObs);
}
}

// advection main calls

//! template function for performing SL advection: specialized version for MAC grids
template<>
void fnAdvectSemiLagrange<MACGrid>(FluidSolver *parent,
                                   const FlagGrid &flags,
                                   const MACGrid &vel,
                                   MACGrid &orig,
                                   int order,
                                   Real strength,
                                   int orderSpace,
                                   int clampMode,
                                   int orderTrace)
{
  Real dt = parent->getDt();

  // forward step
  MACGrid *fwd = new MACGrid(parent, true, false, false);
  SemiLagrangeMAC(flags, vel, *fwd, orig, dt, orderSpace, orderTrace);

  if (orderSpace != 1) {
    debMsg("Warning higher order for MAC grids not yet implemented...", 1);
  }

  if (order == 1) {
    applyOutflowBC(flags, *fwd, orig, dt);
#if OPENMP && OPENMP_OFFLOAD
    orig.copyFrom(*fwd, true, false);
#else
    orig.swap(*fwd);
#endif
  }
  else if (order == 2) {  // MacCormack
    MACGrid bwd(parent);
    MACGrid *newGrid = new MACGrid(parent, true, false, false);

    // bwd <- backwards step
    SemiLagrangeMAC(flags, vel, bwd, *fwd, -dt, orderSpace, orderTrace);

    // newGrid <- compute correction
    MacCormackCorrectMAC<Vec3>(flags, *newGrid, orig, *fwd, bwd, strength, false, true);

    // clamp values
    MacCormackClampMAC(flags, vel, *newGrid, orig, *fwd, dt, clampMode);

    applyOutflowBC(flags, *newGrid, orig, dt);
#if OPENMP && OPENMP_OFFLOAD
    orig.copyFrom(*newGrid, true, false);
#else
    orig.swap(*newGrid);
#endif
    if (newGrid)
      delete newGrid;
  }
  if (fwd)
    delete fwd;
}

//! Perform semi-lagrangian advection of target Real- or Vec3 grid
//! Open boundary handling needs information about width of border
//! Clamping modes: 1 regular clamp leading to more overshoot and sharper results, 2 revert to 1st
//! order slightly smoother less overshoot (enable when 1 gives artifacts)

void advectSemiLagrange(const FlagGrid *flags,
                        const MACGrid *vel,
                        GridBase *grid,
                        int order = 1,
                        Real strength = 1.0,
                        int orderSpace = 1,
                        bool openBounds = false,
                        int boundaryWidth = -1,
                        int clampMode = 2,
                        int orderTrace = 1)
{
  assertMsg(order == 1 || order == 2,
            "AdvectSemiLagrange: Only order 1 (regular SL) and 2 (MacCormack) supported");
  if ((boundaryWidth != -1) || (openBounds)) {
    debMsg(
        "Warning: boundaryWidth and openBounds parameters in AdvectSemiLagrange plugin are "
        "deprecated (and have no more effect), please remove.",
        0);
  }

  // determine type of grid
  if (grid->getType() & GridBase::TypeReal) {
    fnAdvectSemiLagrange<Grid<Real>>(flags->getParent(),
                                     *flags,
                                     *vel,
                                     *((Grid<Real> *)grid),
                                     order,
                                     strength,
                                     orderSpace,
                                     clampMode,
                                     orderTrace);
  }
  else if (grid->getType() & GridBase::TypeMAC) {
    fnAdvectSemiLagrange<MACGrid>(flags->getParent(),
                                  *flags,
                                  *vel,
                                  *((MACGrid *)grid),
                                  order,
                                  strength,
                                  orderSpace,
                                  clampMode,
                                  orderTrace);
  }
  else if (grid->getType() & GridBase::TypeVec3) {
    fnAdvectSemiLagrange<Grid<Vec3>>(flags->getParent(),
                                     *flags,
                                     *vel,
                                     *((Grid<Vec3> *)grid),
                                     order,
                                     strength,
                                     orderSpace,
                                     clampMode,
                                     orderTrace);
  }
  else
    errMsg("AdvectSemiLagrange: Grid Type is not supported (only Real, Vec3, MAC, Levelset)");
}
static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "advectSemiLagrange", !noTiming);
    PyObject *_retval = nullptr;
    {
      ArgLocker _lock;
      const FlagGrid *flags = _args.getPtr<FlagGrid>("flags", 0, &_lock);
      const MACGrid *vel = _args.getPtr<MACGrid>("vel", 1, &_lock);
      GridBase *grid = _args.getPtr<GridBase>("grid", 2, &_lock);
      int order = _args.getOpt<int>("order", 3, 1, &_lock);
      Real strength = _args.getOpt<Real>("strength", 4, 1.0, &_lock);
      int orderSpace = _args.getOpt<int>("orderSpace", 5, 1, &_lock);
      bool openBounds = _args.getOpt<bool>("openBounds", 6, false, &_lock);
      int boundaryWidth = _args.getOpt<int>("boundaryWidth", 7, -1, &_lock);
      int clampMode = _args.getOpt<int>("clampMode", 8, 2, &_lock);
      int orderTrace = _args.getOpt<int>("orderTrace", 9, 1, &_lock);
      _retval = getPyNone();
      advectSemiLagrange(flags,
                         vel,
                         grid,
                         order,
                         strength,
                         orderSpace,
                         openBounds,
                         boundaryWidth,
                         clampMode,
                         orderTrace);
      _args.check();
    }
    pbFinalizePlugin(parent, "advectSemiLagrange", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("advectSemiLagrange", e.what());
    return 0;
  }
}
static const Pb::Register _RP_advectSemiLagrange("", "advectSemiLagrange", _W_1);
extern "C" {
void PbRegister_advectSemiLagrange()
{
  KEEP_UNUSED(_RP_advectSemiLagrange);
}
}

}  // namespace Manta