Switch BKE_smoke from OMP to BLI_task.

Gives a global speedup of about 5% in smoke simulation (as usual, parallelized chunks themselves
are about 15-25% quicker with BLI_task than with OMP), using a simple setup with two generators
(one from mesh, one from particles), an obstacle and a windfield.

1 files changed, 333 insertions, 217 deletions

diff --git a/source/blender/blenkernel/intern/smoke.c b/source/blender/blenkernel/intern/smoke.c
index 09cda606fea..ad2736a6a75 100644
--- a/source/blender/blenkernel/intern/smoke.c
+++ b/source/blender/blenkernel/intern/smoke.c
@@ -44,6 +44,7 @@
 #include "BLI_math.h"
 #include "BLI_kdtree.h"
 #include "BLI_kdopbvh.h"
+#include "BLI_task.h"
 #include "BLI_threads.h"
 #include "BLI_utildefines.h"
 #include "BLI_voxel.h"
@@ -723,7 +724,75 @@ static int get_lamp(Scene *scene, float *light)
  *	Obstacles
  **********************************************************/
 
-static void obstacles_from_derivedmesh(Object *coll_ob, SmokeDomainSettings *sds, SmokeCollSettings *scs, unsigned char *obstacle_map, float *velocityX, float *velocityY, float *velocityZ, float dt)
+typedef struct ObstaclesFromDMData {
+	SmokeDomainSettings *sds;
+	const MVert *mvert;
+	const MLoop *mloop;
+	const MLoopTri *looptri;
+	BVHTreeFromMesh *tree;
+	unsigned char *obstacle_map;
+
+	bool has_velocity;
+	float *vert_vel;
+	float *velocityX, *velocityY, *velocityZ;
+} ObstaclesFromDMData;
+
+static void obstacles_from_derivedmesh_task_cb(void *userdata, void *UNUSED(userdata_chunk), int z)
+{
+	ObstaclesFromDMData *data = userdata;
+	SmokeDomainSettings *sds = data->sds;
+
+	/* slightly rounded-up sqrt(3 * (0.5)^2) == max. distance of cell boundary along the diagonal */
+	const float surface_distance = 0.867f;
+
+	for (int x = sds->res_min[0]; x < sds->res_max[0]; x++) {
+		for (int y = sds->res_min[1]; y < sds->res_max[1]; y++) {
+			const int index = smoke_get_index(x - sds->res_min[0], sds->res[0], y - sds->res_min[1], sds->res[1], z - sds->res_min[2]);
+
+			float ray_start[3] = {(float)x + 0.5f, (float)y + 0.5f, (float)z + 0.5f};
+			BVHTreeNearest nearest = {0};
+			nearest.index = -1;
+			nearest.dist_sq = surface_distance * surface_distance; /* find_nearest uses squared distance */
+
+			/* find the nearest point on the mesh */
+			if (BLI_bvhtree_find_nearest(data->tree->tree, ray_start, &nearest, data->tree->nearest_callback, data->tree) != -1) {
+				const MLoopTri *lt = &data->looptri[nearest.index];
+				float weights[4];
+				int v1, v2, v3;
+
+				/* calculate barycentric weights for nearest point */
+				v1 = data->mloop[lt->tri[0]].v;
+				v2 = data->mloop[lt->tri[1]].v;
+				v3 = data->mloop[lt->tri[2]].v;
+				interp_weights_face_v3(
+				            weights, data->mvert[v1].co, data->mvert[v2].co, data->mvert[v3].co, NULL, nearest.co);
+
+				// DG TODO
+				if (data->has_velocity)
+				{
+					/* apply object velocity */
+					{
+						float hit_vel[3];
+						interp_v3_v3v3v3(hit_vel, &data->vert_vel[v1 * 3], &data->vert_vel[v2 * 3], &data->vert_vel[v3 * 3], weights);
+						data->velocityX[index] += hit_vel[0];
+						data->velocityY[index] += hit_vel[1];
+						data->velocityZ[index] += hit_vel[2];
+					}
+				}
+
+				/* tag obstacle cells */
+				data->obstacle_map[index] = 1;
+
+				if (data->has_velocity)
+					data->obstacle_map[index] |= 8;
+			}
+		}
+	}
+}
+
+static void obstacles_from_derivedmesh(
+        Object *coll_ob, SmokeDomainSettings *sds, SmokeCollSettings *scs,
+        unsigned char *obstacle_map, float *velocityX, float *velocityY, float *velocityZ, float dt)
 {
 	if (!scs->dm) return;
 	{
@@ -732,13 +801,10 @@ static void obstacles_from_derivedmesh(Object *coll_ob, SmokeDomainSettings *sds
 		const MLoopTri *looptri;
 		const MLoop *mloop;
 		BVHTreeFromMesh treeData = {NULL};
-		int numverts, i, z;
-
-		/* slightly rounded-up sqrt(3 * (0.5)^2) == max. distance of cell boundary along the diagonal */
-		const float surface_distance = 0.867f;
+		int numverts, i;
 
 		float *vert_vel = NULL;
-		int has_velocity = 0;
+		bool has_velocity = false;
 
 		tstart();
 
@@ -763,7 +829,7 @@ static void obstacles_from_derivedmesh(Object *coll_ob, SmokeDomainSettings *sds
 				scs->numverts = numverts;
 			}
 			else {
-				has_velocity = 1;
+				has_velocity = true;
 			}
 		}
 
@@ -795,51 +861,14 @@ static void obstacles_from_derivedmesh(Object *coll_ob, SmokeDomainSettings *sds
 		}
 
 		if (bvhtree_from_mesh_looptri(&treeData, dm, 0.0f, 4, 6)) {
-#pragma omp parallel for schedule(static)
-			for (z = sds->res_min[2]; z < sds->res_max[2]; z++) {
-				int x, y;
-				for (x = sds->res_min[0]; x < sds->res_max[0]; x++)
-					for (y = sds->res_min[1]; y < sds->res_max[1]; y++) {
-						int index = smoke_get_index(x - sds->res_min[0], sds->res[0], y - sds->res_min[1], sds->res[1], z - sds->res_min[2]);
-
-						float ray_start[3] = {(float)x + 0.5f, (float)y + 0.5f, (float)z + 0.5f};
-						BVHTreeNearest nearest = {0};
-						nearest.index = -1;
-						nearest.dist_sq = surface_distance * surface_distance; /* find_nearest uses squared distance */
-
-						/* find the nearest point on the mesh */
-						if (BLI_bvhtree_find_nearest(treeData.tree, ray_start, &nearest, treeData.nearest_callback, &treeData) != -1) {
-							const MLoopTri *lt = &looptri[nearest.index];
-							float weights[4];
-							int v1, v2, v3;
-
-							/* calculate barycentric weights for nearest point */
-							v1 = mloop[lt->tri[0]].v;
-							v2 = mloop[lt->tri[1]].v;
-							v3 = mloop[lt->tri[2]].v;
-							interp_weights_face_v3(weights, mvert[v1].co, mvert[v2].co, mvert[v3].co, NULL, nearest.co);
-
-							// DG TODO
-							if (has_velocity)
-							{
-								/* apply object velocity */
-								{
-									float hit_vel[3];
-									interp_v3_v3v3v3(hit_vel, &vert_vel[v1 * 3], &vert_vel[v2 * 3], &vert_vel[v3 * 3], weights);
-									velocityX[index] += hit_vel[0];
-									velocityY[index] += hit_vel[1];
-									velocityZ[index] += hit_vel[2];
-								}
-							}
-
-							/* tag obstacle cells */
-							obstacle_map[index] = 1;
-
-							if (has_velocity)
-								obstacle_map[index] |= 8;
-						}
-					}
-			}
+			ObstaclesFromDMData data = {
+			    .sds = sds, .mvert = mvert, .mloop = mloop, .looptri = looptri,
+			    .tree = &treeData, .obstacle_map = obstacle_map,
+			    .has_velocity = has_velocity, .vert_vel = vert_vel,
+			    .velocityX = velocityX, .velocityY = velocityY, .velocityZ = velocityZ
+			};
+			BLI_task_parallel_range_ex(
+			            sds->res_min[2], sds->res_max[2], &data, NULL, 0, obstacles_from_derivedmesh_task_cb, 0, false);
 		}
 		/* free bvh tree */
 		free_bvhtree_from_mesh(&treeData);
@@ -1220,8 +1249,85 @@ static void em_combineMaps(EmissionMap *output, EmissionMap *em2, int hires_mult
 	em_freeData(&em1);
 }
 
+typedef struct EmitFromParticlesData {
+	SmokeFlowSettings *sfs;
+	KDTree *tree;
+	int hires_multiplier;
+
+	EmissionMap *em;
+	float *particle_vel;
+	float hr;
+
+	int *min, *max, *res;
+
+	float solid;
+	float smooth;
+	float hr_smooth;
+} EmitFromParticlesData;
+
+static void emit_from_particles_task_cb(void *userdata, void *UNUSED(userdata_chunk), int z)
+{
+	EmitFromParticlesData *data = userdata;
+	SmokeFlowSettings *sfs = data->sfs;
+	EmissionMap *em = data->em;
+	const int hires_multiplier = data->hires_multiplier;
+
+	for (int x = data->min[0]; x < data->max[0]; x++) {
+		for (int y = data->min[1]; y < data->max[1]; y++) {
+			/* take low res samples where possible */
+			if (hires_multiplier <= 1 || !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
+				/* get low res space coordinates */
+				const int lx = x / hires_multiplier;
+				const int ly = y / hires_multiplier;
+				const int lz = z / hires_multiplier;
+
+				const int index = smoke_get_index(lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
+				const float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
+
+				/* find particle distance from the kdtree */
+				KDTreeNearest nearest;
+				const float range = data->solid + data->smooth;
+				BLI_kdtree_find_nearest(data->tree, ray_start, &nearest);
+
+				if (nearest.dist < range) {
+					em->influence[index] = (nearest.dist < data->solid) ?
+					                       1.0f : (1.0f - (nearest.dist - data->solid) / data->smooth);
+					/* Uses particle velocity as initial velocity for smoke */
+					if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (sfs->psys->part->phystype != PART_PHYS_NO)) {
+						VECADDFAC(&em->velocity[index * 3], &em->velocity[index * 3],
+						          &data->particle_vel[nearest.index * 3], sfs->vel_multi);
+					}
+				}
+			}
+
+			/* take high res samples if required */
+			if (hires_multiplier > 1) {
+				/* get low res space coordinates */
+				const float lx = ((float)x) * data->hr;
+				const float ly = ((float)y) * data->hr;
+				const float lz = ((float)z) * data->hr;
+
+				const int index = smoke_get_index(
+				                      x - data->min[0], data->res[0], y - data->min[1], data->res[1], z - data->min[2]);
+				const float ray_start[3] = {lx + 0.5f * data->hr, ly + 0.5f * data->hr, lz + 0.5f * data->hr};
+
+				/* find particle distance from the kdtree */
+				KDTreeNearest nearest;
+				const float range = data->solid + data->hr_smooth;
+				BLI_kdtree_find_nearest(data->tree, ray_start, &nearest);
+
+				if (nearest.dist < range) {
+					em->influence_high[index] = (nearest.dist < data->solid) ?
+					                            1.0f : (1.0f - (nearest.dist - data->solid) / data->smooth);
+				}
+			}
+
+		}
+	}
+}
 
-static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, Scene *scene, float dt)
+static void emit_from_particles(
+        Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, Scene *scene, float dt)
 {
 	if (sfs && sfs->psys && sfs->psys->part && ELEM(sfs->psys->part->type, PART_EMITTER, PART_FLUID)) // is particle system selected
 	{
@@ -1235,11 +1341,10 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 		int bounds_margin = 1;
 
 		/* radius based flow */
-		float solid = sfs->particle_size * 0.5f;
-		float smooth = 0.5f; /* add 0.5 cells of linear falloff to reduce aliasing */
+		const float solid = sfs->particle_size * 0.5f;
+		const float smooth = 0.5f; /* add 0.5 cells of linear falloff to reduce aliasing */
 		int hires_multiplier = 1;
-		int z;
-		KDTree *tree;
+		KDTree *tree = NULL;
 
 		sim.scene = scene;
 		sim.ob = flow_ob;
@@ -1350,9 +1455,8 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 			}   // particles loop
 		}
 		else if (valid_particles > 0) { // MOD_SMOKE_FLOW_USE_PART_SIZE
-
 			int min[3], max[3], res[3];
-			float hr = 1.0f / ((float)hires_multiplier);
+			const float hr = 1.0f / ((float)hires_multiplier);
 			/* slightly adjust high res antialias smoothness based on number of divisions
 			 * to allow smaller details but yet not differing too much from the low res size */
 			const float hr_smooth = smooth * powf(hr, 1.0f / 3.0f);
@@ -1366,62 +1470,13 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 
 			BLI_kdtree_balance(tree);
 
-			/* begin thread safe malloc */
-			BLI_begin_threaded_malloc();
-
-#pragma omp parallel for schedule(static)
-			for (z = min[2]; z < max[2]; z++) {
-				int x, y;
-				for (x = min[0]; x < max[0]; x++)
-					for (y = min[1]; y < max[1]; y++) {
-						/* take low res samples where possible */
-						if (hires_multiplier <= 1 || !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
-							/* get low res space coordinates */
-							int lx = x / hires_multiplier;
-							int ly = y / hires_multiplier;
-							int lz = z / hires_multiplier;
-
-							int index = smoke_get_index(lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
-							float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
-
-							/* find particle distance from the kdtree */
-							KDTreeNearest nearest;
-							float range = solid + smooth;
-							BLI_kdtree_find_nearest(tree, ray_start, &nearest);
-
-							if (nearest.dist < range) {
-								em->influence[index] = (nearest.dist < solid) ? 1.0f : (1.0f - (nearest.dist-solid) / smooth);
-								/* Uses particle velocity as initial velocity for smoke */
-								if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO))
-								{
-									VECADDFAC(&em->velocity[index * 3], &em->velocity[index * 3], &particle_vel[nearest.index * 3], sfs->vel_multi);
-								}
-							}
-						}
-
-						/* take high res samples if required */
-						if (hires_multiplier > 1) {
-							/* get low res space coordinates */
-							float lx = ((float)x) * hr;
-							float ly = ((float)y) * hr;
-							float lz = ((float)z) * hr;
+			EmitFromParticlesData data = {
+				.sfs = sfs, .tree = tree, .hires_multiplier = hires_multiplier, .hr = hr,
+			    .em = em, .particle_vel = particle_vel, .min = min, .max = max, .res = res,
+			    .solid = solid, .smooth = smooth, .hr_smooth = hr_smooth,
+			};
 
-							int index = smoke_get_index(x - min[0], res[0], y - min[1], res[1], z - min[2]);
-							float ray_start[3] = {lx + 0.5f*hr, ly + 0.5f*hr, lz + 0.5f*hr};
-
-							/* find particle distance from the kdtree */
-							KDTreeNearest nearest;
-							float range = solid + hr_smooth;
-							BLI_kdtree_find_nearest(tree, ray_start, &nearest);
-
-							if (nearest.dist < range) {
-								em->influence_high[index] = (nearest.dist < solid) ? 1.0f : (1.0f - (nearest.dist-solid) / smooth);
-							}
-						}
-
-					}
-			}
-			BLI_end_threaded_malloc();
+			BLI_task_parallel_range_ex(min[2], max[2], &data, NULL, 0, emit_from_particles_task_cb, 0, false);
 		}
 
 		if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
