Commit 24ad6bb0 authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

Point particles adapted for OpenACC III

parent 4ab6f1aa

myRHEA.cpp

+20 −4
Original line number Diff line number Diff line
@@ -182,6 +182,8 @@ void myRHEA::timeAdvanceVelocityPointParticles() { 
    /// IMPORTANT: This method needs to be modified/overwritten according to the problem under consideration



    /// Explicit Euler time-integration of particles velocity

    int i_local_index, j_local_index, k_local_index;

    double x_position_particle, y_position_particle, z_position_particle;

    double u_velocity_particle, v_velocity_particle, w_velocity_particle;

    double u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, relaxation_time_particle;

    delta_t = this->delta_t;
@@ -189,16 +191,30 @@ void myRHEA::timeAdvanceVelocityPointParticles() { 
    #pragma acc parallel loop collapse (1) private(x_position_particle, y_position_particle, z_position_particle, u_velocity_particle, v_velocity_particle, w_velocity_particle, relaxation_time_particle, dynamic_viscosity_fluid, u_velocity_fluid_particle,w_velocity_fluid_particle,v_velocity_fluid_particle) present(this, mesh, u_field.vector[0:_ls_], v_field.vector[0:_ls_], w_field.vector[0:_ls_], mu_field.vector[0:_ls_], point_particles, point_particles->local_prts_positions_x[0:capacity], point_particles->local_prts_positions_y[0:capacity], point_particles->local_prts_positions_z[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]) copyin(delta_t)

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



        /// Obtain Lagrangian values

        /// Obtain Lagrangian-Eulerian indexes 0

        i_local_index = point_particles->local_prts_indexes_0_i[p];

        j_local_index = point_particles->local_prts_indexes_0_j[p];

        k_local_index = point_particles->local_prts_indexes_0_k[p];



        /// 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];

        z_position_particle = point_particles->local_prts_positions_0_z[p];



        /// Obtain Lagrangian velocities 0

	u_velocity_particle = point_particles->local_prts_velocities_0_x[p];

        v_velocity_particle = point_particles->local_prts_velocities_0_y[p];

        w_velocity_particle = point_particles->local_prts_velocities_0_z[p];



        /// Obtain Lagrangian-Euler values

        this->obtainLagrangianEulerianVelocityDynamicViscosityValues( u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, p );

        /// Interpolate (trilinear) values: u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid 

        u_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

        v_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

        w_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );



        /// Calculate relaxation time particle

        relaxation_time_particle = point_particles->calculate_relaxation_time_prt( p, dynamic_viscosity_fluid );



        /// Update velocity

        /// Update particle velocity

        u_velocity_particle = u_velocity_particle + delta_t*( u_velocity_fluid_particle - u_velocity_particle )/relaxation_time_particle;

        v_velocity_particle = v_velocity_particle + delta_t*( v_velocity_fluid_particle - v_velocity_particle )/relaxation_time_particle;

        w_velocity_particle = w_velocity_particle + delta_t*( w_velocity_fluid_particle - w_velocity_particle )/relaxation_time_particle;

src/FlowSolverRHEA.cpp

+46 −12
Original line number Diff line number Diff line
@@ -2033,7 +2033,10 @@ double FlowSolverRHEA::trilinearInterpolation(const double &x, const double &y, 


};



void FlowSolverRHEA::updateLagrangianEulerianMeshIndexes0() {

void FlowSolverRHEA::updateLagrangianEulerianMeshIndexes0(const int &my_rank) {



    // Get number particles in use

    number_particles_local_in_use = point_particles->get_num_prts_local_in_use( my_rank );



    int i_local_index, j_local_index, k_local_index;

    double x_position_particle, y_position_particle, z_position_particle;
@@ -2140,17 +2143,31 @@ void FlowSolverRHEA::timeAdvancePointParticles(const int &my_rank) { 
    }

   

    /// Update Lagrangian-Eulerian Mesh indexes 0

    this->updateLagrangianEulerianMeshIndexes0();

    this->updateLagrangianEulerianMeshIndexes0( my_rank );

 

    if( activate_pure_tracer_particles ) {



        double u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid;

        int i_local_index, j_local_index, k_local_index;

        double x_position_particle, y_position_particle, z_position_particle;

        double u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle;

	int capacity = point_particles->get_prt_capacity();

        #pragma acc parallel loop collapse (1) private(x_position_particle, y_position_particle, z_position_particle, u_velocity_particle, v_velocity_particle, w_velocity_particle, dynamic_viscosity_fluid, u_velocity_fluid_particle,w_velocity_fluid_particle,v_velocity_fluid_particle, i_local_index, j_local_index, k_local_index) present(this, mesh, point_particles, u_field.vector[0:_ls_], v_field.vector[0:_ls_], w_field.vector[0:_ls_], mu_field.vector[0:_ls_], point_particles->local_prts_positions_x[0:capacity], point_particles->local_prts_positions_y[0:capacity], point_particles->local_prts_positions_z[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]) 

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



            /// Obtain Lagrangian-Euler values

	    this->obtainLagrangianEulerianVelocityDynamicViscosityValues( u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, p );

            /// Obtain Lagrangian-Eulerian indexes 0

            i_local_index = point_particles->local_prts_indexes_0_i[p];

            j_local_index = point_particles->local_prts_indexes_0_j[p];

            k_local_index = point_particles->local_prts_indexes_0_k[p];



            /// 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];

            z_position_particle = point_particles->local_prts_positions_0_z[p];

    

           /// Interpolate (trilinear) values: u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle

           u_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

           v_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

           w_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );



            /// Update velocity

	    point_particles->local_prts_velocities_x[p] = u_velocity_fluid_particle;
@@ -2175,6 +2192,8 @@ void FlowSolverRHEA::timeAdvanceVelocityPointParticles() { 
    /// IMPORTANT: This method needs to be modified/overwritten according to the problem under consideration



    /// Explicit Euler time-integration of particles velocity

    int i_local_index, j_local_index, k_local_index;

    double x_position_particle, y_position_particle, z_position_particle;

    double u_velocity_particle, v_velocity_particle, w_velocity_particle;

    double u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, relaxation_time_particle;

    delta_t = this->delta_t;
@@ -2182,16 +2201,30 @@ void FlowSolverRHEA::timeAdvanceVelocityPointParticles() { 
    #pragma acc parallel loop collapse (1) private(x_position_particle, y_position_particle, z_position_particle, u_velocity_particle, v_velocity_particle, w_velocity_particle, relaxation_time_particle, dynamic_viscosity_fluid, u_velocity_fluid_particle,w_velocity_fluid_particle,v_velocity_fluid_particle) present(this, mesh, u_field.vector[0:_ls_], v_field.vector[0:_ls_], w_field.vector[0:_ls_], mu_field.vector[0:_ls_], point_particles, point_particles->local_prts_positions_x[0:capacity], point_particles->local_prts_positions_y[0:capacity], point_particles->local_prts_positions_z[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]) copyin(delta_t)

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



        /// Obtain Lagrangian values

        /// Obtain Lagrangian-Eulerian indexes 0

        i_local_index = point_particles->local_prts_indexes_0_i[p];

        j_local_index = point_particles->local_prts_indexes_0_j[p];

        k_local_index = point_particles->local_prts_indexes_0_k[p];



        /// 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];

        z_position_particle = point_particles->local_prts_positions_0_z[p];



        /// Obtain Lagrangian velocities 0

	u_velocity_particle = point_particles->local_prts_velocities_0_x[p];

        v_velocity_particle = point_particles->local_prts_velocities_0_y[p];

        w_velocity_particle = point_particles->local_prts_velocities_0_z[p];



        /// Obtain Lagrangian-Euler values

        this->obtainLagrangianEulerianVelocityDynamicViscosityValues( u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, p );

        /// Interpolate (trilinear) values: u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid 

        u_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

        v_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

        w_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );



        /// Calculate relaxation time particle

        relaxation_time_particle = point_particles->calculate_relaxation_time_prt( p, dynamic_viscosity_fluid );



        /// Update velocity

        /// Update particle velocity

        u_velocity_particle = u_velocity_particle + delta_t*( u_velocity_fluid_particle - u_velocity_particle )/relaxation_time_particle;

        v_velocity_particle = v_velocity_particle + delta_t*( v_velocity_fluid_particle - v_velocity_particle )/relaxation_time_particle;

        w_velocity_particle = w_velocity_particle + delta_t*( w_velocity_fluid_particle - w_velocity_particle )/relaxation_time_particle;
@@ -2200,6 +2233,7 @@ void FlowSolverRHEA::timeAdvanceVelocityPointParticles() { 
        point_particles->local_prts_velocities_z[p] = w_velocity_particle;



    }



};



void FlowSolverRHEA::updatePreviousStateConservedVariables() {
@@ -3574,14 +3608,14 @@ void FlowSolverRHEA::execute() { 
    if( use_restart_particles ) {



        point_particles->read_from_file( restart_data_file_particles );

        this->updateLagrangianEulerianMeshIndexes0();

        this->updateLagrangianEulerianMeshIndexes0( my_rank );



    } else {



        point_particles->generate_prts_random( buffer_ratio_particles );

	    this->updateLagrangianEulerianMeshIndexes0();

	this->updateLagrangianEulerianMeshIndexes0( my_rank );

        this->setInitialParticlesPositionsVelocities();

	    this->updateLagrangianEulerianMeshIndexes0();	

	this->updateLagrangianEulerianMeshIndexes0( my_rank );	



    }

    point_particles->copyToDeviceParticles();

src/FlowSolverRHEA.hpp

+1 −1
Original line number Diff line number Diff line
@@ -213,7 +213,7 @@ class FlowSolverRHEA { 
        virtual void setInitialParticlesPositionsVelocities();



	/// Update Lagrangian-Eulerian Mesh indexes 0

	void updateLagrangianEulerianMeshIndexes0();

	void updateLagrangianEulerianMeshIndexes0(const int &my_rank);



	/// Obtain Lagrangian-Eulerian velocity and dynamic viscosity values

        //virtual void obtainLagrangianEulerianVelocityDynamicViscosityValues(double &u_velocity_fluid_point_particle, double &v_velocity_fluid_point_particle, double &w_velocity_fluid_point_particle, double &dynamic_viscosity_fluid, const int &p);

tests/1d_high_pressure_framework/myRHEA.cpp

+21 −4
Original line number Diff line number Diff line
@@ -159,6 +159,8 @@ void myRHEA::timeAdvanceVelocityPointParticles() { 
    /// IMPORTANT: This method needs to be modified/overwritten according to the problem under consideration



    /// Explicit Euler time-integration of particles velocity

    int i_local_index, j_local_index, k_local_index;

    double x_position_particle, y_position_particle, z_position_particle;

    double u_velocity_particle, v_velocity_particle, w_velocity_particle;

    double u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, relaxation_time_particle;

    delta_t = this->delta_t;
@@ -166,16 +168,30 @@ void myRHEA::timeAdvanceVelocityPointParticles() { 
    #pragma acc parallel loop collapse (1) private(x_position_particle, y_position_particle, z_position_particle, u_velocity_particle, v_velocity_particle, w_velocity_particle, relaxation_time_particle, dynamic_viscosity_fluid, u_velocity_fluid_particle,w_velocity_fluid_particle,v_velocity_fluid_particle) present(this, mesh, u_field.vector[0:_ls_], v_field.vector[0:_ls_], w_field.vector[0:_ls_], mu_field.vector[0:_ls_], point_particles, point_particles->local_prts_positions_x[0:capacity], point_particles->local_prts_positions_y[0:capacity], point_particles->local_prts_positions_z[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]) copyin(delta_t)

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



        /// Obtain Lagrangian values

        /// Obtain Lagrangian-Eulerian indexes 0

        i_local_index = point_particles->local_prts_indexes_0_i[p];

        j_local_index = point_particles->local_prts_indexes_0_j[p];

        k_local_index = point_particles->local_prts_indexes_0_k[p];



        /// 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];

        z_position_particle = point_particles->local_prts_positions_0_z[p];



        /// Obtain Lagrangian velocities 0

	u_velocity_particle = point_particles->local_prts_velocities_0_x[p];

        v_velocity_particle = point_particles->local_prts_velocities_0_y[p];

        w_velocity_particle = point_particles->local_prts_velocities_0_z[p];



        /// Obtain Lagrangian-Euler values

        this->obtainLagrangianEulerianVelocityDynamicViscosityValues( u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid, p );

        /// Interpolate (trilinear) values: u_velocity_fluid_particle, v_velocity_fluid_particle, w_velocity_fluid_particle, dynamic_viscosity_fluid 

        u_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], u_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], u_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], u_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

        v_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], v_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], v_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], v_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );

        w_velocity_fluid_particle = this->trilinearInterpolation( x_position_particle, y_position_particle, z_position_particle, mesh->x[i_local_index-1], mesh->x[i_local_index+1], mesh->y[j_local_index-1], mesh->y[j_local_index+1], mesh->z[k_local_index-1], mesh->z[k_local_index+1], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index-1)], w_field[I1D(i_local_index-1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index-1,k_local_index+1)], w_field[I1D(i_local_index-1,j_local_index+1,k_local_index+1)], w_field[I1D(i_local_index+1,j_local_index+1,k_local_index+1)] );



        /// Calculate relaxation time particle

        relaxation_time_particle = point_particles->calculate_relaxation_time_prt( p, dynamic_viscosity_fluid );



        /// Update velocity

        /// Update particle velocity

        u_velocity_particle = u_velocity_particle + delta_t*( u_velocity_fluid_particle - u_velocity_particle )/relaxation_time_particle;

        v_velocity_particle = v_velocity_particle + delta_t*( v_velocity_fluid_particle - v_velocity_particle )/relaxation_time_particle;

        w_velocity_particle = w_velocity_particle + delta_t*( w_velocity_fluid_particle - w_velocity_particle )/relaxation_time_particle;
@@ -187,6 +203,7 @@ void myRHEA::timeAdvanceVelocityPointParticles() { 


};





////////// MAIN //////////

int main(int argc, char** argv) {

tests/1d_sod_shock_tube/myRHEA.cpp

+21 −4

File changed.

Preview size limit exceeded, changes collapsed.