diff --git a/NEWS.md b/NEWS.md index 857a55ef60..5ade32aba8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,9 +4,16 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Version 0.5.2 + +### Features +- Added the ability to restart simulations from VTK solution files generated by `SolutionSavingCallback`. + Users can now pass the `restart_with` keyword argument to `semidiscretize` (#1190). + ## Version 0.5.1 ### Features + - Implemented stage-level coupling for split integration (#1049). - Added `MechanicalWorkCalculatorCallback` (#940). diff --git a/docs/src/callbacks.md b/docs/src/callbacks.md index c96ccc7dba..e356917226 100644 --- a/docs/src/callbacks.md +++ b/docs/src/callbacks.md @@ -14,3 +14,13 @@ The following pre-defined custom quantities can be used with the Modules = [TrixiParticles] Pages = ["general/custom_quantities.jl"] ``` + +# Mechanical Work Calculator + +The `MechanicalWorkCalculator` is a special custom quantity to be used with the +[`PostprocessCallback`](@ref). + +```@autodocs +Modules = [TrixiParticles] +Pages = ["general/mechanical_work_calculator.jl"] +``` 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/fin_2d.jl b/examples/fsi/fin_2d.jl new file mode 100644 index 0000000000..a7359700de --- /dev/null +++ b/examples/fsi/fin_2d.jl @@ -0,0 +1,458 @@ +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.502 +fin_thickness = 30e-3 +real_modulus = 125e9 +poisson_ratio = 0.3 + +# Real blade thickness profile along the flexible blade: +# x = 0 is the attachment to the foot pocket, x = 1 is the blade tip. +function real_thickness(x) + real_thickness_at_attachment = 1.2e-3 + real_thickness_at_tip = 0.7e-3 + + # Clamp to use constant material properties for the clamped region and foot pocket. + x_clamped = clamp(x, 0.0, 1.0) + return real_thickness_at_attachment + + x_clamped * (real_thickness_at_tip - real_thickness_at_attachment) +end + +# The simulated blade is artificially thick for resolution. Scale Young's modulus so +# E_artificial * t_artificial^3 matches the real bending stiffness E_real * t_real^3. +function artificial_modulus(real_modulus, real_thickness, artificial_thickness) + return real_modulus * (real_thickness / artificial_thickness)^3 +end + +# Scale density so rho_artificial * t_artificial keeps the same mass per blade area +# as rho_real * t_real. +function artificial_density(real_density, real_thickness, artificial_thickness) + return real_density * real_thickness / artificial_thickness +end + +fiber_volume_fraction = 0.6 +fiber_density = 1800.0 +epoxy_density = 1250.0 +real_density = fiber_volume_fraction * fiber_density + + (1 - fiber_volume_fraction) * epoxy_density +density = real_density + +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", "hyper_bifins_x.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 + min_corner .-= center + max_corner .-= center + cell_list = FullGridCellList(; min_corner, max_corner) + neighborhood_search = GridNeighborhoodSearch{2}(; 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(); + abstol=1e-8, + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing_foot"), + 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 + + # `union(packed_foot, beam)`, but when particles are too close together, keep the ones + # from `beam` instead of `packed_foot` to ensure that the blade doesn't have holes. + structure = union(setdiff(packed_foot, beam), beam) + 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) + + min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 + cell_list = FullGridCellList(; min_corner, max_corner) + neighborhood_search = GridNeighborhoodSearch{2}(; cell_list, update_strategy=ParallelUpdate()) + + 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()), + abstol=1e-8, + dtmax=1e-2) + + fluid = InitialCondition(sol_packing, fluid_packing_system, semi_packing) + fluid = union(fluid, fixed_fluid) +else + structure = union(fixed_particles, beam) + # Move the fin to the center of the tank + structure.coordinates .+= center + + fluid = setdiff(tank.fluid, structure) +end + +# Convert particle x-positions to the relative blade coordinate used by `real_thickness`. +# A value of 0 corresponds to the blade attachment, and a value of 1 corresponds to the tip. +function normalized_blade_coordinate(coordinates, particle) + return (coordinates[1, particle] - center[1]) / fin_length +end + +real_thickness_structure = [real_thickness(normalized_blade_coordinate(structure.coordinates, + particle)) + for particle in 1:nparticles(structure)] + +modulus = [artificial_modulus(real_modulus, thickness, fin_thickness) + for thickness in real_thickness_structure] + +# Update both density and mass based on the artificial thickness to ensure that +# the structure has the same mass per blade area as the real fin, +# while keeping the same bending stiffness. +structure.density .= [artificial_density(real_density, thickness, fin_thickness) + for thickness in real_thickness_structure] +structure.mass .= structure.density .* particle_spacing^2 + +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) + +viscosity_structure = ArtificialViscosityMonaghan(alpha=1.0) +structure_system = TotalLagrangianSPHSystem(structure; smoothing_kernel, smoothing_length=smoothing_length_structure, + young_modulus=modulus, poisson_ratio, + clamped_particles=1:n_clamped_particles, + clamped_particles_motion=boundary_motion, + boundary_model=boundary_model_structure, + velocity_averaging=TrixiParticles.VelocityAveraging(time_constant=5e-4), + viscosity=viscosity_structure, + penalty_force=PenaltyForceGanzenmueller(alpha=0.1)) + +# ========================================================================================== +# ==== Fluid +fluid_density_calculator = ContinuityDensity() +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(fluid; density_calculator=fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length=smoothing_length_fluid, + viscosity=viscosity_fluid, + 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) + +# ========================================================================================== +# ==== 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.5 +# 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") + +mechanical_work_calculator = MechanicalWorkCalculator(semi.systems[4], semi) +thrust_calculator = ThrustCalculator(semi.systems[4], semi, direction=SVector(1.0, 0.0)) +calculator_cb = PostprocessCallback(; mechanical_work_calculator, thrust_calculator, + interval=100, filename="$(prefix)_efficiency", + write_file_interval=1000) + +callbacks = CallbackSet(info_callback, saving_callback, + stepsize_callback, split_integration, pp_cb, interpolate_cb, + calculator_cb, + UpdateCallback(), SortingCallback(interval=10_000)) + +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/postprocessing/restart_oscillating_beam_2d.jl b/examples/postprocessing/restart_oscillating_beam_2d.jl new file mode 100644 index 0000000000..51080ec265 --- /dev/null +++ b/examples/postprocessing/restart_oscillating_beam_2d.jl @@ -0,0 +1,27 @@ +# ========================================================================================== +# Restart Example: Oscillating Beam 2D +# +# This example demonstrates how to restart a simulation. +# We first run a simulation of oscillating beam up to t=1.0s, then restart from the +# saved state and continue the simulation until t=2.0s. +# ========================================================================================== +using TrixiParticles + +trixi_include(@__MODULE__, + joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), + tspan=(0.0, 1.0)) + +# Get latest iteration +iter = saving_callback.affect!.affect!.latest_saved_iter + +restart_file = joinpath("out", "structure_1_$iter.vtu") + +ode_restart = semidiscretize(semi, (1.0, 2.0); + restart_with=restart_file) + +saving_callback = SolutionSavingCallback(dt=0.02, prefix="restart") + +callbacks = CallbackSet(info_callback, saving_callback) + +sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1e-5, reltol=1e-3, dtmax=1e-2, + save_everystep=false, callback=callbacks) diff --git a/examples/postprocessing/restart_poiseuille_flow_2d.jl b/examples/postprocessing/restart_poiseuille_flow_2d.jl new file mode 100644 index 0000000000..4564220805 --- /dev/null +++ b/examples/postprocessing/restart_poiseuille_flow_2d.jl @@ -0,0 +1,31 @@ +# ========================================================================================== +# Restart Example: Poiseuille Flow 2D +# +# This example demonstrates how to restart a simulation. +# We first run a simulation of 2D Poiseuille flow up to t=0.3s, then restart from the +# saved state and continue the simulation until t=0.6s. +# ========================================================================================== +using TrixiParticles + +trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"), + tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5) + +# Get latest iteration +iter = saving_callback.affect!.affect!.latest_saved_iter + +restart_file_fluid = joinpath("out", "fluid_1_$iter.vtu") +restart_file_open_boundary = joinpath("out", "open_boundary_1_$iter.vtu") +restart_file_boundary = joinpath("out", "boundary_1_$iter.vtu") + +ode_restart = semidiscretize(semi, (0.3, 0.6); + restart_with=(restart_file_fluid, + restart_file_open_boundary, + restart_file_boundary)) + +saving_callback = SolutionSavingCallback(dt=0.02, prefix="restart") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1e-5, reltol=1e-3, dtmax=1e-2, + save_everystep=false, callback=callbacks) diff --git a/examples/preprocessing/data/hyper_bifins_x.dxf b/examples/preprocessing/data/hyper_bifins_x.dxf new file mode 100644 index 0000000000..87ffef8728 --- /dev/null +++ b/examples/preprocessing/data/hyper_bifins_x.dxf @@ -0,0 +1,5218 @@ + 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 +2461208.724513889 + 9 +$TDUCREATE + 40 +2458532.153996898 + 9 +$TDUPDATE + 40 +2461208.724513889 + 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 +94 + 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 +4D9D8938-6997-11F1-A0CF-A088C2115B2C + 9 +$VERSIONGUID + 2 +4D9DEE3C-6997-11F1-A0CF-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 +ACDBPLACEHOLDER + 2 +AcDbPlaceHolder + 3 +ObjectDBX Classes + 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 +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 +93 +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 +0.0 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +30 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.005623843252542472 + 20 +4.8636960090447824e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +31 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.012585488133344366 + 20 +9.214060349372237e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +32 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.02061691015601784 + 20 +0.00012974617782760488 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +33 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.02945008483417531 + 20 +0.00016068893070988926 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +34 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03881698768142893 + 20 +0.00018420410975834403 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +35 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.048449594211391014 + 20 +0.00019952696259076315 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +36 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.05807987993767382 + 20 +0.0002058927368249151 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +37 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06743982037388971 + 20 +0.00020253668007860656 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +38 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0762613910336508 + 20 +0.00018869403996961873 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +39 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0842765674305695 + 20 +0.00016360006411573287 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.091217325078258 + 20 +0.00012649000013471745 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09681563949032863 + 20 +7.659909564437914e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10080348618039366 + 20 +1.3162598262486442e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10367157245296806 + 20 +-1.17394308284186e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10800875657294529 + 20 +-0.00020302270660931183 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11378660650741865 + 20 +-0.0005710572626687967 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +40 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12080334462153497 + 20 +-0.001096656544127838 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +41 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1288571932804413 + 20 +-0.0017606339961074389 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +42 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1377463748492844 + 20 +-0.002543803063728539 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +43 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14726911169321127 + 20 +-0.0034269771921121408 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +44 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.15722362617736882 + 20 +-0.004390969826379197 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +45 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.16740814066690388 + 20 +-0.005416594411650697 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +46 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.17762087752696346 + 20 +-0.006484664393047595 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +47 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18766005912269434 + 20 +-0.00757599321569088 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +48 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19732390781924355 + 20 +-0.008671394324701528 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +49 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20641064598175785 + 20 +-0.009751681165200496 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21471849597538425 + 20 +-0.010797667182308755 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22204568016526965 + 20 +-0.011790165821147313 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2281904209165609 + 20 +-0.01270999052683711 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23295094059440502 + 20 +-0.013537954744499132 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23979996360677527 + 20 +-0.015129520021566228 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.25148260503740166 + 20 +-0.017606857022919985 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +50 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2579546191528017 + 20 +-0.01863173176617797 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +51 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.262226860215378 + 20 +-0.019412214991891653 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +52 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2669374459702987 + 20 +-0.020999484552667723 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +53 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.27060990263285706 + 20 +-0.02168528748379199 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +54 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.27605391151655334 + 20 +-0.02187701263418331 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +55 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28290739393172737 + 20 +-0.021542503621769647 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +56 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29080827118871905 + 20 +-0.02064960406447897 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +57 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2993944645978683 + 20 +-0.019166157580239235 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +58 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3083038954695148 + 20 +-0.017060007786978408 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +59 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3171744851139986 + 20 +-0.014298998302624445 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.32564415484165954 + 20 +-0.01085097274510531 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3333508259628373 + 20 +-0.006683774732348999 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.339932419787872 + 20 +-0.0017652478822834355 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3454881558252771 + 20 +0.004704833462746747 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3483272761239312 + 20 +0.011776135399291367 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3484514218784382 + 20 +0.018861860410078025 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +60 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3463546852158496 + 20 +0.025802732601637842 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +61 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.342531158263217 + 20 +0.03243947608050197 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +62 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.337474933147592 + 20 +0.038612814953201464 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +63 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3316801019960261 + 20 +0.04416347332626759 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +64 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.32564075693557093 + 20 +0.04893217530623134 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +65 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3188458778589197 + 20 +0.05341221504874015 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +66 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3108571163669807 + 20 +0.058297362776538256 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +67 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.3044121089074911 + 20 +0.06182324367359262 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +68 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2985118429121997 + 20 +0.06441062011468514 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +69 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2921573058128556 + 20 +0.06648025447459761 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2865282882030938 + 20 +0.0679285676863985 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2832559809362735 + 20 +0.068596028695872 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.27901929595517305 + 20 +0.06934348373234978 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2738182315326241 + 20 +0.07017093288219028 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.26765278594145847 + 20 +0.07107837623175192 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.26052295745450776 + 20 +0.07206581386739311 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +70 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.25242874434460383 + 20 +0.07313324587547225 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +71 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2433701448845784 + 20 +0.07428067234234775 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +72 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23334715734726305 + 20 +0.07550809335437807 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +73 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22235978000548964 + 20 +0.07681550899792156 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +74 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2104080111320898 + 20 +0.07820291935933668 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +75 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19771046199064185 + 20 +0.0790226280364367 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +76 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.185841662792443 + 20 +0.07936991036791603 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +77 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1750579172700254 + 20 +0.07960915554997243 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +78 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.16535923085392507 + 20 +0.07974036358803642 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +79 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.15674560897467776 + 20 +0.07976353448753852 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14921705706281949 + 20 +0.07967866825390929 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.142773580548886 + 20 +0.07948576489257927 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13741518486341311 + 20 +0.07918482440897898 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13314187543693684 + 20 +0.07877584680853897 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12995365769999287 + 20 +0.07825883209668975 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1254877079691722 + 20 +0.07684097193913388 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +80 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11873677349468834 + 20 +0.07384188739430268 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +81 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11106532790527829 + 20 +0.06986682283410413 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +82 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10309967970240826 + 20 +0.0653738215810658 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +83 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09546613738754436 + 20 +0.06082092695771521 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +84 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08879100946215296 + 20 +0.05666618228657999 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +85 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08370060442770015 + 20 +0.05336763089018767 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +86 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08160927675386886 + 20 +0.05197399227539537 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +87 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07927815509870763 + 20 +0.050450872590407314 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +88 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07613436097711551 + 20 +0.04841720074061562 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +89 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07217789696133127 + 20 +0.04587297621157261 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06740876562359346 + 20 +0.0428181984888305 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06182696953614077 + 20 +0.03925286705794162 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.055432511271211765 + 20 +0.03517698140445818 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04822539340104519 + 20 +0.030590541013932528 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04020561849787964 + 20 +0.025493545371916906 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03137318913395361 + 20 +0.019885993963963525 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +90 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.02172810788150599 + 20 +0.013767886275624787 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +91 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.01127037731277512 + 20 +0.007139221792452783 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +92 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +1.0181425078340335e-16 + 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/src/TrixiParticles.jl b/src/TrixiParticles.jl index 0bc56370f2..ace8012fc2 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -62,6 +62,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("preprocessing/preprocessing.jl") include("io/io.jl") +include("general/restart.jl") include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! @@ -73,8 +74,8 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, - SplitIntegrationCallback, MechanicalWorkCalculatorCallback, - calculated_mechanical_work, + SplitIntegrationCallback, MechanicalWorkCalculator, calculated_mechanical_work, + ThrustCalculator, calculated_thrust, SortingCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 3f295a6bc8..de8b09986a 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -76,5 +76,4 @@ include("stepsize.jl") include("update.jl") include("steady_state_reached.jl") include("split_integration.jl") -include("mechanical_work_calculator.jl") include("sorting.jl") diff --git a/src/callbacks/info.jl b/src/callbacks/info.jl index 45c60a4a79..3134247144 100644 --- a/src/callbacks/info.jl +++ b/src/callbacks/info.jl @@ -1,5 +1,7 @@ mutable struct InfoCallback start_time::Float64 + last_performance_print_time::Float64 + ncalls_at_last_print::Int interval::Int end @@ -33,7 +35,7 @@ beginning of a simulation and then resets the timer. When the returned callback directly, the current timer values are shown. """ function InfoCallback(; interval=0, reset_threads=true) - info_callback = InfoCallback(0.0, interval) + info_callback = InfoCallback(0.0, 0.0, 0, interval) function initialize(cb, u, t, integrator) initialize_info_callback(cb, u, t, integrator; @@ -65,6 +67,17 @@ function (info_callback::InfoCallback)(integrator) integrator.stats.naccept, integrator.dt, t, sim_time_percentage), 71) * @sprintf("│ run time: %.4e s", runtime_absolute)) + + if condition_integrator_interval(integrator, 1 * info_callback.interval) + runtime_since_last_print = (time_ns() - info_callback.last_performance_print_time) + info_callback.last_performance_print_time = time_ns() + n_calls = integrator.stats.nf - info_callback.ncalls_at_last_print + info_callback.ncalls_at_last_print = integrator.stats.nf + pid = runtime_since_last_print / n_calls / sum(nparticles.(integrator.p.semi.systems)) + println("─"^100) + println(@sprintf("time/particle/rhs: %.4e ns", pid)) + println("─"^100) + end end # This callback only reports progress and does not change the result of the right-hand side. @@ -143,6 +156,8 @@ function initialize_info_callback(discrete_callback, u, t, integrator; # Save current time as start_time info_callback.start_time = time_ns() + info_callback.last_performance_print_time = time_ns() + info_callback.ncalls_at_last_print = 0 return nothing end @@ -256,9 +271,16 @@ end function print_summary(integrator) println("─"^100) - println("Trixi simulation finished. Final time: ", integrator.t, - " Time steps: ", integrator.stats.naccept, " (accepted), ", - integrator.iter, " (total)") + # Print the final time as string to let Julia decide on the formatting + @printf(" %-31s %10s\n", "Final time:", string(integrator.t)) + @printf(" %-31s %10d (accepted) %10d (total)\n", + "Time steps:", integrator.stats.naccept, integrator.iter) + if !isnothing(integrator.p.split_integration_data) + @printf(" %-31s %10d (accepted) %10d (total)\n", + "Split integration time steps:", + integrator.p.split_integration_data.integrator.stats.naccept, + integrator.p.split_integration_data.integrator.iter) + end println("─"^100) println() diff --git a/src/callbacks/mechanical_work_calculator.jl b/src/callbacks/mechanical_work_calculator.jl deleted file mode 100644 index f20672af09..0000000000 --- a/src/callbacks/mechanical_work_calculator.jl +++ /dev/null @@ -1,243 +0,0 @@ -""" - MechanicalWorkCalculatorCallback(system::TotalLagrangianSPHSystem, semi; interval=1, - eachparticle=(n_integrated_particles(system) + 1):nparticles(system), - only_compute_force_on_fluid=false) - -Callback that accumulates the work done by a set of particles in a -[`TotalLagrangianSPHSystem`](@ref) by integrating the instantaneous power over time. - -With the default arguments it tracks the work done by the clamped particles -that follow a [`PrescribedMotion`](@ref). By selecting a different particle set, it can also -be used to measure the work done by the structure on the surrounding fluid. - -- **Prescribed/clamped motion work** (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_mechanical_work`](@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 instantaneous power. - It is recommended to keep this at `1` (every time step) or small (≤ 5) - to limit time integration errors in the 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 = TotalLagrangianSPHSystem(RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel=WendlandC2Kernel{2}(), smoothing_length=1.0, young_modulus=1.0, poisson_ratio=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) -semi = Semidiscretization(system) -ode = semidiscretize(semi, (0.0, 1.0)) - -# Note that `Semidiscretization` might create a deep copy of the system, -# which means we have to extract the new system from `semi`. -# When working with GPUs, `semidiscretize` also creates a deep copy of `semi` and another -# copy of the system, so the clean way to get the correct new system is this: -semi_new = ode.p.semi -system_new = semi_new.systems[1] - -# Create a mechanical work calculator callback that is called every 2 time steps -mechanical_work_cb = MechanicalWorkCalculatorCallback(system_new, semi_new; interval=2) - -# After the simulation, retrieve the accumulated mechanical work -mechanical_work = calculated_mechanical_work(mechanical_work_cb) - -# output -[ Info: To create the self-interaction neighborhood search of a `TotalLagrangianSPHSystem`, a deep copy of the system is created inside the `Semidiscretization`. Use `system = semi.systems[i]` to access simulation data. -0.0 -``` -""" -struct MechanicalWorkCalculatorCallback{T, DV, EP} - interval :: Int - t :: T # Time of last call - work :: T - system_index :: Int - dv :: DV - eachparticle :: EP - only_compute_force_on_fluid :: Bool -end - -# This should dispatch on `TotalLagrangianSPHSystem`, but this name is not yet defined -# due to the include order. -function MechanicalWorkCalculatorCallback(system::AbstractStructureSystem, 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, - (v_nvariables(system), nparticles(system))) - - # Note that time and work are initialized in - # `initialize_mechanical_work_calculator_callback`. - cb = MechanicalWorkCalculatorCallback(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), - initialize=initialize_mechanical_work_calculator_callback) -end - -function initialize_mechanical_work_calculator_callback(discrete_callback, u, t, integrator) - work_callback = discrete_callback.affect! - - # Reset time and mechanical work - work_callback.t[] = t - work_callback.work[] = zero(eltype(work_callback.work)) -end - -""" - calculated_mechanical_work(cb::DiscreteCallback{<:Any, <:MechanicalWorkCalculatorCallback}) - -Get the accumulated mechanical work from a [`MechanicalWorkCalculatorCallback`](@ref). - -# Arguments -- `cb`: The `DiscreteCallback` returned by [`MechanicalWorkCalculatorCallback`](@ref). - -# Examples -```jldoctest; output = false, setup = :(system = TotalLagrangianSPHSystem(RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel=WendlandC2Kernel{2}(), smoothing_length=1.0, young_modulus=1.0, poisson_ratio=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) -# Create a mechanical work calculator callback -mechanical_work_cb = MechanicalWorkCalculatorCallback(system, semi) - -# After the simulation, retrieve the accumulated mechanical work -mechanical_work = calculated_mechanical_work(mechanical_work_cb) - -# output -0.0 -``` -""" -function calculated_mechanical_work(cb::DiscreteCallback{<:Any, - <:MechanicalWorkCalculatorCallback}) - return cb.affect!.work[] -end - -# `condition` -function (callback::MechanicalWorkCalculatorCallback)(u, t, integrator) - (; interval) = callback - - return condition_integrator_interval(integrator, interval) -end - -# `affect!` -function (callback::MechanicalWorkCalculatorCallback)(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.semi - v_ode, u_ode = integrator.u.x - work = callback.work - - system = semi.systems[callback.system_index] - update_mechanical_work_calculator!(work, system, callback.eachparticle, - callback.only_compute_force_on_fluid, callback.dv, - v_ode, u_ode, semi, t, dt) - - # This callback only processes results and does not change the result of the right-hand side. - derivative_discontinuity!(integrator, false) - - return integrator -end - -function update_mechanical_work_calculator!(work, system, eachparticle, - only_compute_force_on_fluid, dv, - v_ode, u_ode, semi, t, dt) - @trixi_timeit timer() "calculate mechanical work" 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, system, particle) - 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 mechanical work, 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. - work[] += dot(-F_particle, velocity) * dt - end - end - - return work -end - -function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:MechanicalWorkCalculatorCallback}) - @nospecialize cb # reduce precompilation time - - ELTYPE = eltype(cb.affect!.work) - print(io, "MechanicalWorkCalculatorCallback{$ELTYPE}(interval=", cb.affect!.interval, - ")") -end - -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, <:MechanicalWorkCalculatorCallback}) - @nospecialize cb # reduce precompilation time - - if get(io, :compact, false) - show(io, cb) - else - update_cb = cb.affect! - ELTYPE = eltype(update_cb.work) - setup = [ - "interval" => update_cb.interval - ] - summary_box(io, "MechanicalWorkCalculatorCallback{$ELTYPE}", setup) - end -end diff --git a/src/general/general.jl b/src/general/general.jl index 3e9c868507..0615c8124f 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -8,4 +8,5 @@ include("initial_condition.jl") include("buffer.jl") include("interpolation.jl") include("custom_quantities.jl") +include("mechanical_work_calculator.jl") include("time_integration.jl") diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index b78ae8013f..d49553bf47 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -95,7 +95,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 116388753b..cb24e3079b 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 @@ -177,6 +177,10 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ vtk["density"] = density vtk["velocity"] = velocity vtk["pressure"] = pressure + + if !isnothing(pvd) && t >= 0 + pvd[t] = vtk + end end end diff --git a/src/general/mechanical_work_calculator.jl b/src/general/mechanical_work_calculator.jl new file mode 100644 index 0000000000..8dd78a62ec --- /dev/null +++ b/src/general/mechanical_work_calculator.jl @@ -0,0 +1,335 @@ +""" + MechanicalWorkCalculator(system::TotalLagrangianSPHSystem, semi; + eachparticle=(n_integrated_particles(system) + 1):nparticles(system), + only_compute_force_on_fluid=false) + +Functor that accumulates the work done by a set of particles in a +[`TotalLagrangianSPHSystem`](@ref) by integrating the instantaneous power over time. +It can be passed as a custom quantity to [`PostprocessCallback`](@ref), which controls +when the work is sampled and whether it is written to a file. + +With the default arguments it tracks the work done by the clamped particles +that follow a [`PrescribedMotion`](@ref). By selecting a different particle set, it can also +be used to measure the work done by the structure on the surrounding fluid. + +- **Prescribed/clamped motion work** (default) -- monitor only the clamped particles by + leaving `eachparticle` at its default range + `(n_integrated_particles(system) + 1):nparticles(system)`. +- **Fluid energy transfer** -- 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. + +Internally the calculator 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_mechanical_work`](@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 +- `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. + +# Examples +```jldoctest; output = false, setup = :(system = TotalLagrangianSPHSystem(RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel=WendlandC2Kernel{2}(), smoothing_length=1.0, young_modulus=1.0, poisson_ratio=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) +semi = Semidiscretization(system) +ode = semidiscretize(semi, (0.0, 1.0)) + +# Note that `Semidiscretization` might create a deep copy of the system, +# which means we have to extract the new system from `semi`. +# When working with GPUs, `semidiscretize` also creates a deep copy of `semi` and another +# copy of the system, so the clean way to get the correct new system is this: +semi_new = ode.p.semi +system_new = semi_new.systems[1] + +# Create a mechanical work calculator that is called every 2 time steps. +mechanical_work_calculator = MechanicalWorkCalculator(system_new, semi_new) +postprocess_cb = PostprocessCallback(; interval=2, mechanical_work_calculator) + +# After the simulation, retrieve the accumulated mechanical work. +mechanical_work = calculated_mechanical_work(mechanical_work_calculator) + +# output +[ Info: To create the self-interaction neighborhood search of a `TotalLagrangianSPHSystem`, a deep copy of the system is created inside the `Semidiscretization`. Use `system = semi.systems[i]` to access simulation data. +0.0 +``` +""" +mutable struct MechanicalWorkCalculator{ELTYPE, DV, EP} + t :: ELTYPE # Time of last call + work :: ELTYPE + initialized :: Bool + system_index :: Int + dv :: DV + eachparticle :: EP + only_compute_force_on_fluid :: Bool +end + +# This should dispatch on `TotalLagrangianSPHSystem`, but this name is not yet defined +# due to the include order. +function MechanicalWorkCalculator(system::AbstractStructureSystem, semi; + 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). + # `PostprocessCallback` transfers data to the CPU before calling custom quantities, + # so this can be a regular `Array` even when the simulation is running on a GPU. + dv = Array{ELTYPE, 2}(undef, (v_nvariables(system), nparticles(system))) + + return MechanicalWorkCalculator(zero(ELTYPE), zero(ELTYPE), false, + system_index, dv, eachparticle, + only_compute_force_on_fluid) +end + +function reset!(calculator::MechanicalWorkCalculator) + calculator.t = zero(calculator.t) + calculator.work = zero(calculator.work) + calculator.initialized = false + + return calculator +end + +""" + calculated_mechanical_work(calculator::MechanicalWorkCalculator) + +Get the accumulated mechanical work from a [`MechanicalWorkCalculator`](@ref). + +# Arguments +- `calculator`: The [`MechanicalWorkCalculator`](@ref) functor. + +# Examples +```jldoctest; output = false, setup = :(system = TotalLagrangianSPHSystem(RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel=WendlandC2Kernel{2}(), smoothing_length=1.0, young_modulus=1.0, poisson_ratio=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) +# Create a mechanical work calculator +mechanical_work_calculator = MechanicalWorkCalculator(system, semi) + +# After the simulation, retrieve the accumulated mechanical work +mechanical_work = calculated_mechanical_work(mechanical_work_calculator) + +# output +0.0 +``` +""" +function calculated_mechanical_work(calculator::MechanicalWorkCalculator) + return calculator.work +end + +function (calculator::MechanicalWorkCalculator)(system, dv_ode, du_ode, v_ode, u_ode, + semi, t) + if system_indices(system, semi) != calculator.system_index + return nothing + end + + if !calculator.initialized + calculator.t = t + calculator.work = zero(calculator.work) + calculator.initialized = true + return calculator.work + end + + dt = t - calculator.t + calculator.t = t + + @trixi_timeit timer() "calculate mechanical work" begin + calculator.work = update_mechanical_work(calculator.work, system, + calculator.eachparticle, + calculator.only_compute_force_on_fluid, + calculator.dv, + v_ode, u_ode, semi, t, dt) + end + + return calculator.work +end + +function update_mechanical_work(work, system, eachparticle, + only_compute_force_on_fluid, dv, + v_ode, u_ode, semi, t, dt) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + compute_structure_acceleration!(dv, v, u, system, eachparticle, + only_compute_force_on_fluid, + v_ode, u_ode, semi, t) + + # Note that this is a reduction, so we cannot use `@threaded` here. + 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 mechanical work, 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. + work += dot(-F_particle, velocity) * dt + end + + return work +end + +function compute_structure_acceleration!(dv, v, u, system, eachparticle, + only_compute_force_on_fluid, + v_ode, u_ode, semi, t) + set_zero!(dv) + + foreach_system(semi) do neighbor_system + if only_compute_force_on_fluid && !(neighbor_system isa AbstractFluidSystem) + # Not a fluid system, ignore this system. + return nothing + 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) + + return nothing + end + + if !only_compute_force_on_fluid + @threaded semi for particle in eachparticle + add_acceleration!(dv, system, particle) + add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t) + end + end + + return dv +end + +""" + ThrustCalculator(system::TotalLagrangianSPHSystem, semi; + direction, + eachparticle=eachparticle(system)) + +Functor that computes the instantaneous hydrodynamic force exerted by the fluid on a +[`TotalLagrangianSPHSystem`](@ref), projected onto `direction`. +It can be passed as a custom quantity to [`PostprocessCallback`](@ref). + +For a fixed fluid-interacting structure in a channel flow, choose `direction` as the +direction of useful force and multiply the recorded thrust by the corresponding reference +speed to obtain useful mechanical power. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +# Arguments +- `system`: The [`TotalLagrangianSPHSystem`](@ref) whose hydrodynamic force should be + monitored. +- `semi`: The [`Semidiscretization`](@ref) that contains `system`. + +# Keywords +- `direction`: Direction onto which the hydrodynamic force is projected. The vector is + normalized internally. +- `eachparticle=eachparticle(system)`: Iterator selecting which particles contribute. + +# Examples +```jldoctest; output = false, setup = :(system = TotalLagrangianSPHSystem(RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel=WendlandC2Kernel{2}(), smoothing_length=1.0, young_modulus=1.0, poisson_ratio=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) +# Create a thrust calculator in x-direction +thrust_calculator = ThrustCalculator(system, semi; direction=SVector(1.0, 0.0)) + +# After postprocessing, retrieve the latest thrust value +thrust = calculated_thrust(thrust_calculator) + +# output +0.0 +``` +""" +mutable struct ThrustCalculator{ELTYPE, DV, EP, D} + thrust :: ELTYPE + system_index :: Int + dv :: DV + eachparticle :: EP + direction :: D +end + +# This should dispatch on `TotalLagrangianSPHSystem`, but this name is not yet defined +# due to the include order. +function ThrustCalculator(system::AbstractStructureSystem, semi; + direction, eachparticle=eachparticle(system)) + ELTYPE = eltype(system) + system_index = system_indices(system, semi) + + # Check vector length before converting to `SVector` to avoid extremely long + # compile times when accidentally passing large vectors. + if length(direction) != ndims(system) + throw(ArgumentError("length of `direction` must match the number of dimensions")) + end + direction_ = SVector(Tuple(direction)) + if iszero(direction_) + throw(ArgumentError("`direction` must not be zero")) + end + direction_ = normalize(direction_) + + # Allocate buffer to write hydrodynamic accelerations for all particles. + # `PostprocessCallback` transfers data to the CPU before calling custom quantities, + # so this can be a regular `Array` even when the simulation is running on a GPU. + dv = Array{ELTYPE, 2}(undef, (v_nvariables(system), nparticles(system))) + + return ThrustCalculator(zero(ELTYPE), system_index, dv, eachparticle, direction_) +end + +function reset!(calculator::ThrustCalculator) + calculator.thrust = zero(calculator.thrust) + + return calculator +end + +""" + calculated_thrust(calculator::ThrustCalculator) + +Get the latest projected hydrodynamic force from a [`ThrustCalculator`](@ref). +""" +function calculated_thrust(calculator::ThrustCalculator) + return calculator.thrust +end + +function (calculator::ThrustCalculator)(system, dv_ode, du_ode, v_ode, u_ode, semi, t) + if system_indices(system, semi) != calculator.system_index + return nothing + end + + @trixi_timeit timer() "calculate thrust" begin + calculator.thrust = compute_thrust(system, calculator.eachparticle, + calculator.direction, + calculator.dv, v_ode, u_ode, semi, t) + end + + return calculator.thrust +end + +function compute_thrust(system, eachparticle, direction, dv, v_ode, u_ode, semi, t) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + compute_structure_acceleration!(dv, v, u, system, eachparticle, true, + v_ode, u_ode, semi, t) + + return projected_force(dv, system, eachparticle, direction) +end + +function projected_force(dv, system, eachparticle, direction) + force = zero(eltype(system)) + + # Note that this is a reduction, so we cannot use `@threaded` here. + for particle in eachparticle + dv_particle = extract_svector(dv, system, particle) + force_particle = system.mass[particle] * dv_particle + force += dot(force_particle, direction) + end + + return force +end diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index bc81241e13..4d5fe89555 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -261,6 +261,11 @@ end end # === Initialization === + +function initialize_neighborhood_searches!(semi, u0_ode, restart_with::Nothing) + initialize_neighborhood_searches!(semi) +end + function initialize_neighborhood_searches!(semi) foreach_system(semi) do system foreach_system(semi) do neighbor @@ -275,11 +280,16 @@ function initialize_neighborhood_search!(semi, system, neighbor) # TODO Initialize after adapting to the GPU. # Currently, this cannot use `semi.parallelization_backend` # because data is still on the CPU. + parallelization_backend = if semi.parallelization_backend isa KernelAbstractions.GPU + PolyesterBackend() + else + semi.parallelization_backend + end PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi), initial_coordinates(system), initial_coordinates(neighbor), eachindex_y=each_active_particle(neighbor), - parallelization_backend=PolyesterBackend()) + parallelization_backend=parallelization_backend) return semi end diff --git a/src/general/restart.jl b/src/general/restart.jl new file mode 100644 index 0000000000..33efd006ec --- /dev/null +++ b/src/general/restart.jl @@ -0,0 +1,133 @@ +function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Tuple{Vararg{String}}) + # Check number of systems + if length(semi.systems) != length(restart_with) + throw(ArgumentError("number of systems in `semi` does not match number of restart files provided " * + "in `restart_with`")) + end + + # Check that systems match + expected_system_names = system_names(semi.systems) + foreach_system(semi) do system + system_index = system_indices(system, semi) + filename = restart_with[system_index] + expected_system_name = expected_system_names[system_index] + if !occursin(expected_system_name, basename(splitext(filename)[1])) + throw(ArgumentError("Filename '$filename' for system $system_index does not contain expected name '$expected_system_name'. " * + "Expected a VTK file for system of type $(nameof(typeof(system))).")) + end + end + + # Set initial conditions + foreach_noalloc(semi.systems, restart_with) do (system, restart_file) + v0_system = wrap_v(v0_ode, system, semi) + u0_system = wrap_u(u0_ode, system, semi) + + restart_data = vtk2trixi(restart_file; create_initial_condition=false) + v_restart = restart_v(system, restart_data) + u_restart = restart_u(system, restart_data) + + v0_system .= Adapt.adapt(semi.parallelization_backend, v_restart) + u0_system .= Adapt.adapt(semi.parallelization_backend, u_restart) + + restore_previous_state!(system, restart_file) + end +end + +function time_span(tspan, restart_with::Tuple{Vararg{String}}) + # Read restart times from all files + restart_times = [vtk2trixi(file).time for file in restart_with] + t_restart = convert(eltype(tspan), first(restart_times)) + + # Check if all restart files have the same time + if !all(isapprox(t, t_restart) for t in restart_times) + throw(ArgumentError("All restart files must start from the same time.")) + end + + if !isapprox(tspan[1], t_restart) + @info "Adjusting initial time from $(tspan[1]) to restart time $t_restart" + end + + return (t_restart, tspan[2]) +end + +function write_density_and_pressure!(v_restart, system, density_calculator, + pressure, density) + return v_restart +end + +function write_density_and_pressure!(v_restart, system, + density_calculator::ContinuityDensity, + pressure, density) + v_restart[size(v_restart, 1), :] = density + + return v_restart +end + +function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem, + density_calculator::ContinuityDensity, + pressure, density) + v_restart[size(v_restart, 1), :] = density + v_restart[size(v_restart, 1) - 1, :] = pressure + + return v_restart +end + +function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem, + density_calculator::SummationDensity, + pressure, density) + v_restart[size(v_restart, 1) - 1, :] = pressure + + return v_restart +end + +restore_previous_state!(system, restart_file) = system + +function initialize_neighborhood_searches!(semi, u0_ode, + restart_with::Tuple{Vararg{String}}) + foreach_system(semi) do system + foreach_system(semi) do neighbor + initialize_neighborhood_search!(semi, system, neighbor, u0_ode) + end + end + + return semi +end + +function initialize_neighborhood_search!(semi, system, neighbor, u0_ode) + # TODO Initialize after adapting to the GPU. + # Currently, this cannot use `semi.parallelization_backend` + # because data is still on the CPU. + PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi), + initial_restart_coordinates(system, u0_ode, semi), + initial_restart_coordinates(neighbor, u0_ode, semi), + eachindex_y=each_active_particle(neighbor), + parallelization_backend=PolyesterBackend()) + return semi +end + +function initialize_neighborhood_search!(semi, system::TotalLagrangianSPHSystem, + neighbor::TotalLagrangianSPHSystem, u0_ode) + # For TLSPH, the self-interaction NHS is already initialized in the system constructor + return semi +end + +function initial_restart_coordinates(system, u0_ode, semi) + # Transfer to CPU if data is on the GPU. Do nothing if already on CPU. + return transfer2cpu(wrap_u(u0_ode, system, semi)) +end + +function initial_restart_coordinates(system::Union{WallBoundarySystem, + AbstractStructureSystem}, u0_ode, semi) + return initial_coordinates(system) +end + +function initialize!(semi::Semidiscretization, restart_with::Tuple{Vararg{String}}) + foreach_system(semi) do system + # Initialize this system + initialize_restart!(system, semi) + end + + return semi +end + +initialize_restart!(system, semi) = initialize!(system, semi) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 463b3f0a9e..514f1e6b9f 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -184,7 +184,7 @@ end end """ - semidiscretize(semi, tspan; reset_threads=true) + semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing) Create an `ODEProblem` from the semidiscretization with the specified `tspan`. @@ -193,6 +193,14 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`. - `tspan`: The time span over which the simulation will be run. # Keywords +- `restart_with`: Can be used to restart the simulation from VTK solution files (see [`SolutionSavingCallback`](@ref)). + This can be either `nothing` (default, no restart), a single filename as a `String`, + or a `Tuple` of filenames, one for each system in the [`Semidiscretization`](@ref). + The order of the filenames must match the order of the systems in the [`Semidiscretization`](@ref). + Note that `semidiscretize` replaces the initial time (`tspan[1]`) with the timestamp read + from the VTK files. If the user-provided `tspan[1]` does not match the restart time, + it is adjusted and an info message is logged. If multiple files are provided, their + timestamps must match. - `reset_threads`: A boolean flag to reset Polyester.jl threads before the simulation (default: `true`). After an error within a threaded loop, threading might be disabled. Resetting the threads before the simulation ensures that threading is enabled again for the simulation. @@ -219,9 +227,16 @@ timespan: (0.0, 1.0) u0: ([...], [...]) *this line is ignored by filter* ``` """ -function semidiscretize(semi, tspan; reset_threads=true) +function semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing) (; systems) = semi + if restart_with isa String + restart_with = (restart_with,) + elseif !isnothing(restart_with) && !(restart_with isa NTuple{<:Any, String}) + throw(ArgumentError("`restart_with` must be `nothing`, a string, or a tuple of strings, " * + "got $(typeof(restart_with))")) + end + # Check that all systems have the same eltype first_system = first(systems) if !all(system -> eltype(system) === eltype(first_system), systems) @@ -267,14 +282,11 @@ function semidiscretize(semi, tspan; reset_threads=true) end # Set initial condition - foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system - write_u0!(u0_system, system) - write_v0!(v0_system, system) - end + set_initial_conditions!(v0_ode, u0_ode, semi, restart_with) # TODO initialize after adapting to the GPU. # Requires https://github.com/trixi-framework/PointNeighbors.jl/pull/86. - initialize_neighborhood_searches!(semi) + initialize_neighborhood_searches!(semi, u0_ode, restart_with) if semi.parallelization_backend isa KernelAbstractions.GPU # Convert all arrays in the systems to the correct array type. @@ -304,10 +316,7 @@ function semidiscretize(semi, tspan; reset_threads=true) end # Initialize all particle systems - foreach_system(semi_new) do system - # Initialize this system - initialize!(system, semi_new) - end + initialize!(semi_new, restart_with) # Reset callback flag that will be set by the `UpdateCallback` semi_new.update_callback_used[] = false @@ -318,7 +327,27 @@ function semidiscretize(semi, tspan; reset_threads=true) # here, since we cannot change `p` from within the callback (only its contents). p = @NamedTuple{semi::typeof(semi_new), split_integration_data::Any}((semi_new, nothing)) - return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, p) + + return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, + time_span(tspan, restart_with), p) +end + +function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Nothing) + foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system + write_u0!(u0_system, system) + write_v0!(v0_system, system) + end +end + +time_span(tspan, restart_with::Nothing) = tspan + +function initialize!(semi::Semidiscretization, restart_with::Nothing) + foreach_system(semi) do system + # Initialize this system + initialize!(system, semi) + end + + return semi end """ diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl index 86bae1af77..aefa4a3c72 100644 --- a/src/io/read_vtk.jl +++ b/src/io/read_vtk.jl @@ -51,8 +51,19 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing, point_coords = ReadVTK.get_points(vtk_file) cELTYPE = isnothing(coordinates_eltype) ? eltype(point_coords) : coordinates_eltype - ELTYPE = isnothing(element_type) ? - eltype(first(ReadVTK.get_data(point_data["pressure"]))) : element_type + + if !isnothing(element_type) + ELTYPE = element_type + else + # Try to get element type from pressure or density (whichever exists) + ELTYPE = if "pressure" in keys(point_data) + eltype(first(ReadVTK.get_data(point_data["pressure"]))) + elseif "material_density" in keys(point_data) + eltype(first(ReadVTK.get_data(point_data["material_density"]))) + else + error("Neither 'pressure' nor 'material_density' field found in VTK file") + end + end results = Dict{Symbol, Any}() @@ -87,7 +98,8 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing, results[:particle_spacing] results[:coordinates] = coordinates results[:time] = "time" in keys(field_data) ? - first(ReadVTK.get_data(field_data["time"])) : zero(ELTYPE) + convert.(ELTYPE, first(ReadVTK.get_data(field_data["time"]))) : + zero(ELTYPE) append!(used_keys, ["index", "ndims"]) # Load any custom quantities diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index 0b19f32337..fa3fb80985 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -447,6 +447,16 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem) for particle in eachparticle(system)] vtk["pressure"] = [current_pressure(v, system, particle) for particle in eachparticle(system)] + vtk["zone_id"] = [system.boundary_zone_indices[particle] + for particle in eachparticle(system)] + + if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values) + for (i, pressure_model) in enumerate(system.cache.pressure_reference_values) + if pressure_model isa AbstractPressureModel + vtk["boundary_zone_pressure_$i"] = system.cache.pressure_reference_values[i].pressure[] + end + end + end if system.calculate_flow_rate Q_total = zero(eltype(system)) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 7a56b2cfcc..5e72ae310e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -135,6 +135,10 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; calculate_flow_rate, cache) end +function Base.:(==)(system1::OpenBoundarySystem, system2::OpenBoundarySystem) + return system1.mass === system2.mass +end + function initialize!(system::OpenBoundarySystem, semi) (; boundary_zones) = system @@ -143,6 +147,9 @@ function initialize!(system::OpenBoundarySystem, semi) return system end +# Skip during restart, as boundary zone indices are updated in `restore_previous_state!` +initialize_restart!(system::OpenBoundarySystem, semi) = system + function create_cache_open_boundary(boundary_model, fluid_system, initial_condition, density_diffusion, calculate_flow_rate, boundary_zones) reference_values = map(bz -> bz.reference_values, boundary_zones) @@ -261,6 +268,8 @@ end return system.shifting_technique end +@inline density_calculator(system::OpenBoundarySystem) = nothing + system_sound_speed(system::OpenBoundarySystem) = system_sound_speed(system.fluid_system) @inline hydrodynamic_mass(system::OpenBoundarySystem, particle) = system.mass[particle] @@ -705,6 +714,72 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone, return system end +function restart_u(system::OpenBoundarySystem, data) + coords_total = zeros(coordinates_eltype(system), u_nvariables(system), + n_integrated_particles(system)) + coords_total .= coordinates_eltype(system)(1e16) + + # Since only active particles are written during saving, the loaded `data` contains + # only active particles. These are placed at the beginning of the array, leaving the + # inactive buffer particles at the end. Thus, we can safely activate the first N particles. + coords_active = data.coordinates + for particle in axes(coords_active, 2) + for dim in 1:ndims(system) + coords_total[dim, particle] = coords_active[dim, particle] + end + end + + system.buffer.active_particle .= false + system.buffer.active_particle[1:size(coords_active, 2)] .= true + + update_system_buffer!(system.buffer) + + return coords_total +end + +function restart_v(system::OpenBoundarySystem, data) + v_total = zeros(eltype(system), v_nvariables(system), + n_integrated_particles(system)) + + # Since only active particles are written during saving, the loaded `data` contains + # only active particles. These are placed at the beginning of the array, leaving the + # inactive buffer particles at the end. Thus, we can safely activate the first N particles. + v_active = zeros(eltype(system), v_nvariables(system), size(data.velocity, 2)) + + v_active[1:ndims(system), :] = data.velocity + write_density_and_pressure!(v_active, system.fluid_system, + density_calculator(system), data.pressure, data.density) + + for particle in axes(v_active, 2) + for i in axes(v_active, 1) + v_total[i, particle] = v_active[i, particle] + end + end + + return v_total +end + +function restore_previous_state!(system::OpenBoundarySystem, file) + # We cannot simply use `update_boundary_zone_indices!` because rounding errors during file I/O + # may result in particles being located outside their intended boundary zone, even though they + # were written as active particles. + set_zero!(system.boundary_zone_indices) + + values = vtk2trixi(file; create_initial_condition=false) + system.boundary_zone_indices[each_integrated_particle(system)] .= values.zone_id + + if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values) + for (i, pressure_model) in enumerate(system.cache.pressure_reference_values) + if pressure_model isa AbstractPressureModel + pressure_model.pressure[] = values[Symbol(:boundary_zone_pressure_, i)] + pressure_model.flow_rate[] = values[Symbol(:Q_, i)] + end + end + end + + return system +end + # Open-boundary interpolation should reconstruct the surrounding fluid-like velocity field. # Therefore, only actual fluid systems and other open-boundary particles contribute; # rigid bodies, walls, and other non-fluid systems are intentionally excluded. 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 3f0addc408..c7dfdfe03c 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl @@ -328,6 +328,10 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles) print(io, ")") end +@inline function density_calculator(model::BoundaryModelDummyParticles) + return model.density_calculator +end + # Set the initial pressure to zero for visualization function initial_boundary_pressure(initial_density, density_calculator, _) return zero(initial_density) diff --git a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl index 7300cc03d9..2ea6e82cf6 100644 --- a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl +++ b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl @@ -43,6 +43,10 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar) print(io, ")") end +@inline function density_calculator(model::BoundaryModelMonaghanKajtar) + return nothing +end + @inline function pressure_acceleration(particle_system, neighbor_system::Union{WallBoundarySystem{<:BoundaryModelMonaghanKajtar}, TotalLagrangianSPHSystem{<:BoundaryModelMonaghanKajtar}}, diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl index 3a546f4013..115b6202d3 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -264,6 +264,23 @@ function restart_with!(system::WallBoundarySystem{<:BoundaryModelDummyParticles{ return system end +function restart_u(system::WallBoundarySystem, data) + if n_integrated_particles(system) > 0 + throw(ArgumentError("`WallBoundarySystem` does not support integrated particle coordinates")) + end + + return zeros(eltype(system), u_nvariables(system), n_integrated_particles(system)) +end + +function restart_v(system::WallBoundarySystem, data) + v_ode = zeros(eltype(system), v_nvariables(system), n_integrated_particles(system)) + + write_density_and_pressure!(v_ode, system, density_calculator(system), + data.pressure, data.density) + + return v_ode +end + # To incorporate the effect at boundaries in the viscosity term of the RHS, the neighbor # viscosity model has to be used. @inline function viscosity_model(system::WallBoundarySystem, @@ -315,6 +332,10 @@ function system_correction(system::WallBoundarySystem{<:BoundaryModelDummyPartic return system.boundary_model.correction end +@inline function density_calculator(system::WallBoundarySystem) + return density_calculator(system.boundary_model) +end + function system_data(system::WallBoundarySystem, dv_ode, du_ode, v_ode, u_ode, semi) dv = [current_acceleration(system, particle) for particle in eachparticle(system)] v = wrap_v(v_ode, system, semi) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index eb6338f1cc..97e9bf08b3 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -176,7 +176,7 @@ end # 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 + if particle_system == neighbor_system density_diffusion!(drho_particle, density_diffusion(particle_system), particle_system, particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, grad_kernel) @@ -253,6 +253,50 @@ end return nothing end +function restart_u(system::AbstractFluidSystem, data) + coords_total = zeros(coordinates_eltype(system), u_nvariables(system), + n_integrated_particles(system)) + coords_total .= coordinates_eltype(system)(1e16) + + coords_active = data.coordinates + + for particle in axes(coords_active, 2) + for dim in 1:ndims(system) + coords_total[dim, particle] = coords_active[dim, particle] + end + end + + buf = buffer(system) + if !isnothing(buf) + buf.active_particle .= false + buf.active_particle[1:size(coords_active, 2)] .= true + update_system_buffer!(buf) + end + + return coords_total +end + +function restart_v(system::AbstractFluidSystem, data) + velocity_total = zeros(eltype(system), v_nvariables(system), + n_integrated_particles(system)) + velocity_active = zeros(eltype(system), v_nvariables(system), size(data.velocity, 2)) + + velocity_active[1:ndims(system), :] = data.velocity + + density_calc = hasproperty(system, :density_calculator) ? density_calculator(system) : + nothing + write_density_and_pressure!(velocity_active, system, density_calc, data.pressure, + data.density) + + for particle in axes(velocity_active, 2) + for i in axes(velocity_active, 1) + velocity_total[i, particle] = velocity_active[i, particle] + end + end + + return velocity_total +end + function check_configuration(fluid_system::AbstractFluidSystem, systems, nhs) if !(fluid_system isa ParticlePackingSystem) && !isnothing(fluid_system.surface_tension) foreach_system(systems) do neighbor diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 2e5474a33c..00eb47409a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -173,6 +173,11 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, cache) end +@inline function Base.:(==)(system1::WeaklyCompressibleSPHSystem, + system2::WeaklyCompressibleSPHSystem) + return system1.mass === system2.mass +end + function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl index d3ee370237..232f011838 100644 --- a/src/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl @@ -108,6 +108,8 @@ function interact!(dv, v_particle_system, u_particle_system, neighbor_system::AbstractFluidSystem, semi; 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 diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index fbd4dc0333..9b4f2b2f1f 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -454,6 +454,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 @@ -556,6 +558,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 @@ -634,6 +684,16 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u) restart_with!(system, system.boundary_model, v, u) end +function restart_u(system::TotalLagrangianSPHSystem, data) + # The integration array requires only the particles that need to be integrated + return data.coordinates[:, each_integrated_particle(system)] +end + +function restart_v(system::TotalLagrangianSPHSystem, data) + # The integration array requires only the particles that need to be integrated + return data.velocity[:, each_integrated_particle(system)] +end + # An explanation of these equation can be found in # J. Lubliner, 2008. Plasticity theory. # See here below Equation 5.3.21 for the equation for the equivalent stress. diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl index 6a78b412e2..fd7e821ad1 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) @@ -145,7 +141,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/test/callbacks/mechanical_work_calculator.jl b/test/callbacks/mechanical_work_calculator.jl index 31496fb45d..badc3dfa46 100644 --- a/test/callbacks/mechanical_work_calculator.jl +++ b/test/callbacks/mechanical_work_calculator.jl @@ -1,9 +1,13 @@ -@testset verbose=true "MechanicalWorkCalculatorCallback" begin +@testset verbose=true "MechanicalWorkCalculator" begin # Mock system struct MockSystem <: TrixiParticles.AbstractStructureSystem{2} eltype::Type + mass::Vector + + function MockSystem(ELTYPE) + new(ELTYPE, ELTYPE[1, 2, 3, 4]) + end end - TrixiParticles.nparticles(::MockSystem) = 4 Base.eltype(system::MockSystem) = system.eltype function TrixiParticles.create_neighborhood_search(neighborhood_search, @@ -18,43 +22,120 @@ @testset "Constructor and Basic Properties" begin # Test default constructor - callback = MechanicalWorkCalculatorCallback(system64, semi64) - @test callback.affect!.system_index == 1 - @test callback.affect!.interval == 1 - @test callback.affect!.t[] == 0.0 - @test callback.affect!.work[] == 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_mechanical_work(callback) == 0.0 - - # Test constructor with interval - callback = MechanicalWorkCalculatorCallback(system64, semi64; interval=5) - @test callback.affect!.interval == 5 - @test eltype(callback.affect!.work) == Float64 - @test eltype(callback.affect!.t) == Float64 + calculator = MechanicalWorkCalculator(system64, semi64) + @test calculator.system_index == 1 + @test calculator.t == 0.0 + @test calculator.work == 0.0 + @test !calculator.initialized + @test calculator.dv isa Array{Float64, 2} + @test size(calculator.dv) == (2, 4) + @test calculator.eachparticle == 5:4 + @test calculated_mechanical_work(calculator) == 0.0 + + # Test constructor with explicit particle range + calculator = MechanicalWorkCalculator(system64, semi64; eachparticle=1:2) + @test calculator.eachparticle == 1:2 + @test eltype(calculator.work) == Float64 + @test eltype(calculator.t) == Float64 # Test with specific element type - callback = MechanicalWorkCalculatorCallback(system32, semi32; interval=2) - @test eltype(callback.affect!.work) == Float32 - @test eltype(callback.affect!.t) == Float32 + calculator = MechanicalWorkCalculator(system32, semi32) + @test eltype(calculator.work) == Float32 + @test eltype(calculator.t) == Float32 + end + + @testset "ThrustCalculator constructor and force projection" begin + @test_throws UndefKeywordError ThrustCalculator(system64, semi64) + + calculator = ThrustCalculator(system64, semi64; direction=SVector(1.0, 0.0)) + @test calculator.system_index == 1 + @test calculator.thrust == 0.0 + @test calculator.dv isa Array{Float64, 2} + @test size(calculator.dv) == (2, 4) + @test calculator.eachparticle == eachparticle(system64) + @test calculator.direction == SVector(1.0, 0.0) + @test calculated_thrust(calculator) == 0.0 + + calculator = ThrustCalculator(system32, semi32; direction=(0.0, 2.0), + eachparticle=2:3) + @test eltype(calculator.thrust) == Float32 + @test calculator.direction == SVector(0.0f0, 1.0f0) + @test calculator.eachparticle == 2:3 + + @test_throws ArgumentError ThrustCalculator(system64, semi64; direction=(0.0, 0.0)) + + dv = [2.0 -1.0 0.0 3.0 + 4.0 5.0 -2.0 1.0] + @test TrixiParticles.projected_force(dv, system64, eachparticle(system64), + SVector(1.0, 0.0)) == 12.0 + @test TrixiParticles.projected_force(dv, system64, 2:3, + SVector(0.0, 1.0)) == 4.0 + + TrixiParticles.reset!(calculator) + @test calculated_thrust(calculator) == 0.0f0 end - @testset "show" begin - callback = MechanicalWorkCalculatorCallback(system64, semi64; interval=10) - - # Test compact representation - show_compact = "MechanicalWorkCalculatorCallback{Float64}(interval=10)" - @test repr(callback) == show_compact - - # Test detailed representation - check against expected box format - show_box = """ - ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ - │ MechanicalWorkCalculatorCallback{Float64} │ - │ ═════════════════════════════════════════ │ - │ interval: ……………………………………………………… 10 │ - └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" - @test repr("text/plain", callback) == show_box + @testset "ThrustCalculator FSI force" begin + particle_spacing = 1.0 + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 1.0 + fluid_density = 1000.0 + structure_density = 2000.0 + particle_volume = particle_spacing^2 + + state_equation = StateEquationCole(sound_speed=10.0, + reference_density=fluid_density, + exponent=1.0) + + fluid_ic = InitialCondition(; coordinates=reshape([0.0, 0.0], 2, 1), + velocity=zeros(2, 1), + mass=[particle_volume * fluid_density], + density=[fluid_density], particle_spacing) + + fluid_system = WeaklyCompressibleSPHSystem(fluid_ic; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation, + reference_particle_spacing=particle_spacing) + + structure_coordinates = reshape([1.5, 0.0], 2, 1) + structure_ic = InitialCondition(; coordinates=structure_coordinates, + velocity=zeros(2, 1), + mass=[particle_volume * structure_density], + density=[structure_density], particle_spacing) + + boundary_model = BoundaryModelDummyParticles([fluid_density], + [particle_volume * fluid_density], + AdamiPressureExtrapolation(), + smoothing_kernel, smoothing_length; + state_equation, + reference_particle_spacing=particle_spacing) + + structure_system = TotalLagrangianSPHSystem(structure_ic; smoothing_kernel, + smoothing_length, + young_modulus=1.0e5, + poisson_ratio=0.3, + boundary_model) + + semi_ = Semidiscretization(fluid_system, structure_system) + ode = semidiscretize(semi_, (0.0, 0.01)) + semi = ode.p.semi + + v_ode, u_ode = ode.u0.x + dv_ode = zero(v_ode) + TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0) + + fluid = semi.systems[1] + structure = semi.systems[2] + dv_fluid = TrixiParticles.wrap_v(dv_ode, fluid, semi) + + thrust = ThrustCalculator(structure, semi; direction=SVector(1.0, 0.0)) + thrust(structure, dv_ode, nothing, v_ode, u_ode, semi, 0.0) + + expected_force = -fluid.mass[1] * dv_fluid[1, 1] + @test !iszero(expected_force) + @test isapprox(calculated_thrust(thrust), expected_force; + rtol=sqrt(eps()), atol=sqrt(eps())) end @testset "update_mechanical_work_calculator!" begin @@ -110,56 +191,53 @@ TrixiParticles.update_quantities!(system, v, u, v_ode, u_ode, semi, 0.0) # Set up test parameters - work1 = Ref(0.0) + work1 = 0.0 dt1 = 0.1 # Test that mechanical work is integrated, i.e., values of instantaneous # power are accumulated over time. This initial work should just be an offset. # Also, half the step size means half the work increase. - work2 = Ref(1.0) + work2 = 1.0 dt2 = 0.05 # Test `only_compute_force_on_fluid` - work3 = Ref(0.0) + work3 = 0.0 dt3 = 0.1 eachparticle = (TrixiParticles.n_integrated_particles(system) + 1):nparticles(system) dv = zeros(2, nparticles(system)) - TrixiParticles.update_mechanical_work_calculator!(work1, system, eachparticle, - false, - dv, v_ode, u_ode, semi, 0.0, - dt1) - TrixiParticles.update_mechanical_work_calculator!(work2, system, eachparticle, - false, - dv, v_ode, u_ode, semi, 0.0, - dt2) - TrixiParticles.update_mechanical_work_calculator!(work3, system, eachparticle, - true, - dv, v_ode, u_ode, semi, 0.0, - dt3) + work1 = TrixiParticles.update_mechanical_work!(work1, system, eachparticle, + false, dv, v_ode, u_ode, semi, + 0.0, dt1) + work2 = TrixiParticles.update_mechanical_work!(work2, system, eachparticle, + false, dv, v_ode, u_ode, semi, + 0.0, dt2) + work3 = TrixiParticles.update_mechanical_work!(work3, system, eachparticle, + true, dv, v_ode, u_ode, semi, + 0.0, dt3) if i == 1 - @test isapprox(work1[], 0.8) - @test isapprox(work2[], 1.0 + 0.4) + @test isapprox(work1, 0.8) + @test isapprox(work2, 1.0 + 0.4) elseif i == 2 # For very soft material, we can just pull up the top row of particles # and the work required is almost just the potential energy difference. - @test isapprox(work1[], 0.4080357142857143) - @test isapprox(work2[], 1.0 + 0.5 * 0.4080357142857143) + @test isapprox(work1, 0.4080357142857143) + @test isapprox(work2, 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 work # to pull the material apart. - @test isapprox(work1[], 803.9714285714281) - @test isapprox(work2[], 1.0 + 0.5 * 803.9714285714281) + @test isapprox(work1, 803.9714285714281) + @test isapprox(work2, 1.0 + 0.5 * 803.9714285714281) elseif i == 4 # For a very stiff material, the work is even larger. - @test isapprox(work1[], 80357.5428571428) - @test isapprox(work2[], 1.0 + 0.5 * 80357.5428571428) + @test isapprox(work1, 80357.5428571428) + @test isapprox(work2, 1.0 + 0.5 * 80357.5428571428) end - @test isapprox(work3[], 0.0) + @test isapprox(work3, 0.0) end end end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 7498147374..a98af097ab 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -91,7 +91,7 @@ end end - @trixi_testset "structure/oscillating_beam_2d.jl with MechanicalWorkCalculatorCallback" begin + @trixi_testset "structure/oscillating_beam_2d.jl with MechanicalWorkCalculator" begin # Load variables from the example trixi_include(@__MODULE__, joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), @@ -129,17 +129,18 @@ ode = semidiscretize(semi, (0.0, 1.0)) system = ode.p.semi.systems[1] - mechanical_work_calculator = MechanicalWorkCalculatorCallback(system, semi; - interval=1) + mechanical_work = MechanicalWorkCalculator(system, semi) + postprocess_callback = PostprocessCallback(; interval=1, mechanical_work, + write_file_interval=0) sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL49(), save_everystep=false, - callback=mechanical_work_calculator) + callback=postprocess_callback) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol) == 0 # Potential energy difference should be m * g * h - @test isapprox(calculated_mechanical_work(mechanical_work_calculator), + @test isapprox(calculated_mechanical_work(mechanical_work), sum(system.mass) * gravity * 1, rtol=rtol[name]) end @@ -156,9 +157,9 @@ end @testset verbose=true "FSI" begin - @trixi_testset "fluid/hydrostatic_water_column_2d.jl with MechanicalWorkCalculatorCallback and moving TLSPH walls" begin + @trixi_testset "fluid/hydrostatic_water_column_2d.jl with MechanicalWorkCalculator and moving TLSPH walls" begin # In this test, we move a water-filled tank up against gravity by 1 unit - # and verify that the work accumulated by the `MechanicalWorkCalculatorCallback` + # and verify that the work accumulated by the `MechanicalWorkCalculator` # matches the expected potential energy difference. # Load variables from the example @@ -181,7 +182,7 @@ # Create TLSPH system for the tank walls and clamp all particles. # This is identical to a `WallBoundarySystem`, but now we can - # use the `MechanicalWorkCalculatorCallback` to compute the mechanical work. + # use the `MechanicalWorkCalculator` to compute the mechanical work. boundary_spacing = tank.boundary.particle_spacing tlsph_kernel = WendlandC2Kernel{2}() tlsph_smoothing_length = sqrt(2) * boundary_spacing @@ -203,18 +204,17 @@ tlsph_system_new = ode.p.semi.systems[2] # Mechanical work calculators for fluid + tank and fluid only - mechanical_work_calculator1 = MechanicalWorkCalculatorCallback(tlsph_system_new, - semi; - interval=1) - mechanical_work_calculator2 = MechanicalWorkCalculatorCallback(tlsph_system_new, - semi; - interval=1, - only_compute_force_on_fluid=true) + mechanical_work1 = MechanicalWorkCalculator(tlsph_system_new, semi) + mechanical_work2 = MechanicalWorkCalculator(tlsph_system_new, semi; + only_compute_force_on_fluid=true) + thrust = ThrustCalculator(tlsph_system_new, semi; direction=SVector(0.0, 1.0)) + postprocess_callback = PostprocessCallback(; interval=1, mechanical_work1, + mechanical_work2, thrust, + write_file_interval=0) sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=CallbackSet(info_callback, - mechanical_work_calculator1, - mechanical_work_calculator2)) + postprocess_callback)) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol) == 0 @@ -229,11 +229,12 @@ # 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_mechanical_work(mechanical_work_calculator1), + @test isapprox(calculated_mechanical_work(mechanical_work1), expected_energy_fluid + expected_energy_tank, rtol=5e-4) - @test isapprox(calculated_mechanical_work(mechanical_work_calculator2), + @test isapprox(calculated_mechanical_work(mechanical_work2), expected_energy_fluid, rtol=5e-4) + @test isfinite(calculated_thrust(thrust)) end @trixi_testset "fsi/falling_water_column_2d.jl" begin diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index e3415c0e9d..73e66813fc 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -896,4 +896,55 @@ end end end end + + @testset verbose=true "Restart" begin + trixi_include_changeprecision(Float32, @__MODULE__, + joinpath(examples_dir(), "fluid", + "poiseuille_flow_2d.jl"), + tspan=(0.0f0, 0.6f0), sound_speed_factor=10, + particle_spacing=4.0f-5, + coordinates_eltype=Float32, + parallelization_backend=Main.parallelization_backend) + + # Since this is an open boundary simulation, the number of active particles may + # differ. The results must be interpolated to enable comparison with the restart + # simulation. The fluid domain starts at `x = 10 * particle_spacing`. + n_interpolation_points = 10 + start_point = [0.0f0 + 10 * particle_spacing, domain_size[1] / 2] + end_point = [domain_size[1] - 10 * particle_spacing, domain_size[2] / 2] + result_full = interpolate_line(start_point, end_point, n_interpolation_points, + sol.prob.p.semi, sol.prob.p.semi.systems[1], sol, + cut_off_bnd=false) + + # Run half simulation and safe checkpoint + trixi_include_changeprecision(Float32, @__MODULE__, + joinpath(examples_dir(), "fluid", + "poiseuille_flow_2d.jl"), + tspan=(0.0f0, 0.3f0), sound_speed_factor=10, + particle_spacing=4.0f-5, + coordinates_eltype=Float32, + parallelization_backend=Main.parallelization_backend) + + iter = saving_callback.affect!.affect!.latest_saved_iter + fluid_restart = joinpath("out", "fluid_1_$iter.vtu") + open_boundary_restart = joinpath("out", "open_boundary_1_$iter.vtu") + boundary_restart = joinpath("out", "boundary_1_$iter.vtu") + + ode_restart = semidiscretize(semi, (0.3f0, 0.6f0); + restart_with=(fluid_restart, open_boundary_restart, + boundary_restart)) + + sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1.0f-5, + reltol=1.0f-3, dtmax=1.0f-2, save_everystep=false, + callback=UpdateCallback()) + + result_restart = interpolate_line(start_point, end_point, + n_interpolation_points, sol_restart.prob.p.semi, + sol_restart.prob.p.semi.systems[1], + sol_restart, cut_off_bnd=false) + + @test isapprox(result_full.velocity, result_restart.velocity, rtol=8.0f-3) + @test isapprox(result_full.density, result_restart.density, rtol=8.0f-4) + @test isapprox(result_full.pressure, result_restart.pressure, rtol=8.0f-2) + end end diff --git a/test/general/general.jl b/test/general/general.jl index 2eb24d50a0..acb07de8b9 100644 --- a/test/general/general.jl +++ b/test/general/general.jl @@ -6,4 +6,5 @@ include("interpolation.jl") include("buffer.jl") include("util.jl") include("custom_quantities.jl") +include("restart.jl") include("neighborhood_search.jl") diff --git a/test/general/restart.jl b/test/general/restart.jl new file mode 100644 index 0000000000..21e627e8e6 --- /dev/null +++ b/test/general/restart.jl @@ -0,0 +1,109 @@ +@trixi_testset "Restart" begin + @trixi_testset "Poiseuille Flow Half-Simulation Restart" begin + # Run full simulation + trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"), + tspan=(0.0, 0.6), sound_speed_factor=10, particle_spacing=4e-5, + info_callback=nothing) + + # Since this is an open boundary simulation, the number of active particles may + # differ. The results must be interpolated to enable comparison with the restart + # simulation. The fluid domain starts at `x = 10 * particle_spacing`. + n_interpolation_points = 10 + start_point = [0.0 + 10 * particle_spacing, channel_height / 2] + end_point = [channel_length - 10 * particle_spacing, channel_height / 2] + result_full = interpolate_line(start_point, end_point, n_interpolation_points, + semi, fluid_system, sol, cut_off_bnd=false) + + # Run half simulation + trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"), + tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5, + info_callback=nothing) + + iter = saving_callback.affect!.affect!.latest_saved_iter + fluid_restart = joinpath("out", "fluid_1_$iter.vtu") + open_boundary_restart = joinpath("out", "open_boundary_1_$iter.vtu") + boundary_restart = joinpath("out", "boundary_1_$iter.vtu") + + ode_restart = semidiscretize(semi, (0.3, 0.6); + restart_with=(fluid_restart, open_boundary_restart, + boundary_restart)) + + sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1e-5, reltol=1e-3, + dtmax=1e-2, save_everystep=false, callback=UpdateCallback()) + + result_restart = interpolate_line(start_point, end_point, + n_interpolation_points, sol_restart.prob.p.semi, + sol_restart.prob.p.semi.systems[1], + sol_restart, cut_off_bnd=false) + + @test isapprox(result_full.velocity, result_restart.velocity, rtol=1e-2) + @test isapprox(result_full.density, result_restart.density, rtol=8e-4) + @test isapprox(result_full.pressure, result_restart.pressure, rtol=8e-2) + end + + @trixi_testset "Poiseuille Flow Restore Previous State" begin + R1 = 1.7714 + R2 = 106.66 + C = 1.1808e-2 + pressure_model = RCRWindkesselModel(; peripheral_resistance=R2, + compliance=C, + characteristic_resistance=R1) + # Run half simulation + trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"), + outlet_reference_pressure=pressure_model, info_callback=nothing, + tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5) + + iter = saving_callback.affect!.affect!.latest_saved_iter + fluid_restart = joinpath("out", "fluid_1_$iter.vtu") + open_boundary_restart = joinpath("out", "open_boundary_1_$iter.vtu") + boundary_restart = joinpath("out", "boundary_1_$iter.vtu") + + ode_restart = semidiscretize(semi, (0.3, 0.6); + restart_with=(fluid_restart, open_boundary_restart, + boundary_restart)) + restart_pressure = ode_restart.p.semi.systems[2].cache.pressure_reference_values[2].pressure[] + restart_flow_rate = ode_restart.p.semi.systems[2].cache.pressure_reference_values[2].flow_rate[] + + @test isapprox(restart_pressure, + open_boundary.cache.pressure_reference_values[2].pressure[]) + @test isapprox(restart_flow_rate, + open_boundary.cache.pressure_reference_values[2].flow_rate[]) + end + + @trixi_testset "Oscillating Beam Half-Simulation Restart" begin + # Run full simulation + trixi_include(@__MODULE__, + joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), + tspan=(0.0, 1.0), n_particles_y=5, info_callback=nothing) + + # Store the final solution for comparison + full_sol = sol + positions_full = full_sol.u[end].x[2] + + # Run half simulation + trixi_include(@__MODULE__, + joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), + tspan=(0.0, 0.5), n_particles_y=5) + + # Get latest iteration and create restart + iter = saving_callback.affect!.affect!.latest_saved_iter + restart_file = joinpath("out", "structure_1_$iter.vtu") + + ode_restart = semidiscretize(semi, (0.5, 1.0); + restart_with=restart_file) + + # Create a new saving callback for the restart + saving_callback_restart = SolutionSavingCallback(dt=0.02, prefix="restart"; + deflection_x, deflection_y) + callbacks_restart = CallbackSet(saving_callback_restart) + + sol_restart = solve(ode_restart, RDPK3SpFSAL49(), save_everystep=false, + callback=callbacks_restart) + + # Compare the final particle velocities + @test isapprox(sol_restart.u[end].x[2], positions_full, rtol=1e-6) + end +end