-
Hook PRE_FORCE_ON_ELEMENT_INTERACTIONS
-
Hook PRE_FORCE_ON_TRIANGLE_INTERACTIONS
-
Sum of halo values between GPUs
-
Hook PRE_FORCE_ON_ELEMENTS
-
Copy of halo values between GPUs
-
Hook FORCE_ON_ELEMENT_INTERACTIONS
-
Hook FORCE_ON_TRIANGLE_INTERACTIONS
```
In the `ROCKY_PLUGIN_INITIALIZE_CUDA` hook, the local module data is copied to GPU devices after being allocated the required memory block there.
115 // Pre-post force calls order: on_interactions -> on_elements
116 ROCKY_PLUGIN_PRE_FORCE_ON_SPH_ELEMENT_INTERACTIONS(device_model, sph_model, sph_interaction, _data)
117 {
118 auto data = static_cast<Data*>(_data);
119
120 // Density for interacting elements
121 auto home_element = sph_interaction.get_home_element();
122 auto near_element = sph_interaction.get_near_element();
123 auto home_scalars = home_element.get_scalars();
124 auto near_scalars = near_element.get_scalars();
125 const float home_density = home_element.get_density();
126 const float near_density = near_element.get_density();
127
128 // Distance & velocity difference between elements
129 const float3 d_diff = sph_interaction.calculate_elements_distance();
130 const float3 v_diff = sph_interaction.calculate_elements_relative_velocity();
131 float dd = get_norm(d_diff);
132
133 // Calculating the filtered strain rate per Violeau & Issa (2007)
134 const float kernel_derivative = sph_model.get_kernel_derivative(dd);
135 dd = max(dd, sph_model.get_minimum_distance());
136 const float d_strain_rate = -0.5f * sph_model.get_mass() * (home_density + near_density)
137 / (home_density * near_density) * dot(v_diff) * kernel_derivative / dd;
138
139 // Accumulating strain rate contributions on the turbulent viscosity scalar
140 home_scalars.add_scalar(data->turbulent_viscosity_index, d_strain_rate);
141 near_scalars.add_scalar(data->turbulent_viscosity_index, d_strain_rate);
142 }
143 ROCKY_PLUGIN_PRE_FORCE_ON_SPH_ELEMENT_INTERACTIONS_END()
The block of code above implements the calculation of the rate of strain tensor for the SPH elements that are in contact, which is used to calculate the turbulent viscosity values in the following steps. The rate of strain tensor, `d_shear`, is calculated according to the equation available in section 2.4 of the SPH Technical Manual.
115 ROCKY_PLUGIN_PRE_FORCE_ON_SPH_ELEMENTS(device_model, sph_model, sph_element, _data)
116 {
117 auto data = static_cast<Data*>(_data);
118
119 // Calculating turbulence coefficient
120 float coeff = data->les_smagorinsky_constant * data->les_distance_factor
121 * sph_model.get_initial_element_spacing();
122 coeff *= coeff;
123 coeff *= sph_model.get_fluid_density();
124
125 // Updating element turbulent viscosity
126 auto element_scalars = sph_element.get_scalars();
127 element_scalars.set_scalar(
128 data->turbulent_viscosity_index,
129 coeff * sqrt(2.0f * element_scalars.get_scalar(data->turbulent_viscosity_index)));
130 }
131 ROCKY_PLUGIN_PRE_FORCE_ON_SPH_ELEMENTS_END()
The block of code above implements the turbulence viscosity calculation, including the rate of strain tensor previously calculated and considering the Smagorinsky constant, the LES distance factor parameter and the initial SPH elements spacing. It represents exactly the same equation for the turbulent viscosty term calculation used in the LES model and available in section 2.3 of the SPH Technical Manual.
115 ROCKY_PLUGIN_FORCE_ON_SPH_ELEMENT_INTERACTIONS(
116 rocky_model, sph_model, sph_interaction, _data)
117 {
118 const float mass = sph_model.get_mass();
119 auto data = static_cast<Data*>(_data);
120
121 // Interacting elements
122 auto home_element = sph_interaction.get_home_element();
123 auto near_element = sph_interaction.get_near_element();
124 auto home_scalars = home_element.get_scalars();
125 auto near_scalars = near_element.get_scalars();
126
127 // Distance & relative velocity between interacting elements
128 const float3 d_diff = sph_interaction.calculate_elements_distance();
129 const float3 v_diff = sph_interaction.calculate_elements_relative_velocity();
130 float dd = get_norm(d_diff);
131 const float kernel_derivative = sph_model.get_kernel_derivative(dd);
132
133 // Pressures to density ratio
134 const float home_density = home_element.get_density();
135 const float near_density = near_element.get_density();
136 const float home_pres_dens = home_element.get_pressure() / (home_density * home_density);
137 const float near_pres_dens = near_element.get_pressure() / (near_density * near_density);
138
139 // Acceleration due to pressure forces
140 float acc = kernel_derivative * mass * (home_pres_dens + near_pres_dens);
141
142 dd = max(dd, sph_model.get_minimum_distance());
143 float3 acceleration = -acc * d_diff / dd;
144
145 // Viscous terms + LES turbulent model
146 float home_viscosity = sph_model.get_fluid_viscosity();
147 float near_viscosity = sph_model.get_fluid_viscosity();
148 if (!home_element.is_dem_coupled() && !near_element.is_dem_coupled())
149 {
150 home_viscosity += home_scalars.get_scalar(data->turbulent_viscosity_index);
151 near_viscosity += near_scalars.get_scalar(data->turbulent_viscosity_index);
152 }
153
154 /* Morris et al. (1997) approximation for the viscous term */
155 float viscosity_coefficient = mass * (home_viscosity + near_viscosity) * kernel_derivative
156 / (home_density * near_density) / dd;
157 acceleration += viscosity_coefficient * v_diff;
158
159 // Setting accelerations
160 sph_interaction.add_acceleration(acceleration);
161 }
162 ROCKY_PLUGIN_FORCE_ON_SPH_ELEMENT_INTERACTIONS_END()
The block of code above implements the LES turbulence SPH force model. The elements accelerations due to pressure and viscous forces are calculated, regarding Morris formulation for the laminar viscosity available in section 2.3.3.1 of the SPH Technical Manual. Bear in mind that this hook overrides all internal force calculations. Thus, the pressure forces are calculated besides the viscous ones, which are the only affected by the turbulent viscosity previously calculated.
115 ROCKY_PLUGIN_FORCE_ON_SPH_TRIANGLE_INTERACTIONS(
116 rocky_model, sph_model, sph_interaction, data)
117 {
118 auto home_element = sph_interaction.get_home_element();
119 const float mass = sph_model.get_mass();
120
121 /* Boundary forces on free SPH elements only
122 Viscous forces are not considered (Free Slip Boundary) */
123
124 if (home_element.is_enabled() && !home_element.is_dem_coupled())
125 {
126 // Interaction distance and velocity
127 const float distance = sph_interaction.get_distance_to_home();
128 const float3 unit_vector = sph_interaction.get_unit_vector();
129 const auto boundary_velocity = sph_interaction.get_boundary_velocity();
130
131 // Normal Acceleration
132 float force = 0.0f;
133 float dd = sph_model.get_boundary_distance_normal_factor()
134 * sph_model.get_initial_element_spacing() - distance;
135 float vnorm = 0.0f;
136 const float3 home_velocity = sph_interaction.get_home_element_velocity();
137 if (dd > 0.0f)
138 {
139 // Elastic part
140 force += sph_model.get_stiffness() * dd;
141
142 // Dissipation
143 vnorm = dot(unit_vector, home_velocity - boundary_velocity);
144 force -= sph_model.get_damping_coefficient() * vnorm;
145 force = max(0.0f, force);
146 }
147
148 // Setting forces
149 sph_interaction.add_force(force * unit_vector);
150 }
151 }
152 ROCKY_PLUGIN_FORCE_ON_SPH_TRIANGLE_INTERACTIONS_END()
The block of code above implements the normal boundary forces calculation regarding the free slip boundary condition, that uses the DEM-style interaction between the elements and the boundaries triangles. The normal repulsive forces are calculated for the SPH elements using the spring-dashpot model, with an elastic and a viscous terms, using the equations in section Specification of custom models of the SPH Technical Manual.
115 ROCKY_PLUGIN_SPH_END()
116 ROCKY_PLUGIN_END
After the module is built following the steps given in section Building Procedures, the module will be available to be used in Rocky. When the `SPH LES Turbulence Model` is enabled in a Rocky project, the models of turbulence, viscosity and boundary forces implemented on it will override built-in models in those categories.