Commit 28dfd191 authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

Two-way coupling term in particles fixed

parent d20a7ef3

src/FlowSolverRHEA.cpp

+11 −4
Original line number Diff line number Diff line
@@ -3117,9 +3117,10 @@ void FlowSolverRHEA::timeAdvanceConservedVariables() { 
        double x_position_particle = 0.0, y_position_particle = 0.0, z_position_particle = 0.0;

        double u_velocity_fluid_particle = 0.0, v_velocity_fluid_particle = 0.0, w_velocity_fluid_particle = 0.0;

        double dynamic_viscosity_fluid = 0.0, mass_particle = 0.0, relaxation_time_particle = 0.0;

        double delta_x = 0.0, delta_y = 0.0, delta_z = 0.0, volume = 0.0;

        double f_twc_rhou = 0.0, f_twc_rhov = 0.0, f_twc_rhow = 0.0;

        int capacity = point_particles->get_prt_capacity(); 

        #pragma acc parallel loop collapse (1) private( i_local_index, j_local_index, k_local_index, x_position_particle, y_position_particle, z_position_particle, u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, relaxation_time_particle, f_twc_rhou, f_twc_rhov, f_twc_rhow ) present( this, mesh, point_particles, mesh->x[0:_lNx_], mesh->y[0:_lNy_], mesh->z[0:_lNz_], u_field.vector[0:_ls_], v_field.vector[0:_ls_], w_field.vector[0:_ls_], mu_field.vector[0:_ls_], f_rhou_field.vector[0:_ls_], f_rhov_field.vector[0:_ls_], f_rhow_field.vector[0:_ls_], point_particles->local_prts_densities[0:capacity], point_particles->local_prts_diameters[0:capacity], point_particles->local_prts_positions_0_x[0:capacity], point_particles->local_prts_positions_0_y[0:capacity], point_particles->local_prts_positions_0_z[0:capacity], point_particles->local_prts_velocities_x[0:capacity], point_particles->local_prts_velocities_y[0:capacity], point_particles->local_prts_velocities_z[0:capacity], point_particles->local_prts_velocities_0_x[0:capacity], point_particles->local_prts_velocities_0_y[0:capacity], point_particles->local_prts_velocities_0_z[0:capacity], point_particles->local_prts_indexes_0_i[0:capacity], point_particles->local_prts_indexes_0_j[0:capacity], point_particles->local_prts_indexes_0_k[0:capacity] )         

        #pragma acc parallel loop collapse (1) private( i_local_index, j_local_index, k_local_index, x_position_particle, y_position_particle, z_position_particle, u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, relaxation_time_particle, delta_x, delta_y, delta_z, volume, f_twc_rhou, f_twc_rhov, f_twc_rhow ) present( this, mesh, point_particles, mesh->x[0:_lNx_], mesh->y[0:_lNy_], mesh->z[0:_lNz_], u_field.vector[0:_ls_], v_field.vector[0:_ls_], w_field.vector[0:_ls_], mu_field.vector[0:_ls_], f_rhou_field.vector[0:_ls_], f_rhov_field.vector[0:_ls_], f_rhow_field.vector[0:_ls_], point_particles->local_prts_densities[0:capacity], point_particles->local_prts_diameters[0:capacity], point_particles->local_prts_positions_0_x[0:capacity], point_particles->local_prts_positions_0_y[0:capacity], point_particles->local_prts_positions_0_z[0:capacity], point_particles->local_prts_velocities_x[0:capacity], point_particles->local_prts_velocities_y[0:capacity], point_particles->local_prts_velocities_z[0:capacity], point_particles->local_prts_velocities_0_x[0:capacity], point_particles->local_prts_velocities_0_y[0:capacity], point_particles->local_prts_velocities_0_z[0:capacity], point_particles->local_prts_indexes_0_i[0:capacity], point_particles->local_prts_indexes_0_j[0:capacity], point_particles->local_prts_indexes_0_k[0:capacity] )         

	for( int p = 0; p < this->number_particles_local_in_use; p++ ) {



            /// Obtain Lagrangian-Eulerian indexes 0
@@ -3127,6 +3128,12 @@ void FlowSolverRHEA::timeAdvanceConservedVariables() { 
            j_local_index = point_particles->local_prts_indexes_0_j[p];	/// Local index j

            k_local_index = point_particles->local_prts_indexes_0_k[p];	/// Local index k



            /// Geometric stuff

            delta_x = 0.5*( x_field[I1D(i_local_index+1,j_local_index,k_local_index)] - x_field[I1D(i_local_index-1,j_local_index,k_local_index)] ); 

            delta_y = 0.5*( y_field[I1D(i_local_index,j_local_index+1,k_local_index)] - y_field[I1D(i_local_index,j_local_index-1,k_local_index)] ); 

            delta_z = 0.5*( z_field[I1D(i_local_index,j_local_index,k_local_index+1)] - z_field[I1D(i_local_index,j_local_index,k_local_index-1)] );

	    volume = delta_x*delta_y*delta_z;



            /// Obtain Lagrangian position 0

            x_position_particle = point_particles->local_prts_positions_0_x[p];

            y_position_particle = point_particles->local_prts_positions_0_y[p];
@@ -3144,9 +3151,9 @@ void FlowSolverRHEA::timeAdvanceConservedVariables() { 
	    relaxation_time_particle = ( point_particles->local_prts_densities[p]*pow( point_particles->local_prts_diameters[p], 2.0 ) )/( 18.0*dynamic_viscosity_fluid );



            /// Calculate two-way coupling momentum forces

            f_twc_rhou = ( mass_particle/relaxation_time_particle )*( point_particles->local_prts_velocities_0_x[p]- u_velocity_fluid_particle );

            f_twc_rhov = ( mass_particle/relaxation_time_particle )*( point_particles->local_prts_velocities_0_y[p] - v_velocity_fluid_particle );

            f_twc_rhow = ( mass_particle/relaxation_time_particle )*( point_particles->local_prts_velocities_0_z[p] - w_velocity_fluid_particle );

            f_twc_rhou = ( mass_particle/relaxation_time_particle )*( point_particles->local_prts_velocities_0_x[p] - u_velocity_fluid_particle )/volume;

            f_twc_rhov = ( mass_particle/relaxation_time_particle )*( point_particles->local_prts_velocities_0_y[p] - v_velocity_fluid_particle )/volume;

            f_twc_rhow = ( mass_particle/relaxation_time_particle )*( point_particles->local_prts_velocities_0_z[p] - w_velocity_fluid_particle )/volume;

            

            /// Update momentum source terms

            #pragma acc atomic update