From 1d65f6d4af82c5b6e09941f6ed87a590e7b5adf2 Mon Sep 17 00:00:00 2001 From: yaru0818 Date: Wed, 30 Aug 2023 16:44:45 +0200 Subject: [PATCH] Add rotation Patch and impact patch cases --- .../general_dynamics/general_interaction.cpp | 5 +- .../general_dynamics/general_interaction.h | 1 + .../test_2d_impact_patch/CMakeLists.txt | 25 ++ .../test_2d_impact_patch/impact_patch.cpp | 270 +++++++++++++++++ .../regression_test_tool.py | 44 +++ .../test_2d_ration_patch/CMakeLists.txt | 25 ++ .../test_2d_ration_patch/ratotion_patch.cpp | 283 ++++++++++++++++++ .../regression_test_tool.py | 44 +++ 8 files changed, 696 insertions(+), 1 deletion(-) create mode 100644 tests/user_examples/test_2d_impact_patch/CMakeLists.txt create mode 100644 tests/user_examples/test_2d_impact_patch/impact_patch.cpp create mode 100644 tests/user_examples/test_2d_impact_patch/regression_test_tool/regression_test_tool.py create mode 100644 tests/user_examples/test_2d_ration_patch/CMakeLists.txt create mode 100644 tests/user_examples/test_2d_ration_patch/ratotion_patch.cpp create mode 100644 tests/user_examples/test_2d_ration_patch/regression_test_tool/regression_test_tool.py diff --git a/src/shared/particle_dynamics/general_dynamics/general_interaction.cpp b/src/shared/particle_dynamics/general_dynamics/general_interaction.cpp index 733009d8c8..f4670ca9c3 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_interaction.cpp +++ b/src/shared/particle_dynamics/general_dynamics/general_interaction.cpp @@ -9,6 +9,7 @@ CorrectedConfigurationInner:: : LocalDynamics(inner_relation.getSPHBody()), GeneralDataDelegateInner(inner_relation), beta_(beta), alpha_(alpha), + smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), B_(*particles_->getVariableByName("CorrectionMatrix")) {} //=================================================================================================// void CorrectedConfigurationInner::interaction(size_t index_i, Real dt) @@ -30,7 +31,9 @@ void CorrectedConfigurationInner::update(size_t index_i, Real dt) { Real det_sqr = pow(B_[index_i].determinant(), beta_); Matd inverse = B_[index_i].inverse(); - B_[index_i] = (det_sqr * inverse + alpha_ * Matd::Identity()) / (alpha_ + det_sqr); + Real weight1_ = det_sqr / (alpha_ * pow(smoothing_length_, beta_) + det_sqr); + Real weight2_ = alpha_ * pow(smoothing_length_, beta_) / (alpha_ * pow(smoothing_length_, beta_) + det_sqr); + B_[index_i] = weight1_ * inverse + weight2_ * Matd::Identity(); } //=================================================================================================// CorrectedConfigurationComplex:: diff --git a/src/shared/particle_dynamics/general_dynamics/general_interaction.h b/src/shared/particle_dynamics/general_dynamics/general_interaction.h index f7f34beec0..e8f423ecdf 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_interaction.h +++ b/src/shared/particle_dynamics/general_dynamics/general_interaction.h @@ -42,6 +42,7 @@ class CorrectedConfigurationInner : public LocalDynamics, public GeneralDataDele protected: int beta_; Real alpha_; + Real smoothing_length_; StdLargeVec &B_; void interaction(size_t index_i, Real dt = 0.0); diff --git a/tests/user_examples/test_2d_impact_patch/CMakeLists.txt b/tests/user_examples/test_2d_impact_patch/CMakeLists.txt new file mode 100644 index 0000000000..1e1bc16be3 --- /dev/null +++ b/tests/user_examples/test_2d_impact_patch/CMakeLists.txt @@ -0,0 +1,25 @@ +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_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}) +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} ) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") +target_link_libraries(${PROJECT_NAME} extra_sources_2d) + +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +add_test(NAME ${PROJECT_NAME}_restart COMMAND ${PROJECT_NAME} --restart_step=4000 + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) diff --git a/tests/user_examples/test_2d_impact_patch/impact_patch.cpp b/tests/user_examples/test_2d_impact_patch/impact_patch.cpp new file mode 100644 index 0000000000..dca807e5f6 --- /dev/null +++ b/tests/user_examples/test_2d_impact_patch/impact_patch.cpp @@ -0,0 +1,270 @@ +/** + * @file impactpatch.cpp + * @brief 2D impactpatch example. + * @author Yaru Ren, Chi Zhang and Xiangyu Hu + */ +#include "sphinxsys.h" //SPHinXsys Library. +using namespace SPH; // Namespace cite here. +#define PI 3.1415926 +//---------------------------------------------------------------------- +// Basic geometry parameters and numerical setup. +//---------------------------------------------------------------------- +Real LL = 0.667; /**< Liquid column length. */ +Real LH = 2.0; /**< Liquid column height. */ +Real particle_spacing_ref = LL / 150; /**< Initial reference particle spacing. */ +Real BW = particle_spacing_ref * 4; /**< Extending width for boundary conditions. */ +BoundingBox system_domain_bounds(Vec2d(-LH, -LH), Vec2d(LL + BW, LH + BW)); +//---------------------------------------------------------------------- +// Material parameters. +//---------------------------------------------------------------------- +Real rho0_f = 1.0; /**< Reference density of fluid. */ +Real U_max = 1.0; /**< Characteristic velocity. */ +Real c_f = 100.0; /**< Reference sound speed. */ +//---------------------------------------------------------------------- +// Geometric shapes used in this case. +//---------------------------------------------------------------------- +Vec2d DamP_lb(-LL / 2, -LH / 2); /**< Left bottom. */ +Vec2d DamP_lt(-LL / 2, LH / 2); /**< Left top. */ +Vec2d DamP_rt(LL / 2, LH / 2); /**< Right top. */ +Vec2d DamP_rb(LL / 2, -LH / 2); /**< Right bottom. */ +//---------------------------------------------------------------------- +// Geometry definition. +//---------------------------------------------------------------------- +std::vector createWaterBlockShape() +{ + std::vector water_block_shape; + water_block_shape.push_back(DamP_lb); + water_block_shape.push_back(DamP_lt); + water_block_shape.push_back(DamP_rt); + water_block_shape.push_back(DamP_rb); + water_block_shape.push_back(DamP_lb); + + return water_block_shape; +} +class WaterBlock : public MultiPolygonShape +{ +public: + explicit WaterBlock(const std::string& shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(createWaterBlockShape(), ShapeBooleanOps::add); + } +}; + +/** + * application dependent initial velocity + */ +class InitialVelocity + : public fluid_dynamics::FluidInitialCondition +{ +public: + InitialVelocity(SPHBody& sph_body) + : fluid_dynamics::FluidInitialCondition(sph_body) {}; + + void update(size_t index_i, Real dt) + { + /** initial velocity profile */ + if (pos_[index_i][1]>=0.0) + { + vel_[index_i][1] = -1.0; + } + else + { + vel_[index_i][1] = 1.0; + } + + } +}; +//---------------------------------------------------------------------- +// wave gauge +//---------------------------------------------------------------------- +Real h = 1.3 * particle_spacing_ref; +MultiPolygon createWaveProbeShape() +{ + std::vector pnts; + pnts.push_back(Vecd(1.0 - h, 0.0)); + pnts.push_back(Vecd(1.0 - h, LH)); + pnts.push_back(Vecd(1.0 + h, LH)); + pnts.push_back(Vecd(1.0 + h, 0.0)); + pnts.push_back(Vecd(1.0 - h, 0.0)); + + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add); + return multi_polygon; +} +//---------------------------------------------------------------------- +// Main program starts here. +//---------------------------------------------------------------------- +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up an SPHSystem. + //---------------------------------------------------------------------- + BoundingBox system_domain_bounds(Vec2d(-LH, -LH), Vec2d(LL + BW, LH + BW)); + SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); + /** Tag for run particle relaxation for the initial body fitted distribution. */ + sph_system.setRunParticleRelaxation(false); + /** Tag for computation start with relaxed body fitted particles distribution. */ + sph_system.setReloadParticles(false); + sph_system.handleCommandlineOptions(ac, av); + IOEnvironment io_environment(sph_system); + //---------------------------------------------------------------------- + // Creating bodies with corresponding materials and particles. + //---------------------------------------------------------------------- + FluidBody water_block(sph_system, makeShared("WaterBody")); + water_block.defineAdaptation(1.3, 1.0); + water_block.defineParticlesAndMaterial(rho0_f, c_f); + // Using relaxed particle distribution if needed + (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) + ? water_block.generateParticles(io_environment, water_block.getName()) + : water_block.generateParticles(); + //---------------------------------------------------------------------- + // 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. + //---------------------------------------------------------------------- + InnerRelation water_body_inner(water_block); + //---------------------------------------------------------------------- + //---------------------------------------------------------------------- + // Define the numerical methods used in the simulation. + // Note that there may be data dependence on the sequence of constructions. + //---------------------------------------------------------------------- + SimpleDynamics initial_condition(water_block); + /** time-space method to detect surface particles. */ + InteractionWithUpdate + free_surface_indicator(water_body_inner); + /** Apply transport velocity formulation. */ + InteractionDynamics> transport_velocity_correction(water_body_inner); + + Dynamics1Level fluid_pressure_relaxation_correct(water_body_inner); + Dynamics1Level fluid_density_relaxation(water_body_inner); + InteractionWithUpdate corrected_configuration_fluid(water_body_inner, 2, 1.0); + InteractionWithUpdate fluid_density_by_summation(water_body_inner); + SimpleDynamics fluid_step_initialization(water_block); + ReduceDynamics fluid_advection_time_step(water_block, U_max); + ReduceDynamics fluid_acoustic_time_step(water_block); + /** We can output a method-specific particle data for debug */ + water_block.addBodyStateForRecording("Pressure"); + water_block.addBodyStateForRecording("Indicator"); + water_block.addBodyStateForRecording("PositionDivergence"); + //---------------------------------------------------------------------- + // Define the methods for I/O operations, observations + // and regression tests of the simulation. + //---------------------------------------------------------------------- + BodyStatesRecordingToVtp body_states_recording(io_environment, sph_system.real_bodies_); + RestartIO restart_io(io_environment, sph_system.real_bodies_); + RegressionTestDynamicTimeWarping>> + write_water_mechanical_energy(io_environment, water_block); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + initial_condition.exec(); + //---------------------------------------------------------------------- + // Load restart file if necessary. + //---------------------------------------------------------------------- + if (sph_system.RestartStep() != 0) + { + GlobalStaticVariables::physical_time_ = restart_io.readRestartFiles(sph_system.RestartStep()); + water_block.updateCellLinkedList(); + water_body_inner.updateConfiguration(); + } + //---------------------------------------------------------------------- + // 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; + int restart_output_interval = screen_output_interval * 10; + Real end_time = 1.0; + Real output_interval = 0.005; + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); + TimeInterval interval; + TimeInterval interval_computing_time_step; + TimeInterval interval_computing_fluid_pressure_relaxation; + TimeInterval interval_updating_configuration; + TickCount time_instance; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + body_states_recording.writeToFile(); + write_water_mechanical_energy.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_interval) + { + /** outer loop for dual-time criteria time-stepping. */ + time_instance = TickCount::now(); + fluid_step_initialization.exec(); + Real advection_dt = 0.3 *fluid_advection_time_step.exec(); + free_surface_indicator.exec(); + fluid_density_by_summation.exec(); + corrected_configuration_fluid.exec(); + transport_velocity_correction.exec(); + interval_computing_time_step += TickCount::now() - time_instance; + + time_instance = TickCount::now(); + Real relaxation_time = 0.0; + Real acoustic_dt = 0.0; + while (relaxation_time < advection_dt) + { + /** inner loop for dual-time criteria time-stepping. */ + acoustic_dt = fluid_acoustic_time_step.exec(); + fluid_pressure_relaxation_correct.exec(acoustic_dt); + fluid_density_relaxation.exec(acoustic_dt); + relaxation_time += acoustic_dt; + integration_time += acoustic_dt; + GlobalStaticVariables::physical_time_ += acoustic_dt; + } + interval_computing_fluid_pressure_relaxation += TickCount::now() - time_instance; + + /** screen output, write body reduced values and restart files */ + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ + << " advection_dt = " << advection_dt << " acoustic_dt = " << acoustic_dt << "\n"; + + if (number_of_iterations % restart_output_interval == 0) + restart_io.writeToFile(number_of_iterations); + } + number_of_iterations++; + + write_water_mechanical_energy.writeToFile(number_of_iterations); + + /** Update cell linked list and configuration. */ + time_instance = TickCount::now(); + water_block.updateCellLinkedListWithParticleSort(100); + water_body_inner.updateConfiguration(); + interval_updating_configuration += TickCount::now() - time_instance; + } + + body_states_recording.writeToFile(); + TickCount t2 = TickCount::now(); + 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_fluid_pressure_relaxation = " + << interval_computing_fluid_pressure_relaxation.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_updating_configuration = " + << interval_updating_configuration.seconds() << "\n"; + + return 0; +}; diff --git a/tests/user_examples/test_2d_impact_patch/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_2d_impact_patch/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..c3b7c39dd2 --- /dev/null +++ b/tests/user_examples/test_2d_impact_patch/regression_test_tool/regression_test_tool.py @@ -0,0 +1,44 @@ +# !/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_standing_wave +""" + +case_name = "test_2d_standing_wave" +body_name = "WaterBody" +parameter_name = "TotalMechanicalEnergy" +body_name_1 = "WaveProbe" +parameter_name_1 = "FreeSurfaceHeight" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) +sphinxsys_1 = SphinxsysRegressionTest(case_name, body_name_1, parameter_name_1) + + +while True: + print("Now start a new run......") + sphinxsys.run_case() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + converged_1 = sphinxsys_1.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if (converged == "true") and (converged_1 == "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 + elif converged_1 != "true": + print("The tested parameters of", sphinxsys_1.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break diff --git a/tests/user_examples/test_2d_ration_patch/CMakeLists.txt b/tests/user_examples/test_2d_ration_patch/CMakeLists.txt new file mode 100644 index 0000000000..1e1bc16be3 --- /dev/null +++ b/tests/user_examples/test_2d_ration_patch/CMakeLists.txt @@ -0,0 +1,25 @@ +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_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}) +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} ) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") +target_link_libraries(${PROJECT_NAME} extra_sources_2d) + +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +add_test(NAME ${PROJECT_NAME}_restart COMMAND ${PROJECT_NAME} --restart_step=4000 + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) diff --git a/tests/user_examples/test_2d_ration_patch/ratotion_patch.cpp b/tests/user_examples/test_2d_ration_patch/ratotion_patch.cpp new file mode 100644 index 0000000000..7b107a9741 --- /dev/null +++ b/tests/user_examples/test_2d_ration_patch/ratotion_patch.cpp @@ -0,0 +1,283 @@ +/** + * @file rotationpatch.cpp + * @brief 2D rotationpatch example. + * @author Yaru Ren, Chi Zhang and Xiangyu Hu + */ + +#include "sphinxsys.h" //SPHinXsys Library. +using namespace SPH; // Namespace cite here. +#define PI 3.1415926 +//---------------------------------------------------------------------- +// Basic geometry parameters and numerical setup. +//---------------------------------------------------------------------- +Real LL = 1.0; /**< Liquid column length. */ +Real LH = 1.0; /**< Liquid column height. */ +Real particle_spacing_ref = LL/200; /**< Initial reference particle spacing. */ +Real BW = particle_spacing_ref * 4; /**< Extending width for boundary conditions. */ +BoundingBox system_domain_bounds(Vec2d(-BW, -BW), Vec2d(LL + BW, LH + BW)); +//---------------------------------------------------------------------- +// Material parameters. +//---------------------------------------------------------------------- +Real rho0_f = 1.0; /**< Reference density of fluid. */ +Real U_max = 1.0; /**< Characteristic velocity. */ +Real c_f = 5.0 * sqrt(2.0); /**< Reference sound speed. */ + +Vec2d DamP_lb(-LL / 2, -LH / 2); /**< Left bottom. */ +Vec2d DamP_lt(-LL / 2, LH / 2); /**< Left top. */ +Vec2d DamP_rt(LL / 2, LH / 2); /**< Right top. */ +Vec2d DamP_rb(LL / 2, -LH / 2); /**< Right bottom. */ + +//---------------------------------------------------------------------- +// Geometry definition. +//---------------------------------------------------------------------- +std::vector createWaterBlockShape() +{ + std::vector water_block_shape; + water_block_shape.push_back(DamP_lb); + water_block_shape.push_back(DamP_lt); + water_block_shape.push_back(DamP_rt); + water_block_shape.push_back(DamP_rb); + water_block_shape.push_back(DamP_lb); + + return water_block_shape; +} +class WaterBlock : public MultiPolygonShape +{ +public: + explicit WaterBlock(const std::string& shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(createWaterBlockShape(), ShapeBooleanOps::add); + } +}; + +/** + * application dependent initial velocity + */ +class InitialVelocity + : public fluid_dynamics::FluidInitialCondition +{ +public: + InitialVelocity(SPHBody& sph_body) + : fluid_dynamics::FluidInitialCondition(sph_body), + fluid_particles_(dynamic_cast(&sph_body.getBaseParticles())), + p_(*fluid_particles_->getVariableByName("Pressure")), rho_(fluid_particles_->rho_) {}; + + void update(size_t index_i, Real dt) + { + Real omega = 1.0; + /** initial velocity profile */ + vel_[index_i][0] = omega * pos_[index_i][1]; + vel_[index_i][1] = -omega * pos_[index_i][0]; + + /** initial pressure */ + for (size_t m = 0; m!= 100; ++m) + { + for (size_t n = 0; n!= 100; ++n) + { + if (m % 2 == 1 && n % 2 == 1) + { + Real x_star = pos_[index_i][0] + LL / 2; + Real y_star = pos_[index_i][1] + LL / 2; + Real coefficient1 = m * n * PI * PI * (pow((m * PI / LL), 2) + pow((n * PI / LL), 2)); + p_[index_i] += rho0_f * (-32 * omega * omega) / coefficient1 * sin(m * PI * x_star / LL) + * sin(n * PI * y_star / LL); + } + } + } + + rho_[index_i] = p_[index_i]/ pow(c_f,2) + rho0_f; + } + +protected: + BaseParticles* fluid_particles_; + StdLargeVec&p_,&rho_; +}; +//---------------------------------------------------------------------- +// wave gauge +//---------------------------------------------------------------------- +Real h = 1.3 * particle_spacing_ref; +MultiPolygon createWaveProbeShape() +{ + std::vector pnts; + pnts.push_back(Vecd(0.0, 0.0)); + + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add); + return multi_polygon; +} + +//---------------------------------------------------------------------- +// Main program starts here. +//---------------------------------------------------------------------- +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up an SPHSystem. + //---------------------------------------------------------------------- + BoundingBox system_domain_bounds(Vec2d(-LL, -LL), Vec2d(LL + 10 * BW, LH + 10 * BW)); + SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); + /** Tag for run particle relaxation for the initial body fitted distribution. */ + sph_system.setRunParticleRelaxation(false); + /** Tag for computation start with relaxed body fitted particles distribution. */ + sph_system.setReloadParticles(false); + sph_system.handleCommandlineOptions(ac, av); + IOEnvironment io_environment(sph_system); + //---------------------------------------------------------------------- + // Creating bodies with corresponding materials and particles. + //---------------------------------------------------------------------- + FluidBody water_block(sph_system, makeShared("WaterBody")); + water_block.defineAdaptation(1.3, 1.0); + water_block.defineParticlesAndMaterial(rho0_f, c_f); + // Using relaxed particle distribution if needed + (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) + ? water_block.generateParticles(io_environment, water_block.getName()) + : water_block.generateParticles(); + + ObserverBody fluid_observer(sph_system, "FluidObserver"); + StdVec observation_location = { Vecd(0.0, 0.0) }; + fluid_observer.generateParticles(observation_location); + //------------------------------------------------------------------------------------------------------------------------------------------- + InnerRelation water_body_inner(water_block); + ContactRelation fluid_observer_contact(fluid_observer, { &water_block }); + //---------------------------------------------------------------------- + // Define the numerical methods used in the simulation. + // Note that there may be data dependence on the sequence of constructions. + //---------------------------------------------------------------------- + SimpleDynamics initial_condition(water_block); + /** time-space method to detect surface particles. */ + InteractionWithUpdate + free_surface_indicator(water_body_inner); + /** Apply transport velocity formulation. */ + InteractionDynamics> transport_velocity_correction(water_body_inner); + + Dynamics1Level fluid_pressure_relaxation_correct(water_body_inner); + Dynamics1Level fluid_density_relaxation(water_body_inner); + InteractionWithUpdate corrected_configuration_fluid(water_body_inner, 2, 1.0); + SimpleDynamics fluid_step_initialization(water_block); + ReduceDynamics fluid_advection_time_step(water_block, U_max); + ReduceDynamics fluid_acoustic_time_step(water_block); + /** We can output a method-specific particle data for debug */ + water_block.addBodyStateForRecording("Pressure"); + water_block.addBodyStateForRecording("Indicator"); + //---------------------------------------------------------------------- + // Define the methods for I/O operations, observations + // and regression tests of the simulation. + //---------------------------------------------------------------------- + BodyStatesRecordingToVtp body_states_recording(io_environment, sph_system.real_bodies_); + RestartIO restart_io(io_environment, sph_system.real_bodies_); + RegressionTestDynamicTimeWarping>> + write_water_mechanical_energy(io_environment, water_block); + RegressionTestDynamicTimeWarping> + write_recorded_water_pressure("Pressure", io_environment, fluid_observer_contact); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + initial_condition.exec(); + //---------------------------------------------------------------------- + // Load restart file if necessary. + //---------------------------------------------------------------------- + if (sph_system.RestartStep() != 0) + { + GlobalStaticVariables::physical_time_ = restart_io.readRestartFiles(sph_system.RestartStep()); + water_block.updateCellLinkedList(); + water_body_inner.updateConfiguration(); + } + //---------------------------------------------------------------------- + // 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; + int restart_output_interval = screen_output_interval * 10; + Real end_time = 8.0; + Real output_interval = 0.01; + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); + TimeInterval interval; + TimeInterval interval_computing_time_step; + TimeInterval interval_computing_fluid_pressure_relaxation; + TimeInterval interval_updating_configuration; + TickCount time_instance; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + body_states_recording.writeToFile(); + write_water_mechanical_energy.writeToFile(number_of_iterations); + write_recorded_water_pressure.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_interval) + { + /** outer loop for dual-time criteria time-stepping. */ + time_instance = TickCount::now(); + fluid_step_initialization.exec(); + free_surface_indicator.exec(); + corrected_configuration_fluid.exec(); + transport_velocity_correction.exec(); + interval_computing_time_step += TickCount::now() - time_instance; + + time_instance = TickCount::now(); + Real relaxation_time = 0.0; + Real acoustic_dt = 0.0; + + /** inner loop for dual-time criteria time-stepping. */ + acoustic_dt = fluid_acoustic_time_step.exec(); + fluid_pressure_relaxation_correct.exec(acoustic_dt); + fluid_density_relaxation.exec(acoustic_dt); + relaxation_time += acoustic_dt; + integration_time += acoustic_dt; + GlobalStaticVariables::physical_time_ += acoustic_dt; + + interval_computing_fluid_pressure_relaxation += TickCount::now() - time_instance; + + /** screen output, write body reduced values and restart files */ + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ << " acoustic_dt = " << acoustic_dt << "\n"; + + if (number_of_iterations % restart_output_interval == 0) + restart_io.writeToFile(number_of_iterations); + } + number_of_iterations++; + + write_water_mechanical_energy.writeToFile(number_of_iterations); + write_recorded_water_pressure.writeToFile(number_of_iterations); + + /** Update cell linked list and configuration. */ + time_instance = TickCount::now(); + water_block.updateCellLinkedListWithParticleSort(100); + water_body_inner.updateConfiguration(); + fluid_observer_contact.updateConfiguration(); + interval_updating_configuration += TickCount::now() - time_instance; + } + body_states_recording.writeToFile(); + TickCount t2 = TickCount::now(); + 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_fluid_pressure_relaxation = " + << interval_computing_fluid_pressure_relaxation.seconds() << "\n"; + std::cout << std::fixed << std::setprecision(9) << "interval_updating_configuration = " + << interval_updating_configuration.seconds() << "\n"; + + return 0; +}; diff --git a/tests/user_examples/test_2d_ration_patch/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_2d_ration_patch/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..c3b7c39dd2 --- /dev/null +++ b/tests/user_examples/test_2d_ration_patch/regression_test_tool/regression_test_tool.py @@ -0,0 +1,44 @@ +# !/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_standing_wave +""" + +case_name = "test_2d_standing_wave" +body_name = "WaterBody" +parameter_name = "TotalMechanicalEnergy" +body_name_1 = "WaveProbe" +parameter_name_1 = "FreeSurfaceHeight" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) +sphinxsys_1 = SphinxsysRegressionTest(case_name, body_name_1, parameter_name_1) + + +while True: + print("Now start a new run......") + sphinxsys.run_case() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + converged_1 = sphinxsys_1.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if (converged == "true") and (converged_1 == "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 + elif converged_1 != "true": + print("The tested parameters of", sphinxsys_1.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break