#ifdef __okl__
#endif

void UDF_Setup() {
  auto mesh = nrs->mesh;
  auto [x, y, z] = mesh->xyzHost();
  std::vector<dfloat> U(mesh->dim * nrs->fieldOffset, 0.0);

  const dfloat A{4.5}, B{3.5}, C{1./6};

  // mesh modification
for(int i{0}; i < mesh -> Nlocal; ++i) {
    const dfloat argx{B * (std::abs(x[i] - A) - B)};
    const dfloat A1{C * (1. + std::tanh(argx))};
    y[i] = y[i] + A1 * (3. - y[i]);
  }

  // set initial conditions
  for(int n{0}; n < mesh -> Nlocal; ++n) {
    U[n + 0 * nrs->fieldOffset] = 1.0;
    U[n + 1 * nrs->fieldOffset] = 0.0;
    U[n + 2 * nrs->fieldOffset] = 0.0;
  }
}

void UDF_ExecuteStep(double time, int tstep) {}