diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/CMakeLists.txt b/tests/extra_source_and_tests/test_2d_modified_T_flow/CMakeLists.txt new file mode 100644 index 0000000000..f8ab3ce965 --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/CMakeLists.txt @@ -0,0 +1,26 @@ +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${SPHINXSYS_PROJECT_DIR}/cmake) # main (top) cmake dir + +set(CMAKE_VERBOSE_MAKEFILE on) + +STRING( REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR} ) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) + +aux_source_directory(. DIR_SRCS) +ADD_EXECUTABLE(${PROJECT_NAME} ${EXECUTABLE_OUTPUT_PATH} ${DIR_SRCS}) + +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + +set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "spatial_temporal_identification") + +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") +target_link_libraries(${PROJECT_NAME} extra_sources_2d) diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/modified_T_flow.cpp b/tests/extra_source_and_tests/test_2d_modified_T_flow/modified_T_flow.cpp new file mode 100644 index 0000000000..4bbd88543c --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/modified_T_flow.cpp @@ -0,0 +1,374 @@ +/* +* @file modified_T_shaped_pipe.cpp +* @brief This is the benchmark test of multi -inlet and multi - outlet. +* @details We consider a flow with one inlet and two outlets in a T - shaped pipe in 2D. +* @author Xiangyu Hu,Shuoguo Zhang +*/ + +#include "bidirectional_buffer.h" +#include "density_correciton.h" +#include "density_correciton.hpp" +#include "kernel_summation.h" +#include "kernel_summation.hpp" +#include "pressure_boundary.h" +#include "sphinxsys.h" + +using namespace SPH; +//---------------------------------------------------------------------- +// Basic geometry parameters and numerical setup. +//---------------------------------------------------------------------- +Real DL = 0.2; /**< Reference length. */ +Real DH = 0.1; /**< Reference and the height of main channel. */ +Real DL1 = 0.75 * DL; /**< The length of the main channel. */ +Real resolution_ref = 0.005; /**< Initial reference particle spacing. */ +Real BW = resolution_ref * 4; /**< Reference size of the emitter. */ +Real DL_sponge = resolution_ref * 20; /**< Reference size of the emitter buffer to impose inflow condition. */ +StdVec observer_location = {Vecd(0.5 * DL, 0.5 * DH)}; /**< Displacement observation point. */ + +//---------------------------------------------------------------------- +// Global parameters on the fluid properties +//---------------------------------------------------------------------- +Real Outlet_pressure = 0; +Real rho0_f = 1000.0; /**< Reference density of fluid. */ +Real Re = 100.0; /**< Reynolds number. */ +Real U_f = 1.0; /**< Characteristic velocity. */ +Real mu_f = rho0_f * U_f * DH / Re; /**< Dynamics viscosity. */ +Real c_f = 10.0 * U_f * SMAX(Real(1), DH / (Real(2.0) * (DL - DL1)));/** Reference sound speed needs to consider the flow speed in the narrow channels. */ +Vec2d normal = Vec2d(1.0, 0.0); +//---------------------------------------------------------------------- +// Pressure boundary definition. +//---------------------------------------------------------------------- +struct LeftInflowPressure +{ + template + LeftInflowPressure(BoundaryConditionType &boundary_condition) {} + + Real operator()(Real &p_) + { + return p_; + } +}; + +struct UpOutflowPressure +{ + template + UpOutflowPressure(BoundaryConditionType &boundary_condition) {} + + Real operator()(Real &p_) + { + /*constant pressure*/ + Real pressure = Outlet_pressure; + return pressure; + } +}; + +struct DownOutflowPressure +{ + template + DownOutflowPressure(BoundaryConditionType &boundary_condition) {} + + Real operator()(Real &p_) + { + /*constant pressure*/ + Real pressure = Outlet_pressure; + return pressure; + } +}; +//---------------------------------------------------------------------- +// inflow velocity definition. +//---------------------------------------------------------------------- +struct InflowVelocity +{ + Real u_ref_, t_ref_; + AlignedBoxShape &aligned_box_; + Vecd halfsize_; + + template + InflowVelocity(BoundaryConditionType &boundary_condition) + : u_ref_(U_f), t_ref_(2.0), + aligned_box_(boundary_condition.getAlignedBox()), + halfsize_(aligned_box_.HalfSize()) {} + + Vecd operator()(Vecd &position, Vecd &velocity) + { + Vecd target_velocity = velocity; + Real run_time = GlobalStaticVariables::physical_time_; + Real u_ave = run_time < t_ref_ ? 0.5 * u_ref_ * (1.0 - cos(Pi * run_time / t_ref_)) : u_ref_; + target_velocity[0] = 1.5 * u_ave * SMAX(0.0, 1.0 - position[1] * position[1] / halfsize_[1] / halfsize_[1]); + return target_velocity; + } +}; + +//---------------------------------------------------------------------- +// define geometry of SPH bodies +//---------------------------------------------------------------------- +/** the water block in T shape polygon. */ +std::vector water_block_shape{ + Vecd(-DL_sponge, 0.0), Vecd(-DL_sponge, DH), Vecd(DL1, DH), Vecd(DL1, 2.0 * DH), + Vecd(DL, 2.0 * DH), Vecd(DL, -DH), Vecd(DL1, -DH), Vecd(DL1, 0.0), Vecd(-DL_sponge, 0.0)}; +/** the outer wall polygon. */ +std::vector outer_wall_shape{ + Vecd(-DL_sponge, -BW), Vecd(-DL_sponge, DH + BW), Vecd(DL1 - BW, DH + BW), Vecd(DL1 - BW, 2.0 * DH), + Vecd(DL + BW, 2.0 * DH), Vecd(DL + BW, -DH), Vecd(DL1 - BW, -DH), Vecd(DL1 - BW, -BW), Vecd(-DL_sponge, -BW)}; +/** the inner wall polygon. */ +std::vector inner_wall_shape{ + Vecd(-DL_sponge - BW, 0.0), Vecd(-DL_sponge - BW, DH), Vecd(DL1, DH), Vecd(DL1, 2.0 * DH + BW), + Vecd(DL, 2.0 * DH + BW), Vecd(DL, -DH - BW), Vecd(DL1, -DH - BW), Vecd(DL1, 0.0), Vecd(-DL_sponge - BW, 0.0)}; +//---------------------------------------------------------------------- +// Define case dependent body shapes. +//---------------------------------------------------------------------- +class WaterBlock : public MultiPolygonShape +{ + public: + explicit WaterBlock(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(water_block_shape, ShapeBooleanOps::add); + } +}; + +class WallBoundary : public MultiPolygonShape +{ + public: + explicit WallBoundary(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(outer_wall_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(inner_wall_shape, ShapeBooleanOps::sub); + } +}; + +//---------------------------------------------------------------------- +// Main program starts here. +//---------------------------------------------------------------------- +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up an SPHSystem and IO environment. + //---------------------------------------------------------------------- + BoundingBox system_domain_bounds(Vec2d(-DL_sponge - BW, -DH - BW), Vec2d(DL + BW, 2.0 * DH + BW)); + SPHSystem sph_system(system_domain_bounds, resolution_ref); + sph_system.setGenerateRegressionData(false); + sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment(); + //---------------------------------------------------------------------- + // Creating bodies with corresponding materials and particles. + //---------------------------------------------------------------------- + FluidBody water_block(sph_system, makeShared("WaterBody")); + water_block.defineMaterial(rho0_f, c_f, mu_f); + ParticleBuffer in_outlet_particle_buffer(0.5); + water_block.generateParticlesWithReserve(in_outlet_particle_buffer); + + SolidBody wall_boundary(sph_system, makeShared("WallBoundary")); + wall_boundary.defineMaterial(); + wall_boundary.generateParticles(); + + ObserverBody velocity_observer(sph_system, "VelocityObserver"); + velocity_observer.generateParticles(observer_location); + //---------------------------------------------------------------------- + // Define body relation map. + // The contact map gives the topological connections between the bodies. + // Basically the the range of bodies to build neighbor particle lists. + // Generally, we first define all the inner relations, then the contact relations. + //---------------------------------------------------------------------- + InnerRelation water_block_inner(water_block); + ContactRelation water_block_contact(water_block, {&wall_boundary}); + ContactRelation velocity_observer_contact(velocity_observer, {&water_block}); + //---------------------------------------------------------------------- + // Combined relations built from basic relations + // which is only used for update configuration. + //---------------------------------------------------------------------- + ComplexRelation water_block_complex(water_block_inner, water_block_contact); + //---------------------------------------------------------------------- + // Define the numerical methods used in the simulation. + // Note that there may be data dependence on the sequence of constructions. + // Generally, the geometric models or simple objects without data dependencies, + // such as gravity, should be initiated first. + // Then the major physical particle dynamics model should be introduced. + // Finally, the auxillary models such as time step estimator, initial condition, + // boundary condition and other constraints should be defined. + //---------------------------------------------------------------------- + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + InteractionDynamics kernel_summation(water_block_inner, water_block_contact); + InteractionWithUpdate boundary_indicator(water_block_inner, water_block_contact); + + Dynamics1Level pressure_relaxation(water_block_inner, water_block_contact); + Dynamics1Level density_relaxation(water_block_inner, water_block_contact); + InteractionWithUpdate viscous_acceleration(water_block_inner, water_block_contact); + InteractionWithUpdate> transport_velocity_correction(water_block_inner, water_block_contact); + ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); + ReduceDynamics get_fluid_time_step_size(water_block); + + //---------------------------------------------------------------------- + // Left buffer + //---------------------------------------------------------------------- + Vec2d left_buffer_halfsize = Vec2d(0.5 * BW, 0.5 * DH); + Vec2d left_buffer_translation = Vec2d(-DL_sponge, 0.0) + left_buffer_halfsize; + BodyAlignedBoxByCell left_emitter(water_block, makeShared(Transform(Vec2d(left_buffer_translation)), left_buffer_halfsize)); + fluid_dynamics::NonPrescribedPressureBidirectionalBuffer left_emitter_inflow_injection(left_emitter, in_outlet_particle_buffer, xAxis); + + BodyAlignedBoxByCell left_disposer(water_block, makeShared(Transform(Rotation2d(Pi), Vec2d(left_buffer_translation)), left_buffer_halfsize)); + SimpleDynamics left_disposer_outflow_deletion(left_disposer, xAxis); + //---------------------------------------------------------------------- + // Up buffer + //---------------------------------------------------------------------- + Vec2d up_buffer_halfsize = Vec2d(0.5 * BW, 0.75); + Vec2d up_buffer_translation = Vec2d(0.5 * (DL + DL1), 2.0 * DH - 0.5 * BW); + BodyAlignedBoxByCell up_emitter(water_block, makeShared(Transform(Rotation2d(-0.5 * Pi), Vec2d(up_buffer_translation)), up_buffer_halfsize)); + fluid_dynamics::BidirectionalBuffer up_emitter_inflow_injection(up_emitter, in_outlet_particle_buffer, xAxis); + + BodyAlignedBoxByCell up_disposer(water_block, makeShared(Transform(Rotation2d(0.5 * Pi), Vec2d(up_buffer_translation)), up_buffer_halfsize)); + SimpleDynamics up_disposer_outflow_deletion(up_disposer, xAxis); + //---------------------------------------------------------------------- + // Down buffer + //---------------------------------------------------------------------- + Vec2d down_buffer_halfsize = Vec2d(0.5 * BW, 0.75); + Vec2d down_buffer_translation = Vec2d(0.5 * (DL + DL1), -DH + 0.5 * BW); + BodyAlignedBoxByCell down_emitter(water_block, makeShared(Transform(Rotation2d(0.5 * Pi), Vec2d(down_buffer_translation)), down_buffer_halfsize)); + fluid_dynamics::BidirectionalBuffer down_emitter_inflow_injection(down_emitter, in_outlet_particle_buffer, xAxis); + + BodyAlignedBoxByCell down_disposer(water_block, makeShared(Transform(Rotation2d(-0.5 * Pi), Vec2d(down_buffer_translation)), down_buffer_halfsize)); + SimpleDynamics down_disposer_outflow_deletion(down_disposer, xAxis); + + InteractionWithUpdate update_fluid_density(water_block_inner, water_block_contact); + SimpleDynamics> left_inflow_pressure_condition(left_emitter); + SimpleDynamics> up_inflow_pressure_condition(up_emitter); + SimpleDynamics> down_inflow_pressure_condition(down_emitter); + SimpleDynamics> inflow_velocity_condition(left_emitter); + //---------------------------------------------------------------------- + // Define the methods for I/O operations, observations + // and regression tests of the simulation. + //---------------------------------------------------------------------- + BodyStatesRecordingToVtp body_states_recording(sph_system); + body_states_recording.addVariableRecording(water_block, "Pressure"); + body_states_recording.addVariableRecording(water_block, "Indicator"); + body_states_recording.addVariableRecording(water_block, "Density"); + body_states_recording.addVariableRecording(water_block, "BufferParticleIndicator"); + + RegressionTestDynamicTimeWarping> write_centerline_velocity("Velocity", velocity_observer_contact); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + boundary_indicator.exec(); + left_emitter_inflow_injection.tag_buffer_particles.exec(); + up_emitter_inflow_injection.tag_buffer_particles.exec(); + down_emitter_inflow_injection.tag_buffer_particles.exec(); + wall_boundary_normal_direction.exec(); + //---------------------------------------------------------------------- + // Setup for time-stepping control + //---------------------------------------------------------------------- + size_t number_of_iterations = sph_system.RestartStep(); + int screen_output_interval = 100; + int observation_sample_interval = screen_output_interval * 2; + Real end_time = 30.0; /**< End time. */ + Real Output_Time = end_time/300.0; /**< Time stamps for output of body states. */ + Real dt = 0.0; /**< Default acoustic time step sizes. */ + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); + TimeInterval interval; + + TimeInterval interval_computing_time_step; + TimeInterval interval_computing_pressure_relaxation; + TimeInterval interval_updating_configuration; + TickCount time_instance; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + body_states_recording.writeToFile(); + write_centerline_velocity.writeToFile(number_of_iterations); + //---------------------------------------------------------------------- + // Main loop starts here. + //---------------------------------------------------------------------- + while (GlobalStaticVariables::physical_time_ < end_time) + { + Real integration_time = 0.0; + /** Integrate time (loop) until the next output time. */ + while (integration_time < Output_Time) + { + time_instance = TickCount::now(); + Real Dt = get_fluid_advection_time_step_size.exec(); + update_fluid_density.exec(); + viscous_acceleration.exec(); + transport_velocity_correction.exec(); + interval_computing_time_step += TickCount::now() - time_instance; + + time_instance = TickCount::now(); + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + dt = SMIN(get_fluid_time_step_size.exec(), Dt); + pressure_relaxation.exec(dt); + kernel_summation.exec(); + left_inflow_pressure_condition.exec(dt); + up_inflow_pressure_condition.exec(dt); + down_inflow_pressure_condition.exec(dt); + inflow_velocity_condition.exec(); + density_relaxation.exec(dt); + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + } + interval_computing_pressure_relaxation += TickCount::now() - time_instance; + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ + << " Dt = " << Dt << " dt = " << dt << "\n"; + + if (number_of_iterations % observation_sample_interval == 0 && number_of_iterations != sph_system.RestartStep()) + { + write_centerline_velocity.writeToFile(number_of_iterations); + } + } + number_of_iterations++; + + time_instance = TickCount::now(); + + left_emitter_inflow_injection.injection.exec(); + up_emitter_inflow_injection.injection.exec(); + down_emitter_inflow_injection.injection.exec(); + left_disposer_outflow_deletion.exec(); + up_disposer_outflow_deletion.exec(); + down_disposer_outflow_deletion.exec(); + + water_block.updateCellLinkedListWithParticleSort(100); + water_block_complex.updateConfiguration(); + + interval_updating_configuration += TickCount::now() - time_instance; + boundary_indicator.exec(); + left_emitter_inflow_injection.tag_buffer_particles.exec(); + up_emitter_inflow_injection.tag_buffer_particles.exec(); + down_emitter_inflow_injection.tag_buffer_particles.exec(); + } + TickCount t2 = TickCount::now(); + body_states_recording.writeToFile(); + velocity_observer_contact.updateConfiguration(); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall time for computation: " << tt.seconds() + << " seconds." << std::endl; + std::cout << std::fixed << std::setprecision(9) << "interval_computing_time_step =" + << interval_computing_time_step.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_computing_pressure_relaxation = " + << interval_computing_pressure_relaxation.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_updating_configuration = " + << interval_updating_configuration.seconds() << "\n"; + + if (sph_system.GenerateRegressionData()) + { + write_centerline_velocity.generateDataBase(1.0e-3); + } + else + { + write_centerline_velocity.testResult(); + } + + return 0; +} diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_11_result.xml b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_11_result.xml new file mode 100644 index 0000000000..40ef3b4042 --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_11_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_12_result.xml b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_12_result.xml new file mode 100644 index 0000000000..99007eb935 --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_12_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_13_result.xml b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_13_result.xml new file mode 100644 index 0000000000..c1d2a48f19 --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_Run_13_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_dtwdistance.xml b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_dtwdistance.xml new file mode 100644 index 0000000000..0bb659e3ff --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_runtimes.dat b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_runtimes.dat new file mode 100644 index 0000000000..b3c38ec1f8 --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/VelocityObserver_Velocity_runtimes.dat @@ -0,0 +1,3 @@ +false +16 +0 \ No newline at end of file diff --git a/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/regression_test_tool.py b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..f866842722 --- /dev/null +++ b/tests/extra_source_and_tests/test_2d_modified_T_flow/regression_test_tool/regression_test_tool.py @@ -0,0 +1,38 @@ +# !/usr/bin/env python3 +import os +import sys + +path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') +sys.path.append(path) +from regression_test_base_tool import SphinxsysRegressionTest + +""" +case name: test_2d_modified_T_flow +""" + +case_name = "test_2d_modified_T_flow" +body_name = "VelocityObserver" +parameter_name = "Velocity" + + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) + + +while True: + print("Now start a new run......") + sphinxsys.run_case() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if (converged == "true"): + print("The tested parameters of all variables are converged, and the run will stop here!") + break + elif converged != "true": + print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break