diff --git a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp new file mode 100644 index 0000000000..037b7be579 --- /dev/null +++ b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.cpp @@ -0,0 +1,29 @@ +#include "wetting_coupled_spatial_temporal_method.h" + +namespace SPH +{ +//=====================================================================================================// +namespace fluid_dynamics +{ +//=================================================================================================// +NonWettingSurfaceIndication:: + NonWettingSurfaceIndication(BaseInnerRelation &inner_relation, + BaseContactRelation &contact_relation, Real threshold, Real criterion) + : FreeSurfaceIndicationComplex(inner_relation, contact_relation, threshold), wetting_criterion(criterion) +{ + for (size_t k = 0; k != contact_particles_.size(); ++k) + { + contact_phi_.push_back(this->contact_particles_[k]->template getVariableByName("Phi")); + } +} +//=================================================================================================// +NonWettingSurfaceIndication:: + NonWettingSurfaceIndication(ComplexRelation &complex_relation, Real threshold, Real criterion) + : NonWettingSurfaceIndication(complex_relation.getInnerRelation(), + complex_relation.getContactRelation(), threshold, criterion) {} + +//=================================================================================================// +} // namespace fluid_dynamics + //=================================================================================================// +} // namespace SPH + //=================================================================================================// \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h new file mode 100644 index 0000000000..bed6605bfd --- /dev/null +++ b/tests/user_examples/extra_src/shared/wetting_coupled_spatial_temporal_method.h @@ -0,0 +1,107 @@ +/* ------------------------------------------------------------------------- * + * SPHinXsys * + * ------------------------------------------------------------------------- * + * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle * + * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * + * physical accurate simulation and aims to model coupled industrial dynamic * + * systems including fluid, solid, multi-body dynamics and beyond with SPH * + * (smoothed particle hydrodynamics), a meshless computational method using * + * particle discretization. * + * * + * SPHinXsys is partially funded by German Research Foundation * + * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * + * HU1527/12-1 and HU1527/12-4. * + * * + * Portions copyright (c) 2017-2023 Technical University of Munich and * + * the authors' affiliations. * + * * + * Licensed under the Apache License, Version 2.0 (the "License"); you may * + * not use this file except in compliance with the License. You may obtain a * + * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * + * * + * ------------------------------------------------------------------------- */ +/** + * @file fluid_surface_complex.h + * @brief Here, we define the algorithm classes for fluid dynamics within the body. + * @details We consider here weakly compressible fluids. The algorithms may be + * different for free surface flow and the one without free surface. + * @author Chi Zhang and Xiangyu Hu + */ + +#ifndef WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H +#define WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H + +#include "fluid_surface_complex.h" + +namespace SPH +{ +namespace fluid_dynamics +{ +/** + * @class NonWettingSurfaceIndication + * @brief Non wetting surface particles include free-surface ones and interfacial ones near the non-wetted structure. + * @brief Even the position divergence of interfacial fluid pariticles has satisfied with the threshold of spatial-temporal + identification approach to be identified as internal ones,they will remain as free-surface ones if without + any wetted neighboring solid particles. + */ +class NonWettingSurfaceIndication : public FreeSurfaceIndicationComplex +{ + public: + NonWettingSurfaceIndication(BaseInnerRelation &inner_relation, + BaseContactRelation &contact_relation, Real threshold = 0.75, Real criterion = 0.99); + explicit NonWettingSurfaceIndication(ComplexRelation &complex_relation, Real threshold = 0.75, Real criterion = 0.99); + virtual ~NonWettingSurfaceIndication(){}; + + inline void interaction(size_t index_i, Real dt = 0.0) + { + FreeSurfaceIndicationInner::interaction(index_i, dt); + + Real pos_div = 0.0; + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + pos_div -= contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.r_ij_[n]; + } + } + pos_div_[index_i] += pos_div; + + if (pos_div_[index_i] > this->threshold_by_dimensions_) + { + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t j = contact_neighborhood.j_[n]; + if ((*(contact_phi_[k]))[j] > wetting_criterion) + { + pos_div_[index_i] = 2.0 * this->threshold_by_dimensions_; + break; + } + else + { + pos_div_[index_i] = 0.5 * this->threshold_by_dimensions_; + } + } + if (pos_div_[index_i] == 2.0 * this->threshold_by_dimensions_) + break; + } + } + }; + + protected: + Real wetting_criterion; + StdVec *> contact_phi_; +}; + +using WettingCoupledSpatialTemporalFreeSurfaceIdentificationComplex = + SpatialTemporalFreeSurfaceIdentification; +using SpatialTemporalFreeSurfaceIdentificationComplex = + SpatialTemporalFreeSurfaceIdentification; + + +} // namespace fluid_dynamics +} // namespace SPH +#endif // WETTING_COUPLED_SPATIAL_TEMPORAL_COMPLEX_H \ No newline at end of file diff --git a/tests/user_examples/test_2d_water_entry_exit/CMakeLists.txt b/tests/user_examples/test_2d_water_entry_exit/CMakeLists.txt new file mode 100644 index 0000000000..97eef33f2b --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/CMakeLists.txt @@ -0,0 +1,38 @@ +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} ) + +if(${CMAKE_SYSTEM_NAME} MATCHES "Windows") + add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --r=true + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --r=false --i=true + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +else() + file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/run_test.sh + DESTINATION ${EXECUTABLE_OUTPUT_PATH}) + add_test(NAME ${PROJECT_NAME} COMMAND bash ${EXECUTABLE_OUTPUT_PATH}/run_test.sh + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +endif() + +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/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_0_result.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_0_result.xml new file mode 100644 index 0000000000..8c7499408f --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_10_result.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_10_result.xml new file mode 100644 index 0000000000..30e1e009d6 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_10_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_20_result.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_20_result.xml new file mode 100644 index 0000000000..4e82ba0c88 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_Run_20_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_dtwdistance.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_dtwdistance.xml new file mode 100644 index 0000000000..937895ceed --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_runtimes.dat b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_runtimes.dat new file mode 100644 index 0000000000..fd8ad08880 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/FluidObserver_Pressure_runtimes.dat @@ -0,0 +1,3 @@ +true +21 +4 \ No newline at end of file diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_0_result.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_0_result.xml new file mode 100644 index 0000000000..040faac7d7 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_11_result.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_11_result.xml new file mode 100644 index 0000000000..a4f9d99f27 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_11_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_6_result.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_6_result.xml new file mode 100644 index 0000000000..2b0252d34c --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_Run_6_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_dtwdistance.xml new file mode 100644 index 0000000000..eaf89cd547 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_runtimes.dat b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_runtimes.dat new file mode 100644 index 0000000000..5cc2fa0e29 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/WaterBody_TotalMechanicalEnergy_runtimes.dat @@ -0,0 +1,3 @@ +true +12 +4 \ No newline at end of file diff --git a/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_2d_water_entry_exit/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..8f3d795957 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/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_dambreak +""" + +case_name = "test_2d_dambreak" +body_name = "WaterBody" +parameter_name = "TotalMechanicalEnergy" +body_name_1 = "FluidObserver" +parameter_name_1 = "Pressure" + +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_water_entry_exit/run_test.sh b/tests/user_examples/test_2d_water_entry_exit/run_test.sh new file mode 100644 index 0000000000..bf83adea29 --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/run_test.sh @@ -0,0 +1,2 @@ +./test_2d_water_entry_exit --r=true +./test_2d_water_entry_exit --r=false --i=true diff --git a/tests/user_examples/test_2d_water_entry_exit/test_2d_water_entry_exit.cpp b/tests/user_examples/test_2d_water_entry_exit/test_2d_water_entry_exit.cpp new file mode 100644 index 0000000000..f0e3967c9f --- /dev/null +++ b/tests/user_examples/test_2d_water_entry_exit/test_2d_water_entry_exit.cpp @@ -0,0 +1,539 @@ +/** + * @file water entry and exit.cpp + * @brief 2D water entry and exit example with surface wettability considered. + * @details This is the one of FSI test cases, also one case for + * understanding spatial temporal identification approach, + * especially when coupled with the wetting. + * @author Shuoguo Zhang and Xiangyu Hu + */ +#include "sphinxsys.h" //SPHinXsys Library. +#include "wetting_coupled_spatial_temporal_method.h" + +using namespace SPH; // Namespace cite here. +//---------------------------------------------------------------------- +// Basic geometry parameters and numerical setup. +//---------------------------------------------------------------------- +Real insert_cylinder_radius = 0.055; /**< Cylinder radius. */ +Real DL = 8.0 * insert_cylinder_radius; /**< Water tank length. */ +Real DH = 7.0 * insert_cylinder_radius; /**< Water tank height. */ +Real LL = DL; /**< Water column length. */ +Real LH = 3.0 * insert_cylinder_radius; /**< Water column height. */ +Real particle_spacing_ref = 2.0 * insert_cylinder_radius / 40.0; /**< Initial reference particle spacing. */ +Real BW = particle_spacing_ref * 4; /**< Thickness of tank wall. */ +Vec2d insert_cylinder_center(0.5 * DL, LH + 0.15); /**< Location of the cylinder center. */ +Vecd tethering_point(0.5 * DL, DH); /**< The tethering point. */ +StdVec observer_location = {insert_cylinder_center}; /**< Observation point. */ +//---------------------------------------------------------------------- +// Material parameters. +//---------------------------------------------------------------------- +Real rho0_f = 1.0; /**< Fluid density. */ +Real rho0_s = 0.5; /**< Cylinder density. */ +Real gravity_g = 9.81; /**< Gravity. */ +Real U_max = 2.0 * sqrt(gravity_g * LH); /**< Characteristic velocity. */ +Real c_f = 10.0 * U_max; /**< Reference sound speed. */ +Real mu_f = 8.9e-7; /**< Water dynamics viscosity. */ +//---------------------------------------------------------------------- +// Wetting parameters +//---------------------------------------------------------------------- +Real diffusion_coeff = 330.578 * pow(particle_spacing_ref,2); /**< Wetting coefficient. */ +Real fluid_moisture = 1.0; /**< fluid moisture. */ +Real cylinder_moisture = 0.0; /**< cylinder moisture. */ +Real wall_moisture = 1.0; /**< wall moisture. */ +//---------------------------------------------------------------------- +// Definition for water body +//---------------------------------------------------------------------- +std::vector createWaterBlockShape() +{ + std::vector water_block; + water_block.push_back(Vecd(0.0, 0.0)); + water_block.push_back(Vecd(0.0, LH)); + water_block.push_back(Vecd(LL, LH)); + water_block.push_back(Vecd(LL, 0.0)); + water_block.push_back(Vecd(0.0, 0.0)); + + return water_block; +} +class WettingFluidBody : public MultiPolygonShape +{ + public: + explicit WettingFluidBody(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(createWaterBlockShape(), ShapeBooleanOps::add); + } +}; +class WettingFluidBodyMaterial : public DiffusionReaction +{ + public: + WettingFluidBodyMaterial() + : DiffusionReaction({"Phi"}, SharedPtr(), rho0_f, c_f, mu_f) + { + initializeAnDiffusion("Phi", "Phi"); + }; +}; +using DiffusionFluidParticles = DiffusionReactionParticles; +class WettingFluidBodyInitialCondition + : public DiffusionReactionInitialCondition +{ + protected: + size_t phi_; + + public: + explicit WettingFluidBodyInitialCondition(SPHBody &sph_body) + : DiffusionReactionInitialCondition(sph_body) + { + phi_ = particles_->diffusion_reaction_material_.AllSpeciesIndexMap()["Phi"]; + }; + + void update(size_t index_i, Real dt) + { + all_species_[phi_][index_i] = fluid_moisture; + }; +}; +//---------------------------------------------------------------------- +// Definition for wall body +//---------------------------------------------------------------------- +std::vector createOuterWallShape() +{ + std::vector outer_wall; + outer_wall.push_back(Vecd(-BW, -BW)); + outer_wall.push_back(Vecd(-BW, DH + BW)); + outer_wall.push_back(Vecd(DL + BW, DH + BW)); + outer_wall.push_back(Vecd(DL + BW, -BW)); + outer_wall.push_back(Vecd(-BW, -BW)); + + return outer_wall; +} +std::vector createInnerWallShape() +{ + std::vector inner_wall; + inner_wall.push_back(Vecd(0.0, 0.0)); + inner_wall.push_back(Vecd(0.0, DH)); + inner_wall.push_back(Vecd(DL, DH)); + inner_wall.push_back(Vecd(DL, 0.0)); + inner_wall.push_back(Vecd(0.0, 0.0)); + + return inner_wall; +} +class WettingWallBody : public MultiPolygonShape +{ + public: + explicit WettingWallBody(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(createOuterWallShape(), ShapeBooleanOps::add); + multi_polygon_.addAPolygon(createInnerWallShape(), ShapeBooleanOps::sub); + } +}; +class WettingWallBodyMaterial : public DiffusionReaction +{ + public: + WettingWallBodyMaterial() : DiffusionReaction({"Phi"}, SharedPtr()) + { + initializeAnDiffusion("Phi", "Phi"); + }; +}; +using DiffusionWallParticles = DiffusionReactionParticles; +class WettingWallBodyInitialCondition + : public DiffusionReactionInitialCondition +{ + protected: + size_t phi_; + + public: + explicit WettingWallBodyInitialCondition(SPHBody &sph_body) + : DiffusionReactionInitialCondition(sph_body) + { + phi_ = particles_->diffusion_reaction_material_.AllSpeciesIndexMap()["Phi"]; + }; + + void update(size_t index_i, Real dt) + { + all_species_[phi_][index_i] = wall_moisture; + }; +}; +//---------------------------------------------------------------------- +// Definition for cylinder body +//---------------------------------------------------------------------- +class WettingCylinderBody : public MultiPolygonShape +{ + public: + explicit WettingCylinderBody(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addACircle(insert_cylinder_center, insert_cylinder_radius, 100, ShapeBooleanOps::add); + } +}; +class WettingCylinderBodyMaterial : public DiffusionReaction +{ + public: + WettingCylinderBodyMaterial() : DiffusionReaction({"Phi"}, SharedPtr(), rho0_s) + { + initializeAnDiffusion("Phi", "Phi",diffusion_coeff); + }; +}; +using DiffusionCylinderParticles = DiffusionReactionParticles; +class WettingCylinderBodyInitialCondition + : public DiffusionReactionInitialCondition +{ + protected: + size_t phi_; + + public: + explicit WettingCylinderBodyInitialCondition(SPHBody &sph_body) + : DiffusionReactionInitialCondition(sph_body) + { + phi_ = particles_->diffusion_reaction_material_.AllSpeciesIndexMap()["Phi"]; + }; + + void update(size_t index_i, Real dt) + { + all_species_[phi_][index_i] = cylinder_moisture; + }; +}; + +//---------------------------------------------------------------------- +// Set topology for wetting bodies +//---------------------------------------------------------------------- +using CylinderFluidDiffusionDirichlet = DiffusionRelaxationDirichlet; +class ThermalRelaxationComplex + : public DiffusionRelaxationRK2< + ComplexInteraction> +{ + public: + explicit ThermalRelaxationComplex(BaseContactRelation &body_contact_relation_Dirichlet) + : DiffusionRelaxationRK2>(body_contact_relation_Dirichlet){}; + virtual ~ThermalRelaxationComplex(){}; +}; +//------------------------------------------------------------------------------ +// Constrained part for Simbody +//------------------------------------------------------------------------------ +MultiPolygon createSimbodyConstrainShape(SPHBody &sph_body) +{ + MultiPolygon multi_polygon; + multi_polygon.addACircle(insert_cylinder_center, insert_cylinder_radius, 100, ShapeBooleanOps::add); + return multi_polygon; +}; +//---------------------------------------------------------------------- +// Main program starts here. +//---------------------------------------------------------------------- +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up an SPHSystem. + //---------------------------------------------------------------------- + BoundingBox system_domain_bounds(Vec2d(-BW, -BW), Vec2d(DL + BW, DH + BW)); + SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); + sph_system.setRunParticleRelaxation(false); + sph_system.setReloadParticles(true); + sph_system.handleCommandlineOptions(ac, av); + GlobalStaticVariables::physical_time_=0.0; + IOEnvironment io_environment(sph_system); + //---------------------------------------------------------------------- + // Creating bodies with corresponding materials and particles. + //---------------------------------------------------------------------- + FluidBody water_block(sph_system, makeShared("WaterBody")); + water_block.defineParticlesAndMaterial(); + water_block.generateParticles(); + + SolidBody wall_boundary(sph_system, makeShared("WallBoundary")); + wall_boundary.defineParticlesAndMaterial(); + wall_boundary.generateParticles(); + wall_boundary.addBodyStateForRecording("NormalDirection"); + + SolidBody cylinder(sph_system, makeShared("Cylinder")); + cylinder.defineAdaptationRatios(1.15, 1.0); + cylinder.defineBodyLevelSetShape(); + cylinder.defineParticlesAndMaterial(); + (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles()) + ? cylinder.generateParticles(io_environment, cylinder.getName()) + : cylinder.generateParticles(); + + ObserverBody fluid_observer(sph_system, "FluidObserver"); + fluid_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. + //---------------------------------------------------------------------- + InnerRelation water_block_inner(water_block); + InnerRelation cylinder_inner(cylinder); + ComplexRelation water_block_complex(water_block_inner, {&wall_boundary, &cylinder}); + ContactRelation cylinder_contact(cylinder, {&water_block}); + ContactRelation fluid_observer_contact(fluid_observer, {&cylinder}); + //---------------------------------------------------------------------- + // Run particle relaxation for body-fitted distribution if chosen. + //---------------------------------------------------------------------- + if (sph_system.RunParticleRelaxation()) + { + /** body topology only for particle relaxation */ + InnerRelation cylinder_inner(cylinder); + //---------------------------------------------------------------------- + // Methods used for particle relaxation. + //---------------------------------------------------------------------- + /** Random reset the insert body particle position. */ + SimpleDynamics random_inserted_body_particles(cylinder); + /** Write the body state to Vtp file. */ + BodyStatesRecordingToVtp write_inserted_body_to_vtp(io_environment, {&cylinder}); + /** Write the particle reload files. */ + ReloadParticleIO write_particle_reload_files(io_environment, {&cylinder}); + /** A Physics relaxation step. */ + relax_dynamics::RelaxationStepInner relaxation_step_inner(cylinder_inner); + //---------------------------------------------------------------------- + // Particle relaxation starts here. + //---------------------------------------------------------------------- + random_inserted_body_particles.exec(0.25); + relaxation_step_inner.SurfaceBounding().exec(); + write_inserted_body_to_vtp.writeToFile(0); + //---------------------------------------------------------------------- + // Relax particles of the insert body. + //---------------------------------------------------------------------- + int ite_p = 0; + while (ite_p < 1000) + { + relaxation_step_inner.exec(); + ite_p += 1; + if (ite_p % 200 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the inserted body N = " << ite_p << "\n"; + write_inserted_body_to_vtp.writeToFile(ite_p); + } + } + std::cout << "The physics relaxation process of inserted body finish !" << std::endl; + /** Output results. */ + write_particle_reload_files.writeToFile(0); + return 0; + } + //---------------------------------------------------------------------- + // Define the fluid dynamics used in the simulation. + // Note that there may be data dependence on the sequence of constructions. + //---------------------------------------------------------------------- + SharedPtr gravity_ptr = makeShared(Vecd(0.0, -gravity_g)); + SimpleDynamics fluid_step_initialization(water_block, gravity_ptr); + InteractionWithUpdate + free_stream_surface_indicator(water_block_complex); + InteractionWithUpdate fluid_density_by_summation(water_block_complex); + water_block.addBodyStateForRecording("Pressure"); + water_block.addBodyStateForRecording("Density"); + water_block.addBodyStateForRecording("SurfaceIndicator"); + cylinder.addBodyStateForRecording("Density"); + Dynamics1Level fluid_pressure_relaxation(water_block_complex); + Dynamics1Level fluid_density_relaxation(water_block_complex); + ReduceDynamics fluid_advection_time_step(water_block, U_max); + ReduceDynamics fluid_acoustic_time_step(water_block); + InteractionDynamics viscous_acceleration(water_block_complex); + InteractionDynamics transport_velocity_correction(water_block_complex); + //---------------------------------------------------------------------- + // Define the wetting diffusion dynamics used in the simulation. + //---------------------------------------------------------------------- + SimpleDynamics Wetting_water_initial_condition(water_block); + SimpleDynamics Wetting_wall_initial_condition(wall_boundary); + SimpleDynamics Wetting_cylinder_initial_condition(cylinder); + GetDiffusionTimeStepSize get_thermal_time_step(cylinder); + ThermalRelaxationComplex thermal_relaxation_complex(cylinder_contact); + //---------------------------------------------------------------------- + // Algorithms of FSI. + //---------------------------------------------------------------------- + SimpleDynamics wall_boundary_normal_direction(wall_boundary); + SimpleDynamics cylinder_normal_direction(cylinder); + InteractionDynamics fluid_viscous_force_on_inserted_body(cylinder_contact); + InteractionDynamics + fluid_pressure_force_on_inserted_body(cylinder_contact, fluid_viscous_force_on_inserted_body); + //---------------------------------------------------------------------- + // Building Simbody. + //---------------------------------------------------------------------- + SimTK::MultibodySystem MBsystem; + /** The bodies or matter of the MBsystem. */ + SimTK::SimbodyMatterSubsystem matter(MBsystem); + /** The forces of the MBsystem.*/ + SimTK::GeneralForceSubsystem forces(MBsystem); + /** Mass proeprties of the fixed spot. */ + SimTK::Body::Rigid fixed_spot_info(SimTK::MassProperties(1.0, SimTKVec3(0), SimTK::UnitInertia(1))); + SolidBodyPartForSimbody cylinder_constraint_area(cylinder, makeShared(createSimbodyConstrainShape(cylinder), "cylinder")); + /** Mass properties of the consrained spot. */ + SimTK::Body::Rigid tethered_spot_info(*cylinder_constraint_area.body_part_mass_properties_); + /** Mobility of the fixed spot. */ + SimTK::MobilizedBody::Weld fixed_spot(matter.Ground(), SimTK::Transform(SimTKVec3(tethering_point[0], tethering_point[1], 0.0)), + fixed_spot_info, SimTK::Transform(SimTKVec3(0))); + /** Mobility of the tethered spot. + * Set the mass center as the origin location of the planar mobilizer + */ + Vecd displacement0 = cylinder_constraint_area.initial_mass_center_ - tethering_point; + SimTK::MobilizedBody::Planar tethered_spot(fixed_spot, + SimTK::Transform(SimTKVec3(displacement0[0], displacement0[1], 0.0)), tethered_spot_info, SimTK::Transform(SimTKVec3(0))); + // discreted forces acting on the bodies + SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTK::Vec3(0.0, Real(-9.81), 0.0), 0.0); + SimTK::Force::DiscreteForces force_on_bodies(forces, matter); + fixed_spot_info.addDecoration(SimTK::Transform(), SimTK::DecorativeSphere(0.02)); + tethered_spot_info.addDecoration(SimTK::Transform(), SimTK::DecorativeSphere(0.4)); + SimTK::State state = MBsystem.realizeTopology(); + + /** Time steping method for multibody system.*/ + SimTK::RungeKuttaMersonIntegrator integ(MBsystem); + integ.setAccuracy(1e-3); + integ.setAllowInterpolation(false); + integ.initialize(state); + //---------------------------------------------------------------------- + // Coupling between SimBody and SPH.. + //---------------------------------------------------------------------- + ReduceDynamics + force_on_tethered_spot(cylinder_constraint_area, MBsystem, tethered_spot, integ); + SimpleDynamics + constraint_tethered_spot(cylinder_constraint_area, MBsystem, tethered_spot, integ); + //---------------------------------------------------------------------- + // 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_cylinder_displacement("Position", 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(); + wall_boundary_normal_direction.exec(); + cylinder_normal_direction.exec(); + Wetting_water_initial_condition.exec(); + Wetting_wall_initial_condition.exec(); + Wetting_cylinder_initial_condition.exec(); + Real dt_thermal = get_thermal_time_step.exec(); + free_stream_surface_indicator.exec(); + //---------------------------------------------------------------------- + // Load restart file if necessary. + //---------------------------------------------------------------------- + if (sph_system.RestartStep() != 0) + { + GlobalStaticVariables::physical_time_ = restart_io.readRestartFiles(sph_system.RestartStep()); + water_block.updateCellLinkedList(); + water_block_complex.updateConfiguration(); + fluid_observer_contact.updateConfiguration(); + cylinder.updateCellLinkedList(); + water_block_inner.updateConfiguration(); + cylinder_inner.updateConfiguration(); + cylinder_contact.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 = 0.7; + Real output_interval = end_time/70.0; + //---------------------------------------------------------------------- + // 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_cylinder_displacement.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 Dt = fluid_advection_time_step.exec(); + + fluid_density_by_summation.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; + Real dt = 0.0; + fluid_viscous_force_on_inserted_body.exec(); + while (relaxation_time < Dt) + { + /** inner loop for dual-time criteria time-stepping. */ + dt = SMIN(SMIN(dt_thermal, fluid_acoustic_time_step.exec()), Dt); + fluid_pressure_relaxation.exec(dt); + fluid_pressure_force_on_inserted_body.exec(); + fluid_density_relaxation.exec(dt); + thermal_relaxation_complex.exec(dt); + + integ.stepBy(dt); + SimTK::State &state_for_update = integ.updAdvancedState(); + force_on_bodies.clearAllBodyForces(state_for_update); + force_on_bodies.setOneBodyForce(state_for_update, tethered_spot, force_on_tethered_spot.exec()); + constraint_tethered_spot.exec(); + + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += 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_ + << " Dt = " << Dt << " dt = " << dt << "\n"; + + if (number_of_iterations % observation_sample_interval == 0 && number_of_iterations != sph_system.RestartStep()) + { + write_cylinder_displacement.writeToFile(number_of_iterations); + } + if (number_of_iterations % restart_output_interval == 0) + restart_io.writeToFile(number_of_iterations); + } + number_of_iterations++; + + /** Update cell linked list and configuration. */ + time_instance = TickCount::now(); + water_block.updateCellLinkedListWithParticleSort(100); + cylinder.updateCellLinkedList(); + water_block_inner.updateConfiguration(); + cylinder_inner.updateConfiguration(); + cylinder_contact.updateConfiguration(); + water_block_complex.updateConfiguration(); + fluid_observer_contact.updateConfiguration(); + free_stream_surface_indicator.exec(); + 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"; + + if (sph_system.generate_regression_data_) + { + write_cylinder_displacement.generateDataBase(1.0e-3); + } + else if (sph_system.RestartStep() == 0) + { + write_cylinder_displacement.testResult(); + } + + return 0; +};