@@ -1571,6 +1626,76 @@ static void sample_derivedmesh(
 	influence_map[index] = MAX2(volume_factor, sample_str);
 }
 
+typedef struct EmitFromDMData {
+	SmokeDomainSettings *sds;
+	SmokeFlowSettings *sfs;
+	const MVert *mvert;
+	const MLoop *mloop;
+	const MLoopTri *mlooptri;
+	const MLoopUV *mloopuv;
+	MDeformVert *dvert;
+	int defgrp_index;
+
+	BVHTreeFromMesh *tree;
+	int hires_multiplier;
+	float hr;
+
+	EmissionMap *em;
+	bool has_velocity;
+	float *vert_vel;
+
+	float *flow_center;
+	int *min, *max, *res;
+} EmitFromDMData;
+
+static void emit_from_derivedmesh_task_cb(void *userdata, void *UNUSED(userdata_chunk), int z)
+{
+	EmitFromDMData *data = userdata;
+	EmissionMap *em = data->em;
+	const int hires_multiplier = data->hires_multiplier;
+
+	for (int x = data->min[0]; x < data->max[0]; x++) {
+		for (int y = data->min[1]; y < data->max[1]; y++) {
+			/* take low res samples where possible */
+			if (hires_multiplier <= 1 || !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
+				/* get low res space coordinates */
+				const int lx = x / hires_multiplier;
+				const int ly = y / hires_multiplier;
+				const int lz = z / hires_multiplier;
+
+				const int index = smoke_get_index(
+				                      lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
+				const float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
+
+				sample_derivedmesh(
+				        data->sfs, data->mvert, data->mloop, data->mlooptri, data->mloopuv,
+				        em->influence, em->velocity, index, data->sds->base_res, data->flow_center,
+				        data->tree, ray_start, data->vert_vel, data->has_velocity, data->defgrp_index, data->dvert,
+				        (float)lx, (float)ly, (float)lz);
+			}
+
+			/* take high res samples if required */
+			if (hires_multiplier > 1) {
+				/* get low res space coordinates */
+				const float lx = ((float)x) * data->hr;
+				const float ly = ((float)y) * data->hr;
+				const float lz = ((float)z) * data->hr;
+
+				const int index = smoke_get_index(
+				                      x - data->min[0], data->res[0], y - data->min[1], data->res[1], z - data->min[2]);
+				const float ray_start[3] = {lx + 0.5f * data->hr, ly + 0.5f * data->hr, lz + 0.5f * data->hr};
+
+				sample_derivedmesh(
+				        data->sfs, data->mvert, data->mloop, data->mlooptri, data->mloopuv,
+				        em->influence_high, NULL, index, data->sds->base_res, data->flow_center,
+				        data->tree, ray_start, data->vert_vel, data->has_velocity, data->defgrp_index, data->dvert,
+				        /* x,y,z needs to be always lowres */
+				        lx, ly, lz);
+			}
+		}
+	}
+}
+
 static void emit_from_derivedmesh(Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, float dt)
 {
 	if (sfs->dm) {
@@ -1583,7 +1708,7 @@ static void emit_from_derivedmesh(Object *flow_ob, SmokeDomainSettings *sds, Smo
 		const MLoopUV *mloopuv = NULL;
 		const MLoop *mloop = NULL;
 		BVHTreeFromMesh treeData = {NULL};
-		int numOfVerts, i, z;
+		int numOfVerts, i;
 		float flow_center[3] = {0};
 
 		float *vert_vel = NULL;
@@ -1665,49 +1790,18 @@ static void emit_from_derivedmesh(Object *flow_ob, SmokeDomainSettings *sds, Smo
 		}
 
 		if (bvhtree_from_mesh_looptri(&treeData, dm, 0.0f, 4, 6)) {
-#pragma omp parallel for schedule(static)
-			for (z = min[2]; z < max[2]; z++) {
-				int x, y;
-				for (x = min[0]; x < max[0]; x++)
-					for (y = min[1]; y < max[1]; y++) {
-						/* take low res samples where possible */
-						if (hires_multiplier <= 1 || !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
-							/* get low res space coordinates */
-							int lx = x / hires_multiplier;
-							int ly = y / hires_multiplier;
-							int lz = z / hires_multiplier;
-
-							int index = smoke_get_index(lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
-							float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
-
-							sample_derivedmesh(
-							        sfs, mvert, mloop, mlooptri, mloopuv,
-							        em->influence, em->velocity, index, sds->base_res, flow_center,
-							        &treeData, ray_start,vert_vel, has_velocity, defgrp_index, dvert,
-							        (float)lx, (float)ly, (float)lz);
-						}
+			const float hr = 1.0f / ((float)hires_multiplier);
 
-						/* take high res samples if required */
-						if (hires_multiplier > 1) {
-							/* get low res space coordinates */
-							float hr = 1.0f / ((float)hires_multiplier);
-							float lx = ((float)x) * hr;
-							float ly = ((float)y) * hr;
-							float lz = ((float)z) * hr;
-
-							int index = smoke_get_index(x - min[0], res[0], y - min[1], res[1], z - min[2]);
-							float ray_start[3] = {lx + 0.5f*hr, ly + 0.5f*hr, lz + 0.5f*hr};
-
-							sample_derivedmesh(
-							        sfs, mvert, mloop, mlooptri, mloopuv,
-							        em->influence_high, NULL, index, sds->base_res, flow_center,
-							        &treeData, ray_start, vert_vel, has_velocity, defgrp_index, dvert,
-							        /* x,y,z needs to be always lowres */
-							        lx, ly, lz);
-						}
+			EmitFromDMData data = {
+				.sds = sds, .sfs = sfs,
+			    .mvert = mvert, .mloop = mloop, .mlooptri = mlooptri, .mloopuv = mloopuv,
+			    .dvert = dvert, .defgrp_index = defgrp_index,
+			    .tree = &treeData, .hires_multiplier = hires_multiplier, .hr = hr,
+			    .em = em, .has_velocity = has_velocity, .vert_vel = vert_vel,
+			    .flow_center = flow_center, .min = min, .max = max, .res = res,
+			};
 
-					}
-			}
+			BLI_task_parallel_range_ex(min[2], max[2], &data, NULL, 0, emit_from_derivedmesh_task_cb, 0, false);
 		}
 		/* free bvh tree */
 		free_bvhtree_from_mesh(&treeData);
@@ -2405,6 +2499,73 @@ static void update_flowsfluids(Scene *scene, Object *ob, SmokeDomainSettings *sd
 		MEM_freeN(emaps);
 }
 
+typedef struct UpdateEffectorsData {
+	Scene *scene;
+	SmokeDomainSettings *sds;
+	ListBase *effectors;
+
+	float *density;
+	float *fuel;
+	float *force_x;
+	float *force_y;
+	float *force_z;
+	float *velocity_x;
+	float *velocity_y;
+	float *velocity_z;
+	unsigned char *obstacle;
+} UpdateEffectorsData;
+
+static void update_effectors_task_cb(void *userdata, void *UNUSED(userdata_chunk), int x)
+{
+	UpdateEffectorsData *data = userdata;
+	SmokeDomainSettings *sds = data->sds;
+
+	for (int y = 0; y < sds->res[1]; y++) {
+		for (int z = 0; z < sds->res[2]; z++)
+		{
+			EffectedPoint epoint;
+			float mag;
+			float voxelCenter[3] = {0, 0, 0}, vel[3] = {0, 0, 0}, retvel[3] = {0, 0, 0};
+			const unsigned int index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
+
+			if (((data->fuel ? MAX2(data->density[index], data->fuel[index]) : data->density[index]) < FLT_EPSILON) ||
+			    data->obstacle[index])
+			{
+				continue;
+			}
+
+			vel[0] = data->velocity_x[index];
+			vel[1] = data->velocity_y[index];
+			vel[2] = data->velocity_z[index];
+
+			/* convert vel to global space */
+			mag = len_v3(vel);
+			mul_mat3_m4_v3(sds->obmat, vel);
+			normalize_v3(vel);
+			mul_v3_fl(vel, mag);
+
+			voxelCenter[0] = sds->p0[0] + sds->cell_size[0] * ((float)(x + sds->res_min[0]) + 0.5f);
+			voxelCenter[1] = sds->p0[1] + sds->cell_size[1] * ((float)(y + sds->res_min[1]) + 0.5f);
+			voxelCenter[2] = sds->p0[2] + sds->cell_size[2] * ((float)(z + sds->res_min[2]) + 0.5f);
+			mul_m4_v3(sds->obmat, voxelCenter);
+
+			pd_point_from_loc(data->scene, voxelCenter, vel, index, &epoint);
+			pdDoEffectors(data->effectors, NULL, sds->effector_weights, &epoint, retvel, NULL);
+
+			/* convert retvel to local space */
+			mag = len_v3(retvel);
+			mul_mat3_m4_v3(sds->imat, retvel);
+			normalize_v3(retvel);
+			mul_v3_fl(retvel, mag);
+
+			// TODO dg - do in force!
+			data->force_x[index] = min_ff(max_ff(-1.0f, retvel[0] * 0.2f), 1.0f);
+			data->force_y[index] = min_ff(max_ff(-1.0f, retvel[1] * 0.2f), 1.0f);
+			data->force_z[index] = min_ff(max_ff(-1.0f, retvel[2] * 0.2f), 1.0f);
+		}
+	}
+}
+
 static void update_effectors(Scene *scene, Object *ob, SmokeDomainSettings *sds, float UNUSED(dt))
 {
 	ListBase *effectors;
@@ -2412,65 +2573,23 @@ static void update_effectors(Scene *scene, Object *ob, SmokeDomainSettings *sds,
 	sds->effector_weights->weight[PFIELD_SMOKEFLOW] = 0.0f;
 	effectors = pdInitEffectors(scene, ob, NULL, sds->effector_weights, true);
 
-	if (effectors)
-	{
-		float *density = smoke_get_density(sds->fluid);
-		float *fuel = smoke_get_fuel(sds->fluid);
-		float *force_x = smoke_get_force_x(sds->fluid);
-		float *force_y = smoke_get_force_y(sds->fluid);
-		float *force_z = smoke_get_force_z(sds->fluid);
-		float *velocity_x = smoke_get_velocity_x(sds->fluid);
-		float *velocity_y = smoke_get_velocity_y(sds->fluid);
-		float *velocity_z = smoke_get_velocity_z(sds->fluid);
-		unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
-		int x;
-
+	if (effectors) {
 		// precalculate wind forces
-#pragma omp parallel for schedule(static)
-		for (x = 0; x < sds->res[0]; x++)
-		{
-			int y, z;
-			for (y = 0; y < sds->res[1]; y++)
-				for (z = 0; z < sds->res[2]; z++)
-				{
-					EffectedPoint epoint;
-					float mag;
-					float voxelCenter[3] = {0, 0, 0}, vel[3] = {0, 0, 0}, retvel[3] = {0, 0, 0};
-					unsigned int index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
-
-					if (((fuel ? MAX2(density[index], fuel[index]) : density[index]) < FLT_EPSILON) || obstacle[index])
-						continue;
-
-					vel[0] = velocity_x[index];
-					vel[1] = velocity_y[index];
-					vel[2] = velocity_z[index];
-
-					/* convert vel to global space */
-					mag = len_v3(vel);
-					mul_mat3_m4_v3(sds->obmat, vel);
-					normalize_v3(vel);
-					mul_v3_fl(vel, mag);
-
-					voxelCenter[0] = sds->p0[0] + sds->cell_size[0] * ((float)(x + sds->res_min[0]) + 0.5f);
-					voxelCenter[1] = sds->p0[1] + sds->cell_size[1] * ((float)(y + sds->res_min[1]) + 0.5f);
-					voxelCenter[2] = sds->p0[2] + sds->cell_size[2] * ((float)(z + sds->res_min[2]) + 0.5f);
-					mul_m4_v3(sds->obmat, voxelCenter);
-
-					pd_point_from_loc(scene, voxelCenter, vel, index, &epoint);
-					pdDoEffectors(effectors, NULL, sds->effector_weights, &epoint, retvel, NULL);
-
-					/* convert retvel to local space */
-					mag = len_v3(retvel);
-					mul_mat3_m4_v3(sds->imat, retvel);
-					normalize_v3(retvel);
-					mul_v3_fl(retvel, mag);
-
-					// TODO dg - do in force!
-					force_x[index] = min_ff(max_ff(-1.0f, retvel[0] * 0.2f), 1.0f);
-					force_y[index] = min_ff(max_ff(-1.0f, retvel[1] * 0.2f), 1.0f);
-					force_z[index] = min_ff(max_ff(-1.0f, retvel[2] * 0.2f), 1.0f);
-				}
-		}
+		UpdateEffectorsData data;
+		data.scene = scene;
+		data.sds = sds;
+		data.effectors = effectors;
+		data.density = smoke_get_density(sds->fluid);
+		data.fuel = smoke_get_fuel(sds->fluid);
+		data.force_x = smoke_get_force_x(sds->fluid);
+		data.force_y = smoke_get_force_y(sds->fluid);
+		data.force_z = smoke_get_force_z(sds->fluid);
+		data.velocity_x = smoke_get_velocity_x(sds->fluid);
+		data.velocity_y = smoke_get_velocity_y(sds->fluid);
+		data.velocity_z = smoke_get_velocity_z(sds->fluid);
+		data.obstacle = smoke_get_obstacle(sds->fluid);
+
+		BLI_task_parallel_range_ex(0, sds->res[0], &data, NULL, 0, update_effectors_task_cb, 0, false);
 	}
 
 	pdEndEffectors(&effectors);
@@ -2817,7 +2936,6 @@ static float calc_voxel_transp(float *result, float *input, int res[3], int *pix
 
 	if (result[index] < 0.0f)
 	{
-// #pragma omp critical
 		result[index] = *tRay;
 	}
 
@@ -2932,7 +3050,6 @@ static void smoke_calc_transparency(SmokeDomainSettings *sds, Scene *scene)
 	bv[3] = (float)sds->res[1]; // y
 	bv[5] = (float)sds->res[2]; // z
 
-// #pragma omp parallel for schedule(static, 1)
 	for (z = 0; z < sds->res[2]; z++)
 	{
 		size_t index = z * slabsize;
@@ -2974,7 +3091,6 @@ static void smoke_calc_transparency(SmokeDomainSettings *sds, Scene *scene)
 				bresenham_linie_3D(cell[0], cell[1], cell[2], x, y, z, &tRay, calc_voxel_transp, sds->shadow, density, sds->res, correct);
 
 				// convention -> from a RGBA float array, use G value for tRay
-// #pragma omp critical
 				sds->shadow[index] = tRay;
 			}
 	}