Method step

SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp:42–191

#ifdef USE_AVX

void DragForce_Gissler2017::step()
{
	Simulation *sim = Simulation::getCurrent();
	const Real supportRadius = sim->getSupportRadius();
	const Real radius = sim->getValue<Real>(Simulation::PARTICLE_RADIUS);
	const Real diam = static_cast<Real>(2.0)*radius;
	static const Real pi = static_cast<Real>(M_PI);
	const Real rho_l = m_model->getDensity0();
	const unsigned int fluidModelIndex = m_model->getPointSetIndex();
	const unsigned int nFluids = sim->numberOfFluidModels();
	const unsigned int nBoundaries = sim->numberOfBoundaryModels();
	FluidModel *model = m_model;
	const unsigned int numParticles = m_model->numActiveParticles();
	if (numParticles == 0)
		return;

	const Real fluidParticleVolume = m_model->getVolume(0);

	// Air velocity.
	const Vector3r va(0, 0, 0);

	const Real L = cbrt(static_cast<Real>(0.75) / pi) * diam;

	const Real inv_td = static_cast<Real>(0.5)*C_d * mu_l / (rho_l * L*L);
	const Real td = static_cast<Real>(1.0) / inv_td;
	Real omegaSquare = C_k * sigma / (rho_l * L*L*L) - inv_td*inv_td;
	omegaSquare = std::max(omegaSquare, static_cast<Real>(0.0));
	const Real omega = sqrt(omegaSquare);

	// Equation (6)
	Real val = td*td*omegaSquare;
	val = sqrt(val+ static_cast<Real>(1.0)) + td*omega;
	val = std::max(val, -static_cast<Real>(0.5) * pi);
	val = std::min(val, static_cast<Real>(0.5) * pi);
	const Real t_max = -static_cast<Real>(2.0) * (atan(val) - pi) / omega;

	// Equation (7)
	const Real c_def = static_cast<Real>(1.0) - exp(-t_max / td) * (cos(omega * t_max) + static_cast<Real>(1.0)/(omega*td) * sin(omega * t_max));

	// Weber number without velocity
	const Real We_i_wo_v = rho_a * L / sigma;

	// Equation (8)
	const Real y_coeff = (C_F * We_i_wo_v * c_def) / (C_k * C_b);

	const Real n_full = 38;
	const Real n_full_23 = n_full * static_cast<Real>(2.0/3.0);
	
	#pragma omp parallel default(shared)
	{
		#pragma omp for schedule(static)  
		for (int i = 0; i < (int)numParticles; i++)
		{
			const Vector3r &vi = m_model->getVelocity(i);
			Vector3r v_i_rel = va - vi;
			const Real vi_rel_square = v_i_rel.squaredNorm();
			const Real vi_rel_norm = sqrt(vi_rel_square);
			const Real We_i = We_i_wo_v * vi_rel_square;