diff --git a/docs/src/refs.bib b/docs/src/refs.bib index b514169ceb..c479e59210 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -939,7 +939,3 @@ @article{Zhu2021 doi = {10.1007/s42241-021-0031-y}, publisher = {Springer Science and Business Media LLC}, } - -@Comment{jabref-meta: databaseType:bibtex;} - -@Comment{jabref-meta: saveOrderConfig:specified;bibtexkey;false;date;false;abstract;false;} diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 5dab8cfc4a..a29ecd1e13 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -84,14 +84,7 @@ extra_callback = nothing callbacks = CallbackSet(info_callback, saving_callback, extra_callback) # Use a Runge-Kutta method with automatic (error based) time step size control. -# Limiting of the maximum stepsize is necessary to prevent crashing. -# When particles are approaching a wall in a uniform way, they can be advanced -# with large time steps. Close to the wall, the stepsize has to be reduced drastically. -# Sometimes, the method fails to do so because forces become extremely large when -# fluid particles are very close to boundary particles, and the time integration method -# interprets this as an instability. sol = solve(ode, RDPK3SpFSAL35(), - abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) - reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=1e-2, # Limit stepsize to prevent crashing + abstol=1e-8, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 save_everystep=false, callback=callbacks); diff --git a/examples/fluid/vortex_street.jl b/examples/fluid/vortex_street.jl new file mode 100644 index 0000000000..fc33e9b994 --- /dev/null +++ b/examples/fluid/vortex_street.jl @@ -0,0 +1,193 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.005 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +# Boundary geometry and initial fluid particle positions +tank_size = (1.0, 0.5) +# tank_size = (0.6, 0.6) +initial_fluid_size = tank_size + +Re = 10000 +initial_velocity = (0.0, 0.0) +nu = 1 * 1 / Re + +tspan = (0.0, 2.0) + +fluid_density = 1000.0 +sound_speed = 50.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(false, false, true, true), velocity=initial_velocity) + +center = tank_size ./ 2 +hollow_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, + n_layers=4, sphere_type=RoundSphere()) + +filled_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, + sphere_type=RoundSphere()) + +# hollow_sphere = RectangularShape(fluid_particle_spacing, round.(Int, (0.1, 0.3) ./ fluid_particle_spacing), center .- 0.15, +# density=fluid_density) +# filled_sphere = hollow_sphere + +# using PythonCall +# using CondaPkg + +# CondaPkg.add("svgpathtools") +# CondaPkg.add("ezdxf") +# svgpath = pyimport("svgpathtools") + +# svg_path = "M509.299 100.016C507.966 91.5871 505.915 85.7111 503.145 82.3879C498.991 77.4031 479.883 60.7871 475.521 58.2947C471.16 55.8023 455.167 49.1559 448.936 47.702C442.705 46.2481 355.471 30.463 339.893 28.8014C329.508 27.6937 287.761 23.5397 214.651 16.3394C199.727 15.7768 185.488 15.1656 171.934 14.5056C151.602 13.5157 106.318 5.19544 82.4982 0.830064C58.6781 -3.53532 3.26262 9.37214 0.279574 38.2664C-2.54326 65.6089 16.5085 89.4186 34.9345 99.2843C49.9801 107.34 58.7166 107.403 71.5338 114.275C101.912 130.564 100.849 169.8 123.353 175.939C145.856 182.078 146.637 180.9 155.986 179.279C162.22 178.199 199.267 168.32 267.129 149.644L359.355 117.398L405.801 102.863C446.849 101.008 472.412 100.061 482.489 100.022C497.605 99.9635 507.524 100.062 509.299 100.016Z" + +# path = svgpath.parse_path(svg_path) +# t_range = range(0, 1, length=50 * length(path)) +# points = [(pyconvert(Float64, p.real), -pyconvert(Float64, p.imag)) +# for p in (path.point(t) for t in t_range)] + +# # ezdxf = pyimport("ezdxf") +# # doc = ezdxf.new(dxfversion="R2010") +# # msp = doc.modelspace() +# # msp.add_polyline2d(points) +# # doc.saveas("output.dxf") + +# center = tank_size ./ 2 +# points_matrix = reinterpret(reshape, Float64, points) +# scaling = 0.3 / maximum(points_matrix, dims=2)[1] +# points_matrix .*= scaling +# points_matrix .+= (-0.15, -points_matrix[2, 1]) + +# geometry = TrixiParticles.Polygon(points_matrix) + +# # trixi2vtk(geometry) + +# point_in_geometry_algorithm = WindingNumberJacobson(; geometry, +# winding_number_factor=0.4, +# hierarchical_winding=true) + +# # Returns `InitialCondition` +# shape_sampled = ComplexShape(geometry; particle_spacing=fluid_particle_spacing, density=fluid_density, +# store_winding_number=true, +# point_in_geometry_algorithm) + +# # angle = pi / 4 +# # using StaticArrays +# # rotation_matrix = @SMatrix [cos(angle) -sin(angle) +# # sin(angle) cos(angle)] +# # shape_sampled.initial_condition.coordinates .= rotation_matrix * shape_sampled.initial_condition.coordinates +# shape_sampled.initial_condition.coordinates .+= center + +# hollow_sphere = shape_sampled.initial_condition +# filled_sphere = hollow_sphere + +# n_particles = round(Int, 0.12 / fluid_particle_spacing) +# cylinder = RectangularShape(fluid_particle_spacing, (n_particles, n_particles), (0.2 - 1.0, 0.24 - 1.0), density=fluid_density) + +fluid = setdiff(tank.fluid, filled_sphere) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 2 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +viscosity_fluid = ViscosityAdami(; nu) +# viscosity_fluid = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +viscosity_wall = ViscosityAdami(; nu) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_density_calculator = ContinuityDensity() +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + density_diffusion=density_diffusion, + pressure_acceleration=TrixiParticles.tensile_instability_control, + particle_shifting=TrixiParticles.ParticleShiftingSun2019(5.0), + smoothing_length, viscosity=viscosity_fluid) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation(pressure_offset=10000.0) +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length) + +# Movement function +# https://en.wikipedia.org/wiki/Triangle_wave#Harmonics +# triangle_series(t, N) = 8 / pi^2 * sum(i -> (-1)^i / (2i + 1)^2 * sin(2pi * (2i + 1) * t), 0:(N-1)) +# movement_function(x, t) = x + SVector(0.4 * triangle_series(2 * t, 5), 0.0) +# is_moving(t) = true +# boundary_movement = BoundaryMovement(movement_function, is_moving) + +boundary_movement = TrixiParticles.oscillating_movement(1.0, + SVector(0.4, 0.0), + 0.0, center; + ramp_up=0.5) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model_cylinder = BoundaryModelDummyParticles(hollow_sphere.density, + hollow_sphere.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) + +boundary_system_cylinder = BoundarySPHSystem(hollow_sphere, boundary_model_cylinder, + movement=boundary_movement) + +# boundary_system_cylinder = TotalLagrangianSPHSystem(hollow_sphere, smoothing_kernel, smoothing_length, +# 1e5, 0.0; +# n_fixed_particles=nparticles(hollow_sphere), movement=boundary_movement, +# boundary_model=boundary_model) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(fluid.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(fluid.coordinates, dims=2) .+ fluid_particle_spacing / 2 +periodic_box = PeriodicBox(; min_corner, max_corner) +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list) + +semi = Semidiscretization(fluid_system, boundary_system_cylinder; + neighborhood_search) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="") +shifting_callback = ParticleShiftingCallback() + +stepsize_callback = StepsizeCallback(cfl=1.5) + +callbacks = CallbackSet(info_callback, saving_callback, shifting_callback, + stepsize_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +# fluid_dt = 1e-3 +# sol = solve(ode, RDPK3SpFSAL49(), +# # abstol=1.0e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) +# # reltol=1.0e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) +# adaptive=false, dt=fluid_dt, +# save_everystep=false, callback=callbacks); + +time_step = 1.0 +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=time_step, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); + +# plane = interpolate_plane([0.0, -0.25], [1.0, 0.75], 0.0025, semi, fluid_system, sol) diff --git a/examples/fluid/vortex_street2.jl b/examples/fluid/vortex_street2.jl new file mode 100644 index 0000000000..8f89583693 --- /dev/null +++ b/examples/fluid/vortex_street2.jl @@ -0,0 +1,233 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.002 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +# Boundary geometry and initial fluid particle positions +tank_size = (1.0, 0.5) +# tank_size = (0.6, 0.6) +initial_fluid_size = tank_size + +initial_velocity = (1.0, 0.0) +nu = 1e-4 + +tspan = (0.0, 2.0) + +fluid_density = 1000.0 +sound_speed = 20.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(false, false, true, true), velocity=initial_velocity) + +center = tank_size ./ 2 +# hollow_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, +# n_layers=4, sphere_type=RoundSphere()) + +# filled_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, +# sphere_type=RoundSphere()) + +# fluid = setdiff(tank.fluid, filled_sphere) + + +# ========================================================================================== +# ==== Sample and pack foot pocket +using PythonCall +using CondaPkg + +CondaPkg.add("svgpathtools") +CondaPkg.add("ezdxf") +svgpath = pyimport("svgpathtools") + +svg_path = "M509.299 100.016C507.966 91.5871 505.915 85.7111 503.145 82.3879C498.991 77.4031 479.883 60.7871 475.521 58.2947C471.16 55.8023 455.167 49.1559 448.936 47.702C442.705 46.2481 355.471 30.463 339.893 28.8014C329.508 27.6937 287.761 23.5397 214.651 16.3394C199.727 15.7768 185.488 15.1656 171.934 14.5056C151.602 13.5157 106.318 5.19544 82.4982 0.830064C58.6781 -3.53532 3.26262 9.37214 0.279574 38.2664C-2.54326 65.6089 16.5085 89.4186 34.9345 99.2843C49.9801 107.34 58.7166 107.403 71.5338 114.275C101.912 130.564 100.849 169.8 123.353 175.939C145.856 182.078 146.637 180.9 155.986 179.279C162.22 178.199 199.267 168.32 267.129 149.644L359.355 117.398L405.801 102.863C446.849 101.008 472.412 100.061 482.489 100.022C497.605 99.9635 507.524 100.062 509.299 100.016Z" + +path = svgpath.parse_path(svg_path) +t_range = range(0, 1, length=50 * length(path)) +points = [(pyconvert(Float64, p.real), -pyconvert(Float64, p.imag)) + for p in (path.point(t) for t in t_range)] + +points_matrix = reinterpret(reshape, Float64, points) +scaling = 0.3 / maximum(points_matrix, dims=2)[1] +points_matrix .*= scaling +points_matrix .+= (-0.15, -points_matrix[2, 1]) .+ center .- (0.0, 1e-4) +# # Clamp the blade in one layer of particles by moving the foot down by a particle spacing +# points_matrix .-= (0.0, particle_spacing) +geometry = TrixiParticles.Polygon(points_matrix) + +point_in_geometry_algorithm = WindingNumberJacobson(; geometry, + winding_number_factor=0.4, + hierarchical_winding=true) + +# Returns `InitialCondition` +shape_sampled = ComplexShape(geometry; particle_spacing=fluid_particle_spacing, density=fluid_density, + point_in_geometry_algorithm) + +foot_sdf = SignedDistanceField(geometry, fluid_particle_spacing; + max_signed_distance=4 * fluid_particle_spacing, + use_for_boundary_packing=true) + +boundary_packing = sample_boundary(foot_sdf; boundary_density=fluid_density, + boundary_thickness=4 * fluid_particle_spacing) + +background_pressure = 1.0 +smoothing_length_packing = 0.8 * fluid_particle_spacing +foot_packing_system = ParticlePackingSystem(shape_sampled; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, background_pressure) + +fluid_packing_system = ParticlePackingSystem(boundary_packing; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, is_boundary=true, background_pressure, + boundary_compress_factor=0.8) + +min_corner = minimum(tank.fluid.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(tank.fluid.coordinates, dims=2) .+ fluid_particle_spacing / 2 +periodic_box = PeriodicBox(; min_corner, max_corner) +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list, update_strategy=ParallelUpdate()) + +semi_packing = Semidiscretization(foot_packing_system, fluid_packing_system; + neighborhood_search) + +ode_packing = semidiscretize(semi_packing, (0.0, 10.0)) + +sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing"), + UpdateCallback()), + dtmax=1e-2) + +hollow_sphere = InitialCondition(sol_packing, foot_packing_system, semi_packing) +fluid = setdiff(tank.fluid, hollow_sphere) + +fluid_packing_system = ParticlePackingSystem(fluid; smoothing_length=smoothing_length_packing, + signed_distance_field=nothing, background_pressure) + +boundary_packing_system = ParticlePackingSystem(hollow_sphere; smoothing_length=smoothing_length_packing, + fixed_system=true, signed_distance_field=nothing, background_pressure) + +semi_packing = Semidiscretization(fluid_packing_system, boundary_packing_system; + neighborhood_search) + +ode_packing = semidiscretize(semi_packing, (0.0, 2.0)) + +sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing"), + UpdateCallback()), + dtmax=1e-2) + +fluid = InitialCondition(sol_packing, fluid_packing_system, semi_packing) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 2 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +viscosity_fluid = ViscosityAdami(; nu) +# viscosity_fluid = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +viscosity_wall = ViscosityAdami(; nu) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_density_calculator = ContinuityDensity() +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + density_diffusion=density_diffusion, + pressure_acceleration=TrixiParticles.tensile_instability_control, + particle_shifting=TrixiParticles.ParticleShiftingSun2019(2.0), + smoothing_length, viscosity=viscosity_fluid) + +# ========================================================================================== +# ==== Boundary +# Movement function +frequency = 1.3 # Hz +amplitude = 0.18 # m +rotation_deg = 25 # degrees +rotation_phase_offset = 0.12 # periods +translation_vector = SVector(0.0, amplitude) +rotation_angle = rotation_deg * pi / 180 + +boundary_movement = TrixiParticles.oscillating_movement(frequency, + SVector(0.0, amplitude), + rotation_angle, center; + rotation_phase_offset, ramp_up=0.5) + +boundary_density_calculator = BernoulliPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model_cylinder = BoundaryModelDummyParticles(hollow_sphere.density, + hollow_sphere.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) + +# boundary_system_cylinder = BoundarySPHSystem(hollow_sphere, boundary_model_cylinder, +# movement=boundary_movement) + +boundary_system_cylinder = TotalLagrangianSPHSystem(hollow_sphere, smoothing_kernel, smoothing_length, + 1e9, 0.3; + n_fixed_particles=nparticles(hollow_sphere), movement=boundary_movement, + boundary_model=boundary_model_cylinder, + penalty_force=PenaltyForceGanzenmueller(alpha=0.1)) + +# boundary_system_cylinder = TotalLagrangianSPHSystem(hollow_sphere, smoothing_kernel, smoothing_length, +# 1e5, 0.0; +# n_fixed_particles=nparticles(hollow_sphere), movement=boundary_movement, +# boundary_model=boundary_model) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(tank.fluid.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(tank.fluid.coordinates, dims=2) .+ fluid_particle_spacing / 2 +periodic_box = PeriodicBox(; min_corner, max_corner) +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list) + +semi = Semidiscretization(fluid_system, boundary_system_cylinder; + neighborhood_search) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="") +shifting_callback = ParticleShiftingCallback() + +stepsize_callback = StepsizeCallback(cfl=1.5) + +callbacks = CallbackSet(info_callback, saving_callback, shifting_callback, + stepsize_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +# fluid_dt = 1e-3 +# sol = solve(ode, RDPK3SpFSAL49(), +# # abstol=1.0e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) +# # reltol=1.0e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) +# adaptive=false, dt=fluid_dt, +# save_everystep=false, callback=callbacks); + +time_step = 1.0 +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=time_step, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); + +# plane = interpolate_plane([0.0, -0.25], [1.0, 0.75], 0.0025, semi, fluid_system, sol) diff --git a/examples/fluid/vortex_street_postprocess.jl b/examples/fluid/vortex_street_postprocess.jl new file mode 100644 index 0000000000..155d592907 --- /dev/null +++ b/examples/fluid/vortex_street_postprocess.jl @@ -0,0 +1,168 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +const fluid_particle_spacing = 0.2 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +# Boundary geometry and initial fluid particle positions +tank_size = (65.0, 20.0) +initial_fluid_size = tank_size + +Re = 200 +initial_velocity = (1.0, 0.0) +nu = 1 * 1 / Re + +strouhal_number = 0.198 * (1 - 19.7 / Re) +frequency = strouhal_number * initial_velocity[1] / 1 + +tspan = (0.0, 50.0) + +fluid_density = 1000.0 +sound_speed = 10initial_velocity[1] +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(false, false, true, true), velocity=initial_velocity) + +const radius = 2.0 +const center = SVector(5.0, 10.0) +hollow_sphere = SphereShape(fluid_particle_spacing, radius, Tuple(center), fluid_density, + n_layers=4, sphere_type=RoundSphere()) + +filled_sphere = SphereShape(fluid_particle_spacing, radius, Tuple(center), fluid_density, + sphere_type=RoundSphere()) + +# n_particles = round(Int, 0.12 / fluid_particle_spacing) +# cylinder = RectangularShape(fluid_particle_spacing, (n_particles, n_particles), (0.2 - 1.0, 0.24 - 1.0), density=fluid_density) + +fluid = setdiff(tank.fluid, filled_sphere) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 2 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +viscosity = ViscosityAdami(; nu) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_density_calculator = ContinuityDensity() +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + density_diffusion=density_diffusion, + pressure_acceleration=TrixiParticles.tensile_instability_control, + # transport_velocity=TransportVelocityAdami(50_000.0), + smoothing_length, viscosity=viscosity) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model_cylinder = BoundaryModelDummyParticles(hollow_sphere.density, + hollow_sphere.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system_cylinder = BoundarySPHSystem(hollow_sphere, boundary_model_cylinder) + +# ========================================================================================== +# ==== Simulation +periodic_box = PeriodicBox(min_corner=[0.0, -1.0], max_corner=[65.0, 21.0]) +cell_list = FullGridCellList(min_corner=[0.0, -1.0], max_corner=[65.0, 21.0]) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list) + +semi = Semidiscretization(fluid_system, boundary_system, boundary_system_cylinder; + neighborhood_search) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=1.0, prefix="tvf_old") +shifting_callback = ParticleShiftingCallback() + +# ========================================================================================== +# ==== Postprocessing +circle = SphereShape(fluid_particle_spacing, (2 * radius + fluid_particle_spacing) / 2, + Tuple(center), fluid_density, n_layers=1, + sphere_type=RoundSphere()) + +# Points for pressure interpolation, located at the wall interface +const data_points = copy(circle.coordinates) + +calculate_lift_force(system, v_ode, u_ode, semi, t) = nothing +function calculate_lift_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) + force = zero(SVector{ndims(system), eltype(system)}) + + values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, + clip_negative_pressure=false) + pressure = Array(values.pressure) + + for i in axes(data_points, 2) + point = TrixiParticles.current_coords(data_points, system, i) + + # F = ∑ -p_i * A_i * n_i + force -= pressure[i] * fluid_particle_spacing .* + TrixiParticles.normalize(point - center) + end + + return 2 * force[2] / (fluid_density * 1^2 * 2 * radius) +end + +calculate_drag_force(system, v_ode, u_ode, semi, t) = nothing +function calculate_drag_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) + force = zero(SVector{ndims(system), eltype(system)}) + + values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, + clip_negative_pressure=false) + pressure = Array(values.pressure) + + for i in axes(data_points, 2) + point = TrixiParticles.current_coords(data_points, system, i) + + # F = ∑ -p_i * A_i * n_i + force -= pressure[i] * fluid_particle_spacing .* + TrixiParticles.normalize(point - center) + end + + return 2 * force[1] / (fluid_density * 1^2 * 2 * radius) +end + +pp_callback = PostprocessCallback(; dt=0.5, + f_l=calculate_lift_force, f_d=calculate_drag_force, + filename="resulting_force_pst", + write_csv=true, write_file_interval=10) + +callbacks = CallbackSet(info_callback, saving_callback, shifting_callback, + UpdateCallback(), pp_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL49(), + abstol=1.0e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1.0e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + save_everystep=false, callback=callbacks); + +# sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), +# dt=1.0, # This is overwritten by the stepsize callback +# save_everystep=false, callback=callbacks); + +# plane = interpolate_plane([0.0, -0.25], [1.0, 0.75], 0.0025, semi, fluid_system, sol) diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl index 44d86c5317..4b1bae4084 100644 --- a/examples/fsi/dam_break_plate_2d.jl +++ b/examples/fsi/dam_break_plate_2d.jl @@ -39,7 +39,8 @@ state_equation = StateEquationCole(; sound_speed, reference_density=fluid_densit tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, - acceleration=(0.0, -gravity), state_equation=state_equation) + acceleration=(0.0, -gravity), state_equation=state_equation, + coordinates_eltype=Float64) # Elastic plate/beam length_beam = 0.08 @@ -65,10 +66,12 @@ non_fixed_position = (plate_position[1], plate = RectangularShape(structure_particle_spacing, (n_particles_x, n_particles_y - 1), non_fixed_position, - density=structure_density, place_on_shell=true) + density=structure_density, place_on_shell=true, + coordinates_eltype=Float64) clamped_particles = RectangularShape(structure_particle_spacing, (n_particles_x, 1), plate_position, - density=structure_density, place_on_shell=true) + density=structure_density, place_on_shell=true, + coordinates_eltype=Float64) structure = union(clamped_particles, plate) diff --git a/examples/fsi/fin_2d.jl b/examples/fsi/fin_2d.jl new file mode 100644 index 0000000000..bc3a925ce8 --- /dev/null +++ b/examples/fsi/fin_2d.jl @@ -0,0 +1,408 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK +using OrdinaryDiffEqSymplecticRK + +function convert_ic(ic, T) + return InitialCondition{ndims(ic)}(ic.coordinates, ic.velocity, ic.mass, ic.density, + ic.pressure, T(ic.particle_spacing)) +end + +# ========================================================================================== +# ==== Resolution +n_particles_y = 4 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +fin_length = 0.6 +fin_thickness = 30e-3 +real_thickness = 1e-3 +real_modulus = 125e9 +poisson_ratio = 0.3 +flexural_rigidity = real_modulus * real_thickness^3 / (1 - poisson_ratio^2) / 12 +modulus = 12 * (1 - poisson_ratio^2) * flexural_rigidity / fin_thickness^3 + +fiber_volume_fraction = 0.6 +fiber_density = 1800.0 +epoxy_density = 1250.0 +density = fiber_volume_fraction * fiber_density + + (1 - fiber_volume_fraction) * epoxy_density + +clamp_radius = 0.05 + +tank_size = (2.0, 1.0) +center = (tank_size[2] / 2, tank_size[2] / 2) +initial_fluid_size = tank_size +initial_velocity = (1.0, 0.0) + +# The structure starts at the position of the first particle and ends +# at the position of the last particle. +particle_spacing = fin_thickness / (n_particles_y - 1) +fluid_particle_spacing = particle_spacing + +smoothing_length_structure = sqrt(2) * particle_spacing +smoothing_length_fluid = 1.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +file = joinpath(examples_dir(), "preprocessing", "data", "fin.dxf") +geometry = load_geometry(file) + +# trixi2vtk(geometry) + +point_in_geometry_algorithm = WindingNumberJacobson(; geometry, + winding_number_factor=0.4, + hierarchical_winding=true) + +# Returns `InitialCondition` +shape_sampled = ComplexShape(geometry; particle_spacing, density=density, + grid_offset=center, point_in_geometry_algorithm) +shape_sampled = TrixiParticles.@set shape_sampled.coordinates = Float64.(shape_sampled.coordinates) + +# Beam and clamped particles +length_clamp = round(Int, 0.15 / particle_spacing) * particle_spacing # m +n_particles_per_dimension = (round(Int, (fin_length + length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x, + n_particles_y) + +# Note that the `RectangularShape` puts the first particle half a particle spacing away +# from the boundary, which is correct for fluids, but not for structures. +# We therefore need to pass `place_on_shell=true`. +beam = RectangularShape(particle_spacing, n_particles_per_dimension, + (-length_clamp, 0.0), density=density, place_on_shell=true) + +fixed_particles = setdiff(shape_sampled, beam) + +# structure = union(beam, fixed_particles) + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 10 + +fluid_density = 1000.0 +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, + faces=(false, false, true, true), velocity=initial_velocity) +# fluid = setdiff(tank.fluid, structure) + +open_boundary_size = (fluid_particle_spacing * open_boundary_layers, tank_size[2]) + +min_coords_inlet = (-open_boundary_layers * fluid_particle_spacing, 0.0) +inlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_inlet, + velocity=initial_velocity, + faces=(false, false, true, true)) + +min_coords_outlet = (tank.fluid_size[1], 0.0) +outlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_outlet, + velocity=initial_velocity, + faces=(false, false, true, true)) + + +NDIMS = ndims(tank.fluid) +n_buffer_particles = 20 * tank.n_particles_per_dimension[2]^(NDIMS - 1) + +# ========================================================================================== +# ==== Packing +packing = false +if packing + foot_sdf = SignedDistanceField(geometry, particle_spacing; + max_signed_distance=4 * particle_spacing, + use_for_boundary_packing=true) + + boundary_packing = sample_boundary(foot_sdf; boundary_density=density, + boundary_thickness=4 * particle_spacing) + boundary_packing = TrixiParticles.@set boundary_packing.coordinates = Float64.(boundary_packing.coordinates) + boundary_packing = setdiff(boundary_packing, beam) + + background_pressure = 1.0 + smoothing_length_packing = 0.8 * particle_spacing + foot_packing_system = ParticlePackingSystem(fixed_particles; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, background_pressure) + + fluid_packing_system = ParticlePackingSystem(boundary_packing; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, is_boundary=true, background_pressure, + boundary_compress_factor=0.8) + + blade_packing_system = ParticlePackingSystem(beam; smoothing_length=smoothing_length_packing, + fixed_system=true, signed_distance_field=nothing, background_pressure) + + min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 + periodic_box = PeriodicBox(; min_corner, max_corner) + cell_list = FullGridCellList(; min_corner, max_corner) + neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list, update_strategy=ParallelUpdate()) + + semi_packing = Semidiscretization(foot_packing_system, fluid_packing_system, + blade_packing_system; neighborhood_search) + + ode_packing = semidiscretize(semi_packing, (0.0, 10.0)) + + sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing"), + UpdateCallback()), + dtmax=1e-2) + + packed_foot = InitialCondition(sol_packing, foot_packing_system, semi_packing) + + # Move the fin to the center of the tank + packed_foot.coordinates .+= center + beam.coordinates .+= center + + structure = union(beam, packed_foot) + fluid = setdiff(tank.fluid, structure) + + # Pack the fluid against the fin and the tank boundary + pack_window = TrixiParticles.Polygon(stack([ + [0.15, 0.42], + [0.3, 0.42], + [0.44, 0.48], + [1.12, 0.48], + [1.12, 0.52], + [0.55, 0.52], + [0.5, 0.56], + [0.24, 0.6], + [0.15, 0.6], + [0.15, 0.42] + ])) + + # Then, we extract the particles that fall inside this window + pack_fluid = intersect(fluid, pack_window) + # and those outside the window + fixed_fluid = setdiff(fluid, pack_fluid) + fixed_union = union(fixed_fluid, structure) + + fluid_packing_system = ParticlePackingSystem(pack_fluid; smoothing_length=smoothing_length_packing, + signed_distance_field=nothing, background_pressure) + + fixed_packing_system = ParticlePackingSystem(fixed_union; smoothing_length=smoothing_length_packing, + fixed_system=true, signed_distance_field=nothing, background_pressure) + + semi_packing = Semidiscretization(fluid_packing_system, fixed_packing_system; + neighborhood_search) + + ode_packing = semidiscretize(semi_packing, (0.0, 2.0)) + + sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing"), + UpdateCallback()), + dtmax=1e-2) + + fluid = InitialCondition(sol_packing, fluid_packing_system, semi_packing) + fluid = union(fluid, fixed_fluid) +else + structure = union(beam, fixed_particles) + # Move the fin to the center of the tank + structure.coordinates .+= center + + fluid = setdiff(tank.fluid, structure) +end + +n_clamped_particles = nparticles(structure) - nparticles(beam) + +# Movement function +frequency = 1.3 # Hz +amplitude = 0.18 # m +rotation_deg = 25 # degrees +rotation_phase_offset = 0.12 # periods +translation_vector = SVector(0.0, amplitude) +rotation_angle = rotation_deg * pi / 180 + +boundary_motion = OscillatingMotion2D(; frequency, + translation_vector=SVector(0.0, amplitude), + rotation_angle, rotation_center=center, + rotation_phase_offset, ramp_up_tspan=(0.0, 0.5)) + +sound_speed = 40.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +# ========================================================================================== +# ==== Structure +boundary_density_calculator = AdamiPressureExtrapolation() +viscosity_fluid = ViscosityAdami(nu=1e-4) +viscosity_fin = ViscosityAdami(nu=1e-4) + +# For the FSI we need the hydrodynamic masses and densities in the structure boundary model +hydrodynamic_densites = fluid_density * ones(size(structure.density)) +hydrodynamic_masses = hydrodynamic_densites * particle_spacing^2 + +boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densites, + hydrodynamic_masses, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid, + viscosity=viscosity_fin) + +# k_structure = 1.0 +# beta_structure = fluid_particle_spacing / particle_spacing +# boundary_model_structure = BoundaryModelMonaghanKajtar(k_structure, beta_structure, +# particle_spacing, +# hydrodynamic_masses) + +structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length_structure, + modulus, poisson_ratio; + n_clamped_particles, clamped_particles_motion=boundary_motion, + boundary_model=boundary_model_structure, + velocity_averaging=TrixiParticles.VelocityAveraging(5e-4), + viscosity=ArtificialViscosityMonaghan(alpha=0.1), + penalty_force=PenaltyForceGanzenmueller(alpha=0.1)) + +# ========================================================================================== +# ==== Fluid +fluid_density_calculator = ContinuityDensity() +# density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length_fluid, viscosity=viscosity_fluid, + density_diffusion=density_diffusion, + shifting_technique=ParticleShiftingTechnique(sound_speed_factor=0.2, v_max_factor=0.0), + pressure_acceleration=tensile_instability_control, + buffer_size=n_buffer_particles) +# fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, +# sound_speed, viscosity=ViscosityAdami(; nu), +# transport_velocity=TransportVelocityAdami(10 * sound_speed^2 * fluid_density)) + +# ========================================================================================== +# ==== Open Boundaries +periodic = false +if periodic + min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 + min_corner = convert.(typeof(fluid_particle_spacing), min_corner) + max_corner = convert.(typeof(fluid_particle_spacing), max_corner) + periodic_box = PeriodicBox(; min_corner, max_corner) + open_boundary_system = nothing + wall = tank.boundary +else + periodic_box = nothing + + open_boundary_model = BoundaryModelDynamicalPressureZhang() + # open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring()) + reference_velocity_in = SVector(1.0, 0.0) + reference_pressure_in = 0.0 + reference_density_in = nothing + boundary_type_in = InFlow() + face_in = ([0.0, 0.0], [0.0, tank_size[2]]) + flow_direction = [1.0, 0.0] + inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing, + reference_density=reference_density_in, + reference_pressure=reference_pressure_in, + reference_velocity=reference_velocity_in, + initial_condition=inlet.fluid, boundary_type=boundary_type_in) + + reference_velocity_out = SVector(1.0, 0.0) + reference_pressure_out = nothing + reference_density_out = nothing + boundary_type_out = OutFlow() + face_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], tank_size[2]]) + outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction), + open_boundary_layers, density=fluid_density, particle_spacing, + reference_density=reference_density_out, + reference_pressure=reference_pressure_out, + reference_velocity=reference_velocity_out, + initial_condition=outlet.fluid, boundary_type=boundary_type_out) + + open_boundary_system = OpenBoundarySystem(inflow, outflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles) + + wall = union(tank.boundary, inlet.boundary, outlet.boundary) + min_corner = minimum(wall.coordinates, dims=2) .- 5 * fluid_particle_spacing + max_corner = maximum(wall.coordinates, dims=2) .+ 5 * fluid_particle_spacing +end + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid) + +boundary_system = WallBoundarySystem(wall, boundary_model) + +# ========================================================================================== +# ==== Simulation +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list, + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, boundary_system, open_boundary_system, structure_system; neighborhood_search, + parallelization_backend=PolyesterBackend()) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +prefix = "" +saving_callback = SolutionSavingCallback(dt=0.01; prefix) + +split_cfl = 1.6 +# SSPRK104 CFL = 2.5, 15k RHS evaluations +# CarpenterKennedy2N54 CFL = 1.6, 11k RHS evaluations +# RK4 CFL = 1.2, 12k RHS evaluations +# VerletLeapfrog CFL = 0.5, 6.75k RHS evaluations +# VelocityVerlet CFL = 0.5, 6.75k RHS evaluations +# DPRKN4 CFL = 1.7, 9k RHS evaluations + +# function tip_velocity(system::TotalLagrangianSPHSystem, data, t) +# return data.velocity[2254] +# end +# pp_tip = PostprocessCallback(; tip_velocity, interval=1, +# filename="$(prefix)_tip_velocity", write_file_interval=10_000) +split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), adaptive=false, + stage_coupling=true, + dt=1e-5, # This is overwritten by the stepsize callback + callback=StepsizeCallback(cfl=split_cfl), + maxiters=10^8) + +fluid_cfl = 1.2 +stepsize_callback = StepsizeCallback(cfl=fluid_cfl) + +function total_volume(system::WeaklyCompressibleSPHSystem, data, t) + return sum(data.mass ./ data.density) +end +function total_volume(system, data, t) + return nothing +end +pp_cb = PostprocessCallback(; total_volume, interval=100, + filename="$(prefix)_total_volume", write_file_interval=50) + +function plane_vtk(system, dv_ode, du_ode, v_ode, u_ode, semi, t) + return nothing +end +function plane_vtk(system::WeaklyCompressibleSPHSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) + resolution = fluid_particle_spacing / 6 + pvd = TrixiParticles.paraview_collection("out/$(prefix)_plane"; append=t > 0) + interpolate_plane_2d_vtk(min_corner, max_corner, resolution, + semi, semi.systems[1], v_ode, u_ode, include_wall_velocity=true, + filename="$(prefix)_plane_$(round(Int, t * 1000))", pvd=pvd, t=t) + TrixiParticles.vtk_save(pvd) + return nothing +end +interpolate_cb = PostprocessCallback(; plane_vtk, dt=0.01, filename="plane") + +callbacks = CallbackSet(info_callback, saving_callback, + stepsize_callback, split_integration, pp_cb, interpolate_cb, + UpdateCallback()) + +dt_fluid = 1.25e-4 +sol = solve(ode, + # RDPK3SpFSAL35(), + CarpenterKennedy2N54(williamson_condition=false), + dt=dt_fluid, # This is overwritten by the stepsize callback + # reltol=1e-5, abstol=1e-7, + save_everystep=false, callback=callbacks, maxiters=10^8); diff --git a/examples/fsi/fin_2d_simplified.jl b/examples/fsi/fin_2d_simplified.jl new file mode 100644 index 0000000000..2234fe1df8 --- /dev/null +++ b/examples/fsi/fin_2d_simplified.jl @@ -0,0 +1,253 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK +using OrdinaryDiffEqSymplecticRK + +# ========================================================================================== +# ==== Resolution +n_particles_y = 4 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +fin_length = 0.6 +fin_thickness = 2e-3 +real_thickness = 1e-3 +real_modulus = 125e9 +poisson_ratio = 0.3 +flexural_rigidity = real_modulus * real_thickness^3 / (1 - poisson_ratio^2) / 12 +modulus = 12 * (1 - poisson_ratio^2) * flexural_rigidity / fin_thickness^3 + +fiber_volume_fraction = 0.6 +fiber_density = 1800.0 +epoxy_density = 1250.0 +density = fiber_volume_fraction * fiber_density + + (1 - fiber_volume_fraction) * epoxy_density + +clamp_radius = 0.05 + +tank_size = (0.8, 0.2) +center = (0.05, 0.1) +initial_fluid_size = tank_size +initial_velocity = (1.0, 0.0) + +# The structure starts at the position of the first particle and ends +# at the position of the last particle. +particle_spacing = fin_thickness / (n_particles_y - 1) +fluid_particle_spacing = particle_spacing + +smoothing_length_structure = sqrt(2) * particle_spacing +smoothing_length_fluid = 1.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +# Beam and clamped particles +length_clamp = round(Int, 0.01 / particle_spacing) * particle_spacing # m +n_particles_per_dimension = (round(Int, (length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x, + n_particles_y) +shape_sampled = RectangularShape(particle_spacing, n_particles_per_dimension, + (-length_clamp, 0.0), density=density, place_on_shell=true) +length_clamp = 0.0 +n_particles_per_dimension = (round(Int, (fin_length + length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x, + n_particles_y) + +# Note that the `RectangularShape` puts the first particle half a particle spacing away +# from the boundary, which is correct for fluids, but not for structures. +# We therefore need to pass `place_on_shell=true`. +beam = RectangularShape(particle_spacing, n_particles_per_dimension, + (-length_clamp, 0.0), density=density, place_on_shell=true) + +fixed_particles = setdiff(shape_sampled, beam) + +# structure = union(beam, fixed_particles) + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 6 + +fluid_density = 1000.0 +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, + faces=(false, false, true, true), velocity=initial_velocity) +# fluid = setdiff(tank.fluid, structure) + +open_boundary_size = (fluid_particle_spacing * open_boundary_layers, tank_size[2]) + +min_coords_inlet = (-open_boundary_layers * fluid_particle_spacing, 0.0) +inlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_inlet, + faces=(false, false, true, true)) + +min_coords_outlet = (tank.fluid_size[1], 0.0) +outlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_outlet, + faces=(false, false, true, true)) + + +NDIMS = ndims(tank.fluid) +n_buffer_particles = 10 * tank.n_particles_per_dimension[2]^(NDIMS - 1) + + +structure = union(beam, fixed_particles) +# Move the fin to the center of the tank +structure.coordinates .+= center .+ (fluid_particle_spacing / 2, fluid_particle_spacing / 2) + +fluid = setdiff(tank.fluid, structure) + +n_clamped_particles = nparticles(structure) - nparticles(beam) + +# Movement function +frequency = 1.3 # Hz +amplitude = 0.18 # m +rotation_deg = 25 # degrees +rotation_phase_offset = 0.12 # periods +translation_vector = SVector(0.0, amplitude) +rotation_angle = rotation_deg * pi / 180 + +boundary_motion = OscillatingMotion2D(; frequency, + translation_vector=SVector(0.0, amplitude), + rotation_angle, rotation_center=center, + rotation_phase_offset, ramp_up_tspan=(0.0, 0.5)) + +sound_speed = 40.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +# ========================================================================================== +# ==== Structure +boundary_density_calculator = AdamiPressureExtrapolation() +viscosity_fluid = ViscosityAdami(nu=1e-4) +viscosity_fin = ViscosityAdami(nu=1e-4) + +# For the FSI we need the hydrodynamic masses and densities in the structure boundary model +hydrodynamic_densites = fluid_density * ones(size(structure.density)) +hydrodynamic_masses = hydrodynamic_densites * particle_spacing^2 + +boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densites, + hydrodynamic_masses, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid, + viscosity=viscosity_fin) + +# k_structure = 1.0 +# beta_structure = fluid_particle_spacing / particle_spacing +# boundary_model_structure = BoundaryModelMonaghanKajtar(k_structure, beta_structure, +# particle_spacing, +# hydrodynamic_masses) + +structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length_structure, + modulus, poisson_ratio; + n_clamped_particles, #clamped_particles_motion=boundary_motion, + velocity_averaging=nothing, + boundary_model=boundary_model_structure, + viscosity=ArtificialViscosityMonaghan(alpha=0.1), + penalty_force=PenaltyForceGanzenmueller(alpha=0.1)) + +# ========================================================================================== +# ==== Fluid +fluid_density_calculator = ContinuityDensity() +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +# density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length_fluid, viscosity=viscosity_fluid, + density_diffusion=density_diffusion, + shifting_technique=ParticleShiftingTechnique(sound_speed_factor=0.2, v_max_factor=0.0), + pressure_acceleration=tensile_instability_control, + buffer_size=n_buffer_particles) +# fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, +# sound_speed, viscosity=ViscosityAdami(; nu), +# transport_velocity=TransportVelocityAdami(10 * sound_speed^2 * fluid_density)) + +min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 +min_corner = convert.(typeof(fluid_particle_spacing), min_corner) +max_corner = convert.(typeof(fluid_particle_spacing), max_corner) +periodic_box = PeriodicBox(; min_corner, max_corner) +open_boundary_system = nothing +wall = tank.boundary + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid) + +boundary_system = WallBoundarySystem(wall, boundary_model) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(wall.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(wall.coordinates, dims=2) .+ fluid_particle_spacing / 2 +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list, + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, boundary_system, open_boundary_system, structure_system; neighborhood_search, + parallelization_backend=PolyesterBackend()) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +prefix = "simplified" +saving_callback = SolutionSavingCallback(dt=0.01, prefix=prefix) + +split_cfl = 1.6 +# SSPRK104 CFL = 2.5, 15k RHS evaluations +# CarpenterKennedy2N54 CFL = 1.6, 11k RHS evaluations +# RK4 CFL = 1.2, 12k RHS evaluations +# VerletLeapfrog CFL = 0.5, 6.75k RHS evaluations +# VelocityVerlet CFL = 0.5, 6.75k RHS evaluations +# DPRKN4 CFL = 1.7, 9k RHS evaluations + +split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), adaptive=false, + stage_coupling=true, + dt=1e-5, # This is overwritten by the stepsize callback + callback=StepsizeCallback(cfl=split_cfl), + maxiters=10^8) + +fluid_cfl = 1.2 +stepsize_callback = StepsizeCallback(cfl=fluid_cfl) + +function total_volume(system::WeaklyCompressibleSPHSystem, data, t) + return sum(data.mass ./ data.density) +end +function total_volume(system, data, t) + return nothing +end +pp_cb = PostprocessCallback(; total_volume, interval=100, + filename="$(prefix)_total_volume", write_file_interval=50) + +function plane_vtk(system, dv_ode, du_ode, v_ode, u_ode, semi, t) + return nothing +end +function plane_vtk(system::WeaklyCompressibleSPHSystem, dv_ode, du_ode, v_ode, u_ode, semi_, t) + resolution = fluid_particle_spacing / 6 + pvd = TrixiParticles.paraview_collection("out/$(prefix)_plane"; append=t > 0) + interpolate_plane_2d_vtk(min_corner, max_corner, resolution, + semi_, semi_.systems[1], v_ode, u_ode, include_wall_velocity=true, + filename="$(prefix)_plane_$(round(Int, t * 1000))", pvd=pvd, t=t) + TrixiParticles.vtk_save(pvd) + return nothing +end +interpolate_cb = PostprocessCallback(; plane_vtk, dt=0.01, filename="plane") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), + stepsize_callback, split_integration, pp_cb, interpolate_cb) + +dt_fluid = 1.25e-4 +sol = solve(ode, + # RDPK3SpFSAL35(), + CarpenterKennedy2N54(williamson_condition=false), + dt=dt_fluid, # This is overwritten by the stepsize callback + # reltol=1e-5, abstol=1e-7, + save_everystep=false, callback=callbacks, maxiters=10^8); diff --git a/examples/n_body/n_body_benchmark_trixi.jl b/examples/n_body/n_body_benchmark_trixi.jl index 828f046bf2..a9e0843328 100644 --- a/examples/n_body/n_body_benchmark_trixi.jl +++ b/examples/n_body/n_body_benchmark_trixi.jl @@ -85,15 +85,16 @@ function symplectic_euler!(velocity, coordinates, semi) u = vec(coordinates) dv = copy(v) du = copy(u) + p = (; semi, split_integration_data=nothing) @time for _ in 1:50_000_000 - TrixiParticles.kick!(dv, v, u, semi, 0.0) + TrixiParticles.kick!(dv, v, u, p, 0.0) @inbounds for i in eachindex(v) v[i] += 0.01 * dv[i] end - TrixiParticles.drift!(du, v, u, semi, 0.0) + TrixiParticles.drift!(du, v, u, p, 0.0) @inbounds for i in eachindex(u) u[i] += 0.01 * du[i] diff --git a/examples/preprocessing/data/fin.dxf b/examples/preprocessing/data/fin.dxf new file mode 100644 index 0000000000..3f56143d5a --- /dev/null +++ b/examples/preprocessing/data/fin.dxf @@ -0,0 +1,6758 @@ + 0 +SECTION + 2 +HEADER + 9 +$ACADVER + 1 +AC1024 + 9 +$ACADMAINTVER + 70 +6 + 9 +$DWGCODEPAGE + 3 +ANSI_1252 + 9 +$LASTSAVEDBY + 1 +ezdxf + 9 +$INSBASE + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$EXTMIN + 10 +-100 + 20 +-100 + 30 +-100 + 9 +$EXTMAX + 10 +100 + 20 +100 + 30 +100 + 9 +$LIMMIN + 10 +0.0 + 20 +0.0 + 9 +$LIMMAX + 10 +420.0 + 20 +297.0 + 9 +$ORTHOMODE + 70 +0 + 9 +$REGENMODE + 70 +1 + 9 +$FILLMODE + 70 +1 + 9 +$QTEXTMODE + 70 +0 + 9 +$MIRRTEXT + 70 +1 + 9 +$LTSCALE + 40 +1.0 + 9 +$ATTMODE + 70 +1 + 9 +$TEXTSIZE + 40 +2.5 + 9 +$TRACEWID + 40 +1.0 + 9 +$TEXTSTYLE + 7 +Standard + 9 +$CLAYER + 8 +0 + 9 +$CELTYPE + 6 +ByLayer + 9 +$CECOLOR + 62 +256 + 9 +$CELTSCALE + 40 +1.0 + 9 +$DISPSILH + 70 +0 + 9 +$DIMSCALE + 40 +1.0 + 9 +$DIMASZ + 40 +2.5 + 9 +$DIMEXO + 40 +0.625 + 9 +$DIMDLI + 40 +3.75 + 9 +$DIMRND + 40 +0.0 + 9 +$DIMDLE + 40 +0.0 + 9 +$DIMEXE + 40 +1.25 + 9 +$DIMTP + 40 +0.0 + 9 +$DIMTM + 40 +0.0 + 9 +$DIMTXT + 40 +2.5 + 9 +$DIMCEN + 40 +2.5 + 9 +$DIMTSZ + 40 +0.0 + 9 +$DIMTOL + 70 +0 + 9 +$DIMLIM + 70 +0 + 9 +$DIMTIH + 70 +0 + 9 +$DIMTOH + 70 +0 + 9 +$DIMSE1 + 70 +0 + 9 +$DIMSE2 + 70 +0 + 9 +$DIMTAD + 70 +1 + 9 +$DIMZIN + 70 +8 + 9 +$DIMBLK + 1 + + 9 +$DIMASO + 70 +1 + 9 +$DIMSHO + 70 +1 + 9 +$DIMPOST + 1 + + 9 +$DIMAPOST + 1 + + 9 +$DIMALT + 70 +0 + 9 +$DIMALTD + 70 +3 + 9 +$DIMALTF + 40 +0.03937007874 + 9 +$DIMLFAC + 40 +1.0 + 9 +$DIMTOFL + 70 +1 + 9 +$DIMTVP + 40 +0.0 + 9 +$DIMTIX + 70 +0 + 9 +$DIMSOXD + 70 +0 + 9 +$DIMSAH + 70 +0 + 9 +$DIMBLK1 + 1 + + 9 +$DIMBLK2 + 1 + + 9 +$DIMSTYLE + 2 +ISO-25 + 9 +$DIMCLRD + 70 +0 + 9 +$DIMCLRE + 70 +0 + 9 +$DIMCLRT + 70 +0 + 9 +$DIMTFAC + 40 +1.0 + 9 +$DIMGAP + 40 +0.625 + 9 +$DIMJUST + 70 +0 + 9 +$DIMSD1 + 70 +0 + 9 +$DIMSD2 + 70 +0 + 9 +$DIMTOLJ + 70 +0 + 9 +$DIMTZIN + 70 +8 + 9 +$DIMALTZ + 70 +0 + 9 +$DIMALTTZ + 70 +0 + 9 +$DIMUPT + 70 +0 + 9 +$DIMDEC + 70 +2 + 9 +$DIMTDEC + 70 +2 + 9 +$DIMALTU + 70 +2 + 9 +$DIMALTTD + 70 +3 + 9 +$DIMTXSTY + 7 +Standard + 9 +$DIMAUNIT + 70 +0 + 9 +$DIMADEC + 70 +0 + 9 +$DIMALTRND + 40 +0.0 + 9 +$DIMAZIN + 70 +0 + 9 +$DIMDSEP + 70 +44 + 9 +$DIMATFIT + 70 +3 + 9 +$DIMFRAC + 70 +0 + 9 +$DIMLDRBLK + 1 + + 9 +$DIMLUNIT + 70 +2 + 9 +$DIMLWD + 70 +-2 + 9 +$DIMLWE + 70 +-2 + 9 +$DIMTMOVE + 70 +0 + 9 +$DIMFXL + 40 +1.0 + 9 +$DIMFXLON + 70 +0 + 9 +$DIMJOGANG + 40 +0.785398163397 + 9 +$DIMTFILL + 70 +0 + 9 +$DIMTFILLCLR + 70 +0 + 9 +$DIMARCSYM + 70 +0 + 9 +$DIMLTYPE + 6 + + 9 +$DIMLTEX1 + 6 + + 9 +$DIMLTEX2 + 6 + + 9 +$DIMTXTDIRECTION + 70 +0 + 9 +$LUNITS + 70 +2 + 9 +$LUPREC + 70 +4 + 9 +$SKETCHINC + 40 +1.0 + 9 +$FILLETRAD + 40 +10.0 + 9 +$AUNITS + 70 +0 + 9 +$AUPREC + 70 +2 + 9 +$MENU + 1 +. + 9 +$ELEVATION + 40 +0.0 + 9 +$PELEVATION + 40 +0.0 + 9 +$THICKNESS + 40 +0.0 + 9 +$LIMCHECK + 70 +0 + 9 +$CHAMFERA + 40 +0.0 + 9 +$CHAMFERB + 40 +0.0 + 9 +$CHAMFERC + 40 +0.0 + 9 +$CHAMFERD + 40 +0.0 + 9 +$SKPOLY + 70 +0 + 9 +$TDCREATE + 40 +2461046.8383217594 + 9 +$TDUCREATE + 40 +2458532.153996898 + 9 +$TDUPDATE + 40 +2461046.8383217594 + 9 +$TDUUPDATE + 40 +2458532.1544311 + 9 +$TDINDWG + 40 +0.0 + 9 +$TDUSRTIMER + 40 +0.0 + 9 +$USRTIMER + 70 +1 + 9 +$ANGBASE + 50 +0.0 + 9 +$ANGDIR + 70 +0 + 9 +$PDMODE + 70 +0 + 9 +$PDSIZE + 40 +0.0 + 9 +$PLINEWID + 40 +0.0 + 9 +$SPLFRAME + 70 +0 + 9 +$SPLINETYPE + 70 +6 + 9 +$SPLINESEGS + 70 +8 + 9 +$HANDSEED + 5 +DA + 9 +$SURFTAB1 + 70 +6 + 9 +$SURFTAB2 + 70 +6 + 9 +$SURFTYPE + 70 +6 + 9 +$SURFU + 70 +6 + 9 +$SURFV + 70 +6 + 9 +$UCSBASE + 2 + + 9 +$UCSNAME + 2 + + 9 +$UCSORG + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSXDIR + 10 +1.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSYDIR + 10 +0.0 + 20 +1.0 + 30 +0.0 + 9 +$UCSORTHOREF + 2 + + 9 +$UCSORTHOVIEW + 70 +0 + 9 +$UCSORGTOP + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGBOTTOM + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGLEFT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGRIGHT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGFRONT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGBACK + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSBASE + 2 + + 9 +$PUCSNAME + 2 + + 9 +$PUCSORG + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSXDIR + 10 +1.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSYDIR + 10 +0.0 + 20 +1.0 + 30 +0.0 + 9 +$PUCSORTHOREF + 2 + + 9 +$PUCSORTHOVIEW + 70 +0 + 9 +$PUCSORGTOP + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGBOTTOM + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGLEFT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGRIGHT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGFRONT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGBACK + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$USERI1 + 70 +0 + 9 +$USERI2 + 70 +0 + 9 +$USERI3 + 70 +0 + 9 +$USERI4 + 70 +0 + 9 +$USERI5 + 70 +0 + 9 +$USERR1 + 40 +0.0 + 9 +$USERR2 + 40 +0.0 + 9 +$USERR3 + 40 +0.0 + 9 +$USERR4 + 40 +0.0 + 9 +$USERR5 + 40 +0.0 + 9 +$WORLDVIEW + 70 +1 + 9 +$SHADEDGE + 70 +3 + 9 +$SHADEDIF + 70 +70 + 9 +$TILEMODE + 70 +1 + 9 +$MAXACTVP + 70 +64 + 9 +$PINSBASE + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PLIMCHECK + 70 +0 + 9 +$PEXTMIN + 10 +1e+20 + 20 +1e+20 + 30 +1e+20 + 9 +$PEXTMAX + 10 +-1e+20 + 20 +-1e+20 + 30 +-1e+20 + 9 +$PLIMMIN + 10 +0.0 + 20 +0.0 + 9 +$PLIMMAX + 10 +420.0 + 20 +297.0 + 9 +$UNITMODE + 70 +0 + 9 +$VISRETAIN + 70 +1 + 9 +$PLINEGEN + 70 +0 + 9 +$PSLTSCALE + 70 +1 + 9 +$TREEDEPTH + 70 +3020 + 9 +$CMLSTYLE + 2 +Standard + 9 +$CMLJUST + 70 +0 + 9 +$CMLSCALE + 40 +20.0 + 9 +$PROXYGRAPHICS + 70 +1 + 9 +$MEASUREMENT + 70 +1 + 9 +$CELWEIGHT +370 +-1 + 9 +$ENDCAPS +280 +0 + 9 +$JOINSTYLE +280 +0 + 9 +$LWDISPLAY +290 +0 + 9 +$INSUNITS + 70 +6 + 9 +$HYPERLINKBASE + 1 + + 9 +$STYLESHEET + 1 + + 9 +$XEDIT +290 +1 + 9 +$CEPSNTYPE +380 +0 + 9 +$PSTYLEMODE +290 +1 + 9 +$FINGERPRINTGUID + 2 +BD1C87EE-EA69-11F0-BD04-A088C2115B2C + 9 +$VERSIONGUID + 2 +BD1D174A-EA69-11F0-BD04-A088C2115B2C + 9 +$EXTNAMES +290 +1 + 9 +$PSVPSCALE + 40 +0.0 + 9 +$OLESTARTUP +290 +0 + 9 +$SORTENTS +280 +127 + 9 +$INDEXCTL +280 +0 + 9 +$HIDETEXT +280 +1 + 9 +$XCLIPFRAME +280 +2 + 9 +$HALOGAP +280 +0 + 9 +$OBSCOLOR + 70 +257 + 9 +$OBSLTYPE +280 +0 + 9 +$INTERSECTIONDISPLAY +280 +0 + 9 +$INTERSECTIONCOLOR + 70 +257 + 9 +$DIMASSOC +280 +2 + 9 +$PROJECTNAME + 1 + + 9 +$CAMERADISPLAY +290 +0 + 9 +$LENSLENGTH + 40 +50.0 + 9 +$CAMERAHEIGHT + 40 +0.0 + 9 +$STEPSPERSEC + 40 +24.0 + 9 +$STEPSIZE + 40 +100.0 + 9 +$3DDWFPREC + 40 +2.0 + 9 +$PSOLWIDTH + 40 +0.005 + 9 +$PSOLHEIGHT + 40 +0.08 + 9 +$LOFTANG1 + 40 +1.570796326795 + 9 +$LOFTANG2 + 40 +1.570796326795 + 9 +$LOFTMAG1 + 40 +0.0 + 9 +$LOFTMAG2 + 40 +0.0 + 9 +$LOFTPARAM + 70 +7 + 9 +$LOFTNORMALS +280 +1 + 9 +$LATITUDE + 40 +37.795 + 9 +$LONGITUDE + 40 +-122.394 + 9 +$NORTHDIRECTION + 40 +0.0 + 9 +$TIMEZONE + 70 +-8000 + 9 +$LIGHTGLYPHDISPLAY +280 +1 + 9 +$TILEMODELIGHTSYNCH +280 +1 + 9 +$CMATERIAL +347 +20 + 9 +$SOLIDHIST +280 +0 + 9 +$SHOWHIST +280 +1 + 9 +$DWFFRAME +280 +2 + 9 +$DGNFRAME +280 +2 + 9 +$REALWORLDSCALE +290 +1 + 9 +$INTERFERECOLOR + 62 +256 + 9 +$CSHADOW +280 +0 + 9 +$SHADOWPLANELOCATION + 40 +0.0 + 0 +ENDSEC + 0 +SECTION + 2 +CLASSES + 0 +CLASS + 1 +ACDBDICTIONARYWDFLT + 2 +AcDbDictionaryWithDefault + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +SUN + 2 +AcDbSun + 3 +SCENEOE + 90 +1153 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +VISUALSTYLE + 2 +AcDbVisualStyle + 3 +ObjectDBX Classes + 90 +4095 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +MATERIAL + 2 +AcDbMaterial + 3 +ObjectDBX Classes + 90 +1153 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +SCALE + 2 +AcDbScale + 3 +ObjectDBX Classes + 90 +1153 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +TABLESTYLE + 2 +AcDbTableStyle + 3 +ObjectDBX Classes + 90 +4095 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +MLEADERSTYLE + 2 +AcDbMLeaderStyle + 3 +ACDB_MLEADERSTYLE_CLASS + 90 +4095 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +DICTIONARYVAR + 2 +AcDbDictionaryVar + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +CELLSTYLEMAP + 2 +AcDbCellStyleMap + 3 +ObjectDBX Classes + 90 +1152 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +MENTALRAYRENDERSETTINGS + 2 +AcDbMentalRayRenderSettings + 3 +SCENEOE + 90 +1024 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +ACDBDETAILVIEWSTYLE + 2 +AcDbDetailViewStyle + 3 +ObjectDBX Classes + 90 +1025 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +ACDBSECTIONVIEWSTYLE + 2 +AcDbSectionViewStyle + 3 +ObjectDBX Classes + 90 +1025 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +RASTERVARIABLES + 2 +AcDbRasterVariables + 3 +ISM + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +LAYOUT + 2 +AcDbLayout + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +ACDBPLACEHOLDER + 2 +AcDbPlaceHolder + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +ENDSEC + 0 +SECTION + 2 +TABLES + 0 +TABLE + 2 +VPORT + 5 +8 +330 +0 +100 +AcDbSymbolTable + 70 +1 + 0 +VPORT + 5 +23 +330 +8 +100 +AcDbSymbolTableRecord +100 +AcDbViewportTableRecord + 2 +*Active + 70 +0 + 10 +0.0 + 20 +0.0 + 11 +1.0 + 21 +1.0 + 12 +70.0 + 22 +50.0 + 13 +0.0 + 23 +0.0 + 14 +0.5 + 24 +0.5 + 15 +0.5 + 25 +0.5 + 16 +0.0 + 26 +0.0 + 36 +1.0 + 17 +0.0 + 27 +0.0 + 37 +0.0 + 40 +1.0 + 41 +1.34 + 42 +50.0 + 43 +0.0 + 44 +0.0 + 50 +0.0 + 51 +0.0 + 71 +0 + 72 +1000 + 73 +1 + 74 +3 + 75 +0 + 76 +0 + 77 +0 + 78 +0 +281 +0 + 65 +0 +146 +0.0 + 0 +ENDTAB + 0 +TABLE + 2 +LTYPE + 5 +2 +330 +0 +100 +AcDbSymbolTable + 70 +3 + 0 +LTYPE + 5 +24 +330 +2 +100 +AcDbSymbolTableRecord +100 +AcDbLinetypeTableRecord + 2 +ByBlock + 70 +0 + 3 + + 72 +65 + 73 +0 + 40 +0.0 + 0 +LTYPE + 5 +25 +330 +2 +100 +AcDbSymbolTableRecord +100 +AcDbLinetypeTableRecord + 2 +ByLayer + 70 +0 + 3 + + 72 +65 + 73 +0 + 40 +0.0 + 0 +LTYPE + 5 +26 +330 +2 +100 +AcDbSymbolTableRecord +100 +AcDbLinetypeTableRecord + 2 +Continuous + 70 +0 + 3 + + 72 +65 + 73 +0 + 40 +0.0 + 0 +ENDTAB + 0 +TABLE + 2 +LAYER + 5 +1 +330 +0 +100 +AcDbSymbolTable + 70 +2 + 0 +LAYER + 5 +27 +330 +1 +100 +AcDbSymbolTableRecord +100 +AcDbLayerTableRecord + 2 +0 + 70 +0 + 62 +7 + 6 +Continuous +370 +-3 +390 +13 +347 +21 + 0 +LAYER + 5 +28 +330 +1 +100 +AcDbSymbolTableRecord +100 +AcDbLayerTableRecord + 2 +Defpoints + 70 +0 + 62 +7 + 6 +Continuous +290 +0 +370 +-3 +390 +13 +347 +21 + 0 +ENDTAB + 0 +TABLE + 2 +STYLE + 5 +5 +330 +0 +100 +AcDbSymbolTable + 70 +1 + 0 +STYLE + 5 +29 +330 +5 +100 +AcDbSymbolTableRecord +100 +AcDbTextStyleTableRecord + 2 +Standard + 70 +0 + 40 +0.0 + 41 +1.0 + 50 +0.0 + 71 +0 + 42 +2.5 + 3 +txt + 4 + + 0 +ENDTAB + 0 +TABLE + 2 +VIEW + 5 +7 +330 +0 +100 +AcDbSymbolTable + 70 +0 + 0 +ENDTAB + 0 +TABLE + 2 +UCS + 5 +6 +330 +0 +100 +AcDbSymbolTable + 70 +0 + 0 +ENDTAB + 0 +TABLE + 2 +APPID + 5 +3 +330 +0 +100 +AcDbSymbolTable + 70 +2 + 0 +APPID + 5 +2A +330 +3 +100 +AcDbSymbolTableRecord +100 +AcDbRegAppTableRecord + 2 +ACAD + 70 +0 + 0 +APPID + 5 +D9 +330 +3 +100 +AcDbSymbolTableRecord +100 +AcDbRegAppTableRecord + 2 +HATCHBACKGROUNDCOLOR + 70 +0 + 0 +ENDTAB + 0 +TABLE + 2 +DIMSTYLE + 5 +4 +330 +0 +100 +AcDbSymbolTable + 70 +1 +100 +AcDbDimStyleTable + 0 +DIMSTYLE +105 +2B +330 +4 +100 +AcDbSymbolTableRecord +100 +AcDbDimStyleTableRecord + 2 +Standard + 70 +0 + 40 +1.0 + 41 +2.5 + 42 +0.625 + 43 +3.75 + 44 +1.25 + 45 +0.0 + 46 +0.0 + 47 +0.0 + 48 +0.0 + 49 +2.5 +140 +2.5 +141 +2.5 +142 +0.0 +143 +0.03937007874 +144 +1.0 +145 +0.0 +146 +1.0 +147 +0.625 +148 +0.0 + 69 +0 + 70 +0 + 71 +0 + 72 +0 + 73 +0 + 74 +0 + 75 +0 + 76 +0 + 77 +1 + 78 +8 + 79 +0 +170 +0 +171 +3 +172 +1 +173 +0 +174 +0 +175 +0 +176 +0 +177 +0 +178 +0 +179 +0 +271 +0 +272 +2 +273 +2 +274 +3 +275 +0 +276 +0 +277 +2 +278 +44 +279 +0 +280 +0 +281 +0 +282 +0 +283 +0 +284 +8 +285 +0 +286 +0 +288 +0 +289 +3 +290 +0 +371 +-2 +372 +-2 + 0 +ENDTAB + 0 +TABLE + 2 +BLOCK_RECORD + 5 +9 +330 +0 +100 +AcDbSymbolTable + 70 +2 + 0 +BLOCK_RECORD + 5 +17 +330 +9 +100 +AcDbSymbolTableRecord +100 +AcDbBlockTableRecord + 2 +*Model_Space +340 +1A + 70 +0 +280 +1 +281 +0 + 0 +BLOCK_RECORD + 5 +1B +330 +9 +100 +AcDbSymbolTableRecord +100 +AcDbBlockTableRecord + 2 +*Paper_Space +340 +1E + 70 +0 +280 +1 +281 +0 + 0 +ENDTAB + 0 +ENDSEC + 0 +SECTION + 2 +BLOCKS + 0 +BLOCK + 5 +18 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbBlockBegin + 2 +*Model_Space + 70 +0 + 10 +0.0 + 20 +0.0 + 30 +0.0 + 3 +*Model_Space + 1 + + 0 +ENDBLK + 5 +19 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbBlockEnd + 0 +BLOCK + 5 +1C +330 +1B +100 +AcDbEntity + 8 +0 +100 +AcDbBlockBegin + 2 +*Paper_Space + 70 +0 + 10 +0.0 + 20 +0.0 + 30 +0.0 + 3 +*Paper_Space + 1 + + 0 +ENDBLK + 5 +1D +330 +1B +100 +AcDbEntity + 8 +0 +100 +AcDbBlockEnd + 0 +ENDSEC + 0 +SECTION + 2 +ENTITIES + 0 +POLYLINE + 5 +2D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDb2dPolyline + 66 +1 + 10 +0.0 + 20 +0.0 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +2F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +0.0 + 20 +5.551115123125783e-17 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +30 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0010647730582502013 + 20 +-2.8770736736127844e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +31 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0034971934090743395 + 20 +-4.953574980059994e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +32 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.007027977637432459 + 20 +-6.326524879191053e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +33 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.011387842328284381 + 20 +-7.092944330866491e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +34 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.01630750406658993 + 20 +-7.34985429493018e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +35 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.02151767943730898 + 20 +-7.19427573124265e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +36 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.026749085025401576 + 20 +-6.723229599658875e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +37 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.031732437415827486 + 20 +-6.0337368600282826e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +38 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03619845319354664 + 20 +-5.222818472205848e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +39 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03987784894351887 + 20 +-4.387495396040997e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04378446456114665 + 20 +-0.0008783586002438781 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.048103522430103185 + 20 +-0.0025049262890193824 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.051718721680916624 + 20 +-0.003948651160244676 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0549385925607315 + 20 +-0.0052322536218777915 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.05807166531669217 + 20 +-0.006378454081876428 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.061532460806013545 + 20 +-0.007673859629111657 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +40 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06516986230184571 + 20 +-0.009171395348674838 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +41 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06880726379767793 + 20 +-0.010668931068237963 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +42 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07244466529351018 + 20 +-0.0121664667878012 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +43 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07608206678934243 + 20 +-0.013664002507364326 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +44 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0797194682851746 + 20 +-0.015161538226927451 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +45 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08335686978100681 + 20 +-0.016659073946490688 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +46 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08697507316482567 + 20 +-0.01820113484491409 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +47 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09056097837690974 + 20 +-0.019818103206756144 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +48 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09414688358899384 + 20 +-0.02143507156859814 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +49 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09773278880107797 + 20 +-0.023052039930440138 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1013186940131621 + 20 +-0.02466900829228219 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10490459922524614 + 20 +-0.026285976654124188 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10849050443733027 + 20 +-0.027902945015966185 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1120764096494144 + 20 +-0.029519913377808238 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11566231486149847 + 20 +-0.031136881739650235 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11924822007358254 + 20 +-0.03275385010149223 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +50 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12283412528566659 + 20 +-0.03437081846333434 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +51 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12642003049775075 + 20 +-0.03598778682517634 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +52 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1300059357098348 + 20 +-0.037604755187018335 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +53 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1335918409219189 + 20 +-0.03922172354886039 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +54 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13717774613400305 + 20 +-0.040838691910702385 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +55 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14372856522610838 + 20 +-0.04327496838266798 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +56 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14995755911305256 + 20 +-0.04558405899092027 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +57 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.15582399389417462 + 20 +-0.04775445959307795 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +58 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.16132787111901833 + 20 +-0.04978616888891935 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +59 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.16646919233712754 + 20 +-0.05167918557822265 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.17124795909804597 + 20 +-0.053433508360766435 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1756641729513177 + 20 +-0.05504913593632882 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.17971783544648637 + 20 +-0.05652606700468815 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18340894813309583 + 20 +-0.05786430026562289 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18673751256069004 + 20 +-0.05906383441891133 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18970353027881282 + 20 +-0.06012466816433182 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +60 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1923070028370079 + 20 +-0.06104680020166264 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +61 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19454793178481922 + 20 +-0.0618302292306821 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +62 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19642631867179056 + 20 +-0.06247495395116859 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +63 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1979421650474658 + 20 +-0.06298097306290046 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +64 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1990954724613888 + 20 +-0.063348285265656 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +65 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2000151253740754 + 20 +-0.06361126741164846 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +66 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20276874829217467 + 20 +-0.06435484737989106 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +67 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.205246815645514 + 20 +-0.06479549956536129 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +68 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20836029970179235 + 20 +-0.06480866119589457 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +69 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21302017272870882 + 20 +-0.06426976949932639 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21934980512311386 + 20 +-0.06318850973792922 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2223415227029405 + 20 +-0.062085340964722235 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22495973231109506 + 20 +-0.060292609672437714 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2273012549555563 + 20 +-0.057927912877090926 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22946291164430282 + 20 +-0.05510884759469742 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2315415233853133 + 20 +-0.05195301084127257 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +70 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2336339111865664 + 20 +-0.04857799963283188 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +71 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23583689605604086 + 20 +-0.04510141098539078 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +72 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23824729900171532 + 20 +-0.04164084191496459 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +73 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24096194103156843 + 20 +-0.038313889437569026 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +74 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24407764315357894 + 20 +-0.03523815056921925 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +75 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24769122637572547 + 20 +-0.032531222325930864 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +76 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.25189527154487523 + 20 +-0.03031250748991371 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +77 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2556710587823801 + 20 +-0.028897726461885542 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +78 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.259176350329518 + 20 +-0.027895338568229233 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +79 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.26261507934256245 + 20 +-0.027062875912520212 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2661911789777869 + 20 +-0.026157870598334243 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2701085823914646 + 20 +-0.02493785472924681 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2742767322269781 + 20 +-0.023287611479191006 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2773537617614356 + 20 +-0.02171288824969564 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28041119952074683 + 20 +-0.019761313469147057 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.283386801947182 + 20 +-0.01745533373208813 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +80 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28621832548301135 + 20 +-0.014817395633061614 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +81 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28884352657050505 + 20 +-0.011869945766610546 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +82 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29120016165193335 + 20 +-0.008635430727277682 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +83 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29322598716956644 + 20 +-0.005136297109605947 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +84 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2948587595656746 + 20 +-0.0013949915081380992 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +85 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.296036235282528 + 20 +0.0025660394825827715 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +86 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2966961707623968 + 20 +0.006724349268014074 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +87 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29677131237060905 + 20 +0.010823895979494569 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +88 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29606911246228473 + 20 +0.014181500952964832 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +89 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2945903225789705 + 20 +0.01733619705562711 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29243224340788126 + 20 +0.020277968983556938 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28969217563623184 + 20 +0.022996801432829728 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28646741995123726 + 20 +0.025482679099520955 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2828552770401123 + 20 +0.02772558667970615 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.278953047590072 + 20 +0.029715508869460727 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.27485803228833117 + 20 +0.03144243036486016 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +90 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.27066753182210485 + 20 +0.032896335861979986 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +91 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2664788468786078 + 20 +0.034067210056895614 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +92 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.26238927814505497 + 20 +0.034945037645682575 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +93 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2584961263086613 + 20 +0.03551980332441629 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +94 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2548966920566418 + 20 +0.035781491789172226 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +95 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2516882760762112 + 20 +0.03572008773602592 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +96 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24842500237321694 + 20 +0.03541667366752854 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +97 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24477200340806995 + 20 +0.035075636354958084 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +98 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24082262719541409 + 20 +0.03471373073340123 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +99 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23664618716419447 + 20 +0.034341520250685376 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23231199674335595 + 20 +0.03396956835463805 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22788936936184365 + 20 +0.033608438493086756 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22344761844860245 + 20 +0.03326869411385891 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21905605743257744 + 20 +0.03296089866478208 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21478399974271362 + 20 +0.03269561559368367 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2107007588079559 + 20 +0.03248340834839114 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20687564805724937 + 20 +0.03233484037673212 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2033779809195389 + 20 +0.032260475126533905 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20027707082376967 + 20 +0.03227087604562412 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19698654858711234 + 20 +0.03238592413156666 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19318285904709753 + 20 +0.0325402622832367 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18931970887041638 + 20 +0.03270409348936071 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18539709781429586 + 20 +0.032877414975028385 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18141502563596268 + 20 +0.033060223965329594 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.17737349209264375 + 20 +0.03325251768535398 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A9 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1726381225467503 + 20 +0.03338911516149823 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AA +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1659840849214246 + 20 +0.033316250989470886 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AB +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.15964888467818186 + 20 +0.03324643511774267 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AC +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1536325212730734 + 20 +0.033179667593903006 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AD +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14793499416215058 + 20 +0.03311594846554122 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AE +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14255630280146458 + 20 +0.03305527778024664 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AF +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13749644664706695 + 20 +0.032997655585608576 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13275542515500882 + 20 +0.03294308192921647 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12833323778134154 + 20 +0.032891556858659576 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12422988398211648 + 20 +0.032843080421527227 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12044536321338492 + 20 +0.03279765266540885 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11697967493119824 + 20 +0.03275527363789377 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1138328185916076 + 20 +0.032715943386571245 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11100479365066446 + 20 +0.03267966195903066 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10849559956442015 + 20 +0.03264642940286139 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1063052357889259 + 20 +0.03261624576565281 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B9 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10443370178023303 + 20 +0.03258911109499418 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BA +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10288099699439296 + 20 +0.03256502543847484 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BB +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10164712088745684 + 20 +0.03254398884368426 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BC +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09994337892822275 + 20 +0.03249160766593884 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BD +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09722703434661115 + 20 +0.03234615182366568 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BE +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.093771592897088 + 20 +0.032118852135510856 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BF +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08970115316225025 + 20 +0.03182178565127325 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0851398137246947 + 20 +0.03146702942075169 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08021167316701835 + 20 +0.031066660493744958 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07504083007181808 + 20 +0.03063275592005199 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06975138302169073 + 20 +0.03017739274947162 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06446743059923335 + 20 +0.029712648031802624 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0593130713870427 + 20 +0.029250598816843887 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.05441240396771574 + 20 +0.028803322154394295 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04988952692384935 + 20 +0.028382895094252625 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.045868538838040485 + 20 +0.028001394686217707 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C9 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04247353829288597 + 20 +0.027670897980088482 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CA +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03982862387098274 + 20 +0.027403482025663728 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CB +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03805789415492772 + 20 +0.02721122387274233 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CC +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.036264793274075446 + 20 +0.0269216897066466 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CD +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.032322371975447184 + 20 +0.025941611713292345 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CE +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.027830560327799014 + 20 +0.0246016350498387 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CF +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.023812359863461197 + 20 +0.02324602396599168 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.021294698693667236 + 20 +0.022218126153409845 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.01872618634040374 + 20 +0.02064227908050753 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.015093256261121235 + 20 +0.018146755728090103 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.011096281688223297 + 20 +0.01523911360159308 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.007435635854113387 + 20 +0.012426910206452146 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.004811691991195022 + 20 +0.010217703048102877 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.003023760328472236 + 20 +0.008011242975078714 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0014008764041681387 + 20 +0.004592179313590128 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +0.0 + 20 +0.0 + 30 +0.0 + 70 +0 + 0 +SEQEND + 5 +2E +330 +17 +100 +AcDbEntity + 8 +0 + 0 +ENDSEC + 0 +SECTION + 2 +OBJECTS + 0 +DICTIONARY + 5 +A +330 +0 +100 +AcDbDictionary +281 +1 + 3 +ACAD_COLOR +350 +B + 3 +ACAD_GROUP +350 +C + 3 +ACAD_LAYOUT +350 +D + 3 +ACAD_MATERIAL +350 +E + 3 +ACAD_MLEADERSTYLE +350 +F + 3 +ACAD_MLINESTYLE +350 +10 + 3 +ACAD_PLOTSETTINGS +350 +11 + 3 +ACAD_PLOTSTYLENAME +350 +12 + 3 +ACAD_SCALELIST +350 +14 + 3 +ACAD_TABLESTYLE +350 +15 + 3 +ACAD_VISUALSTYLE +350 +16 + 0 +DICTIONARY + 5 +B +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +C +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +D +330 +A +100 +AcDbDictionary +281 +1 + 3 +Model +350 +1A + 3 +Layout1 +350 +1E + 0 +DICTIONARY + 5 +E +330 +A +100 +AcDbDictionary +281 +1 + 3 +ByBlock +350 +1F + 3 +ByLayer +350 +20 + 3 +Global +350 +21 + 0 +DICTIONARY + 5 +F +330 +A +100 +AcDbDictionary +281 +1 + 3 +Standard +350 +2C + 0 +DICTIONARY + 5 +10 +330 +A +100 +AcDbDictionary +281 +1 + 3 +Standard +350 +22 + 0 +DICTIONARY + 5 +11 +330 +A +100 +AcDbDictionary +281 +1 + 0 +ACDBDICTIONARYWDFLT + 5 +12 +330 +A +100 +AcDbDictionary +281 +1 + 3 +Normal +350 +13 +100 +AcDbDictionaryWithDefault +340 +13 + 0 +ACDBPLACEHOLDER + 5 +13 +330 +12 + 0 +DICTIONARY + 5 +14 +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +15 +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +16 +330 +A +100 +AcDbDictionary +281 +1 + 0 +LAYOUT + 5 +1A +330 +D +100 +AcDbPlotSettings + 1 + + 2 +Adobe PDF + 4 +A4 + 6 + + 40 +3.175 + 41 +3.175 + 42 +3.175 + 43 +3.175 + 44 +209.91 + 45 +297.03 + 46 +0.0 + 47 +0.0 + 48 +0.0 + 49 +0.0 +140 +0.0 +141 +0.0 +142 +1.0 +143 +1.0 + 70 +1024 + 72 +0 + 73 +1 + 74 +5 + 7 + + 75 +16 + 76 +0 + 77 +2 + 78 +300 +147 +1.0 +148 +0.0 +149 +0.0 +100 +AcDbLayout + 1 +Model + 70 +1 + 71 +0 + 10 +-3.175 + 20 +-3.175 + 11 +293.857 + 21 +206.735 + 12 +0.0 + 22 +0.0 + 32 +0.0 + 14 +29.068 + 24 +20.356 + 34 +0.0 + 15 +261.614 + 25 +183.204 + 35 +0.0 +146 +0.0 + 13 +0.0 + 23 +0.0 + 33 +0.0 + 16 +1.0 + 26 +0.0 + 36 +0.0 + 17 +0.0 + 27 +1.0 + 37 +0.0 + 76 +1 +330 +17 + 0 +LAYOUT + 5 +1E +330 +D +100 +AcDbPlotSettings + 1 + + 2 +Adobe PDF + 4 +A4 + 6 + + 40 +3.175 + 41 +3.175 + 42 +3.175 + 43 +3.175 + 44 +209.91 + 45 +297.03 + 46 +0.0 + 47 +0.0 + 48 +0.0 + 49 +0.0 +140 +0.0 +141 +0.0 +142 +1.0 +143 +1.0 + 70 +0 + 72 +0 + 73 +1 + 74 +5 + 7 + + 75 +16 + 76 +0 + 77 +2 + 78 +300 +147 +1.0 +148 +0.0 +149 +0.0 +100 +AcDbLayout + 1 +Layout1 + 70 +1 + 71 +1 + 10 +-3.175 + 20 +-3.175 + 11 +293.857 + 21 +206.735 + 12 +0.0 + 22 +0.0 + 32 +0.0 + 14 +29.068 + 24 +20.356 + 34 +0.0 + 15 +261.614 + 25 +183.204 + 35 +0.0 +146 +0.0 + 13 +0.0 + 23 +0.0 + 33 +0.0 + 16 +1.0 + 26 +0.0 + 36 +0.0 + 17 +0.0 + 27 +1.0 + 37 +0.0 + 76 +1 +330 +1B + 0 +MATERIAL + 5 +1F +102 +{ACAD_REACTORS +330 +E +102 +} +330 +E +100 +AcDbMaterial + 1 +ByBlock + 2 + + 70 +0 + 40 +1.0 + 71 +1 + 41 +1.0 + 91 +-1023410177 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 44 +0.5 + 73 +0 + 45 +1.0 + 46 +1.0 + 77 +1 + 4 + + 78 +1 + 79 +1 +170 +1 + 48 +1.0 +171 +1 + 6 + +172 +1 +173 +1 +174 +1 +140 +1.0 +141 +1.0 +175 +1 + 7 + +176 +1 +177 +1 +178 +1 +143 +1.0 +179 +1 + 8 + +270 +1 +271 +1 +272 +1 +145 +1.0 +146 +1.0 +273 +1 + 9 + +274 +1 +275 +1 +276 +1 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 94 +63 + 0 +MATERIAL + 5 +20 +102 +{ACAD_REACTORS +330 +E +102 +} +330 +E +100 +AcDbMaterial + 1 +ByLayer + 2 + + 70 +0 + 40 +1.0 + 71 +1 + 41 +1.0 + 91 +-1023410177 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 44 +0.5 + 73 +0 + 45 +1.0 + 46 +1.0 + 77 +1 + 4 + + 78 +1 + 79 +1 +170 +1 + 48 +1.0 +171 +1 + 6 + +172 +1 +173 +1 +174 +1 +140 +1.0 +141 +1.0 +175 +1 + 7 + +176 +1 +177 +1 +178 +1 +143 +1.0 +179 +1 + 8 + +270 +1 +271 +1 +272 +1 +145 +1.0 +146 +1.0 +273 +1 + 9 + +274 +1 +275 +1 +276 +1 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 94 +63 + 0 +MATERIAL + 5 +21 +102 +{ACAD_REACTORS +330 +E +102 +} +330 +E +100 +AcDbMaterial + 1 +Global + 2 + + 70 +0 + 40 +1.0 + 71 +1 + 41 +1.0 + 91 +-1023410177 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 44 +0.5 + 73 +0 + 45 +1.0 + 46 +1.0 + 77 +1 + 4 + + 78 +1 + 79 +1 +170 +1 + 48 +1.0 +171 +1 + 6 + +172 +1 +173 +1 +174 +1 +140 +1.0 +141 +1.0 +175 +1 + 7 + +176 +1 +177 +1 +178 +1 +143 +1.0 +179 +1 + 8 + +270 +1 +271 +1 +272 +1 +145 +1.0 +146 +1.0 +273 +1 + 9 + +274 +1 +275 +1 +276 +1 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 94 +63 + 0 +MLINESTYLE + 5 +22 +102 +{ACAD_REACTORS +330 +10 +102 +} +330 +10 +100 +AcDbMlineStyle + 2 +Standard + 70 +0 + 3 + + 62 +256 + 51 +90.0 + 52 +90.0 + 71 +2 + 49 +0.5 + 62 +256 + 6 +BYLAYER + 49 +-0.5 + 62 +256 + 6 +BYLAYER + 0 +MLEADERSTYLE + 5 +2C +102 +{ACAD_REACTORS +330 +F +102 +} +330 +F +100 +AcDbMLeaderStyle +179 +2 +170 +2 +171 +1 +172 +0 + 90 +2 + 40 +0.0 + 41 +0.0 +173 +1 + 91 +-1056964608 + 92 +-2 +290 +1 + 42 +2.0 +291 +1 + 43 +8.0 + 3 +Standard + 44 +4.0 +300 + +342 +29 +174 +1 +175 +1 +176 +0 +178 +1 + 93 +-1056964608 + 45 +4.0 +292 +0 +297 +0 + 46 +4.0 + 94 +-1056964608 + 47 +1.0 + 49 +0.0 +140 +1.0 +294 +1 +141 +0.0 +177 +0 +142 +1.0 +295 +0 +296 +0 +143 +3.75 +271 +0 +272 +9 +272 +9 + 0 +ENDSEC + 0 +EOF diff --git a/ext/TrixiParticlesOrdinaryDiffEqExt.jl b/ext/TrixiParticlesOrdinaryDiffEqExt.jl index fa9eb5fe56..a8d257a7f3 100644 --- a/ext/TrixiParticlesOrdinaryDiffEqExt.jl +++ b/ext/TrixiParticlesOrdinaryDiffEqExt.jl @@ -105,7 +105,7 @@ end # The following is equivalent to `du = duprev + dt * kdu` for the velocity, but when # the density is integrated, a different update is used for the density. - semi = p + semi = p.semi TrixiParticles.foreach_system(semi) do system kdu_system = TrixiParticles.wrap_v(kdu, system, semi) du_system = TrixiParticles.wrap_v(du, system, semi) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bc3cbca3a2..d928d30608 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -72,7 +72,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, - SplitIntegrationCallback + SplitIntegrationCallback, EnergyCalculatorCallback, calculated_energy export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique, ParticleShiftingTechniqueSun2017, ConsistentShiftingSun2019, @@ -103,7 +103,7 @@ export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, export SourceTermDamping export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection -export nparticles +export nparticles, eachparticle export available_data, kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density export interpolate_line, interpolate_points, interpolate_plane_3d, interpolate_plane_2d, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 70dfef31a2..006fc9eaf7 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -71,3 +71,4 @@ include("stepsize.jl") include("update.jl") include("steady_state_reached.jl") include("split_integration.jl") +include("energy_calculator.jl") diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index ded47149f9..e1d565b154 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -67,7 +67,7 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr # Reinitialize initial solution if cb.reinit_initial_solution # Update systems to compute quantities like density and pressure. - semi = integrator.p + semi = integrator.p.semi v_ode, u_ode = u.x update_systems_and_nhs(v_ode, u_ode, semi, t) @@ -97,7 +97,7 @@ end # affect! function (reinit_callback::DensityReinitializationCallback)(integrator) vu_ode = integrator.u - semi = integrator.p + semi = integrator.p.semi @trixi_timeit timer() "reinit density" reinit_density!(vu_ode, semi) diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl new file mode 100644 index 0000000000..d463f4ec81 --- /dev/null +++ b/src/callbacks/energy_calculator.jl @@ -0,0 +1,216 @@ +""" + EnergyCalculatorCallback(system, semi; interval=1, + eachparticle=(n_integrated_particles(system) + 1):nparticles(system), + only_compute_force_on_fluid=false) + +Callback to integrate the work done by particles in a +[`TotalLagrangianSPHSystem`](@ref). + +With the default arguments it tracks the energy contribution of the clamped particles +that follow a [`PrescribedMotion`](@ref). By selecting different particles it can also be +used to measure the work done by the structure on the surrounding fluid. + +- **Prescribed/clamped motion energy** (default) -- monitor only the clamped particles by + leaving `eachparticle` at its default range + `(n_integrated_particles(system) + 1):nparticles(system)`. +- **Fluid load measurement** -- set `eachparticle=eachparticle(system)` together with + `only_compute_force_on_fluid=true` to accumulate the work that the entire structure + exerts on the surrounding fluid (useful for drag or lift estimates). + +Internally the callback integrates the instantaneous power, i.e. the dot product between +the force exerted by the particle and its prescribed velocity, using an explicit Euler +time integration scheme. + +The accumulated value can be retrieved via [`calculated_energy`](@ref). + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +# Arguments +- `system`: The [`TotalLagrangianSPHSystem`](@ref) whose particles should be monitored. +- `semi`: The [`Semidiscretization`](@ref) that contains `system`. + +# Keywords +- `interval=1`: Interval (in number of time steps) at which to compute the energy. + It is recommended to keep this at `1` (every time step) or small (≤ 5) + to limit time integration errors in the energy integral. +- `eachparticle=(n_integrated_particles(system) + 1):nparticles(system)`: Iterator + selecting which particles contribute. The default includes all clamped + particles in the system; pass `eachparticle(system)` to include every particle. +- `only_compute_force_on_fluid=false`: When `true`, only interactions with + fluid systems are accounted for. Combined with + `eachparticle=eachparticle(system)`, this accumulates the work that the + entire structure exerts on the fluid, which is useful for drag or lift + estimates. + +# Examples +```jldoctest; output = false, setup = :(system = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) +# Create an energy calculator callback that is called every 2 time steps +energy_cb = EnergyCalculatorCallback(system, semi; interval=2) + +# After the simulation, retrieve the calculated energy +total_energy = calculated_energy(energy_cb) + +# output +0.0 +``` +""" +struct EnergyCalculatorCallback{T, DV, EP} + interval :: Int + t :: T # Time of last call + energy :: T + system_index :: Int + dv :: DV + eachparticle :: EP + only_compute_force_on_fluid :: Bool +end + +function EnergyCalculatorCallback(system, semi; interval=1, + eachparticle=(n_integrated_particles(system) + 1):nparticles(system), + only_compute_force_on_fluid=false) + ELTYPE = eltype(system) + system_index = system_indices(system, semi) + + # Allocate buffer to write accelerations for all particles (including clamped ones) + dv = allocate(semi.parallelization_backend, ELTYPE, (ndims(system), nparticles(system))) + + cb = EnergyCalculatorCallback(interval, Ref(zero(ELTYPE)), Ref(zero(ELTYPE)), + system_index, dv, eachparticle, + only_compute_force_on_fluid) + + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(cb, cb, save_positions=(false, false)) +end + +""" + calculated_energy(cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + +Get the current calculated energy from an [`EnergyCalculatorCallback`](@ref). + +# Arguments +- `cb`: The `DiscreteCallback` returned by [`EnergyCalculatorCallback`](@ref). + +# Examples +```jldoctest; output = false, setup = :(system = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) +# Create an energy calculator callback +energy_cb = EnergyCalculatorCallback(system, semi) + +# After the simulation, retrieve the calculated energy +total_energy = calculated_energy(energy_cb) + +# output +0.0 +``` +""" +function calculated_energy(cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + return cb.affect!.energy[] +end + +# `condition` +function (callback::EnergyCalculatorCallback)(u, t, integrator) + (; interval) = callback + + return condition_integrator_interval(integrator, interval) +end + +# `affect!` +function (callback::EnergyCalculatorCallback)(integrator) + # Determine time step size as difference to last time this callback was called + t = integrator.t + dt = t - callback.t[] + + # Update time of last call + callback.t[] = t + + semi = integrator.p + v_ode, u_ode = integrator.u.x + energy = callback.energy + + system = semi.systems[callback.system_index] + update_energy_calculator!(energy, system, callback.eachparticle, + callback.only_compute_force_on_fluid, callback.dv, + v_ode, u_ode, semi, t, dt) + + # Tell OrdinaryDiffEq that `u` has not been modified + u_modified!(integrator, false) + + return integrator +end + +function update_energy_calculator!(energy, system, eachparticle, + only_compute_force_on_fluid, dv, + v_ode, u_ode, semi, t, dt) + @trixi_timeit timer() "calculate energy" begin + # Update quantities that are stored in the systems. These quantities (e.g. pressure) + # still have the values from the last stage of the previous step if not updated here. + @trixi_timeit timer() "update systems and nhs" begin + # Don't create sub-timers here to avoid cluttering the timer output + @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) + end + + set_zero!(dv) + + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + foreach_system(semi) do neighbor_system + if only_compute_force_on_fluid && !(neighbor_system isa AbstractFluidSystem) + # Not a fluid system, ignore this system + return + end + + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + + interact!(dv, v, u, v_neighbor, u_neighbor, + system, neighbor_system, semi, + integrate_tlsph=true, # Required when using split integration + eachparticle=eachparticle) + end + + if !only_compute_force_on_fluid + @threaded semi for particle in eachparticle + add_acceleration!(dv, particle, system) + add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t) + end + end + + for particle in eachparticle + velocity = current_velocity(v, system, particle) + dv_particle = extract_svector(dv, system, particle) + + # The force on the clamped particle is mass times acceleration + F_particle = system.mass[particle] * dv_particle + + # To obtain energy, we need to integrate the instantaneous power. + # Instantaneous power is force applied BY the particle times its velocity. + # The force applied BY the particle is the negative of the force applied ON it. + energy[] += dot(-F_particle, velocity) * dt + end + end + + return energy +end + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + @nospecialize cb # reduce precompilation time + + ELTYPE = eltype(cb.affect!.energy) + print(io, "EnergyCalculatorCallback{$ELTYPE}(interval=", cb.affect!.interval, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + update_cb = cb.affect! + ELTYPE = eltype(update_cb.energy) + setup = [ + "interval" => update_cb.interval + ] + summary_box(io, "EnergyCalculatorCallback{$ELTYPE}", setup) + end +end diff --git a/src/callbacks/info.jl b/src/callbacks/info.jl index 59b6d3454e..5a846878ce 100644 --- a/src/callbacks/info.jl +++ b/src/callbacks/info.jl @@ -98,7 +98,7 @@ function initialize_info_callback(discrete_callback, u, t, integrator; :total_width => 100, :indentation_level => 0) - semi = integrator.p + semi = integrator.p.semi show(io_context, MIME"text/plain"(), semi) println(io, "\n") foreach_system(semi) do system diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 42b85fd669..ba6c180474 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -122,7 +122,7 @@ function initialize_postprocess_callback!(cb, u, t, integrator) end function initialize_postprocess_callback!(cb::PostprocessCallback, u, t, integrator) - semi = integrator.p + semi = integrator.p.semi set_callbacks_used!(semi, integrator) cb.git_hash[] = compute_git_hash() @@ -149,7 +149,7 @@ function (pp::PostprocessCallback)(integrator) dv_ode, du_ode = dvdu_ode.x end - semi = integrator.p + semi = integrator.p.semi t = integrator.t v_ode, u_ode = integrator.u.x filenames = system_names(semi.systems) diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index a694b3f5ae..847e54ada4 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -127,7 +127,7 @@ function initialize_save_cb!(cb, u, t, integrator) end function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, integrator) - semi = integrator.p + semi = integrator.p.semi set_callbacks_used!(semi, integrator) # Reset `latest_saved_iter` @@ -164,7 +164,7 @@ function (solution_callback::SolutionSavingCallback)(integrator) end vu_ode = integrator.u - semi = integrator.p + semi = integrator.p.semi iter = get_iter(interval, integrator) if iter == solution_callback.latest_saved_iter diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl index 372ec752f0..68b029b1b6 100644 --- a/src/callbacks/split_integration.jl +++ b/src/callbacks/split_integration.jl @@ -1,17 +1,16 @@ -mutable struct SplitIntegrationCallback - # Type-stability does not matter here, and it is impossible to implement because - # the integration will be created during the initialization. - integrator :: Any - alg :: Any - kwargs :: Any +mutable struct SplitIntegrationCallback{A, K} + alg :: A + stage_coupling :: Bool + predict_positions :: Bool + kwargs :: K end @doc raw""" - SplitIntegrationCallback(alg; kwargs...) + SplitIntegrationCallback(alg; stage_coupling=false, predict_positions=true, kwargs...) Callback to integrate the `TotalLagrangianSPHSystem`s in a `Semidiscretization` separately from the other systems. -After each time step of the main integrator (in which TLSPH systems are ignored), +For each time step of the main integrator (in which TLSPH systems are ignored), the TLSPH systems are integrated for multiple smaller time steps with their own integrator. This is useful if the TLSPH systems require much smaller time steps than the fluid systems, @@ -28,31 +27,78 @@ of fluid to solid particles is large enough (e.g. 100:1 or more). - `alg`: The time integration algorithm to use for the TLSPH systems. # Keywords -- `kwargs...`: Additional keyword arguments passed to the integrator of the TLSPH systems. +- `stage_coupling=false`: If `false`, the TLSPH systems are only updated between full + time steps of the main integrator. + If `true`, the TLSPH systems are integrated to the intermediate + stage times of the main integrator as well. The sub-integrator + integrates from the previous fluid stage time to the next stage + time, using the intermediate stage predictions for the fluid + state. This strategy is highly efficient (no sub-steps have to be + repeated) but less accurate than repeating the sub-integration + with the final (as opposed to predicted) fluid state. + Note that this type of stage-level coupling is still more accurate + than step-level coupling (`stage_coupling=false`). + For large time step size ratios, `stage_coupling=false` might + require a significantly (often 2x) smaller fluid time step size + for stability at the FSI interface. + For small time step size ratios, `stage_coupling=false` might be + sufficiently stable and more efficient than `stage_coupling=true`. + Note that `stage_coupling=true` is only compatible with fluid + time integration schemes that have monotonically increasing + stage times and no stage time smaller than the time + of the previous full step. +- `predict_positions=true`: The force on the structure due to the fluid is kept constant + during one sub-integration call. When computing this force, + the new fluid state and the old structure state are available. + To avoid inconsistencies and improve accuracy (not stability), + we can predict the structure positions at the new time with + a simple Euler step, + ``u \leftarrow u + v\,(t_{\mathrm{new}} - t_{\mathrm{previous}})``, + only for this fluid force calculation. + If `false`, use the old structure state together with the new + fluid state. If `true`, predict the structure positions for the + fluid force calculation. +- `kwargs...`: Additional keyword arguments passed to the integrator + of the TLSPH systems. Use this for callbacks like the + [`StepsizeCallback`](@ref) for choosing the sub-integration + time step. + # Examples ```jldoctest; output=false using OrdinaryDiffEq -# Low-storage RK method with fixed step size +# Low-storage RK method with CFL condition for time step size callback = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), - dt=1e-5) - -# RK method with automatic error-based step size control -callback = SplitIntegrationCallback(RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4) + dt=1.0, # This is overwritten by the stepsize callback + callback=StepsizeCallback(cfl=1.6), + stage_coupling=true) # output ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ │ SplitIntegrationCallback │ │ ════════════════════════ │ -│ alg: …………………………………………………………………… RDPK3SpFSAL49 │ -│ abstol: …………………………………………………………… 1.0e-6 │ -│ reltol: …………………………………………………………… 0.0001 │ +│ alg: …………………………………………………………………… CarpenterKennedy2N54 │ +│ stage_coupling: ……………………………………… true │ +│ predict_positions: ……………………………… true │ +│ dt: ……………………………………………………………………… 1.0 │ +│ callback: ……………………………………………………… StepsizeCallback(is_constant=true, cfl_number=1.6) │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` """ -function SplitIntegrationCallback(alg; kwargs...) - split_integration_callback = SplitIntegrationCallback(nothing, alg, kwargs) +function SplitIntegrationCallback(alg; stage_coupling=false, predict_positions=true, + kwargs...) + # Add lightweight callback to (potentially) update the averaged velocity + # during the split integration. + if haskey(kwargs, :callback) + # Note that `CallbackSet`s can be nested + kwargs = (; kwargs..., callback=CallbackSet(values(kwargs).callback, + UpdateAveragedVelocityCallback())) + else + kwargs = (; kwargs..., callback=UpdateAveragedVelocityCallback()) + end + split_integration_callback = SplitIntegrationCallback(alg, stage_coupling, + predict_positions, pairs(kwargs)) # The first one is the `condition`, the second the `affect!` return DiscreteCallback(split_integration_callback, split_integration_callback, @@ -60,10 +106,10 @@ function SplitIntegrationCallback(alg; kwargs...) save_positions=(false, false)) end -function initialize_split_integration!(cb, u, t, integrator) - semi = integrator.p +function initialize_split_integration!(cb, vu_ode, t, integrator) + semi = integrator.p.semi split_integration_callback = cb.affect! - (; alg, kwargs) = split_integration_callback + (; alg, stage_coupling, predict_positions, kwargs) = split_integration_callback # Disable TLSPH integration in the original integrator semi.integrate_tlsph[] = false @@ -80,45 +126,29 @@ function initialize_split_integration!(cb, u, t, integrator) periodic_box = extract_periodic_box(semi.neighborhood_searches[1][1]) neighborhood_search = TrivialNeighborhoodSearch{ndims(first(systems))}(; periodic_box) semi_split = Semidiscretization(systems..., - neighborhood_search=TrivialNeighborhoodSearch{ndims(first(systems))}(), + neighborhood_search=neighborhood_search, parallelization_backend=semi.parallelization_backend) - # Verify that a NHS implementation is used that does not require updates - # for tlsph-neighbor interaction when neighbor is static. - foreach_system(semi_split) do system - foreach_system(semi) do neighbor - if system === neighbor - # TLSPH self-interaction is using its own NHS - return - end - neighborhood_search = get_neighborhood_search(system, neighbor, semi) - # The first element indicates if the NHS requires an update when the first - # system (the TLSPH system) changed. - # If `false`, the NHS does not need to be updated in the substeps - # where the neighbor is static. - if PointNeighbors.requires_update(neighborhood_search)[1] - throw(ArgumentError("Split integration is only supported (and useful) " * - "when the neighborhood search does not require updates " * - "for static neighbors (e.g. GridNeighborhoodSearch).")) - end - end - end - sizes_u = (u_nvariables(system) * n_integrated_particles(system) for system in systems) sizes_v = (v_nvariables(system) * n_integrated_particles(system) for system in systems) - v_ode, u_ode = integrator.u.x + v_ode, u_ode = vu_ode.x v0_ode_split = similar(v_ode, sum(sizes_v)) u0_ode_split = similar(u_ode, sum(sizes_u)) + # A buffer to store interaction forces between structure and fluid + dv_ode_split = similar(v0_ode_split) + # Copy the initial conditions to the split integrator - copy_to_split!(v_ode, u_ode, v0_ode_split, u0_ode_split, semi, semi_split) + copy_to_split!(v0_ode_split, u0_ode_split, semi_split, v_ode, u_ode, semi) # A zero `tspan` sets `tdir` to zero, which breaks adding tstops. # Therefore, we have to use an arbitrary non-zero `tspan` here, # but we remove the final tstop later. - tspan = (integrator.t, integrator.t + 1) - p = (; v_ode, u_ode, semi, semi_split) + tspan = (t, t + 1) + + # Payload that we need inside the split `kick!` function + p = (; semi=semi_split, semi_large=semi, dv_ode_split) ode_split = DynamicalODEProblem(kick_split!, drift_split!, v0_ode_split, u0_ode_split, tspan, p) @@ -127,17 +157,28 @@ function initialize_split_integration!(cb, u, t, integrator) @trixi_timeit timer() "split integration" begin @trixi_timeit timer() "init" begin TimerOutputs.@notimeit timer() begin + # Use `save_everystep=false` to avoid saving multiple copies + # of `v_ode` and `u_ode`. + # We set the final time and solve the integrator to that time in each split + # integration call, which will save the final solution if we don't set + # `save_end=false`, leading to one copy of the split state per fluid stage. split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, - kwargs...) + save_end=false, kwargs...) end + + # Remove the `tstop` for the final time (equivalent to zero `tspan`) + SciMLBase.pop_tstop!(split_integrator) end end - # Remove the `tstop` (equivalent to zero `tspan`) - SciMLBase.pop_tstop!(split_integrator) - - # Store the integrator in the callback - split_integration_callback.integrator = split_integrator + # Store the required data as payload in `integrator.p`. + # `integrator.p` is a `NamedTuple` that contains `p.semi` + # and `p.split_integration_data`. The latter is usually `nothing`, but we can use it + # to store another `NamedTuple` containing all split integration data. + vu_ode_split = copy(split_integrator.u) + payload = (; stage_coupling, predict_positions, integrator=split_integrator, + large_integrator=integrator, vu_ode_split, t_ref=Ref(t), n_reject=Ref(0)) + integrator.p = (; semi, split_integration_data=payload) return cb end @@ -150,118 +191,294 @@ end # `affect!` function (split_integration_callback::SplitIntegrationCallback)(integrator) - # Function barrier for type stability (the callback struct is untyped) - affect_inner!(integrator, split_integration_callback.integrator) + new_t = integrator.t + v_ode, u_ode = integrator.u.x + data = integrator.p.split_integration_data + data.n_reject[] = integrator.stats.nreject + + # Advance the split integrator to the new time + split_integrate!(v_ode, u_ode, new_t, data) + + # Update split solution stored in `integrator.p` + data.vu_ode_split .= data.integrator.u + data.t_ref[] = new_t + + if data.stage_coupling + # Tell OrdinaryDiffEq that `u` has NOT been modified. + # Theoretically, the TLSPH part has been modified, but in the FSAL case, + # the time at the last stage is the same as the step time, so the split integration + # above is skipped and `u` is not modified at all. + # Therefore, the derivative at the last stage can be reused for the next step. + # TODO is `u_modified` ever relevant for non-FSAL methods? + u_modified!(integrator, false) + else + # Tell OrdinaryDiffEq that `u` has been modified + u_modified!(integrator, true) + end + + return integrator end -function affect_inner!(integrator, split_integrator) - semi_split = split_integrator.p.semi_split +# No `SplitIntegrationCallback` used +function split_integrate_stage!(v_ode, u_ode, t, split_integration_data::Nothing) + return v_ode +end - semi = integrator.p - new_t = integrator.t +function split_integrate_stage!(v_ode, u_ode, t, split_integration_data) + (; stage_coupling) = split_integration_data - v_ode, u_ode = integrator.u.x + if stage_coupling + split_integrate!(v_ode, u_ode, t, split_integration_data) + end - @assert semi == split_integrator.p.semi - split_integrator.p = (; v_ode, u_ode, semi, semi_split) + return v_ode +end + +function split_integrate!(v_ode, u_ode, t_new, data) + split_integrator = data.integrator + + t_previous = split_integrator.t + + restart = false + rejected = data.large_integrator.stats.nreject > data.n_reject[] + data.n_reject[] = data.large_integrator.stats.nreject + if rejected + # The previous time step was rejected. + # We have to restart the split integration from the last full time step. + # Note that we even have to do this if `t_new == t_previous` because + # the split integrator was advanced to `t_new` with a rejected fluid state. + restart = true + + # Make sure we don't have to go back before the last time step + @assert t_new >= data.t_ref[] + elseif t_new < t_previous - eps(t_previous) + # @info "" t_new - t_previous data.large_integrator.dt t_new - data.t_ref[] + # The stage time is smaller than the previous stage/step time, but the last step + # was not rejected. This means that either the RK scheme contains a negative + # node value (requesting a stage time before the time of the last full step) + # or the node values are non-monotonic (requesting a stage time before the previous + # stage time). In both cases, the scheme is not suitable for stage-level coupling + # in the form implemented here. + msg = "stage-level coupling with `SplitIntegrationCallback` requires that stage " * + "times are monotonically increasing and that no stage time is smaller " * + "than the time of the last full step. This time integration scheme can " * + "only be used with `stage_coupling=false` when creating the " * + "`SplitIntegrationCallback`. It is recommended to use a different time " * + "integration scheme with monotonic stage times." + throw(ArgumentError(msg)) + end + + if restart + # Restart the split integration from the last full time step + t_previous = data.t_ref[] + vu_ode_split = data.vu_ode_split + else + # Continue the split integration from the previous stage + vu_ode_split = split_integrator.u + end + + semi_split = split_integrator.p.semi + dv_ode_split = split_integrator.p.dv_ode_split + semi_large = split_integrator.p.semi_large @trixi_timeit timer() "split integration" begin - # Update quantities that are stored in the systems. These quantities (e.g. pressure) - # still have the values from the last stage of the previous step if not updated here. + # Copy the solutions back to the original integrator. + # We modify `v_ode` and `u_ode`, which is technically not allowed during stages, + # hence there are no guarantees about the structure part of `v_ode` and `u_ode`. + # By copying the current split integration values, we make sure that it's correct. + predict_positions = Val(data.predict_positions) + @trixi_timeit timer() "copy back" copy_from_split!(v_ode, u_ode, + vu_ode_split.x..., + semi_large, semi_split, + t_new, t_previous; + predict_positions) + end + + if !rejected && isapprox(t_new, t_previous) + # This stage time is the same as the previous stage time and we don't need to + # re-initialize the split integrator due to a rejected step. + # There is nothing to do here. + # IMPORTANT: This has to be after copying the solution to the large integrator. + # Otherwise, the large integrator might contain arbitrary values since we + # are modifying `v_ode` and `u_ode` during stages, which is not allowed. + # IMPORTANT: If we try to `return` from within the `@trixi_timeit` block, + # the timer output will be messed up for some reason. + return true + end + + @trixi_timeit timer() "split integration" begin + if restart + # Reset the split integrator to the state at the last full time step + @trixi_timeit timer() "init" begin + TimerOutputs.@notimeit timer() begin + SciMLBase.reinit!(split_integrator, vu_ode_split; + t0=t_previous, tf=t_new) + end + end + else + # Continue from the previous state + add_tstop!(split_integrator, t_new) + end + + # Update the large semidiscretization with the predicted structure positions @trixi_timeit timer() "update systems and nhs" begin - update_systems_and_nhs(v_ode, u_ode, semi, new_t) + update_systems_and_nhs(v_ode, u_ode, semi_large, t_new) + end + + # Compute structure-fluid interaction forces that are kept constant + # when applied during the split integration. + @trixi_timeit timer() "system interaction" begin + set_zero!(dv_ode_split) + + other_interaction_split!(dv_ode_split, semi_large, v_ode, u_ode, semi_split) end # Integrate the split integrator up to the new time - add_tstop!(split_integrator, new_t) - @trixi_timeit timer() "solve" SciMLBase.solve!(split_integrator) + sol = SciMLBase.solve!(split_integrator) + if sol.retcode != SciMLBase.ReturnCode.Success + error("`SplitIntegrationCallback` failed with return code $(sol.retcode). " * + "Try reducing the split integration time step size.") + end - v_ode_split, u_ode_split = split_integrator.u.x + # compute_average_velocity(split_integrator, semi_split, t_new, t_previous, u_ode, semi_large) # Copy the solutions back to the original integrator - @trixi_timeit timer() "copy back" copy_from_split!(v_ode, u_ode, v_ode_split, - u_ode_split, semi, semi_split) + @trixi_timeit timer() "copy back" copy_from_split!(v_ode, u_ode, + split_integrator.u.x..., + semi_large, semi_split, + t_new, t_previous) + + # foreach_system(semi_split) do system + # # Compute current acceleration + # kick_split!(dv_ode_split, split_integrator.u.x..., split_integrator.p, t_new) + + # # Save the acceleration in the system + # dv = wrap_v(dv_ode_split, system, semi_split) + + # @threaded semi_split for particle in each_integrated_particle(system) + # for i in 1:ndims(system) + # system.cache.acceleration[i, particle] = dv[i, particle] + # end + # end + # end end - # Tell OrdinaryDiffEq that `u` has been modified - u_modified!(integrator, true) + return true +end - return integrator +function compute_average_velocity(split_integrator, semi_split, t_new, t_previous, u_ode, semi_large) + foreach_system(semi_split) do system + # Save the acceleration in the system + u_after = wrap_u(split_integrator.u.x[2], system, semi_split) + u_before = wrap_u(u_ode, system, semi_large) + v_avg = system.cache.average_velocity + + @threaded semi_split for particle in each_integrated_particle(system) + for i in 1:ndims(system) + v_avg[i, particle] = 0 * v_avg[i, particle] + 1 * (u_after[i, particle] - u_before[i, particle]) / (t_new - t_previous) + end + end + end end function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t) - (; v_ode, u_ode, semi, semi_split) = p + semi_large = p.semi_large + semi_split = p.semi @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode_split) # Update the TLSPH systems - @trixi_timeit timer() "update systems" update_systems_split!(semi_split, v_ode_split, - u_ode_split, t) + @trixi_timeit timer() "update systems and nhs" begin + update_systems_split!(semi_split, v_ode_split, u_ode_split, t) + end + # Only compute structure-structure self-interaction. + # structure-fluid interaction forces are computed once + # before the split time integration loop and are applied below. @trixi_timeit timer() "system interaction" begin - system_interaction_split!(dv_ode_split, v_ode, u_ode, semi, - v_ode_split, u_ode_split, semi_split) + self_interaction_split!(dv_ode_split, v_ode_split, u_ode_split, + semi_split, semi_large) end - @trixi_timeit timer() "source terms" add_source_terms!(dv_ode_split, v_ode_split, - u_ode_split, semi, t; - semi_wrap=semi_split) + # Add structure-fluid interaction forces + dv_ode_split .+= p.dv_ode_split + + add_source_terms!(dv_ode_split, v_ode_split, u_ode_split, semi_large, t; + semi_wrap=semi_split) end function drift_split!(du_ode, v_ode, u_ode, p, t) - drift!(du_ode, v_ode, u_ode, p.semi_split, t) + @trixi_timeit timer() "drift!" begin + # Avoid cluttering the timer output with sub-timers of `drift!` + TimerOutputs.@notimeit timer() begin + drift!(du_ode, v_ode, u_ode, p, t) + end + end end # Update the systems before calling `interact!` to compute forces -function update_systems_split!(semi, v_ode, u_ode, t) - # First update step before updating the NHS - # (for example for writing the current coordinates in the solid system) - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) +function update_systems_split!(semi_split, v_ode_split, u_ode_split, t) + # First update step before updating the NHS. + # This is used for writing the current coordinates into the TLSPH system. + foreach_system(semi_split) do system + v = wrap_v(v_ode_split, system, semi_split) + u = wrap_u(u_ode_split, system, semi_split) - update_positions!(system, v, u, v_ode, u_ode, semi, t) + update_positions!(system, v, u, v_ode_split, u_ode_split, semi_split, t) end # Second update step. - # This is used to calculate density and pressure of the fluid systems - # before updating the boundary systems, - # since the fluid pressure is needed by the Adami interpolation. - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) + # This is used to calculate the deformation gradient and stress tensor. + foreach_system(semi_split) do system + v = wrap_v(v_ode_split, system, semi_split) + u = wrap_u(u_ode_split, system, semi_split) - update_quantities!(system, v, u, v_ode, u_ode, semi, t) + update_quantities!(system, v, u, v_ode_split, u_ode_split, semi_split, t) end - # Perform correction and pressure calculation - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - update_pressure!(system, v, u, v_ode, u_ode, semi, t) - end + # The `TotalLagrangianSPHSystem` doesn't have an `update_pressure!` method # No `update_boundary_interpolation!` for performance reasons, or we will lose # a lot of the speedup that we can gain with split integration. # We assume that the TLSPH particles move so little during the substeps # that the extrapolated pressure/density values can be treated as constant. - # Final update step for all remaining systems - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) + # The `TotalLagrangianSPHSystem` doesn't have an `update_final!` method +end + +function self_interaction_split!(dv_ode_split, v_ode_split, u_ode_split, semi_split, semi) + # Only loop over (TLSPH) systems in the split integrator + foreach_system(semi_split) do system + # Construct string for the interactions timer. + # Avoid allocations from string construction when no timers are used. + # TODO do we need to disable timers in split integration? + if timeit_debug_enabled() + system_index = system_indices(system, semi) + timer_str = "$(timer_name(system))$system_index-$(timer_name(system))$system_index" + else + timer_str = "" + end + + dv = wrap_v(dv_ode_split, system, semi_split) + v = wrap_v(v_ode_split, system, semi_split) + u = wrap_u(u_ode_split, system, semi_split) - update_final!(system, v, u, v_ode, u_ode, semi, t) + @trixi_timeit timer() timer_str begin + interact!(dv, v, u, v, u, system, system, semi; integrate_tlsph=true) + end end end -function system_interaction_split!(dv_ode_split, v_ode, u_ode, semi, - v_ode_split, u_ode_split, semi_split) - # Only loop over systems in the split integrator +function other_interaction_split!(dv_ode_split, semi, v_ode, u_ode, semi_split) + # Only loop over (TLSPH) systems in the split integrator foreach_system(semi_split) do system # Loop over all neighbors in the big integrator foreach_system(semi) do neighbor + if system === neighbor + # Only compute interaction with other systems + return + end + # Construct string for the interactions timer. # Avoid allocations from string construction when no timers are used. if timeit_debug_enabled() @@ -273,8 +490,8 @@ function system_interaction_split!(dv_ode_split, v_ode, u_ode, semi, end dv = wrap_v(dv_ode_split, system, semi_split) - v_system = wrap_v(v_ode_split, system, semi_split) - u_system = wrap_u(u_ode_split, system, semi_split) + v_system = wrap_v(v_ode, system, semi) + u_system = wrap_u(u_ode, system, semi) v_neighbor = wrap_v(v_ode, neighbor, semi) u_neighbor = wrap_u(u_ode, neighbor, semi) @@ -285,10 +502,12 @@ function system_interaction_split!(dv_ode_split, v_ode, u_ode, semi, end end end + + return dv_ode_split end # Copy the solution from the large integrator to the split integrator -@inline function copy_to_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split) +@inline function copy_to_split!(v_ode_split, u_ode_split, semi_split, v_ode, u_ode, semi) foreach_system(semi_split) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) @@ -308,7 +527,9 @@ end end # Copy the solution from the split integrator to the large integrator -@inline function copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split) +@inline function copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split, + t_new, t_previous; + predict_positions::Val{PREDICT}=Val(false)) where {PREDICT} foreach_system(semi_split) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) @@ -322,19 +543,20 @@ end for i in axes(u, 1) u[i, particle] = u_split[i, particle] + if PREDICT + # Predict positions at `t_new` with a simple Euler step + u[i, particle] += v[i, particle] * (t_new - t_previous) + end end end end end -function calculate_dt(v_ode, u_ode, cfl_number, p::NamedTuple) - # The split integrator contains a `NamedTuple` - return calculate_dt(v_ode, u_ode, cfl_number, p.semi_split) -end - function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SplitIntegrationCallback}) @nospecialize cb # reduce precompilation time - print(io, "SplitIntegrationCallback(alg=", cb.affect!.alg, ")") + print(io, "SplitIntegrationCallback(alg=", cb.affect!.alg, + ", stage_coupling=", cb.affect!.stage_coupling, + ", predict_positions=", cb.affect!.predict_positions, ")") end function Base.show(io::IO, ::MIME"text/plain", @@ -344,15 +566,61 @@ function Base.show(io::IO, ::MIME"text/plain", if get(io, :compact, false) show(io, cb) else - update_cb = cb.affect! + split_cb = cb.affect! setup = [ - "alg" => update_cb.alg |> typeof |> nameof |> string + "alg" => split_cb.alg |> typeof |> nameof |> string, + "stage_coupling" => split_cb.stage_coupling, + "predict_positions" => split_cb.predict_positions ] - for (key, value) in update_cb.kwargs + for (key, value) in split_cb.kwargs push!(setup, string(key) => string(value)) end summary_box(io, "SplitIntegrationCallback", setup) end end + +# === Non-public callback for updating the averaged velocity === +# When no split integration is used, this is done from the `UpdateCallback`. +# With split integration, we use this lightweight callback to avoid updating the systems. +function UpdateAveragedVelocityCallback() + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(update_averaged_velocity_callback!, + update_averaged_velocity_callback!, + initialize=(initialize_averaged_velocity_callback!), + save_positions=(false, false)) +end + +# `initialize` +function initialize_averaged_velocity_callback!(cb, vu_ode, t, integrator) + v_ode, u_ode = vu_ode.x + semi = integrator.p.semi + + foreach_system(semi) do system + initialize_averaged_velocity!(system, v_ode, semi, t) + end + + return cb +end + +# `condition` +function update_averaged_velocity_callback!(u, t, integrator) + return true +end + +# `affect!` +function update_averaged_velocity_callback!(integrator) + t_new = integrator.t + semi = integrator.p.semi + v_ode, u_ode = integrator.u.x + + foreach_system(semi) do system + compute_averaged_velocity!(system, v_ode, semi, t_new) + end + + # Tell OrdinaryDiffEq that `integrator.u` has not been modified + u_modified!(integrator, false) + + return integrator +end diff --git a/src/callbacks/steady_state_reached.jl b/src/callbacks/steady_state_reached.jl index dfbb8ec75f..bca0faafea 100644 --- a/src/callbacks/steady_state_reached.jl +++ b/src/callbacks/steady_state_reached.jl @@ -81,7 +81,7 @@ end end v_ode, u_ode = integrator.u.x - semi = integrator.p + semi = integrator.p.semi # Calculate kinetic energy ekin = sum(semi.systems) do system diff --git a/src/callbacks/stepsize.jl b/src/callbacks/stepsize.jl index bcf1b23e2e..2c958b39bc 100644 --- a/src/callbacks/stepsize.jl +++ b/src/callbacks/stepsize.jl @@ -46,7 +46,7 @@ end function initialize_stepsize_callback(discrete_callback, u, t, integrator) stepsize_callback = discrete_callback.affect! - semi = integrator.p + semi = integrator.p.semi set_callbacks_used!(semi, integrator) stepsize_callback(integrator) @@ -64,7 +64,7 @@ function (stepsize_callback::StepsizeCallback)(integrator) (; cfl_number) = stepsize_callback v_ode, u_ode = integrator.u.x - semi = integrator.p + semi = integrator.p.semi dt = @trixi_timeit timer() "calculate dt" calculate_dt(v_ode, u_ode, cfl_number, semi) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index ab76dfacc2..b4951dadaf 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -31,7 +31,7 @@ function UpdateCallback(; interval::Integer=-1, dt=0.0) update_callback! = UpdateCallback(interval) if dt > 0 - # Add a `tstop` every `dt`, and save the final solution. + # Add a `tstop` every `dt` return PeriodicCallback(update_callback!, dt, initialize=(initial_update!), save_positions=(false, false)) @@ -52,12 +52,17 @@ function initial_update!(cb, u, t, integrator) initial_update!(cb.affect!, u, t, integrator) end -function initial_update!(cb::UpdateCallback, u, t, integrator) - semi = integrator.p +function initial_update!(cb::UpdateCallback, vu_ode, t, integrator) + v_ode, u_ode = vu_ode.x + semi = integrator.p.semi # Tell the semidiscretization that the `UpdateCallback` is used semi.update_callback_used[] = true + foreach_system(semi) do system + initialize_averaged_velocity!(system, v_ode, semi, t) + end + return cb(integrator) end @@ -71,7 +76,7 @@ end # `affect!` function (update_callback!::UpdateCallback)(integrator) t = integrator.t - semi = integrator.p + semi = integrator.p.semi v_ode, u_ode = integrator.u.x # Tell OrdinaryDiffEq that `integrator.u` has NOT been modified. @@ -105,6 +110,10 @@ function (update_callback!::UpdateCallback)(integrator) particle_shifting_from_callback!(u_ode, shifting_technique(system), system, v_ode, semi, integrator) end + + foreach_system(semi) do system + compute_averaged_velocity!(system, v_ode, semi, t) + end end return integrator diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 1a39f594e8..ead895ce6d 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -91,7 +91,7 @@ initial_condition = InitialCondition(; coordinates, velocity=x -> 2x, mass=1.0, └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` """ -struct InitialCondition{ELTYPE, C, MATRIX, VECTOR} +struct InitialCondition{ELTYPE, C, MATRIX <: AbstractArray{ELTYPE}, VECTOR <: AbstractVector{ELTYPE}} particle_spacing :: ELTYPE coordinates :: C # Array{coordinates_eltype, 2} velocity :: MATRIX # Array{ELTYPE, 2} diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index fa77dab0e0..713342f00d 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -155,7 +155,7 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ v_ode, u_ode; include_wall_velocity=false, smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true, clip_negative_pressure=false, - output_directory="out", filename="plane") + output_directory="out", filename="plane", pvd=nothing, t=-1) # Don't filter out particles without neighbors to keep 2D grid structure filter_no_neighbors = false @trixi_timeit timer() "interpolate plane" begin @@ -172,11 +172,20 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ length(y_range)) pressure = reshape(results.pressure, length(x_range), length(y_range)) + # If data is on the GPU, first transfer it to CPU before writing to VTK + density = transfer2cpu(density) + velocity = transfer2cpu(velocity) + pressure = transfer2cpu(pressure) + @trixi_timeit timer() "write to vtk" vtk_grid(joinpath(output_directory, filename), x_range, y_range) do vtk vtk["density"] = density vtk["velocity"] = velocity vtk["pressure"] = pressure + + if !isnothing(pvd) && t >= 0 + pvd[t] = vtk + end end end @@ -560,6 +569,7 @@ end foreach_system(systems) do neighbor_system system_id = system_indices(neighbor_system, semi) + is_ref_system = Val(system_id == ref_id) nhs = neighborhood_searches[system_id] v = wrap_v(v_ode, neighbor_system, semi) @@ -587,8 +597,15 @@ end # According to: # u(r_a) = (∑_b u(r_b) ⋅ V_b ⋅ W(r_a-r_b)) / (∑_b V_b ⋅ W(r_a-r_b)), # where V_b = m_b / ρ_b. + # + # The equality check above is dynamic, so all methods + # of `interpolate_system!` are compiled on the GPU, even though most + # are skipped. This will fail if we try to access non-existing + # fields in the cache (because the cache belongs to `ref_system`). + # `is_ref_system` allows us to make this check static. interpolate_system!(cache, v, neighbor_system, - point, neighbor, volume_b, W_ab, clip_negative_pressure) + point, neighbor, volume_b, W_ab, + clip_negative_pressure, is_ref_system) else other_density[point] += m_b * W_ab @@ -674,13 +691,20 @@ end end function interpolate_system!(cache, v, neighbor_system, - point, neighbor, volume_b, W_ab, clip_negative_pressure) + point, neighbor, volume_b, W_ab, + clip_negative_pressure, is_ref_system) return cache end -@inline function interpolate_system!(cache, v, system::AbstractFluidSystem, +# `ref_system` and `neighbor_system` are always the same system. But this is behind +# a dynamic `==` on the system index, so all methods of `interpolate_system!` are compiled +# even though most are skipped. This will fail on the GPU if we try to access non-existing +# fields in the cache (because the cache belongs to `ref_system`). +# A simple fix is to pass both systems and skip distinct system types by dispatching. +@inline function interpolate_system!(cache, v, ref_system::AbstractFluidSystem, + system::AbstractFluidSystem, point, neighbor, volume_b, W_ab, - clip_negative_pressure) + clip_negative_pressure, ::Val{true}) velocity = current_velocity(v, system, neighbor) for i in axes(cache.velocity, 1) cache.velocity[i, point] += velocity[i] * volume_b * W_ab @@ -698,9 +722,15 @@ end return cache end -@inline function interpolate_system!(cache, v, system::AbstractStructureSystem, +# `ref_system` and `neighbor_system` are always the same system. But this is behind +# a dynamic `==` on the system index, so all methods of `interpolate_system!` are compiled +# even though most are skipped. This will fail on the GPU if we try to access non-existing +# fields in the cache (because the cache belongs to `ref_system`). +# A simple fix is to pass both systems and skip distinct system types by dispatching. +@inline function interpolate_system!(cache, v, ref_system::AbstractStructureSystem, + system::AbstractStructureSystem, point, neighbor, volume_b, W_ab, - clip_negative_pressure) + clip_negative_pressure, ::Val{true}) velocity = current_velocity(v, system, neighbor) for i in axes(cache.velocity, 1) cache.velocity[i, point] += velocity[i] * volume_b * W_ab diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 38531b37ff..90c97efe2e 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -159,7 +159,7 @@ end @inline function system_indices(system, semi) # Note that this takes only about 5 ns, while mapping systems to indices with a `Dict` # is ~30x slower because `hash(::System)` is very slow. - index = findfirst(==(system), semi.systems) + index = findfirst(s -> s === system, semi.systems) if isnothing(index) throw(ArgumentError("system is not in the semidiscretization")) @@ -286,7 +286,7 @@ function semidiscretize(semi, tspan; reset_threads=true) semi_new = @set semi_.systems = set_system_links.(semi_.systems, Ref(semi_)) @info "To move data to the GPU, `semidiscretize` creates a deep copy of the passed " * - "`Semidiscretization`. Use `semi = ode.p` to access simulation data." + "`Semidiscretization`. Use `semi = ode.p.semi` to access simulation data." else semi_new = semi end @@ -300,7 +300,9 @@ function semidiscretize(semi, tspan; reset_threads=true) # Reset callback flag that will be set by the `UpdateCallback` semi_new.update_callback_used[] = false - return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new) + p = @NamedTuple{semi::typeof(semi_new), split_integration_data::Any}((semi_new, + nothing)) + return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, p) end """ @@ -404,64 +406,81 @@ function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) end end -function drift!(du_ode, v_ode, u_ode, semi, t) +function drift!(du_ode, v_ode, u_ode, p, t) + (; semi) = p + @trixi_timeit timer() "drift!" begin - @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du_ode) - - @trixi_timeit timer() "velocity" begin - # Set velocity and add acceleration for each system - foreach_system(semi) do system - du = wrap_u(du_ode, system, semi) - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - @threaded semi for particle in each_integrated_particle(system) - # This can be dispatched per system - add_velocity!(du, v, u, particle, system, semi, t) - end - end + # Set velocity and add acceleration for each system + foreach_system(semi) do system + du = wrap_u(du_ode, system, semi) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + set_velocity!(du, v, u, system, semi, t) end end return du_ode end -@inline function add_velocity!(du, v, u, particle, system, semi::Semidiscretization, t) - add_velocity!(du, v, u, particle, system, t) +# Generic fallback for all systems that don't define this function +function set_velocity!(du, v, u, system, semi, t) + set_velocity_default!(du, v, u, system, semi, t) end -@inline function add_velocity!(du, v, u, particle, system::TotalLagrangianSPHSystem, - semi::Semidiscretization, t) - # Only add velocity for TLSPH systems if they are integrated +# Only set velocity for TLSPH systems if they are integrated +function set_velocity!(du, v, u, system::TotalLagrangianSPHSystem, semi, t) if semi.integrate_tlsph[] - add_velocity!(du, v, u, particle, system, t) - end -end - -@inline function add_velocity!(du, v, u, particle, system, t) - # Generic fallback for all systems that don't define this function - for i in 1:ndims(system) - @inbounds du[i, particle] = v[i, particle] + set_velocity_default!(du, v, u, system, semi, t) + else + set_zero!(du) end return du end # Solid wall boundary system doesn't integrate the particle positions -@inline add_velocity!(du, v, u, particle, system::WallBoundarySystem, t) = du +function set_velocity!(du, v, u, system::WallBoundarySystem, semi, t) + # Note that `du` is of length zero, so we don't have to set it to zero + return du +end -@inline function add_velocity!(du, v, u, particle, system::AbstractFluidSystem, t) - # This is zero unless a shifting technique is used - delta_v_ = delta_v(system, particle) +# Fluid systems integrate the particle positions and can have a shifting velocity +function set_velocity!(du, v, u, system::AbstractFluidSystem, semi, t) + @threaded semi for particle in each_integrated_particle(system) + delta_v_ = @inbounds delta_v(system, particle) - for i in 1:ndims(system) - @inbounds du[i, particle] = v[i, particle] + delta_v_[i] + for i in 1:ndims(system) + @inbounds du[i, particle] = v[i, particle] + delta_v_[i] + end + end + + return du +end + +function set_velocity_default!(du, v, u, system, semi, t) + @threaded semi for particle in each_integrated_particle(system) + for i in 1:ndims(system) + @inbounds du[i, particle] = v[i, particle] + end end return du end -function kick!(dv_ode, v_ode, u_ode, semi, t) +# This defaults to optimized GPU copy that is about 4x faster than the threaded version above +function set_velocity_default!(du::AbstractGPUArray, v, u, system, semi, t) + indices = CartesianIndices(du) + copyto!(du, indices, v, indices) +end + +function kick!(dv_ode, v_ode, u_ode, p, t) + (; semi, split_integration_data) = p + + # This is a no-op if no split integration + # or split integration without stage-coupling is used. + split_integrate_stage!(v_ode, u_ode, t, split_integration_data) + @trixi_timeit timer() "kick!" begin # Check that the `UpdateCallback` is used if required check_update_callback(semi) @@ -474,13 +493,29 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode, semi) - @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, - semi, t) + add_source_terms!(dv_ode, v_ode, u_ode, semi, t) end return dv_ode end +@inline save_acceleration!(system, dv_ode, semi) = system + +@inline function save_acceleration!(system::TotalLagrangianSPHSystem, dv_ode, semi) + if semi.integrate_tlsph[] + # Save the acceleration in the system + dv = wrap_v(dv_ode, system, semi) + + @threaded semi for particle in each_moving_particle(system) + for i in 1:ndims(system) + system.cache.acceleration[i, particle] = dv[i, particle] + end + end + end + + return system +end + # Update the systems and neighborhood searches (NHS) for a simulation # before calling `interact!` to compute forces. function update_systems_and_nhs(v_ode, u_ode, semi, t) @@ -537,6 +572,7 @@ end # The `SplitIntegrationCallback` overwrites `semi_wrap` to use a different # semidiscretization for wrapping arrays. +# `semi_wrap` is the small semidiscretization, `semi` is the large semidiscretization. # TODO `semi` is not used yet, but will be used when the source terms API is modified # to match the custom quantities API. function add_source_terms!(dv_ode, v_ode, u_ode, semi, t; semi_wrap=semi) @@ -545,54 +581,65 @@ function add_source_terms!(dv_ode, v_ode, u_ode, semi, t; semi_wrap=semi) v = wrap_v(v_ode, system, semi_wrap) u = wrap_u(u_ode, system, semi_wrap) - @threaded semi for particle in each_integrated_particle(system) - # Dispatch by system type to exclude boundary systems. - # `integrate_tlsph` is extracted from the `semi_wrap`, so that this function - # can be used in the `SplitIntegrationCallback` as well. - add_acceleration!(dv, particle, system, semi_wrap.integrate_tlsph[]) - add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t, - semi_wrap.integrate_tlsph[]) - end + # `integrate_tlsph` is extracted from the `semi_wrap`, so that this function + # can be used in the `SplitIntegrationCallback` as well. + # In this case, `semi_wrap` will be the small sub-integration semidiscretization. + add_source_terms!(dv, v, u, system, semi, t, semi_wrap.integrate_tlsph[]) end return dv_ode end -@inline source_terms(system) = nothing -@inline source_terms(system::Union{AbstractFluidSystem, AbstractStructureSystem}) = system.source_terms - -@inline function add_acceleration!(dv, particle, system, integrate_tlsph) - add_acceleration!(dv, particle, system) +# This is a no-op by default but can be dispatched by system type +function add_source_terms!(dv, v, u, system, semi, t, integrate_tlsph) + return dv end -@inline function add_acceleration!(dv, particle, system::TotalLagrangianSPHSystem, - integrate_tlsph) - integrate_tlsph && add_acceleration!(dv, particle, system) +function add_source_terms!(dv, v, u, + system::Union{AbstractFluidSystem, AbstractStructureSystem}, + semi, t, integrate_tlsph) + add_source_terms_inner!(dv, v, u, system, semi, t) end -@inline add_acceleration!(dv, particle, system) = dv +function add_source_terms!(dv, v, u, system::TotalLagrangianSPHSystem, + semi, t, integrate_tlsph) + if integrate_tlsph + add_source_terms_inner!(dv, v, u, system, semi, t) + end -@inline function add_acceleration!(dv, particle, - system::Union{AbstractFluidSystem, - AbstractStructureSystem}) - (; acceleration) = system + return dv +end - for i in 1:ndims(system) - dv[i, particle] += acceleration[i] +function add_source_terms_inner!(dv, v, u, + system::Union{AbstractFluidSystem, + AbstractStructureSystem}, + semi, t) + if iszero(system.acceleration) && isnothing(source_terms(system)) + # Nothing to do + return dv + end + + @trixi_timeit timer() "source terms" begin + @threaded semi for particle in each_integrated_particle(system) + add_acceleration!(dv, system, particle) + add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t) + end end return dv end -@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t, - integrate_tlsph) - add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) -end +@inline source_terms(system) = nothing +@inline source_terms(system::Union{AbstractFluidSystem, AbstractStructureSystem}) = system.source_terms + +@inline function add_acceleration!(dv, system, particle) + (; acceleration) = system -@inline function add_source_terms_inner!(dv, v, u, particle, - system::TotalLagrangianSPHSystem, - source_terms_, t, integrate_tlsph) - integrate_tlsph && add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) + for i in 1:ndims(system) + @inbounds dv[i, particle] += acceleration[i] + end + + return dv end @inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) @@ -652,18 +699,25 @@ end function system_interaction!(dv_ode, v_ode, u_ode, semi) # Call `interact!` for each pair of systems foreach_system(semi) do system - foreach_system(semi) do neighbor - # Construct string for the interactions timer. - # Avoid allocations from string construction when no timers are used. - if timeit_debug_enabled() - system_index = system_indices(system, semi) - neighbor_index = system_indices(neighbor, semi) - timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" - else - timer_str = "" - end + if system isa WeaklyCompressibleSPHSystem && false + system_index = system_indices(system, semi) + timer_str = "$(timer_name(system))$system_index-*" + + @trixi_timeit timer() timer_str interact_all!(dv_ode, v_ode, u_ode, system, semi) + else + foreach_system(semi) do neighbor + # Construct string for the interactions timer. + # Avoid allocations from string construction when no timers are used. + if timeit_debug_enabled() + system_index = system_indices(system, semi) + neighbor_index = system_indices(neighbor, semi) + timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" + else + timer_str = "" + end - interact!(dv_ode, v_ode, u_ode, system, neighbor, semi, timer_str=timer_str) + interact!(dv_ode, v_ode, u_ode, system, neighbor, semi, timer_str=timer_str) + end end end @@ -690,7 +744,7 @@ end function check_update_callback(semi) foreach_system(semi) do system # This check will be optimized away if the system does not require the callback - if requires_update_callback(system) && !semi.update_callback_used[] + if requires_update_callback(system, semi) && !semi.update_callback_used[] system_name = system |> typeof |> nameof throw(ArgumentError("`UpdateCallback` is required for `$system_name`")) end diff --git a/src/io/io.jl b/src/io/io.jl index 4e9469a637..39a76f1ed0 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -25,7 +25,7 @@ end function create_meta_data_dict(callback, integrator) git_hash = callback.git_hash prefix = hasproperty(callback, :prefix) ? callback.prefix : "" - semi = integrator.p + semi = integrator.p.semi names = system_names(semi.systems) meta_data = Dict{String, Any}() @@ -69,7 +69,7 @@ function add_simulation_info!(info, git_hash, integrator) end info["technical_setup"] = Dict{String, Any}() - info["technical_setup"]["parallelization_backend"] = type2string(integrator.p.parallelization_backend) + info["technical_setup"]["parallelization_backend"] = type2string(integrator.p.semi.parallelization_backend) info["technical_setup"]["number_of_threads"] = Threads.nthreads() end diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index d6439741ac..bc72e155b3 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -214,7 +214,7 @@ end return ndims(system) end -@inline requires_update_callback(system::ParticlePackingSystem) = true +@inline requires_update_callback(system::ParticlePackingSystem, semi) = true function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem) vtk["velocity"] = [advection_velocity(v, system, particle) @@ -259,7 +259,13 @@ function kinetic_energy(system::ParticlePackingSystem, v_ode, u_ode, semi, t) end @inline source_terms(system::ParticlePackingSystem) = nothing -@inline add_acceleration!(dv, particle, system::ParticlePackingSystem) = dv + +# Special case for `ParticlePackingSystem`, which is an `AbstractFluidSystem` +# but doesn't have source terms or gravity/acceleration. +function add_source_terms!(dv, v, u, system::ParticlePackingSystem, + semi, t, integrate_tlsph) + return dv +end # Update from `UpdateCallback` (between time steps) update_particle_packing(system, v_ode, u_ode, semi, integrator) = system @@ -390,15 +396,19 @@ end end # Skip for fixed systems -@inline function add_velocity!(du, v, u, particle, - system::ParticlePackingSystem{<:Any, true}, t) +@inline function set_velocity!(du, v, u, + system::ParticlePackingSystem{<:Any, true}, + semi, t) + # Note that `du` is of size `(0, n_particles)`, so we don't have to set it to zero return du end # Add advection velocity -@inline function add_velocity!(du, v, u, particle, system::ParticlePackingSystem, t) - for i in 1:ndims(system) - du[i, particle] = system.advection_velocity[i, particle] +@inline function set_velocity!(du, v, u, system::ParticlePackingSystem, semi, t) + @threaded semi for particle in each_integrated_particle(system) + for i in 1:ndims(system) + @inbounds du[i, particle] = system.advection_velocity[i, particle] + end end return du diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl index 186a9e7b7d..fff6294b63 100644 --- a/src/schemes/boundary/open_boundary/rhs.jl +++ b/src/schemes/boundary/open_boundary/rhs.jl @@ -65,7 +65,7 @@ function interact!(dv, v_particle_system, u_particle_system, density_diffusion!(dv, density_diffusion(particle_system), v_particle_system, particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, - particle_system, grad_kernel) + particle_system, grad_kernel, Val(true)) # Open boundary pressure evolution matches the corresponding fluid system: # - EDAC: Compute pressure evolution like the fluid system diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index d74a69b34b..07114b6b17 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -244,7 +244,7 @@ end @inline buffer(system::OpenBoundarySystem) = system.buffer # The `UpdateCallback` is required to update particle positions between time steps -@inline requires_update_callback(system::OpenBoundarySystem) = true +@inline requires_update_callback(system::OpenBoundarySystem, semi) = true function smoothing_length(system::OpenBoundarySystem, particle) return system.smoothing_length @@ -289,17 +289,20 @@ function calculate_dt(v_ode, u_ode, cfl_number, system::OpenBoundarySystem, semi return Inf end -@inline function add_velocity!(du, v, u, particle, system::OpenBoundarySystem, t) - boundary_zone = current_boundary_zone(system, particle) +@inline function set_velocity!(du, v, u, system::OpenBoundarySystem, semi, t) + @threaded semi for particle in each_integrated_particle(system) + boundary_zone = current_boundary_zone(system, particle) - pos = current_coords(u, system, particle) - v_particle = reference_velocity(boundary_zone, v, system, particle, pos, t) + pos = @inbounds current_coords(u, system, particle) + v_particle = @inbounds reference_velocity(boundary_zone, v, system, + particle, pos, t) - # This is zero unless a shifting technique is used - delta_v_ = delta_v(system, particle) + # This is zero unless a shifting technique is used + delta_v_ = @inbounds delta_v(system, particle) - for i in 1:ndims(system) - @inbounds du[i, particle] = v_particle[i] + delta_v_[i] + for i in 1:ndims(system) + @inbounds du[i, particle] = v_particle[i] + delta_v_[i] + end end return du diff --git a/src/schemes/boundary/prescribed_motion.jl b/src/schemes/boundary/prescribed_motion.jl index 6e5fff7272..4ea5b19ef2 100644 --- a/src/schemes/boundary/prescribed_motion.jl +++ b/src/schemes/boundary/prescribed_motion.jl @@ -83,9 +83,10 @@ function initialize_prescribed_motion!(::Nothing, initial_condition, end function (prescribed_motion::PrescribedMotion)(coordinates, velocity, acceleration, - ismoving, system, semi, t) + ismoving, system, semi, t_) (; movement_function, is_moving, moving_particles) = prescribed_motion + t = convert(Float64, t_) ismoving[] = is_moving(t) is_moving(t) || return nothing @@ -160,19 +161,22 @@ PrescribedMotion{...} function OscillatingMotion2D(; frequency, translation_vector, rotation_angle, rotation_center, rotation_phase_offset=0, tspan=(-Inf, Inf), ramp_up_tspan=(0.0, 0.0), moving_particles=nothing) - translation_vector_ = SVector{2}(translation_vector) - rotation_center_ = SVector{2}(rotation_center) + translation_vector_ = SVector{2, Float64}(translation_vector) + rotation_center_ = SVector{2, Float64}(rotation_center) + frequency_ = convert(Float64, frequency) + rotation_angle_ = convert(Float64, rotation_angle) + rotation_phase_offset_ = convert(Float64, rotation_phase_offset) @inline function movement_function(x, t) if isfinite(tspan[1]) t = t - tspan[1] end - sin_scaled = sin(frequency * 2pi * t) + sin_scaled = sin(frequency_ * 2pi * t) translation = sin_scaled * translation_vector_ x_centered = x .- rotation_center_ - sin_scaled_offset = sin(2pi * (frequency * t - rotation_phase_offset)) - angle = rotation_angle * sin_scaled_offset + sin_scaled_offset = sin(2pi * (frequency_ * t - rotation_phase_offset_)) + angle = rotation_angle_ * sin_scaled_offset rotated = SVector(x_centered[1] * cos(angle) - x_centered[2] * sin(angle), x_centered[1] * sin(angle) + x_centered[2] * cos(angle)) diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl index 9f5d8926a1..db02afd18f 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl @@ -438,7 +438,7 @@ end # Otherwise, `@threaded` does not work here with Julia ARM on macOS. # See https://github.com/JuliaSIMD/Polyester.jl/issues/88. @inline function apply_state_equation!(boundary_model, density, particle) - boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0) + boundary_model.pressure[particle] = max(boundary_model.state_equation(density), -Inf) end function compute_pressure!(boundary_model, @@ -509,7 +509,7 @@ function compute_adami_density!(boundary_model, system, v, particle) # Limit pressure to be non-negative to avoid attractive forces between fluid and # boundary particles at free surfaces (sticking artifacts). - @inbounds pressure[particle] = max(pressure[particle], 0) + # @inbounds pressure[particle] = max(pressure[particle], 0) # Apply inverse state equation to compute density (not used with EDAC) inverse_state_equation!(density, state_equation, pressure, particle) @@ -681,8 +681,10 @@ end (; volume, wall_velocity) = cache # Prescribed velocity of the boundary particle. - # This velocity is zero when not using moving boundaries. - v_boundary = current_velocity(v, system, particle) + # This velocity is zero when not using moving boundaries or TLSPH. + # If not using TLSPH with velocity averaging, this function simply + # forwards to `current_velocity`. + v_boundary = velocity_for_viscosity(v, system, particle) for dim in eachindex(v_boundary) # The second term is the precalculated smoothed velocity field of the fluid diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 8725190506..c60741d99c 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -132,7 +132,8 @@ end particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, grad_kernel) + m_b, rho_a, rho_b, grad_kernel; + apply_density_diffusion=Val(false)) return dv end @@ -142,7 +143,8 @@ end neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, grad_kernel) + m_b, rho_a, rho_b, grad_kernel; + apply_density_diffusion=Val(false)) vdiff = current_velocity(v_particle_system, particle_system, particle) - current_velocity(v_neighbor_system, neighbor_system, neighbor) @@ -152,15 +154,10 @@ end dv[end, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel) - # Artificial density diffusion should only be applied to systems representing a fluid - # with the same physical properties i.e. density and viscosity. - # TODO: shouldn't be applied to particles on the interface (depends on PR #539) - if particle_system === neighbor_system - density_diffusion!(dv, density_diffusion(particle_system), - v_particle_system, particle, neighbor, - pos_diff, distance, m_b, rho_a, rho_b, particle_system, - grad_kernel) - end + density_diffusion!(dv, density_diffusion(particle_system), + v_particle_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, particle_system, + grad_kernel, apply_density_diffusion) end function calculate_dt(v_ode, u_ode, cfl_number, system::AbstractFluidSystem, semi) diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index ceaff12a21..3aecb4ea45 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -5,7 +5,9 @@ abstract type AbstractShiftingTechnique end # WARNING: Be careful if defining this function for a specific system type. # The version for a specific system type will override this generic version. -requires_update_callback(system) = requires_update_callback(shifting_technique(system)) +function requires_update_callback(system, semi) + return requires_update_callback(shifting_technique(system)) +end requires_update_callback(::Nothing) = false requires_update_callback(::AbstractShiftingTechnique) = false diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index ff8d8e079e..325e79b3b9 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -212,7 +212,7 @@ end pos_diff, distance, m_b, rho_a, rho_b, particle_system::Union{AbstractFluidSystem, OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}}, - grad_kernel) + grad_kernel, ::Val{true}) # Density diffusion terms are all zero for distance zero. # See `src/general/smoothing_kernels.jl` for more details. distance^2 < eps(initial_smoothing_length(particle_system)^2) && return @@ -235,6 +235,6 @@ end # Density diffusion `nothing` or interaction other than fluid-fluid @inline function density_diffusion!(dv, density_diffusion, v_particle_system, particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, - particle_system, grad_kernel) + particle_system, grad_kernel, apply_density_diffusion) return dv end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index acb0f06894..abc5b62de4 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -7,6 +7,11 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::WeaklyCompressibleSPHSystem, neighbor_system, semi) (; density_calculator, correction) = particle_system + # Artificial density diffusion should only be applied to systems representing a fluid + # with the same physical properties i.e. density and viscosity. + # TODO: shouldn't be applied to particles on the interface (depends on PR #539) + apply_density_diffusion = Val(particle_system === neighbor_system) + sound_speed = system_sound_speed(particle_system) surface_tension_a = surface_tension_model(particle_system) @@ -100,7 +105,8 @@ function interact!(dv, v_particle_system, u_particle_system, @inbounds continuity_equation!(dv, density_calculator, particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, - pos_diff, distance, m_b, rho_a, rho_b, grad_kernel) + pos_diff, distance, m_b, rho_a, rho_b, grad_kernel; + apply_density_diffusion) end # Debug example # periodic_box = neighborhood_search.periodic_box @@ -112,6 +118,146 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end +function interact_all!(dv_ode, v_ode, u_ode, + particle_system::WeaklyCompressibleSPHSystem, semi) + (; density_calculator, correction) = particle_system + + dv = wrap_v(dv_ode, particle_system, semi) + v_particle_system = wrap_v(v_ode, particle_system, semi) + u_particle_system = wrap_u(u_ode, particle_system, semi) + + vs = Tuple(wrap_v(v_ode, system, semi) for system in semi.systems) + us = Tuple(wrap_u(u_ode, system, semi) for system in semi.systems) + + neighborhood_seaches = semi.neighborhood_searches[system_indices(particle_system, semi)] + + @threaded semi for particle in each_integrated_particle(particle_system) + foreach_noalloc(semi.systems, neighborhood_seaches, vs, us) do (neighbor_system, nhs, v_neighbor_system, u_neighbor_system) + sound_speed = system_sound_speed(particle_system) + + surface_tension_a = surface_tension_model(particle_system) + surface_tension_b = surface_tension_model(neighbor_system) + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff + PointNeighbors.foreach_neighbor(system_coords, neighbor_system_coords, + nhs, particle) do particle, neighbor, pos_diff, distance + # `foreach_point_neighbor` makes sure that `particle` and `neighbor` are + # in bounds of the respective system. For performance reasons, we use `@inbounds` + # in this hot loop to avoid bounds checking when extracting particle quantities. + rho_a = @inbounds current_density(v_particle_system, particle_system, particle) + rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor) + rho_mean = (rho_a + rho_b) / 2 + + # Determine correction factors. + # This can be ignored, as these are all 1 when no correction is used. + (viscosity_correction, pressure_correction, + surface_tension_correction) = free_surface_correction(correction, particle_system, + rho_mean) + + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + + m_a = @inbounds hydrodynamic_mass(particle_system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + + # The following call is equivalent to + # `p_a = current_pressure(v_particle_system, particle_system, particle)` + # `p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)` + # Only when the neighbor system is a `WallBoundarySystem` or a `TotalLagrangianSPHSystem` + # with the boundary model `PressureMirroring`, this will return `p_b = p_a`, which is + # the pressure of the fluid particle. + p_a, + p_b = @inbounds particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) + + dv_pressure = pressure_correction * + pressure_acceleration(particle_system, neighbor_system, + particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, grad_kernel, correction) + + # Propagate `@inbounds` to the viscosity function, which accesses particle data + dv_viscosity_ = viscosity_correction * + @inbounds dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, + grad_kernel) + + # Extra terms in the momentum equation when using a shifting technique + dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), + particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, m_a, m_b, rho_a, rho_b, + pos_diff, distance, grad_kernel, correction) + + dv_surface_tension = surface_tension_correction * + surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel) + + dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, + particle, neighbor, pos_diff, distance) + + dv_particle = dv_pressure + dv_viscosity_ + dv_tvf + dv_surface_tension + + dv_adhesion + + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += dv_particle[i] + # Debug example + # debug_array[i, particle] += dv_pressure[i] + end + + # TODO If variable smoothing_length is used, this should use the neighbor smoothing length + # Propagate `@inbounds` to the continuity equation, which accesses particle data + @inbounds continuity_equation2!(dv, density_calculator, particle_system, + neighbor_system, v_particle_system, + v_neighbor_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, grad_kernel) + end + end + end + # Debug example + # periodic_box = neighborhood_search.periodic_box + # Note: this saves a file in every stage of the integrator + # if !@isdefined iter; iter = 0; end + # TODO: This call should use public API. This requires some additional changes to simplify the calls. + # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter += 1) + + return dv +end + +@inline function continuity_equation2!(dv, density_calculator, + particle_system, + neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, grad_kernel) + continuity_equation!(dv, density_calculator, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, grad_kernel; + apply_density_diffusion=Val(false)) +end + +@inline function continuity_equation2!(dv, density_calculator, + particle_system::WeaklyCompressibleSPHSystem, + neighbor_system::WeaklyCompressibleSPHSystem, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, grad_kernel) + continuity_equation!(dv, density_calculator, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, grad_kernel; + apply_density_diffusion=Val(true)) +end + @propagate_inbounds function particle_neighbor_pressure(v_particle_system, v_neighbor_system, particle_system, neighbor_system, diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl index c2b63c0eab..7a33bcf431 100644 --- a/src/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl @@ -3,18 +3,21 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::TotalLagrangianSPHSystem, semi; - integrate_tlsph=semi.integrate_tlsph[]) + integrate_tlsph=semi.integrate_tlsph[], + eachparticle=each_integrated_particle(particle_system)) # Different structures do not interact with each other (yet) particle_system === neighbor_system || return dv # Skip interaction if TLSPH systems are integrated separately integrate_tlsph || return dv - interact_structure_structure!(dv, v_particle_system, particle_system, semi) + interact_structure_structure!(dv, v_particle_system, particle_system, semi; + eachparticle) end # Function barrier without dispatch for unit testing -@inline function interact_structure_structure!(dv, v_system, system, semi) +@inline function interact_structure_structure!(dv, v_system, system, semi; + eachparticle=each_integrated_particle(system)) (; penalty_force) = system # Everything here is done in the initial coordinates @@ -23,9 +26,8 @@ end # Loop over all pairs of particles and neighbors within the kernel cutoff. # For structure-structure interaction, this has to happen in the initial coordinates. foreach_point_neighbor(system, system, system_coords, system_coords, semi; - points=each_integrated_particle(system)) do particle, neighbor, - initial_pos_diff, - initial_distance + points=eachparticle) do particle, neighbor, + initial_pos_diff, initial_distance # Only consider particles with a distance > 0. # See `src/general/smoothing_kernels.jl` for more details. initial_distance^2 < eps(initial_smoothing_length(system)^2) && return @@ -77,7 +79,10 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::AbstractFluidSystem, semi; - integrate_tlsph=semi.integrate_tlsph[]) + integrate_tlsph=semi.integrate_tlsph[], + eachparticle=each_integrated_particle(particle_system)) + (; boundary_model) = particle_system + # Skip interaction if TLSPH systems are integrated separately integrate_tlsph || return dv @@ -89,12 +94,9 @@ function interact!(dv, v_particle_system, u_particle_system, # Loop over all pairs of particles and neighbors within the kernel cutoff foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, semi; - points=each_integrated_particle(particle_system)) do particle, - neighbor, - pos_diff, - distance - # Only consider particles with a distance > 0. - # See `src/general/smoothing_kernels.jl` for more details. + points=eachparticle) do particle, neighbor, + pos_diff, distance + # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details. distance^2 < eps(initial_smoothing_length(particle_system)^2) && return # Apply the same force to the structure particle @@ -110,10 +112,12 @@ function interact!(dv, v_particle_system, u_particle_system, rho_a = current_density(v_particle_system, particle_system, particle) rho_b = current_density(v_neighbor_system, neighbor_system, neighbor) - # Use kernel from the fluid system in order to get the same force here in - # structure-fluid interaction as for fluid-structure interaction. + # Use kernel from the boundary model. + # This should generally be the same as the kernel and smoothing length + # of the fluid in order to get the same force here in solid-fluid interaction + # as for fluid-solid interaction. # TODO this will not use corrections if the fluid uses corrections. - grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance, particle) + grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance, particle) # In fluid-structure interaction, use the "hydrodynamic pressure" of the structure particles # corresponding to the chosen boundary model. @@ -184,7 +188,8 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::Union{WallBoundarySystem, OpenBoundarySystem}, semi; - integrate_tlsph=semi.integrate_tlsph[]) + integrate_tlsph=semi.integrate_tlsph[], + eachparticle=each_integrated_particle(particle_system)) # TODO continuity equation? return dv end diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 2d2792cd21..b166d44496 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -7,7 +7,8 @@ acceleration=ntuple(_ -> 0.0, NDIMS), penalty_force=nothing, viscosity=nothing, source_terms=nothing, boundary_model=nothing, - self_interaction_nhs=:default) + self_interaction_nhs=:default, + velocity_averaging=nothing) System for particles of an elastic structure. @@ -61,6 +62,9 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. `transpose_backend` of the [`PrecomputedNeighborhoodSearch`](@ref) is set to `true` when running on GPUs and `false` on CPUs. Alternatively, a user-defined neighborhood search can be passed here. +- `velocity_averaging`: Velocity averaging technique to be applied on the velocity field + to obtain an averaged velocity to be used in fluid-structure interaction. + See [`VelocityAveraging`](@ref) for details. !!! note If specifying the clamped particles manually (via `n_clamped_particles`), @@ -83,7 +87,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. where `beam` and `clamped_particles` are of type [`InitialCondition`](@ref). """ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D, - YM, PR, LL, LM, K, PF, V, ST, M, IM, NHS, + YM, PR, LL, LM, K, PF, V, ST, M, IM, NHS, VA, C} <: AbstractStructureSystem{NDIMS} initial_condition :: IC initial_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle] @@ -109,6 +113,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, clamped_particles_motion :: M clamped_particles_moving :: IM self_interaction_nhs :: NHS + velocity_averaging :: VA cache :: C end @@ -121,7 +126,8 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing ndims(smoothing_kernel)), penalty_force=nothing, viscosity=nothing, source_terms=nothing, boundary_model=nothing, - self_interaction_nhs=:default) + self_interaction_nhs=:default, + velocity_averaging=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -186,7 +192,8 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing initialize_prescribed_motion!(clamped_particles_motion, initial_condition_sorted, n_clamped_particles) - cache = create_cache_tlsph(clamped_particles_motion, initial_condition_sorted) + cache = (; create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)..., + create_cache_tlsph(velocity_averaging, initial_condition_sorted)...) return TotalLagrangianSPHSystem(initial_condition_sorted, initial_coordinates, current_coordinates, mass, correction_matrix, @@ -197,7 +204,7 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing smoothing_length, acceleration_, boundary_model, penalty_force, viscosity, source_terms, clamped_particles_motion, ismoving, - self_interaction_nhs, cache) + self_interaction_nhs, velocity_averaging, cache) end # Initialize self-interaction neighborhood search if not provided by the user @@ -259,7 +266,8 @@ function initialize_self_interaction_nhs(system::TotalLagrangianSPHSystem, system.viscosity, system.source_terms, system.clamped_particles_motion, system.clamped_particles_moving, - self_interaction_nhs, system.cache) + self_interaction_nhs, system.velocity_averaging, + system.cache) end extract_periodic_box(::Nothing) = nothing @@ -274,6 +282,22 @@ function create_cache_tlsph(::PrescribedMotion, initial_condition) return (; velocity, acceleration) end +function create_cache_tlsph(::VelocityAveraging, initial_condition) + averaged_velocity = zero(initial_condition.velocity) + t_last_averaging = Ref(zero(eltype(initial_condition))) + + return (; averaged_velocity, t_last_averaging) +end + +@inline function requires_update_callback(system::TotalLagrangianSPHSystem, semi) + # Velocity averaging used and `integrate_tlsph == false` means that split integration + # is used and the averaged velocity is updated there. + # Velocity averaging used, `integrate_tlsph == true`, and no fluid system in the semi + # means we are inside the split integration + return !isnothing(system.velocity_averaging) && semi.integrate_tlsph[] && + any(system -> system isa AbstractFluidSystem, semi.systems) +end + @inline function Base.eltype(::TotalLagrangianSPHSystem{<:Any, <:Any, ELTYPE}) where {ELTYPE} return ELTYPE end @@ -384,6 +408,10 @@ end return poisson_ratio[particle] end +@inline function velocity_averaging(system::TotalLagrangianSPHSystem) + return system.velocity_averaging +end + function initialize!(system::TotalLagrangianSPHSystem, semi) (; correction_matrix) = system @@ -418,6 +446,17 @@ function update_tlsph_positions!(system::TotalLagrangianSPHSystem, u, semi) return system end +# This defaults to optimized GPU copy that is about 4x faster than the threaded version above +function update_tlsph_positions!(system::TotalLagrangianSPHSystem, + u::AbstractGPUArray, semi) + (; current_coordinates) = system + + indices = CartesianIndices(u) + copyto!(current_coordinates, indices, u, indices) + + return system +end + function apply_prescribed_motion!(system::TotalLagrangianSPHSystem, prescribed_motion::PrescribedMotion, semi, t) (; clamped_particles_moving, current_coordinates, cache) = system @@ -439,6 +478,8 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode # Precompute PK1 stress tensor @trixi_timeit timer() "stress tensor" compute_pk1_corrected!(system, semi) + # @trixi_timeit timer() "KV damping" apply_kelvin_voigt_damping!(system.pk1_rho2, v, system, semi) + return system end @@ -504,6 +545,54 @@ end return deformation_grad end +@inline function apply_kelvin_voigt_damping!(pk1_rho2, v, system, semi) + (; mass, material_density) = system + + # Compute bulk modulus from Young's modulus and Poisson's ratio. + # See the table at the end of https://en.wikipedia.org/wiki/Lam%C3%A9_parameters + # TODO Should we compute the sound speed per particle and then use the maximum? + E = maximum(system.young_modulus) + K = E / (ndims(system) * (1 - 2 * maximum(system.poisson_ratio))) + + # Newton–Laplace equation + sound_speed = sqrt(K / minimum(system.material_density)) + + # Loop over all pairs of particles and neighbors within the kernel cutoff + initial_coords = initial_coordinates(system) + foreach_point_neighbor(system, system, initial_coords, initial_coords, + semi) do particle, neighbor, initial_pos_diff, initial_distance + # Only consider particles with a distance > 0. + # See `src/general/smoothing_kernels.jl` for more details. + initial_distance^2 < eps(initial_smoothing_length(system)^2) && return + + volume = @inbounds mass[neighbor] / material_density[neighbor] + v_diff = @inbounds current_velocity(v, system, particle) - + current_velocity(v, system, neighbor) + + grad_kernel = smoothing_kernel_grad(system, initial_pos_diff, + initial_distance, particle) + + # Multiply by L_{0a} + L = @inbounds correction_matrix(system, particle) + + dF = -volume * v_diff * grad_kernel' * L' + + F = deformation_gradient(system, particle) + + alpha = 1.0 + rho = material_density[particle] + + h = smoothing_length(system, particle) + P_rho2 = alpha / rho * sound_speed * h / 2 * F * (dF' * F + F' * dF) + + for j in 1:ndims(system), i in 1:ndims(system) + @inbounds pk1_rho2[i, j, particle] += P_rho2[i, j] + end + end + + return pk1_rho2 +end + # First Piola-Kirchhoff stress tensor @propagate_inbounds function pk1_stress_tensor(system, particle) (; lame_lambda, lame_mu) = system diff --git a/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl b/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl index 168d6b634e..a9b859b661 100644 --- a/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl +++ b/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl @@ -1,5 +1,6 @@ # Penalty force needs to be included first, so that `dv_penalty_force` is available # in the closure of `foreach_point_neighbor`. include("penalty_force.jl") +include("velocity_averaging.jl") include("system.jl") include("viscosity.jl") diff --git a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl new file mode 100644 index 0000000000..79cb99f01f --- /dev/null +++ b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl @@ -0,0 +1,79 @@ +struct VelocityAveraging{ELTYPE} + time_constant::ELTYPE +end + +# No velocity averaging for a system by default +@inline function velocity_averaging(system) + return nothing +end + +@inline function velocity_for_viscosity(v, system, particle) + return velocity_for_viscosity(v, velocity_averaging(system), system, particle) +end + +@inline function velocity_for_viscosity(v, ::Nothing, system, particle) + return current_velocity(v, system, particle) +end + +@inline function velocity_for_viscosity(v, ::VelocityAveraging, system, particle) + if particle <= system.n_integrated_particles + return extract_svector(system.cache.averaged_velocity, system, particle) + end + + return current_clamped_velocity(v, system, system.clamped_particles_motion, particle) +end + +function initialize_averaged_velocity!(system, v_ode, semi, t) + initialize_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t) +end + +initialize_averaged_velocity!(system, velocity_averaging, v_ode, semi, t) = system + +function initialize_averaged_velocity!(system, ::VelocityAveraging, v_ode, semi, t) + v = wrap_v(v_ode, system, semi) + + # Make sure the averaged velocity is zero for non-integrated particles + set_zero!(system.cache.averaged_velocity) + + @threaded semi for particle in each_integrated_particle(system) + v_particle = current_velocity(v, system, particle) + + for i in eachindex(v_particle) + system.cache.averaged_velocity[i, particle] = v_particle[i] + end + end + system.cache.t_last_averaging[] = t + + return system +end + +function compute_averaged_velocity!(system, v_ode, semi, t_new) + compute_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t_new) +end + +compute_averaged_velocity!(system, velocity_averaging, v_ode, semi, t_new) = system + +function compute_averaged_velocity!(system, velocity_averaging::VelocityAveraging, + v_ode, semi, t_new) + time_constant = velocity_averaging.time_constant + averaged_velocity = system.cache.averaged_velocity + dt = t_new - system.cache.t_last_averaging[] + system.cache.t_last_averaging[] = t_new + + @trixi_timeit timer() "compute averaged velocity" begin + v = wrap_v(v_ode, system, semi) + + # Compute an exponential moving average of the velocity with the given time constant + alpha = 1 - exp(-dt / time_constant) + + @threaded semi for particle in each_integrated_particle(system) + v_particle = current_velocity(v, system, particle) + + for i in eachindex(v_particle) + averaged_velocity[i, particle] = (1 - alpha) * averaged_velocity[i, particle] + alpha * v_particle[i] + end + end + end + + return system +end diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl index 55d3c996af..068db8d31a 100644 --- a/src/setups/complex_shape.jl +++ b/src/setups/complex_shape.jl @@ -46,16 +46,12 @@ function ComplexShape(geometry; particle_spacing, density, point_in_geometry_algorithm=WindingNumberJacobson(; geometry, hierarchical_winding=true, winding_number_factor=sqrt(eps())), - store_winding_number=false, grid_offset::Real=0.0, + store_winding_number=false, grid_offset=0.0, max_nparticles=10^7, pad_initial_particle_grid=2particle_spacing) if ndims(geometry) == 3 && point_in_geometry_algorithm isa WindingNumberHormann throw(ArgumentError("`WindingNumberHormann` only supports 2D geometries")) end - if grid_offset < 0.0 - throw(ArgumentError("only a positive `grid_offset` is supported")) - end - grid = particle_grid(geometry, particle_spacing; padding=pad_initial_particle_grid, grid_offset, max_nparticles) @@ -146,7 +142,12 @@ function particle_grid(geometry, particle_spacing; padding=2particle_spacing, grid_offset=0.0, max_nparticles=10^7) (; max_corner) = geometry + # First subtract the grid offset, then add it again after rounding. + # This is making sure that `min_corner` is `n * particle_spacing` + # away from `grid_offset`, so that a particle is placed at `grid_offset`. min_corner = geometry.min_corner .- grid_offset .- padding + min_corner = floor.(Int, min_corner ./ particle_spacing) .* particle_spacing + min_corner = min_corner .+ grid_offset n_particles_per_dimension = Tuple(ceil.(Int, (max_corner .- min_corner .+ 2padding) ./ diff --git a/src/util.jl b/src/util.jl index c1c5a8aa41..90972e418e 100644 --- a/src/util.jl +++ b/src/util.jl @@ -25,6 +25,24 @@ end @inline foreach_noalloc(func, collection1::Tuple{}, collection2::Tuple{}) = nothing +@inline function foreach_noalloc(func, collection1, collection2, collection3, collection4) + element1 = first(collection1) + remaining_collection1 = Base.tail(collection1) + element2 = first(collection2) + remaining_collection2 = Base.tail(collection2) + element3 = first(collection3) + remaining_collection3 = Base.tail(collection3) + element4 = first(collection4) + remaining_collection4 = Base.tail(collection4) + + func((element1, element2, element3, element4)) + + # Process remaining collection + foreach_noalloc(func, remaining_collection1, remaining_collection2, remaining_collection3, remaining_collection4) +end + +@inline foreach_noalloc(func, collection1::Tuple{}, collection2::Tuple{}, collection3::Tuple{}, collection4::Tuple{}) = nothing + # Returns `functions[index](args...)`, but in a type-stable way for a heterogeneous tuple `functions` @inline function apply_ith_function(functions, index, args...) if index == 1 @@ -181,8 +199,8 @@ Base.pointer(A::ThreadedBroadcastArray) = pointer(parent(A)) Base.size(A::ThreadedBroadcastArray) = size(parent(A)) Base.IndexStyle(::Type{<:ThreadedBroadcastArray}) = IndexLinear() -function Base.similar(A::ThreadedBroadcastArray, ::Type{T}) where {T} - return ThreadedBroadcastArray(similar(A.array, T); +function Base.similar(A::ThreadedBroadcastArray, ::Type{T}, dims::Base.Dims) where {T} + return ThreadedBroadcastArray(similar(A.array, T, dims); parallelization_backend=A.parallelization_backend) end diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl index ad70124b44..97c69229fd 100644 --- a/src/visualization/recipes_plots.jl +++ b/src/visualization/recipes_plots.jl @@ -2,12 +2,15 @@ const TrixiParticlesODESolution = ODESolution{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ODEProblem{<:Any, <:Any, <:Any, - <:Semidiscretization}} + <:NamedTuple{(:semi, + :split_integration_data), + <:Tuple{Semidiscretization, + Any}}}} # This is the main recipe RecipesBase.@recipe function f(sol::TrixiParticlesODESolution) # Redirect everything to the next recipe - return sol.u[end].x..., sol.prob.p + return sol.u[end].x..., sol.prob.p.semi end # GPU version diff --git a/test/callbacks/callbacks.jl b/test/callbacks/callbacks.jl index b54fd19c3f..40d75053f9 100644 --- a/test/callbacks/callbacks.jl +++ b/test/callbacks/callbacks.jl @@ -11,7 +11,7 @@ # and `semi.update_callback_used[]`. semi = (; integrate_tlsph=Ref(true), update_callback_used=Ref(false)) - return (; opts=(; callback=(; discrete_callbacks)), p=semi) + return (; opts=(; callback=(; discrete_callbacks)), p=(; semi)) end stepsize_cb = StepsizeCallback(cfl=1.2) @@ -25,11 +25,11 @@ integrator_with_update_periodic = make_integrator([stepsize_cb, update_cb_periodic]) integrator_with_both = make_integrator([stepsize_cb, split_cb, update_cb]) - semi_without_split = integrator_without_split.p - semi_with_split = integrator_with_split.p - semi_with_update = integrator_with_update.p - semi_with_update_periodic = integrator_with_update_periodic.p - semi_with_both = integrator_with_both.p + semi_without_split = integrator_without_split.p.semi + semi_with_split = integrator_with_split.p.semi + semi_with_update = integrator_with_update.p.semi + semi_with_update_periodic = integrator_with_update_periodic.p.semi + semi_with_both = integrator_with_both.p.semi TrixiParticles.set_callbacks_used!(semi_without_split, integrator_without_split) TrixiParticles.set_callbacks_used!(semi_with_split, integrator_with_split) @@ -56,4 +56,5 @@ include("update.jl") include("solution_saving.jl") include("steady_state_reached.jl") + include("energy_calculator.jl") end diff --git a/test/callbacks/energy_calculator.jl b/test/callbacks/energy_calculator.jl new file mode 100644 index 0000000000..e3cc73fb06 --- /dev/null +++ b/test/callbacks/energy_calculator.jl @@ -0,0 +1,158 @@ +@testset verbose=true "EnergyCalculatorCallback" begin + # Mock system + struct MockSystem <: TrixiParticles.AbstractSystem{2} + eltype::Type + end + TrixiParticles.nparticles(::MockSystem) = 4 + Base.eltype(system::MockSystem) = system.eltype + + function TrixiParticles.create_neighborhood_search(neighborhood_search, + system::MockSystem, neighbor) + return nothing + end + + system32 = MockSystem(Float32) + semi32 = Semidiscretization(system32) + system64 = MockSystem(Float64) + semi64 = Semidiscretization(system64) + + @testset "Constructor and Basic Properties" begin + # Test default constructor + callback = EnergyCalculatorCallback(system64, semi64) + @test callback.affect!.system_index == 1 + @test callback.affect!.interval == 1 + @test callback.affect!.t[] == 0.0 + @test callback.affect!.energy[] == 0.0 + @test callback.affect!.dv isa Array{Float64, 2} + @test size(callback.affect!.dv) == (2, 4) + @test callback.affect!.eachparticle == 5:4 + @test calculated_energy(callback) == 0.0 + + # Test constructor with interval + callback = EnergyCalculatorCallback(system64, semi64; interval=5) + @test callback.affect!.interval == 5 + @test eltype(callback.affect!.energy) == Float64 + @test eltype(callback.affect!.t) == Float64 + + # Test with specific element type + callback = EnergyCalculatorCallback(system32, semi32; interval=2) + @test eltype(callback.affect!.energy) == Float32 + @test eltype(callback.affect!.t) == Float32 + end + + @testset "show" begin + callback = EnergyCalculatorCallback(system64, semi64; interval=10) + + # Test compact representation + show_compact = "EnergyCalculatorCallback{Float64}(interval=10)" + @test repr(callback) == show_compact + + # Test detailed representation - check against expected box format + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ EnergyCalculatorCallback{Float64} │ + │ ═════════════════════════════════ │ + │ interval: ……………………………………………………… 10 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + end + + @testset "update_energy_calculator!" begin + # In the first test, we just move the 2x2 grid of particles up against gravity + # and test that the energy calculated is just the potential energy difference. + # In the other tests, we clamp the top row of particles and offset them to create + # stress. We can then test how much energy is required to pull the particles further + # apart against the elastic forces. + n_clamped_particles = [4, 2, 2, 2] + E = [1e6, 1e1, 1e6, 1e8] + @testset "n_clamped_particles = $(n_clamped_particles[i]), E = $(E[i])" for i in + eachindex(E) + # Create a simple 2D system with 2x2 particles + coordinates = [0.0 1.0 0.0 1.0; 0.0 0.0 1.0 1.0] + mass = [1.0, 1.0, 1.0, 1.0] + density = [1000.0, 1000.0, 1000.0, 1000.0] + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = sqrt(2) + young_modulus = E[i] + poisson_ratio = 0.4 + + initial_condition = InitialCondition(; coordinates, mass, density) + + function movement_function(initial_pos, t) + # Offset clamped particles to create stress + return initial_pos + SVector(0.0, t + 0.5) + end + is_moving(t) = true + prescribed_motion = PrescribedMotion(movement_function, is_moving) + + # Create TLSPH system with energy calculator support + system = TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, young_modulus, + poisson_ratio, + clamped_particles_motion=prescribed_motion, + n_clamped_particles=n_clamped_particles[i], + acceleration=(0.0, -2.0)) + + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 1.0)) + + # Create dummy ODE state vectors + v = zeros(2, TrixiParticles.n_integrated_particles(system)) + u = coordinates[:, 1:TrixiParticles.n_integrated_particles(system)] + v_ode = vec(v) + u_ode = vec(u) + + # Update system + TrixiParticles.initialize!(system, semi) + TrixiParticles.update_positions!(system, v, u, v_ode, u_ode, semi, 0.0) + TrixiParticles.update_quantities!(system, v, u, v_ode, u_ode, semi, 0.0) + + # Set up test parameters + energy1 = Ref(0.0) + dt1 = 0.1 + + # Test that energy is integrated, i.e., values of instantaneous power + # are accumulated over time. This initial energy should just be an offset. + # Also, half the step size means half the energy increase. + energy2 = Ref(1.0) + dt2 = 0.05 + + # Test `only_compute_force_on_fluid` + energy3 = Ref(0.0) + dt3 = 0.1 + + eachparticle = (TrixiParticles.n_integrated_particles(system) + 1):nparticles(system) + dv = zeros(2, nparticles(system)) + + TrixiParticles.update_energy_calculator!(energy1, system, eachparticle, false, + dv, v_ode, u_ode, semi, 0.0, dt1) + TrixiParticles.update_energy_calculator!(energy2, system, eachparticle, false, + dv, v_ode, u_ode, semi, 0.0, dt2) + TrixiParticles.update_energy_calculator!(energy3, system, eachparticle, true, + dv, v_ode, u_ode, semi, 0.0, dt3) + + if i == 1 + @test isapprox(energy1[], 0.8) + @test isapprox(energy2[], 1.0 + 0.4) + elseif i == 2 + # For very soft material, we can just pull up the top row of particles + # and the energy required is almost just the potential energy difference. + @test isapprox(energy1[], 0.4080357142857143) + @test isapprox(energy2[], 1.0 + 0.5 * 0.4080357142857143) + elseif i == 3 + # For a stiffer material, the stress from the offset creates larger forces + # pulling the clamped particles back down, so we need a lot of energy + # to pull the material apart. + @test isapprox(energy1[], 803.9714285714281) + @test isapprox(energy2[], 1.0 + 0.5 * 803.9714285714281) + elseif i == 4 + # For a very stiff material, the energy is even larger. + @test isapprox(energy1[], 80357.5428571428) + @test isapprox(energy2[], 1.0 + 0.5 * 80357.5428571428) + end + + @test isapprox(energy3[], 0.0) + end + end +end diff --git a/test/callbacks/info.jl b/test/callbacks/info.jl index 7c8cf881d0..489ed9a55b 100644 --- a/test/callbacks/info.jl +++ b/test/callbacks/info.jl @@ -24,7 +24,7 @@ semi = (; systems=(:system1, :system2), parallelization_backend=PolyesterBackend()) - integrator = (; p=semi, + integrator = (; p=(; semi), opts=(; callback=(; continuous_callbacks, discrete_callbacks), adaptive=true, abstol=1e-2, reltol=1e-1, diff --git a/test/count_allocations.jl b/test/count_allocations.jl index 95c7bffbe5..392ac136ca 100644 --- a/test/count_allocations.jl +++ b/test/count_allocations.jl @@ -37,7 +37,7 @@ end end # Count allocations of one call to the right-hand side (`kick!` + `drift!`) -function count_rhs_allocations(sol, semi) +function count_rhs_allocations(sol) t = sol.t[end] v_ode_, u_ode_ = sol.u[end].x @@ -47,8 +47,15 @@ function count_rhs_allocations(sol, semi) dv_ode = similar(v_ode) du_ode = similar(u_ode) - # Wrap neighborhood searches to avoid counting alloctations in the NHS update - semi_no_nhs_update = copy_semi_with_no_update_nhs(semi) + # Wrap neighborhood searches to avoid counting allocations in the NHS update. + # Note that `p.split_integration_data.split_integrator.p.semi_large` is still the old + # semidiscretization with the original NHS. + # However, we call `kick!` twice with the same `t` below, so in the second call, + # where the allocations are counted, the split integration does nothing because + # the split integrator is already at the time `t`. + p = sol.prob.p + semi_no_nhs_update = copy_semi_with_no_update_nhs(p.semi) + p_no_update = TrixiParticles.@set p.semi = semi_no_nhs_update try # Disable timers, which cause extra allocations @@ -58,26 +65,26 @@ function count_rhs_allocations(sol, semi) # `TrixiParticles.timeit_debug_enabled()` is called, which is redefined in # `disable_debug_timings` above. return @invokelatest count_rhs_allocations_inner(dv_ode, du_ode, v_ode, u_ode, - semi_no_nhs_update, t) + p_no_update, t) finally # Enable timers again @invokelatest TrixiParticles.enable_debug_timings() end end -# Function barrier to avoid type instabilites with `semi_no_nhs_update`, which will +# Function barrier to avoid type instabilites with `p_no_update`, which will # cause extra allocations. @inline function count_rhs_allocations_inner(dv_ode, du_ode, v_ode, u_ode, - semi_no_nhs_update, t) + p_no_update, t) # Run RHS once to avoid counting allocations from compilation - TrixiParticles.kick!(dv_ode, v_ode, u_ode, semi_no_nhs_update, t) - TrixiParticles.drift!(du_ode, v_ode, u_ode, semi_no_nhs_update, t) + TrixiParticles.kick!(dv_ode, v_ode, u_ode, p_no_update, t) + TrixiParticles.drift!(du_ode, v_ode, u_ode, p_no_update, t) # Count allocations allocations_kick = @allocated TrixiParticles.kick!(dv_ode, v_ode, u_ode, - semi_no_nhs_update, t) + p_no_update, t) allocations_drift = @allocated TrixiParticles.drift!(du_ode, v_ode, u_ode, - semi_no_nhs_update, t) + p_no_update, t) return allocations_kick + allocations_drift end diff --git a/test/examples/dam_break_2d_corrections.jl b/test/examples/dam_break_2d_corrections.jl index 9c3e29fbe7..8f0aaba45e 100644 --- a/test/examples/dam_break_2d_corrections.jl +++ b/test/examples/dam_break_2d_corrections.jl @@ -74,7 +74,7 @@ density_diffusion=nothing) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @testset verbose=true "$correction_name" for correction_name in keys(correction_dict) @@ -107,6 +107,6 @@ sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index c7d9fadb92..b6178330d3 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -15,9 +15,9 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -34,9 +34,9 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -60,14 +60,137 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else + @test count_rhs_allocations(sol) == 0 + end + end + + @trixi_testset "structure/oscillating_beam_2d.jl with EnergyCalculatorCallback" begin + # Load variables from the example + trixi_include(@__MODULE__, + joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), + ode=nothing, sol=nothing, E=1e5) + + # We simply clamp all particles, move them up against gravity, and verify that + # the energy calculated is just the potential energy difference. + movement_function(x, t) = x + SVector(0.0, t) + is_moving(t) = true + + tests = Dict( + "all particles clamped" => eachparticle(structure), + # Clamp everything but the top two layers of the beam + "some particles clamped" => 1:(nparticles(structure) - 162) + ) + rtol = Dict( + "all particles clamped" => sqrt(eps()), + # Clamp everything but the top two layers of the beam. + # We don't expect very accurate results here, this is more of a smoke test. + "some particles clamped" => 0.1 + ) + @testset "$name" for (name, clamped_particles) in tests + # We need a new `PrescribedMotion` for each test due to + # https://github.com/trixi-framework/TrixiParticles.jl/issues/1020 + prescribed_motion = PrescribedMotion(movement_function, is_moving) + structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, + smoothing_length, + material.E, material.nu, + clamped_particles=clamped_particles, + acceleration=(0.0, -gravity), + clamped_particles_motion=prescribed_motion) + + semi = Semidiscretization(structure_system) + ode = semidiscretize(semi, (0.0, 1.0)) + + energy_calculator = EnergyCalculatorCallback(structure_system, semi; + interval=1) + + sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL49(), save_everystep=false, + callback=energy_calculator) + + @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 + + # Potential energy difference should be m * g * h + @test isapprox(calculated_energy(energy_calculator), + sum(structure_system.mass) * gravity * 1, + rtol=rtol[name]) end end end @testset verbose=true "FSI" begin + @trixi_testset "fluid/hydrostatic_water_column_2d.jl with moving TLSPH walls" begin + # In this test, we move a water-filled tank up against gravity by 1 unit + # and verify that the energy calculated by the `EnergyCalculatorCallback` + # matches the expected potential energy. + + # Load variables from the example + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "hydrostatic_water_column_2d.jl"), + # Use tank without airspace + tank_size=(1.0, 0.9), + # Higher speed of sound for smaller error + # due to compressibility. + sound_speed=50.0, + ode=nothing, sol=nothing) + + # Move the tank up against gravity smoothly by 1 unit over 1 second + function movement_function(x, t) + return x + SVector(0.0, 0.5 * sin(pi * (t - 0.5)) + 0.5) + end + is_moving(t) = true + prescribed_motion = PrescribedMotion(movement_function, is_moving) + + # Create TLSPH system for the tank walls and clamp all particles. + # This is identical to a `WallBoundarySystem`, but now we can + # use the `EnergyCalculatorCallback` to compute the energy. + boundary_spacing = tank.boundary.particle_spacing + tlsph_kernel = WendlandC2Kernel{2}() + tlsph_smoothing_length = sqrt(2) * boundary_spacing + + tlsph_system = TotalLagrangianSPHSystem(tank.boundary, tlsph_kernel, + tlsph_smoothing_length, + 1.0e6, 0.3, + clamped_particles=eachparticle(tank.boundary), + acceleration=(0.0, -gravity), + clamped_particles_motion=prescribed_motion, + boundary_model=boundary_model) + + semi = Semidiscretization(fluid_system, tlsph_system, + parallelization_backend=PolyesterBackend()) + ode = semidiscretize(semi, (0.0, 1.0)) + + # Energy calculators for fluid + tank and fluid only + energy_calculator1 = EnergyCalculatorCallback(tlsph_system, semi; interval=1) + energy_calculator2 = EnergyCalculatorCallback(tlsph_system, semi; interval=1, + only_compute_force_on_fluid=true) + + sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL35(), save_everystep=false, + callback=CallbackSet(info_callback, + energy_calculator1, + energy_calculator2)) + + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + + # Potential energy difference should be m * g * h and h = 1. + # Since all particles are clamped, the energy added through clamped particles + # should be the same as that added to the fluid. + expected_energy_fluid = sum(fluid_system.mass) * gravity * 1 + expected_energy_tank = sum(tlsph_system.mass) * gravity * 1 + + # We don't expect very accurate results here because the fluid is weakly + # compressible and is deformed during the simulation. + # A slower prescribed motion (e.g., over 2 seconds instead of 1) or a higher + # speed of sound in the fluid would improve accuracy (and increase runtime). + @test isapprox(calculated_energy(energy_calculator1), + expected_energy_fluid + expected_energy_tank, rtol=5e-4) + @test isapprox(calculated_energy(energy_calculator2), expected_energy_fluid, + rtol=5e-4) + end + @trixi_testset "fsi/falling_water_column_2d.jl" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fsi", @@ -79,9 +202,9 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -98,9 +221,9 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -125,66 +248,70 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end - # Now use split integration and verify that we need less than 400 iterations - split_integration = SplitIntegrationCallback(RDPK3SpFSAL35(), adaptive=false, + # Use split integration and verify that we need fewer than 400 iterations + split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), dt=5e-5) - @trixi_test_nowarn trixi_include(@__MODULE__, - joinpath(examples_dir(), "fsi", - "dam_break_plate_2d.jl"), - # Use rounded dimensions to avoid warnings - initial_fluid_size=(0.15, 0.29), - # Move plate closer to be able to use a shorter - # tspan and make CI faster. - plate_position=(0.2, 0.0), - tspan=(0.0, 0.2), - E=1e7, # Stiffer plate - maxiters=400, - extra_callback=split_integration) [ - r"\[ Info: To create the self-interaction neighborhood search.*\n" - ] - @test sol.retcode == ReturnCode.Success + callbacks = CallbackSet(info_callback, saving_callback, split_integration) + # Don't use `sol` here, as this local variable would shadow the global `sol` + # from the `trixi_include` above and break the `sol.retcode` above. + sol2 = @trixi_test_nowarn solve(ode, RDPK3SpFSAL49(), maxiters=400, + save_everystep=false, callback=callbacks) + + @test sol2.retcode == ReturnCode.Success if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol2) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol2) == 0 end - # Now use split integration and verify that it is actually used for TLSPH - # by using a time step that is too large and verifying that it is crashing. - split_integration = SplitIntegrationCallback(RDPK3SpFSAL35(), adaptive=false, - dt=2e-4) - @trixi_test_nowarn trixi_include(@__MODULE__, - joinpath(examples_dir(), "fsi", - "dam_break_plate_2d.jl"), - # Use rounded dimensions to avoid warnings - initial_fluid_size=(0.15, 0.29), - # Move plate closer to be able to use a shorter - # tspan and make CI faster. - plate_position=(0.2, 0.0), - tspan=(0.0, 0.2), - E=1e7, # Stiffer plate - maxiters=500, - extra_callback=split_integration) [ - r"\[ Info: To create the self-interaction neighborhood search.*\n", - "┌ Warning: Instability detected. Aborting\n", - r".*dt was forced below floating point epsilon.*\n", - r"└ @ SciMLBase.*\n" - ] - @test sol.retcode == ReturnCode.Unstable + # Use stage-level coupling and verify that it is not compatible with + # the fluid time integration scheme `RDPK3SpFSAL49`. + split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), + dt=5e-5, stage_coupling=true) + callbacks = CallbackSet(info_callback, saving_callback, split_integration) + + msg = "stage-level coupling with `SplitIntegrationCallback` requires that" + @test_throws msg solve(ode, RDPK3SpFSAL49(), maxiters=400, + save_everystep=false, callback=callbacks) + + # Use stage-level coupling with a compatible fluid time integration scheme + split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), + stage_coupling=true, dt=5e-5) + stepsize_callback = StepsizeCallback(cfl=1.2) + callbacks = CallbackSet(info_callback, saving_callback, split_integration, + stepsize_callback) + sol2 = @trixi_test_nowarn solve(ode, + CarpenterKennedy2N54(williamson_condition=false), + maxiters=400, dt=1.0, + save_everystep=false, callback=callbacks) + + @test sol2.retcode == ReturnCode.Success if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol2) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol2) == 0 end + + # Use split integration and verify that it is actually used for TLSPH + # by using a time step that is too large and verifying that it is crashing. + split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), + stage_coupling=true, dt=2e-4) + callbacks = CallbackSet(info_callback, saving_callback, split_integration, + stepsize_callback) + + msg = "`SplitIntegrationCallback` failed with return code Unstable." + @test_throws msg solve(ode, CarpenterKennedy2N54(williamson_condition=false), + maxiters=400, dt=1.0, + save_everystep=false, callback=callbacks) end @trixi_testset "fsi/dam_break_gate_2d.jl" begin @@ -199,9 +326,9 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -216,9 +343,9 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 500 + @test count_rhs_allocations(sol) < 500 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end end @@ -229,7 +356,7 @@ joinpath(examples_dir(), "n_body", "n_body_solar_system.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "n_body/n_body_benchmark_trixi.jl" begin @@ -314,7 +441,7 @@ "rectangular_tank_2d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @trixi_testset "dem/collapsing_sand_pile_3d.jl" begin @@ -323,6 +450,6 @@ "collapsing_sand_pile_3d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index e0f7eed13d..ee90b7fa4b 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -96,7 +96,7 @@ kwargs...) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end end @@ -108,7 +108,7 @@ @test sol.retcode == ReturnCode.Success # This error varies between serial and multithreaded runs @test isapprox(error_A, 0, atol=2e-4) - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/hydrostatic_water_column_3d.jl" begin @@ -117,7 +117,7 @@ "hydrostatic_water_column_3d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/hydrostatic_water_column_3d.jl with SummationDensity" begin @@ -128,7 +128,7 @@ fluid_density_calculator=SummationDensity(), clip_negative_pressure=true) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/accelerated_tank_2d.jl" begin @@ -136,7 +136,7 @@ joinpath(examples_dir(), "fluid", "accelerated_tank_2d.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/dam_break_2d.jl" begin @@ -184,7 +184,7 @@ ] @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -199,7 +199,7 @@ r"└ New tank length in y-direction.*\n" ] @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 v_ode, u_ode = sol[end].x @test eltype(v_ode) == Float32 @test eltype(u_ode) == Float64 @@ -217,9 +217,9 @@ @test sol.retcode == ReturnCode.Success if VERSION < v"1.11" # For some reason, 1.10 produces allocations here - @test count_rhs_allocations(sol, semi) <= 32 + @test count_rhs_allocations(sol) <= 32 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -235,9 +235,9 @@ @test sol.retcode == ReturnCode.Success if VERSION < v"1.11" # For some reason, 1.10 produces allocations here - @test count_rhs_allocations(sol, semi) <= 32 + @test count_rhs_allocations(sol) <= 32 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -253,9 +253,9 @@ @test sol.retcode == ReturnCode.Success if VERSION < v"1.11" # For some reason, 1.10 produces allocations here - @test count_rhs_allocations(sol, semi) <= 32 + @test count_rhs_allocations(sol) <= 32 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -271,9 +271,9 @@ @test sol.retcode == ReturnCode.Success if VERSION < v"1.11" # For some reason, 1.10 produces allocations here - @test count_rhs_allocations(sol, semi) <= 32 + @test count_rhs_allocations(sol) <= 32 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -288,9 +288,9 @@ @test sol.retcode == ReturnCode.Success if VERSION < v"1.11" # For some reason, 1.10 produces allocations here - @test count_rhs_allocations(sol, semi) <= 32 + @test count_rhs_allocations(sol) <= 32 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end @@ -304,7 +304,7 @@ ] @test semi.neighborhood_searches[1][1].cell_list isa FullGridCellList @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/dam_break_oil_film_2d.jl" begin @@ -316,7 +316,7 @@ r"└ New tank length in y-direction.*\n" ] @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/dam_break_2phase_2d.jl" begin @@ -328,7 +328,7 @@ r"└ New tank length in y-direction.*\n" ] @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/dam_break_3d.jl" begin @@ -337,7 +337,7 @@ "dam_break_3d.jl"), tspan=(0.0, 0.1), fluid_particle_spacing=0.1) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/falling_water_column_2d.jl" begin @@ -346,7 +346,7 @@ "falling_water_column_2d.jl"), tspan=(0.0, 0.4)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/periodic_channel_2d.jl" begin @@ -355,7 +355,7 @@ "periodic_channel_2d.jl"), tspan=(0.0, 0.4)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/periodic_channel_2d.jl with PST" begin @@ -366,7 +366,7 @@ shifting_technique=ParticleShiftingTechniqueSun2017(), extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/periodic_channel_2d.jl with TVF" begin @@ -377,7 +377,7 @@ shifting_technique=TransportVelocityAdami(background_pressure=50_000.0), extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/periodic_channel_2d.jl with PST and TIC" begin @@ -389,7 +389,7 @@ pressure_acceleration=tensile_instability_control, extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/periodic_channel_2d.jl with consistent PST and TIC" begin @@ -401,7 +401,7 @@ pressure_acceleration=tensile_instability_control, extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelCharacteristicsLastiwka (WCSPH)" begin @@ -410,7 +410,7 @@ joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelCharacteristicsLastiwka (EDAC)" begin @@ -419,7 +419,7 @@ joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelMirroringTafuni (EDAC)" begin @@ -429,7 +429,7 @@ boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelMirroringTafuni (WCSPH)" begin @@ -439,7 +439,7 @@ boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`dt`)" begin @@ -463,7 +463,7 @@ joinpath(examples_dir(), "fluid", "pipe_flow_3d.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/poiseuille_flow_2d.jl (WCSPH)" begin @@ -473,7 +473,7 @@ tspan=(0.0, 0.02)) @test fluid_system isa WeaklyCompressibleSPHSystem @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/poiseuille_flow_2d.jl (EDAC)" begin @@ -483,7 +483,7 @@ tspan=(0.0, 0.02), wcsph=false) @test fluid_system isa EntropicallyDampedSPHSystem @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/lid_driven_cavity_2d.jl (EDAC)" begin @@ -492,7 +492,7 @@ "lid_driven_cavity_2d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/lid_driven_cavity_2d.jl (WCSPH)" begin @@ -501,7 +501,7 @@ "lid_driven_cavity_2d.jl"), tspan=(0.0, 0.1), wcsph=true) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/taylor_green_vortex_2d.jl (EDAC)" begin @@ -510,7 +510,7 @@ "taylor_green_vortex_2d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/taylor_green_vortex_2d.jl (WCSPH)" begin @@ -519,7 +519,7 @@ "taylor_green_vortex_2d.jl"), tspan=(0.0, 0.1), wcsph=true) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/sphere_surface_tension_2d.jl" begin @@ -527,7 +527,7 @@ joinpath(examples_dir(), "fluid", "sphere_surface_tension_2d.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/periodic_array_of_cylinders_2d.jl" begin @@ -536,7 +536,7 @@ "periodic_array_of_cylinders_2d.jl"), tspan=(0.0, 20.0)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/sphere_surface_tension_3d.jl" begin @@ -544,7 +544,7 @@ joinpath(examples_dir(), "fluid", "sphere_surface_tension_3d.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "fluid/falling_water_spheres_2d.jl" begin @@ -575,7 +575,7 @@ @test sol.retcode == ReturnCode.Success # Optionally, verify no unexpected RHS allocations - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end end @@ -618,7 +618,7 @@ @test sol.retcode == ReturnCode.Success # Optionally, verify no unexpected RHS allocations - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end end @@ -634,7 +634,7 @@ joinpath(examples_dir(), "fluid", "moving_wall_2d.jl")) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end include("dam_break_2d_corrections.jl") @@ -652,7 +652,7 @@ dt=1.0, # This is overwritten by the stepsize callback save_everystep=false, callback=callbacks) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 # Unstable with this CFL @test maximum(abs, sol.u[end]) > 2^15 @@ -671,7 +671,7 @@ dt=1.0, # This is overwritten by the stepsize callback save_everystep=false, callback=callbacks) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 # Stable with this CFL @test maximum(abs, sol.u[end]) < 2^15 @@ -690,7 +690,7 @@ dt=1.0, # This is overwritten by the stepsize callback save_everystep=false, callback=callbacks) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 @test maximum(abs, sol.u[end]) < 2^15 end end diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index fe0b3c23a4..c076156874 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -566,8 +566,6 @@ end min_corner = minimum(tank.boundary.coordinates, dims=2) max_corner = maximum(tank.boundary.coordinates, dims=2) max_corner[2] = gate_height + movement_function([0, 0], 0.1f0)[2] - # We need a very high `max_points_per_cell` because the plate resolution - # is much finer than the fluid resolution. cell_list = FullGridCellList(; min_corner, max_corner) semi_fullgrid = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_gate, structure_system, @@ -590,6 +588,51 @@ end backend = TrixiParticles.KernelAbstractions.get_backend(sol.u[end].x[1]) @test backend == Main.parallelization_backend end + + @trixi_testset "fsi/dam_break_plate_2d.jl split integration" begin + # Use split integration and verify that we need fewer than 400 iterations. + # See the CPU test for more details. + + # Import variables into scope + trixi_include_changeprecision(Float32, @__MODULE__, + joinpath(examples_dir(), "fsi", + "dam_break_plate_2d.jl"), + coordinates_eltype=Float32, + # Use rounded dimensions to avoid warnings + initial_fluid_size=(0.15f0, 0.29f0), + # Move plate closer to be able to use a shorter + # tspan and make CI faster. + plate_position=(0.2f0, 0.0f0), + E=1.0f7, # Stiffer plate + sol=nothing, ode=nothing) + + # Neighborhood search with `FullGridCellList` for GPU compatibility + min_corner = minimum(tank.boundary.coordinates, dims=2) + max_corner = maximum(tank.boundary.coordinates, dims=2) + cell_list = FullGridCellList(; min_corner, max_corner) + semi = Semidiscretization(fluid_system, boundary_system, structure_system, + neighborhood_search=GridNeighborhoodSearch{2}(; + cell_list), + parallelization_backend=Main.parallelization_backend) + ode = semidiscretize(semi, (0.0f0, 0.2f0)) + + # Set up callbacks + split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), + stage_coupling=true, dt=5.0f-5) + stepsize_callback = StepsizeCallback(cfl=1.2f0) + callbacks = CallbackSet(info_callback, saving_callback, split_integration, + stepsize_callback) + + # Run the simulation + sol = @trixi_test_nowarn solve(ode, + CarpenterKennedy2N54(williamson_condition=false), + maxiters=400, dt=1.0f0, + save_everystep=false, callback=callbacks) + + @test sol.retcode == ReturnCode.Success + backend = TrixiParticles.KernelAbstractions.get_backend(sol.u[end].x[1]) + @test backend == Main.parallelization_backend + end end @testset verbose=true "DEM" begin @@ -632,7 +675,7 @@ end tspan=(0.0f0, 0.01f0), parallelization_backend=Main.parallelization_backend) - semi_new = sol.prob.p + semi_new = sol.prob.p.semi @testset verbose=true "Line" begin # Interpolation parameters diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index cd482ec257..ea676cf923 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -1,8 +1,8 @@ # Use `@trixi_testset` to isolate the mock functions in a separate namespace @trixi_testset "Semidiscretization" begin # Mock systems - struct System1 <: TrixiParticles.AbstractSystem{3} end - struct System2 <: TrixiParticles.AbstractSystem{3} end + struct System1 <: TrixiParticles.AbstractStructureSystem{3} end + struct System2 <: TrixiParticles.AbstractStructureSystem{3} end system1 = System1() system2 = System2() @@ -10,9 +10,9 @@ Base.eltype(::System1) = Float64 TrixiParticles.coordinates_eltype(::System1) = Float32 TrixiParticles.u_nvariables(::System1) = 3 - TrixiParticles.u_nvariables(::System2) = 4 + TrixiParticles.u_nvariables(::System2) = 3 TrixiParticles.v_nvariables(::System1) = 3 - TrixiParticles.v_nvariables(::System2) = 2 + TrixiParticles.v_nvariables(::System2) = 4 TrixiParticles.nparticles(::System1) = 2 TrixiParticles.nparticles(::System2) = 3 TrixiParticles.n_integrated_particles(::System1) = 2 @@ -25,8 +25,8 @@ semi = Semidiscretization(system1, system2, neighborhood_search=nothing) # Verification - @test semi.ranges_u == (1:6, 7:18) - @test semi.ranges_v == (1:6, 7:12) + @test semi.ranges_u == (1:6, 7:15) + @test semi.ranges_v == (1:6, 7:18) nhs = ((TrixiParticles.TrivialNeighborhoodSearch{3}(search_radius=0.2, eachpoint=1:2), @@ -145,14 +145,27 @@ end @testset verbose=true "Source Terms" begin + function Base.getproperty(::System1, name::Symbol) + if name == :acceleration + return SVector(0.0, 0.0, 0.0) + end + error("property $name not defined for `System1`") + end + function Base.getproperty(::System2, name::Symbol) + if name == :acceleration + return SVector(0.0, 0.0, 1.0) + end + error("property $name not defined for `System2`") + end TrixiParticles.source_terms(::System1) = SourceTermDamping(damping_coefficient=0.1) + TrixiParticles.source_terms(::System2) = nothing TrixiParticles.current_density(v, system::System1, particle) = 0.0 TrixiParticles.current_pressure(v, system::System1, particle) = 0.0 semi = Semidiscretization(system1, system2, neighborhood_search=nothing) - dv_ode = zeros(3 * 2 + 2 * 3) - du_ode = zeros(3 * 2 + 4 * 3) + dv_ode = zeros(3 * 2 + 4 * 3) + du_ode = zeros(3 * 2 + 3 * 3) u_ode = zero(du_ode) v1 = [1.0 2.0 @@ -167,6 +180,31 @@ @test dv1 == -0.1 * v1 dv2 = TrixiParticles.wrap_v(dv_ode, system2, semi) - @test iszero(dv2) + @test dv2 == vcat(zeros(2, 3), ones(1, 3), zeros(1, 3)) + end + + @testset verbose=true "drift!" begin + semi = Semidiscretization(system1, system2, neighborhood_search=nothing) + + du_ode = fill(NaN, 3 * 2 + 3 * 3) + u_ode = zeros(3 * 2 + 3 * 3) + + v1 = [1.0 2.0 + 3.0 4.0 + 5.0 6.0] + v2 = [10.0 11.0 12.0 + 20.0 21.0 22.0 + 30.0 31.0 32.0 + -999.0 -999.0 -999.0] + v_ode = vcat(vec(v1), vec(v2)) + + returned = TrixiParticles.drift!(du_ode, v_ode, u_ode, semi, 0.0) + @test returned === du_ode + + du1 = TrixiParticles.wrap_u(du_ode, system1, semi) + @test du1 == v1 + + du2 = TrixiParticles.wrap_u(du_ode, system2, semi) + @test du2 == v2[1:3, :] end end diff --git a/test/general/util.jl b/test/general/util.jl index a7e11da4a1..4aefddd4e3 100644 --- a/test/general/util.jl +++ b/test/general/util.jl @@ -19,6 +19,16 @@ @test typeof(A .+ B) == typeof(A) @test typeof(B .+ A) == typeof(A) + # Test that the resulting type of `similar` is correct + C = similar(A, Float64, (2, 2)) + @test typeof(C) == typeof(A) + @test size(C) == (2, 2) + @test typeof(similar(A, Float64)) == typeof(A) + C = similar(A, (2, 2)) + @test typeof(C) == typeof(A) + @test size(C) == (2, 2) + @test typeof(similar(A)) == typeof(A) + # Test that these operations all use the correct backend struct FailingBackend end diff --git a/test/schemes/structure/total_lagrangian_sph/rhs.jl b/test/schemes/structure/total_lagrangian_sph/rhs.jl index 69f4a920c4..6383f6977e 100644 --- a/test/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/test/schemes/structure/total_lagrangian_sph/rhs.jl @@ -93,9 +93,6 @@ TrixiParticles.each_integrated_particle(::MockSystem) = each_integrated_particle TrixiParticles.smoothing_length(::MockSystem, _) = eps() - function TrixiParticles.add_acceleration!(_, _, ::MockSystem) - return nothing - end TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) = kernel_deriv TrixiParticles.compact_support(::MockSystem, ::MockSystem) = 1000.0 function TrixiParticles.current_coords(system::MockSystem, particle) diff --git a/test/validation/validation.jl b/test/validation/validation.jl index f295e6e1b2..dad439143c 100644 --- a/test/validation/validation.jl +++ b/test/validation/validation.jl @@ -5,7 +5,7 @@ "investigate_relaxation.jl"), tspan=(0.0, 1.0)) @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 # Verify number of plots @test plot1.n == 4 end @@ -21,9 +21,9 @@ if VERSION < v"1.12" # Older Julia versions produce allocations because `get_neighborhood_search` # is not type-stable with TLSPH. - @test count_rhs_allocations(sol, semi) < 200 + @test count_rhs_allocations(sol) < 200 else - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @test isapprox(error_deflection_x, 0, atol=eps()) @test isapprox(error_deflection_y, 0, atol=eps()) @@ -59,7 +59,7 @@ r"WARNING: Method definition interpolated_pressure.*\n" ] @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 # Note that pressure values are in the order of 1e5 @test isapprox(error_wcsph_P1, 0, atol=eps(1e5)) @@ -104,7 +104,7 @@ r"WARNING: Method definition initial_velocity_function.*\n" ] @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end @trixi_testset "LDC_2D" begin @@ -118,6 +118,6 @@ r"WARNING: Method definition is_moving.*\n" ] @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + @test count_rhs_allocations(sol) == 0 end end