Skip to content
Snippets Groups Projects
Commit d587c8b5 authored by Samuel Kemmler's avatar Samuel Kemmler
Browse files

Merge branch 'fix_ci_tests' into 'master'

Fix race condition when building with OpenMP

See merge request walberla/walberla!726
parents 0ed52b31 41921914
Branches
No related merge requests found
...@@ -188,35 +188,56 @@ class ReduceParticleForcesSweep ...@@ -188,35 +188,56 @@ class ReduceParticleForcesSweep
// For every cell, reduce the hydrodynamic forces and torques of the overlapping particles // For every cell, reduce the hydrodynamic forces and torques of the overlapping particles
const real_t dx = block->getAABB().xSize() / real_t(nOverlappingParticlesField->xSize()); const real_t dx = block->getAABB().xSize() / real_t(nOverlappingParticlesField->xSize());
WALBERLA_FOR_ALL_CELLS_XYZ( // Do not use WALBERLA_FOR_ALL_CELLS_XYZ to avoid race condition if waLBerla is built with OpenMP, i.e., this loop
particleForcesField, const Vector3< real_t > cellCenter = // is not OpenMP parallel
Vector3< real_t >(real_t(x) + real_t(0.5) * dx, real_t(y) + real_t(0.5) * dx, #ifdef WALBERLA_BUILD_WITH_OPENMP
real_t(z) + real_t(0.5) * dx) + # pragma omp parallel for schedule(static)
block->getAABB().minCorner(); #endif
for (uint_t p = 0; p < nOverlappingParticlesField->get(x, y, z); p++) { for (cell_idx_t z = cell_idx_t(0); z < cell_idx_t(particleForcesField->zSize()); ++z)
Vector3< real_t > forceOnParticle(particleForcesField->get(x, y, z, p * 3 + 0), {
particleForcesField->get(x, y, z, p * 3 + 1), for (cell_idx_t y = cell_idx_t(0); y < cell_idx_t(particleForcesField->ySize()); ++y)
particleForcesField->get(x, y, z, p * 3 + 2)); {
forceOnParticle[0] *= forceScalingFactor; for (cell_idx_t x = cell_idx_t(0); x < cell_idx_t(particleForcesField->xSize()); ++x)
forceOnParticle[1] *= forceScalingFactor; {
forceOnParticle[2] *= forceScalingFactor; const Vector3< real_t > cellCenter =
Vector3< real_t >(real_t(x) + real_t(0.5) * dx, real_t(y) + real_t(0.5) * dx,
hydrodynamicForces[idxField->get(x, y, z, p) * 3 + 0] += forceOnParticle[0]; real_t(z) + real_t(0.5) * dx) +
hydrodynamicForces[idxField->get(x, y, z, p) * 3 + 1] += forceOnParticle[1]; block->getAABB().minCorner();
hydrodynamicForces[idxField->get(x, y, z, p) * 3 + 2] += forceOnParticle[2]; for (uint_t p = 0; p < nOverlappingParticlesField->get(x, y, z); p++)
Vector3< real_t > torqueOnParticle = {
cross(Vector3< real_t >( Vector3< real_t > forceOnParticle(particleForcesField->get(x, y, z, p * 3 + 0),
particleForcesField->get(x, y, z, p * 3 + 1),
particleForcesField->get(x, y, z, p * 3 + 2));
forceOnParticle[0] *= forceScalingFactor;
forceOnParticle[1] *= forceScalingFactor;
forceOnParticle[2] *= forceScalingFactor;
#ifdef WALBERLA_BUILD_WITH_OPENMP
# pragma omp critical
#endif
{
hydrodynamicForces[idxField->get(x, y, z, p) * 3 + 0] += forceOnParticle[0];
hydrodynamicForces[idxField->get(x, y, z, p) * 3 + 1] += forceOnParticle[1];
hydrodynamicForces[idxField->get(x, y, z, p) * 3 + 2] += forceOnParticle[2];
}
Vector3< real_t > torqueOnParticle = cross(
Vector3< real_t >(
cellCenter[0] - particleAndVolumeFractionSoA_.positions[idxField->get(x, y, z, p) * 3 + 0], cellCenter[0] - particleAndVolumeFractionSoA_.positions[idxField->get(x, y, z, p) * 3 + 0],
cellCenter[1] - particleAndVolumeFractionSoA_.positions[idxField->get(x, y, z, p) * 3 + 1], cellCenter[1] - particleAndVolumeFractionSoA_.positions[idxField->get(x, y, z, p) * 3 + 1],
cellCenter[2] - particleAndVolumeFractionSoA_.positions[idxField->get(x, y, z, p) * 3 + 2]), cellCenter[2] - particleAndVolumeFractionSoA_.positions[idxField->get(x, y, z, p) * 3 + 2]),
forceOnParticle); forceOnParticle);
hydrodynamicTorques[idxField->get(x, y, z, p) * 3 + 0] += torqueOnParticle[0]; #ifdef WALBERLA_BUILD_WITH_OPENMP
hydrodynamicTorques[idxField->get(x, y, z, p) * 3 + 1] += torqueOnParticle[1]; # pragma omp critical
hydrodynamicTorques[idxField->get(x, y, z, p) * 3 + 2] += torqueOnParticle[2]; #endif
{
hydrodynamicTorques[idxField->get(x, y, z, p) * 3 + 0] += torqueOnParticle[0];
hydrodynamicTorques[idxField->get(x, y, z, p) * 3 + 1] += torqueOnParticle[1];
hydrodynamicTorques[idxField->get(x, y, z, p) * 3 + 2] += torqueOnParticle[2];
}
}
}
} }
}
)
// Copy forces and torques of particles // Copy forces and torques of particles
size_t idxMapped = 0; size_t idxMapped = 0;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment