Commit f0565bb8 authored by ProjectRHEA's avatar ProjectRHEA
Browse files

Solver modified to accept non-uniform Cartesian meshes ... compiles, not tested

parent 973682eb

src/FlowSolverRHEA.cpp

+4 −4
Original line number Diff line number Diff line
@@ -1577,7 +1577,7 @@ void FlowSolverRHEA::updateBoundaries() { 
    for(int i = topo->iter_bound[_NORTH_][_INIX_]; i <= topo->iter_bound[_NORTH_][_ENDX_]; i++) {

        for(int j = topo->iter_bound[_NORTH_][_INIY_]; j <= topo->iter_bound[_NORTH_][_ENDY_]; j++) {

            for(int k = topo->iter_bound[_NORTH_][_INIZ_]; k <= topo->iter_bound[_NORTH_][_ENDZ_]; k++) {

                if( ( bocos_type[_NORTH_] == _DIRICHLET_ ) or ( bocos_type[_NORTH_] == _SUBSONIC_INFLOW_ ) or ( bocos_type[_NORTH_] == _SUPERSONIC_INFLOW_ ) ) 

                if( ( bocos_type[_NORTH_] == _DIRICHLET_ ) or ( bocos_type[_NORTH_] == _SUBSONIC_INFLOW_ ) or ( bocos_type[_NORTH_] == _SUPERSONIC_INFLOW_ ) ) { 

                    wg_g  = 1.0 - ( y_field[I1D(i,j+1,k)] - 0.5*( y_field[I1D(i,j,k)] + y_field[I1D(i,j+1,k)] ) )/( y_field[I1D(i,j+1,k)] - y_field[I1D(i,j,k)] );

                    wg_in = 1.0 - ( 0.5*( y_field[I1D(i,j,k)] + y_field[I1D(i,j+1,k)] ) - y_field[I1D(i,j,k)] )/( y_field[I1D(i,j+1,k)] - y_field[I1D(i,j,k)] );

                }
@@ -2250,7 +2250,7 @@ void FlowSolverRHEA::updateLagrangianEulerianMeshIndexes0(const int &my_rank) { 
 

        /// Search in x-direction

        for( int i = topo->iter_common[_INNER_][_INIX_]; i <= topo->iter_common[_INNER_][_ENDX_]; ++i ) {

            if( x_position_particle >= 0.5*( x_field[I1D(i-1,j,k)] + x_field[I1D(i,j,k)] ) && x_position_particle < 0.5*( x_field[I1D(i,j,k)] + x_field[I1D(i+1,j,k)] ) ) {

            if( x_position_particle >= 0.5*( x_field[I1D(i-1,0,0)] + x_field[I1D(i,0,0)] ) && x_position_particle < 0.5*( x_field[I1D(i,0,0)] + x_field[I1D(i+1,0,0)] ) ) {

                i_local_index = i;

                break;

            }
@@ -2258,7 +2258,7 @@ void FlowSolverRHEA::updateLagrangianEulerianMeshIndexes0(const int &my_rank) { 


        /// Search in y-direction

        for( int j = topo->iter_common[_INNER_][_INIY_]; j <= topo->iter_common[_INNER_][_ENDY_]; ++j ) {

            if( y_position_particle >= 0.5*( y_field[I1D(i,j-1,k)] + y_field[I1D(i,j,k)] ) && y_position_particle < 0.5*( y_field[I1D(i,j,k)] + y_field[I1D(i,j+1,k)] ) ) {

            if( y_position_particle >= 0.5*( y_field[I1D(0,j-1,0)] + y_field[I1D(0,j,0)] ) && y_position_particle < 0.5*( y_field[I1D(0,j,0)] + y_field[I1D(0,j+1,0)] ) ) {

                j_local_index = j;

                break;

            }
@@ -2266,7 +2266,7 @@ void FlowSolverRHEA::updateLagrangianEulerianMeshIndexes0(const int &my_rank) { 


        /// Search in z-direction

        for( int k = topo->iter_common[_INNER_][_INIZ_]; k <= topo->iter_common[_INNER_][_ENDZ_]; ++k ) {

            if( z_position_particle >= 0.5*( z_field[I1D(i,j,k-1)] + z_field[I1D(i,j,k)] ) && z_position_particle < 0.5*( z_field[I1D(i,j,k)] + z_field[I1D(i,j,k+1)] ) ) {

            if( z_position_particle >= 0.5*( z_field[I1D(0,0,k-1)] + z_field[I1D(0,0,k)] ) && z_position_particle < 0.5*( z_field[I1D(0,0,k)] + z_field[I1D(0,0,k+1)] ) ) {

                k_local_index = k;

                break;

            }

src/InputOutputManager.cpp

+10 −6
Original line number Diff line number Diff line
@@ -213,7 +213,7 @@ void WriteReadHDF5::write(int it ) //, double time, bool xdmf_file) 


}



void WriteReadHDF5::write2dDataOutputSlice(const DistributedArray &x_field, const DistributedArray &y_field, const DistributedArray &z_field, ComputationalDomain* mesh, ParallelTopology* topo, char const* outname_dos, char const* dos_normal, double dos_x_position_adjusted, double dos_y_position_adjusted, double dos_z_position_adjusted, bool dos_gen_xdmf, int it) {

void WriteReadHDF5::write2dDataOutputSlice(DistributedArray &x_field, DistributedArray &y_field, DistributedArray &z_field, ComputationalDomain* mesh, ParallelTopology* topo, char const* outname_dos, char const* dos_normal, double dos_x_position_adjusted, double dos_y_position_adjusted, double dos_z_position_adjusted, bool dos_gen_xdmf, int it) {



    char filename[100];
@@ -252,7 +252,7 @@ void WriteReadHDF5::write2dDataOutputSlice(const DistributedArray &x_field, cons 
        dim_gsize_2dPlane[2] = 1;

        bool plane_contained = false;

        for(int i =  myTopo->iter_slab[_INIX_]; i <=  myTopo->iter_slab[_ENDX_]; i++){

            if( abs( x_field[I1D(i,j,k)] - dos_x_position_adjusted ) < ( 1.1*numeric_limits<double>::min() ) ) {

            if( abs( x_field[I1D(i,0,0)] - dos_x_position_adjusted ) < ( 1.1*numeric_limits<double>::min() ) ) {

                plane_contained = true;

            }

        }
@@ -275,7 +275,7 @@ void WriteReadHDF5::write2dDataOutputSlice(const DistributedArray &x_field, cons 
        dim_gsize_2dPlane[2] = (topo->getMesh())->getGNx() + 2;

        bool plane_contained = false;

        for(int j =  myTopo->iter_slab[_INIY_]; j <=  myTopo->iter_slab[_ENDY_]; j++){

            if( abs( y_field[I1D(i,j,k)] - dos_y_position_adjusted ) < ( 1.1*numeric_limits<double>::min() ) ) {

            if( abs( y_field[I1D(0,j,0)] - dos_y_position_adjusted ) < ( 1.1*numeric_limits<double>::min() ) ) {

                plane_contained = true;

            }

        }
@@ -298,7 +298,7 @@ void WriteReadHDF5::write2dDataOutputSlice(const DistributedArray &x_field, cons 
        dim_gsize_2dPlane[2] = (topo->getMesh())->getGNx() + 2;

        bool plane_contained = false;

        for(int k =  myTopo->iter_slab[_INIZ_]; k <=  myTopo->iter_slab[_ENDZ_]; k++){

            if( abs( z_field[I1D(i,j,k)] - dos_z_position_adjusted ) < ( 1.1*numeric_limits<double>::min() ) ) {

            if( abs( z_field[I1D(0,0,k)] - dos_z_position_adjusted ) < ( 1.1*numeric_limits<double>::min() ) ) {

                plane_contained = true;

            }

        }
@@ -553,7 +553,7 @@ TemporalPointProbe::TemporalPointProbe(ComputationalDomain* mesh_, ParallelTopol 


};



TemporalPointProbe::TemporalPointProbe(const double &x_position_, const double &y_position_, const double &z_position_, const DistributedArray &x_field, const DistributedArray &y_field, const DistributedArray &z_field, const string output_file_name_, ComputationalDomain* mesh_, ParallelTopology* topo_) {

TemporalPointProbe::TemporalPointProbe(const double &x_position_, const double &y_position_, const double &z_position_, DistributedArray &x_field, DistributedArray &y_field, DistributedArray &z_field, const string output_file_name_, ComputationalDomain* mesh_, ParallelTopology* topo_) {



    /// Set position, output file name, mesh & topo

    x_position = x_position_;
@@ -563,6 +563,10 @@ TemporalPointProbe::TemporalPointProbe(const double &x_position_, const double & 
    mesh = mesh_;

    topo = topo_;



    _lNx_ = topo->getlNx();

    _lNy_ = topo->getlNy();

    _lNz_ = topo->getlNz();



    /// Locate closest grid point to probe

    this->locateClosestGridPointToProbe(x_field, y_field, z_field);
@@ -576,7 +580,7 @@ TemporalPointProbe::~TemporalPointProbe() { 


};



void TemporalPointProbe::locateClosestGridPointToProbe(const DistributedArray &x_field, const DistributedArray &y_field, const DistributedArray &z_field) {

void TemporalPointProbe::locateClosestGridPointToProbe(DistributedArray &x_field, DistributedArray &y_field, DistributedArray &z_field) {



    /// Initialize MPI

    int my_rank, world_size;

src/InputOutputManager.hpp

+8 −3
Original line number Diff line number Diff line
@@ -22,7 +22,7 @@ class WriteReadHDF5 { 
        void addField(DistributedArray*);

        void printOnScreen();

        void write(int);

        void write2dDataOutputSlice(const DistributedArray &x_field, const DistributedArray &y_field, const DistributedArray &z_field, ComputationalDomain*, ParallelTopology*, char const*, char const*, double, double, double, bool, int);

        void write2dDataOutputSlice(DistributedArray &x_field, DistributedArray &y_field, DistributedArray &z_field, ComputationalDomain*, ParallelTopology*, char const*, char const*, double, double, double, bool, int);

        void read(char const*);



        void addAttributeDouble(std::string str){ double dval=0.0; dattrib[str]=dval;};
@@ -79,7 +79,7 @@ class TemporalPointProbe { 
        ////////// CONSTRUCTORS & DESTRUCTOR //////////

        TemporalPointProbe();								/// Default constructor

        TemporalPointProbe(ComputationalDomain* mesh_, ParallelTopology* topo_);	/// Parametrized constructor

        TemporalPointProbe(const double &x_position_, const double &y_position_, const double &z_position_, const DistributedArray &x_field, const DistributedArray &y_field, const DistributedArray &z_field, const std::string output_file_name_, ComputationalDomain* mesh_, ParallelTopology* topo_);	/// Parametrized constructor

        TemporalPointProbe(const double &x_position_, const double &y_position_, const double &z_position_, DistributedArray &x_field, DistributedArray &y_field, DistributedArray &z_field, const std::string output_file_name_, ComputationalDomain* mesh_, ParallelTopology* topo_);	/// Parametrized constructor

        virtual ~TemporalPointProbe();							/// Destructor



	////////// GET FUNCTIONS //////////
@@ -105,7 +105,7 @@ class TemporalPointProbe { 
	////////// PROBE METHODS //////////

	

	/// Locate closest grid point to probe

        virtual void locateClosestGridPointToProbe(const DistributedArray &x_field, const DistributedArray &y_field, const DistributedArray &z_field);

        virtual void locateClosestGridPointToProbe(DistributedArray &x_field, DistributedArray &y_field, DistributedArray &z_field);

        

	/// Write data string to output file

	virtual void writeDataStringToOutputFile(const std::string output_header_string, const std::string output_data_string);
@@ -132,6 +132,11 @@ class TemporalPointProbe { 


    private:



        /// Local mesh values for I1D macro

        int _lNx_;

        int _lNy_;

        int _lNz_;



};



#endif /*_INPUT_OUTPUT_MANAGER_*/

tests/3d_high_pressure_jet/myRHEA.cpp

+1 −1
Original line number Diff line number Diff line
@@ -820,7 +820,7 @@ void myRHEA::updateBoundaries() { 
    for(int i = topo->iter_bound[_NORTH_][_INIX_]; i <= topo->iter_bound[_NORTH_][_ENDX_]; i++) {

        for(int j = topo->iter_bound[_NORTH_][_INIY_]; j <= topo->iter_bound[_NORTH_][_ENDY_]; j++) {

            for(int k = topo->iter_bound[_NORTH_][_INIZ_]; k <= topo->iter_bound[_NORTH_][_ENDZ_]; k++) {

                if( ( bocos_type[_NORTH_] == _DIRICHLET_ ) or ( bocos_type[_NORTH_] == _SUBSONIC_INFLOW_ ) or ( bocos_type[_NORTH_] == _SUPERSONIC_INFLOW_ ) ) 

                if( ( bocos_type[_NORTH_] == _DIRICHLET_ ) or ( bocos_type[_NORTH_] == _SUBSONIC_INFLOW_ ) or ( bocos_type[_NORTH_] == _SUPERSONIC_INFLOW_ ) ) { 

                    wg_g  = 1.0 - ( y_field[I1D(i,j+1,k)] - 0.5*( y_field[I1D(i,j,k)] + y_field[I1D(i,j+1,k)] ) )/( y_field[I1D(i,j+1,k)] - y_field[I1D(i,j,k)] );

                    wg_in = 1.0 - ( 0.5*( y_field[I1D(i,j,k)] + y_field[I1D(i,j+1,k)] ) - y_field[I1D(i,j,k)] )/( y_field[I1D(i,j+1,k)] - y_field[I1D(i,j,k)] );

                }