336 {
337 using V = typename CELL::value_t;
338 using DESCRIPTOR = typename CELL::descriptor_t;
339
340 const auto seeds =
parameters.template get<SEEDS>();
341 const auto seedsVorticity =
parameters.template get<SEEDS_VORTICITY>();
342 const auto axisDirection =
parameters.template get<AXIS_DIRECTION>();
343 const auto sigma =
parameters.template get<SIGMA>();
344
345 Vector<V,DESCRIPTOR::d> x = cell.template getField<descriptors::LOCATION>();
346 Vector<V,DESCRIPTOR::d> output = cell.template getField<U_PROFILE>();
347 const V conversionFactorVelocity =
parameters.template get<CONVERSION_FACTOR_VELOCITY>();
348 output /= conversionFactorVelocity;
349
350
351 Vector<V,3> uVortex {0, 0, 0};
352 for (std::size_t j=0; j <
parameters.template get<SEEDS_COUNT>(); ++j) {
353 Vector<V,DESCRIPTOR::d> diffX = {0,0,0};
354 for (int iD = 0; iD < DESCRIPTOR::d; ++iD) {
355 diffX[iD] = seeds[iD][j] - x[iD];
356 }
358 V
norm2 = diffX[0]*diffX[0] + diffX[1]*diffX[1] + diffX[2]*diffX[2];
360 uVortex += (V{0.5}/
M_PI * seedsVorticity[j] * crossP /
norm2 * (V{1} - expTerm)*expTerm) / conversionFactorVelocity;
361 }
362
363
366 if (normGrad == V{0}) {
367 normGrad = 1;
368 }
369 V streamFluct = -(uVortex*grad);
370 streamFluct /= normGrad;
371
372 output += uVortex + streamFluct * axisDirection;
373
374 cell.defineU(output.data());
375 }
T norm2(const std::vector< T > &a)
l2 norm to the power of 2 of a vector of arbitrary length
ADf< T, DIM > exp(const ADf< T, DIM > &a)
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
constexpr Vector< T, 3 > crossProduct3D(const ScalarVector< T, 3, IMPL > &a, const ScalarVector< T, 3, IMPL_ > &b) any_platform
Vector< V, 3 > computeVelocityGradient(CELL &cell, PARAMETERS ¶meters) any_platform
meta::list< SEEDS_COUNT, SEEDS, SEEDS_VORTICITY, CONVERSION_FACTOR_LENGTH, CONVERSION_FACTOR_VELOCITY, AXIS_DIRECTION, SIGMA > parameters