diff --git a/NEWS.md b/NEWS.md index f004f1410..fc9b832e3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### Features +- New viscosity models `ViscosityMorrisSGS` and `ViscosityAdamiSGS` were added, which use a simplified Smagorinsky-type SGS (#753). + - With all CPU backends, a new array type is used for the integration array, which defines broadcasting to be multithreaded, leading to speedups of up to 5x with large thread counts when combined with thread pinning (#722). diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 10c1b2ef8..b857c7a12 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -715,7 +715,7 @@ @article{MeloEspinosa2014 journal = {Energy Conversion and Management}, volume = {84}, pages = {50--60}, - doi = {https://doi.org/10.1016/j.enconman.2014.08.004}, + doi = {10.1016/j.enconman.2014.08.004}, year = {2014} } @@ -736,6 +736,28 @@ @book{Poling2001 address = {New York} } +@article{Smagorinsky1963, + author = {Smagorinsky, Joseph}, + title = {General Circulation Experiments with the Primitive Equations. I. The Basic Experiment}, + journal = {Monthly Weather Review}, + volume = {91}, + number = {3}, + pages = {99--164}, + year = {1963}, + doi = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2} +} + +@inproceedings{Lilly1967, + author = {Lilly, Douglas K.}, + title = {The Representation of Small‑Scale Turbulence in Numerical Simulation Experiments}, + booktitle = {Proceedings of the {IBM} Scientific Computing Symposium on Environmental Sciences}, + editor = {Goldstine, Herman H.}, + address = {Yorktown Heights, New York}, + pages = {195--210}, + year = {1967}, + doi = {10.5065/D62R3PMM} +} + @article{Zhu2021, author = {Zhu, Yujie and Zhang, Chi and Yu, Yongchuan and Hu, Xiangyu}, journal = {Journal of Hydrodynamics}, diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 35493a8c2..950f8c806 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -47,10 +47,21 @@ smoothing_length = 1.75 * fluid_particle_spacing smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() -viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) -# nu = 0.02 * smoothing_length * sound_speed/8 -# viscosity = ViscosityMorris(nu=nu) -# viscosity = ViscosityAdami(nu=nu) +alpha = 0.02 +viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +# A typical formula to convert Artificial viscosity to a +# kinematic viscosity is provided by Monaghan as +# nu = alpha * smoothing_length * sound_speed/8 + +# Alternatively a kinematic viscosity for water can be set +# nu = 1.0e-6 + +# This allows the use of a physical viscosity model like: +# viscosity_fluid = ViscosityAdami(nu=nu) +# or with additional dissipation through a Smagorinsky model +# viscosity_fluid = ViscosityAdamiSGS(nu=nu) +# For more details see the documentation "Viscosity model overview". + # Alternatively the density diffusion model by Molteni & Colagrossi can be used, # which will run faster. # density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) @@ -58,7 +69,7 @@ density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, - smoothing_length, viscosity=viscosity, + smoothing_length, viscosity=viscosity_fluid, density_diffusion=density_diffusion, acceleration=(0.0, -gravity), correction=nothing, surface_tension=nothing, @@ -67,12 +78,17 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() +viscosity_wall = nothing +# For a no-slip condition the corresponding wall viscosity without SGS can be set +# viscosity_wall = ViscosityAdami(nu=nu) +# viscosity_wall = ViscosityMorris(nu=nu) boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, smoothing_kernel, smoothing_length, correction=nothing, - reference_particle_spacing=0) + reference_particle_spacing=0, + viscosity=viscosity_wall) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coefficient=0.0) diff --git a/examples/fluid/dam_break_2phase_2d.jl b/examples/fluid/dam_break_2phase_2d.jl index de99a5f39..247462b52 100644 --- a/examples/fluid/dam_break_2phase_2d.jl +++ b/examples/fluid/dam_break_2phase_2d.jl @@ -37,7 +37,7 @@ water_viscosity = ViscosityMorris(nu=nu_sim_water) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, - viscosity=water_viscosity, smoothing_length=smoothing_length, + viscosity_fluid=water_viscosity, smoothing_length=smoothing_length, gravity=gravity, tspan=tspan, density_diffusion=nothing, sound_speed=sound_speed, exponent=7, tank_size=(floor(5.366 * H / fluid_particle_spacing) * fluid_particle_spacing, diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index 320f18509..f301f9964 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -30,7 +30,8 @@ oil_viscosity = ViscosityMorris(nu=nu_sim_oil) # TODO: broken if both systems use surface tension trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, tspan=tspan, - viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length, + viscosity_fluid=ViscosityMorris(nu=nu_sim_water), + smoothing_length=smoothing_length, gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed, prefix="", reference_particle_spacing=fluid_particle_spacing) diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index a59d25e25..0eaa5d4e7 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -87,8 +87,7 @@ semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=50) -saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", - prefix="", write_meta_data=true) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="") callbacks = CallbackSet(info_callback, saving_callback) diff --git a/examples/fluid/hydrostatic_water_column_2d.jl b/examples/fluid/hydrostatic_water_column_2d.jl index b49229075..f60a76f96 100644 --- a/examples/fluid/hydrostatic_water_column_2d.jl +++ b/examples/fluid/hydrostatic_water_column_2d.jl @@ -32,7 +32,7 @@ smoothing_length = 1.2 * fluid_particle_spacing smoothing_kernel = SchoenbergCubicSplineKernel{2}() alpha = 0.02 -viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) fluid_density_calculator = ContinuityDensity() @@ -40,7 +40,7 @@ fluid_density_calculator = ContinuityDensity() system_acceleration = (0.0, -gravity) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, - smoothing_length, viscosity=viscosity, + smoothing_length, viscosity=viscosity_fluid, acceleration=system_acceleration, source_terms=nothing) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index c1f0f72fd..75e638718 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -3,7 +3,7 @@ # P. Ramachandran, K. Puri # "Entropically damped artificial compressibility for SPH". # In: Computers and Fluids, Volume 179 (2019), pages 579-594. -# https://doi.org/10.1016/j.compfluid.2018.11.023 +# https://doi.org/10.1016/j.compfluid.2018.11.023 using TrixiParticles using OrdinaryDiffEq diff --git a/examples/fsi/falling_sphere_3d.jl b/examples/fsi/falling_sphere_3d.jl index d44911409..990a8bdf5 100644 --- a/examples/fsi/falling_sphere_3d.jl +++ b/examples/fsi/falling_sphere_3d.jl @@ -10,5 +10,4 @@ trixi_include(@__MODULE__, solid_smoothing_kernel=WendlandC2Kernel{3}(), sphere_type=RoundSphere(), output_directory="out", prefix="", - write_meta_data=false, # Files with meta data can't be read by meshio tspan=(0.0, 1.0), abstol=1e-6, reltol=1e-3) diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index 0c3213c98..1e3f06843 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/examples/fsi/falling_spheres_2d.jl @@ -117,8 +117,7 @@ semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_s ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) -saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="", - write_meta_data=true) +saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="") callbacks = CallbackSet(info_callback, saving_callback) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index 06e760fbe..87660fe81 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -105,7 +105,7 @@ end TrixiParticles.vtkname(system::NBodySystem) = "n-body" -function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem; write_meta_data=true) +function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem) (; mass) = system vtk["velocity"] = v diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bab7e7a81..ff23f9490 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -73,7 +73,8 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel export StateEquationCole, StateEquationIdealGas -export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris +export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, + ViscosityMorrisSGS export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono export tensile_instability_control diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 9370bcdad..7fdce08bf 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -216,6 +216,8 @@ function initialize_postprocess_callback!(cb::PostprocessCallback, u, t, integra # Apply the callback cb(integrator) + write_meta_data(cb, integrator) + return cb end diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index f9d8244f2..deab0a0d1 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -2,8 +2,7 @@ SolutionSavingCallback(; interval::Integer=0, dt=0.0, save_times=Array{Float64, 1}([]), save_initial_solution=true, save_final_solution=true, output_directory="out", append_timestamp=false, prefix="", - verbose=false, write_meta_data=true, max_coordinates=2^15, - custom_quantities...) + verbose=false, max_coordinates=2^15, custom_quantities...) Callback to save the current numerical solution in VTK format in regular intervals. @@ -28,7 +27,6 @@ To ignore a custom quantity for a specific system, return `nothing`. - `append_timestamp=false`: Append current timestamp to the output directory. - 'prefix=""': Prefix added to the filename. - `custom_quantities...`: Additional user-defined quantities. -- `write_meta_data=true`: Write meta data. - `verbose=false`: Print to standard IO when a file is written. - `max_coordinates=2^15`: The coordinates of particles will be clipped if their absolute values exceed this threshold. @@ -70,7 +68,6 @@ mutable struct SolutionSavingCallback{I, CQ} save_times :: Vector{Float64} save_initial_solution :: Bool save_final_solution :: Bool - write_meta_data :: Bool verbose :: Bool output_directory :: String prefix :: String @@ -84,8 +81,8 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0, save_times=Float64[], save_initial_solution=true, save_final_solution=true, output_directory="out", append_timestamp=false, - prefix="", verbose=false, write_meta_data=true, - max_coordinates=Float64(2^15), custom_quantities...) + prefix="", verbose=false, max_coordinates=Float64(2^15), + custom_quantities...) if (dt > 0 && interval > 0) || (length(save_times) > 0 && (dt > 0 || interval > 0)) throw(ArgumentError("Setting multiple save times for the same solution " * "callback is not possible. Use either `dt`, `interval` or `save_times`.")) @@ -101,8 +98,8 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0, solution_callback = SolutionSavingCallback(interval, Float64.(save_times), save_initial_solution, save_final_solution, - write_meta_data, verbose, output_directory, - prefix, max_coordinates, custom_quantities, + verbose, output_directory, prefix, + max_coordinates, custom_quantities, -1, Ref("UnknownVersion")) if length(save_times) > 0 @@ -133,6 +130,8 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in solution_callback.latest_saved_iter = -1 solution_callback.git_hash[] = compute_git_hash() + write_meta_data(solution_callback, integrator) + # Save initial solution if solution_callback.save_initial_solution solution_callback(integrator) @@ -151,8 +150,8 @@ end # `affect!` function (solution_callback::SolutionSavingCallback)(integrator) - (; interval, output_directory, custom_quantities, write_meta_data, git_hash, - verbose, prefix, latest_saved_iter, max_coordinates) = solution_callback + (; interval, output_directory, custom_quantities, git_hash, verbose, + prefix, latest_saved_iter, max_coordinates) = solution_callback vu_ode = integrator.u semi = integrator.p @@ -175,8 +174,8 @@ function (solution_callback::SolutionSavingCallback)(integrator) @trixi_timeit timer() "save solution" trixi2vtk(vu_ode, semi, integrator.t; iter, output_directory, prefix, - write_meta_data, git_hash=git_hash[], - max_coordinates, custom_quantities...) + git_hash=git_hash[], max_coordinates, + custom_quantities...) # Tell OrdinaryDiffEq that `u` has not been modified u_modified!(integrator, false) diff --git a/src/io/io.jl b/src/io/io.jl index 380151bcf..4a0f7e51c 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -1,2 +1,157 @@ include("write_vtk.jl") include("read_vtk.jl") + +function write_meta_data(callback::Union{SolutionSavingCallback, PostprocessCallback}, + integrator) + # handle "_" on optional prefix strings + add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") + + git_hash = callback.git_hash + prefix = hasproperty(callback, :prefix) ? callback.prefix : "" + + semi = integrator.p + names = system_names(semi.systems) + + # fill `systems` with metadata for each system + systems = Dict{String, Any}() + foreach_system(semi) do system + idx = system_indices(system, semi) + name = add_opt_str_pre(prefix) * "$(names[idx])" + + meta_data = Dict{String, Any}() + add_meta_data!(meta_data, system) + + systems[name] = meta_data + end + + # initialize `simulation_meta_data` and add `systems` + simulation_meta_data = Dict{String, Any}( + "info" => Dict( + "solver_name" => "TrixiParticles.jl", + "solver_version" => git_hash, + "julia_version" => string(VERSION) + ), + "systems" => systems + ) + + # write JSON-file + output_directory = callback.output_directory + mkpath(output_directory) + + json_file = joinpath(output_directory, "simulation_metadata.json") + + open(json_file, "w") do file + JSON.print(file, simulation_meta_data, 2) + end +end + +function add_meta_data!(meta_data, system) + return meta_data +end + +function add_meta_data!(meta_data, system::DEMSystem) + return meta_data +end + +function add_meta_data!(meta_data, system::WeaklyCompressibleSPHSystem) + # general `FLuidSystem`-Metadata + meta_data["acceleration"] = system.acceleration + meta_data["viscosity"] = type2string(system.viscosity) + add_meta_data!(meta_data, system.viscosity) + meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) + meta_data["smoothing_length_factor"] = system.cache.smoothing_length_factor + meta_data["density_calculator"] = type2string(system.density_calculator) + + # specific `WCSPH`-Metadata + meta_data["solver"] = "WCSPH" + meta_data["state_equation"] = type2string(system.state_equation) + meta_data["state_equation_rho0"] = system.state_equation.reference_density + meta_data["state_equation_pa"] = system.state_equation.background_pressure + meta_data["state_equation_c"] = system.state_equation.sound_speed + meta_data["correction_method"] = type2string(system.correction) + if system.correction isa AkinciFreeSurfaceCorrection + meta_data["correction_rho0"] = system.correction.rho0 + end + if system.state_equation isa StateEquationCole + meta_data["state_equation_exponent"] = system.state_equation.exponent + end + if system.state_equation isa StateEquationIdealGas + meta_data["state_equation_gamma"] = system.state_equation.gamma + end +end + +function add_meta_data!(meta_data, system::EntropicallyDampedSPHSystem) + # general `FLuidSystem`-Metadata + meta_data["acceleration"] = system.acceleration + meta_data["viscosity"] = type2string(system.viscosity) + add_meta_data!(meta_data, system.viscosity) + meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) + meta_data["smoothing_length_factor"] = system.cache.smoothing_length_factor + meta_data["density_calculator"] = type2string(system.density_calculator) + + # specific `EDAC`-Metadata + meta_data["solver"] = "EDAC" + meta_data["sound_speed"] = system.sound_speed + meta_data["background_pressure_TVF"] = system.transport_velocity isa Nothing ? "-" : + system.transport_velocity.background_pressure +end + +add_meta_data!(meta_data, viscosity::Nothing) = meta_data + +function add_meta_data!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) + meta_data["viscosity_nu"] = viscosity.nu + meta_data["viscosity_epsilon"] = viscosity.epsilon +end + +function add_meta_data!(meta_data, viscosity::ArtificialViscosityMonaghan) + meta_data["viscosity_alpha"] = viscosity.alpha + meta_data["viscosity_beta"] = viscosity.beta + meta_data["viscosity_epsilon"] = viscosity.epsilon +end + +function add_meta_data!(meta_data, system::TotalLagrangianSPHSystem) + meta_data["young_modulus"] = system.young_modulus + meta_data["poisson_ratio"] = system.poisson_ratio + meta_data["lame_lambda"] = system.lame_lambda + meta_data["lame_mu"] = system.lame_mu + meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) + meta_data["smoothing_length_factor"] = initial_smoothing_length(system) / + particle_spacing(system, 1) + + add_meta_data!(meta_data, system.boundary_model, system) +end + +function add_meta_data!(meta_data, system::OpenBoundarySPHSystem) + meta_data["boundary_zone"] = type2string(system.boundary_zone) + meta_data["width"] = round(system.boundary_zone.zone_width, digits=3) + meta_data["velocity_function"] = type2string(system.reference_velocity) + meta_data["pressure_function"] = type2string(system.reference_pressure) + meta_data["density_function"] = type2string(system.reference_density) +end + +function add_meta_data!(meta_data, system::BoundarySPHSystem) + add_meta_data!(meta_data, system.boundary_model, system) +end + +function add_meta_data!(meta_data, model::Nothing, system) + return meta_data +end + +function add_meta_data!(meta_data, model::BoundaryModelMonaghanKajtar, system) + meta_data["boundary_model"] = "BoundaryModelMonaghanKajtar" + meta_data["boundary_spacing_ratio"] = model.beta + meta_data["boundary_K"] = model.K +end + +function add_meta_data!(meta_data, model::BoundaryModelDummyParticles, system) + meta_data["boundary_model"] = "BoundaryModelDummyParticles" + meta_data["smoothing_kernel"] = type2string(model.smoothing_kernel) + meta_data["smoothing_length"] = model.smoothing_length + meta_data["density_calculator"] = type2string(model.density_calculator) + meta_data["state_equation"] = type2string(model.state_equation) + meta_data["viscosity_model"] = type2string(model.viscosity) +end + +function add_meta_data!(meta_data, system::BoundaryDEMSystem) + return meta_data +end diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index 0e5aaf4c1..729baf8dd 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -10,7 +10,7 @@ end """ trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", - write_meta_data=true, max_coordinates=Inf, custom_quantities...) + max_coordinates=Inf, custom_quantities...) Convert Trixi simulation data to VTK format. @@ -25,7 +25,6 @@ Convert Trixi simulation data to VTK format. separate files. This number is just appended to the filename. - `output_directory="out"`: Output directory path. - `prefix=""`: Prefix for output files. -- `write_meta_data=true`: Write meta data. - `max_coordinates=Inf` The coordinates of particles will be clipped if their absolute values exceed this threshold. - `custom_quantities...`: Additional custom quantities to include in the VTK output. @@ -49,8 +48,7 @@ trixi2vtk(sol.u[end], semi, 0.0, iter=1, my_custom_quantity=kinetic_energy) ``` """ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", - write_meta_data=true, git_hash=compute_git_hash(), - max_coordinates=Inf, custom_quantities...) + git_hash=compute_git_hash(), max_coordinates=Inf, custom_quantities...) (; systems) = semi v_ode, u_ode = vu_ode.x @@ -67,14 +65,14 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix trixi2vtk(system, v_ode, u_ode, semi, t, periodic_box; system_name=filenames[system_index], output_directory, iter, prefix, - write_meta_data, git_hash, max_coordinates, custom_quantities...) + git_hash, max_coordinates, custom_quantities...) end end # Convert data for a single TrixiParticle system to VTK format function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_directory="out", prefix="", iter=nothing, system_name=vtkname(system_), - write_meta_data=true, max_coordinates=Inf, git_hash=compute_git_hash(), + max_coordinates=Inf, git_hash=compute_git_hash(), custom_quantities...) mkpath(output_directory) @@ -118,7 +116,7 @@ function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_direc @trixi_timeit timer() "write to vtk" vtk_grid(file, points, cells) do vtk # Dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem - write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) + write2vtk!(vtk, v, u, t, system) # Store particle index vtk["index"] = active_particles(system) @@ -128,11 +126,6 @@ function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_direc vtk["particle_spacing"] = [particle_spacing(system, particle) for particle in active_particles(system)] - if write_meta_data - vtk["solver_version"] = git_hash - vtk["julia_version"] = string(VERSION) - end - # Extract custom quantities for this system for (key, quantity) in custom_quantities value = custom_quantity(quantity, system, v_ode, u_ode, semi, t) @@ -251,13 +244,13 @@ function trixi2vtk(initial_condition::InitialCondition; output_directory="out", pressure=pressure, custom_quantities...) end -function write2vtk!(vtk, v, u, t, system; write_meta_data=true) +function write2vtk!(vtk, v, u, t, system) vtk["velocity"] = view(v, 1:ndims(system), :) return vtk end -function write2vtk!(vtk, v, u, t, system::DEMSystem; write_meta_data=true) +function write2vtk!(vtk, v, u, t, system::DEMSystem) vtk["velocity"] = view(v, 1:ndims(system), :) vtk["mass"] = [hydrodynamic_mass(system, particle) for particle in active_particles(system)] @@ -266,7 +259,7 @@ function write2vtk!(vtk, v, u, t, system::DEMSystem; write_meta_data=true) return vtk end -function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) +function write2vtk!(vtk, v, u, t, system::FluidSystem) vtk["velocity"] = [current_velocity(v, system, particle) for particle in active_particles(system)] vtk["density"] = current_density(v, system) @@ -314,50 +307,14 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) end end - if write_meta_data - vtk["acceleration"] = system.acceleration - vtk["viscosity"] = type2string(system.viscosity) - write2vtk!(vtk, system.viscosity) - vtk["smoothing_kernel"] = type2string(system.smoothing_kernel) - vtk["smoothing_length_factor"] = system.cache.smoothing_length_factor - vtk["density_calculator"] = type2string(system.density_calculator) - - if system isa WeaklyCompressibleSPHSystem - vtk["solver"] = "WCSPH" - - vtk["correction_method"] = type2string(system.correction) - if system.correction isa AkinciFreeSurfaceCorrection - vtk["correction_rho0"] = system.correction.rho0 - end - - if system.state_equation isa StateEquationCole - vtk["state_equation_exponent"] = system.state_equation.exponent - end - - if system.state_equation isa StateEquationIdealGas - vtk["state_equation_gamma"] = system.state_equation.gamma - end - - vtk["state_equation"] = type2string(system.state_equation) - vtk["state_equation_rho0"] = system.state_equation.reference_density - vtk["state_equation_pa"] = system.state_equation.background_pressure - vtk["state_equation_c"] = system.state_equation.sound_speed - vtk["solver"] = "WCSPH" - else - vtk["solver"] = "EDAC" - vtk["sound_speed"] = system.sound_speed - vtk["background_pressure_TVF"] = system.transport_velocity isa Nothing ? - "-" : - system.transport_velocity.background_pressure - end - end - return vtk end write2vtk!(vtk, viscosity::Nothing) = vtk -function write2vtk!(vtk, viscosity::Union{ViscosityAdami, ViscosityMorris}) +function write2vtk!(vtk, + viscosity::Union{ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, + ViscosityMorrisSGS}) vtk["viscosity_nu"] = viscosity.nu vtk["viscosity_epsilon"] = viscosity.epsilon end @@ -368,7 +325,7 @@ function write2vtk!(vtk, viscosity::ArtificialViscosityMonaghan) vtk["viscosity_epsilon"] = viscosity.epsilon end -function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_data=true) +function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem) n_fixed_particles = nparticles(system) - n_moving_particles(system) vtk["velocity"] = hcat(view(v, 1:ndims(system), :), @@ -387,62 +344,31 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d vtk["material_density"] = system.material_density - if write_meta_data - vtk["lame_lambda"] = system.lame_lambda - vtk["lame_mu"] = system.lame_mu - vtk["smoothing_kernel"] = type2string(system.smoothing_kernel) - vtk["smoothing_length_factor"] = initial_smoothing_length(system) / - particle_spacing(system, 1) - end - - write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) + write2vtk!(vtk, v, u, t, system.boundary_model, system) end -function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data=true) +function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem) vtk["velocity"] = [current_velocity(v, system, particle) for particle in active_particles(system)] vtk["density"] = current_density(v, system) vtk["pressure"] = current_pressure(v, system) - if write_meta_data - vtk["boundary_zone"] = type2string(first(typeof(system.boundary_zone).parameters)) - vtk["width"] = round(system.boundary_zone.zone_width, digits=3) - vtk["velocity_function"] = type2string(system.reference_velocity) - vtk["pressure_function"] = type2string(system.reference_pressure) - vtk["density_function"] = type2string(system.reference_density) - end - return vtk end -function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem; write_meta_data=true) - write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) +function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem) + write2vtk!(vtk, v, u, t, system.boundary_model, system) end -function write2vtk!(vtk, v, u, t, model::Nothing, system; write_meta_data=true) +function write2vtk!(vtk, v, u, t, model::Nothing, system) return vtk end -function write2vtk!(vtk, v, u, t, model::BoundaryModelMonaghanKajtar, system; - write_meta_data=true) - if write_meta_data - vtk["boundary_model"] = "BoundaryModelMonaghanKajtar" - vtk["boundary_spacing_ratio"] = model.beta - vtk["boundary_K"] = model.K - end +function write2vtk!(vtk, v, u, t, model::BoundaryModelMonaghanKajtar, system) + return vtk end -function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system; - write_meta_data=true) - if write_meta_data - vtk["boundary_model"] = "BoundaryModelDummyParticles" - vtk["smoothing_kernel"] = type2string(model.smoothing_kernel) - vtk["smoothing_length"] = model.smoothing_length - vtk["density_calculator"] = type2string(model.density_calculator) - vtk["state_equation"] = type2string(model.state_equation) - vtk["viscosity_model"] = type2string(model.viscosity) - end - +function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system) vtk["hydrodynamic_density"] = current_density(v, system) vtk["pressure"] = model.pressure @@ -457,6 +383,6 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system; end end -function write2vtk!(vtk, v, u, t, system::BoundaryDEMSystem; write_meta_data=true) +function write2vtk!(vtk, v, u, t, system::BoundaryDEMSystem) return vtk end diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index be6697fb1..8b713bd40 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -223,12 +223,14 @@ end update_callback_used!(system::ParticlePackingSystem) = system.update_callback_used[] = true -function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true) +function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem) vtk["velocity"] = [advection_velocity(v, system, particle) for particle in active_particles(system)] - if write_meta_data - vtk["signed_distances"] = system.signed_distances - end + vtk["signed_distances"] = system.signed_distances +end + +function add_meta_data!(meta_data, system::ParticlePackingSystem) + return meta_data end # Skip for fixed systems diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index e39e7205d..ddd74a748 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -176,6 +176,32 @@ struct ViscosityAdami{ELTYPE} end end +function adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, + m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) + eta_a = nu_a * rho_a + eta_b = nu_b * rho_b + + eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) + + tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) + + volume_a = m_a / rho_a + volume_b = m_b / rho_b + + # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 + # They argued that the formulation is more flexible because of the possibility to formulate + # different inter-particle averages or to assume different inter-particle distributions. + # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration. + # + # TODO: Is there a better formulation to discretize the Laplace operator? + # Because when using this formulation for the pressure acceleration, it is not + # energy conserving. + # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 + visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a + + return visc .* v_diff +end + @inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, @@ -185,6 +211,7 @@ end smoothing_length_particle = smoothing_length(particle_system, particle) smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 nu_a = kinematic_viscosity(particle_system, viscosity_model(neighbor_system, particle_system), @@ -197,36 +224,225 @@ end v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b - eta_a = nu_a * rho_a - eta_b = nu_b * rho_b + return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, + m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) +end - eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) +function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length, + sound_speed) + return viscosity.nu +end +@propagate_inbounds function viscous_velocity(v, system, particle) + return current_velocity(v, system, particle) +end + +@doc raw""" + ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) + +Viscosity model that extends the standard [Adami formulation](@ref ViscosityAdami) +by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure. +The effective kinematic viscosity is defined as + +```math +\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}}, +``` + +with + +```math +\nu_{\text{SGS}} = (C_S h)^2 |S|, +``` + +and an approximation for the strain rate magnitude given by + +```math +|S| \approx \frac{\|v_{ab}\|}{\|r_{ab}\| + \epsilon}, +``` + +where: +- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2), +- ``h`` is the local smoothing length. + +The effective dynamic viscosities are then computed as +```math +\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, +``` +and averaged as +```math +\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}. +``` + +This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. + +# Keywords +- `nu`: Standard kinematic viscosity. +- `C_S`: Smagorinsky constant. +- `epsilon=0.01`: Parameter to prevent singularities +""" +struct ViscosityAdamiSGS{ELTYPE} + nu :: ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] + C_S :: ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] + epsilon :: ELTYPE # Epsilon for singularity prevention [e.g., 0.001] +end + +ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon) + +@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system, + neighbor_system, + v_particle_system, + v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + epsilon = viscosity.epsilon + + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 - tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) - volume_a = m_a / rho_a - volume_b = m_b / rho_b + nu_a = kinematic_viscosity(particle_system, + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle, sound_speed) + nu_b = kinematic_viscosity(neighbor_system, + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor, sound_speed) - # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 - # They argued that the formulation is more flexible because of the possibility to formulate - # different inter-particle averages or to assume different inter-particle distributions. - # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration. + v_a = viscous_velocity(v_particle_system, particle_system, particle) + v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) + v_diff = v_a - v_b + + # ------------------------------------------------------------------------------ + # SGS part: Compute the subgrid-scale eddy viscosity. + # ------------------------------------------------------------------------------ + # In classical LES [Lilly (1967)](@cite Lilly1967) the Smagorinsky model defines: + # ν_SGS = (C_S Δ)^2 |S|, + # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ). # - # TODO: Is there a better formulation to discretize the Laplace operator? - # Because when using this formulation for the pressure acceleration, it is not - # energy conserving. - # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 - visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a + # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly. + # A common low-order surrogate is to approximate the strain‐rate magnitude by a + # finite difference along each particle pair: + # + # |S| ≈ ‖v_ab‖ / (‖r_ab‖ + δ), + # + # where δ regularizes the denominator to avoid singularities when particles are very close. + # + # This yields: + # S_mag = norm(v_diff) / (distance + ε) + # + # and then the Smagorinsky eddy viscosity: + # ν_SGS = (C_S * h̄)^2 * S_mag. + # + S_mag = norm(v_diff) / (distance + epsilon) + nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag - return visc .* v_diff + # Effective kinematic viscosity is the sum of the standard and SGS parts. + nu_a = nu_a + nu_SGS + nu_b = nu_b + nu_SGS + + return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, + m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) end -function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length, +function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length, sound_speed) return viscosity.nu end -@propagate_inbounds function viscous_velocity(v, system, particle) - return current_velocity(v, system, particle) +@doc raw""" + ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) + +Subgrid-scale (SGS) viscosity model based on the formulation by [Morris (1997)](@cite Morris1997), +by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure. +The effective kinematic viscosity is defined as + +```math +\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}}, +``` + +with + +```math +\nu_{\text{SGS}} = (C_S h)^2 |S|, +``` + +and an approximation for the strain rate magnitude given by + +```math +|S| \approx \frac{\|v_{ab}\|}{\|r_{ab}\| + \epsilon}, +``` + +where: +- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2), +- ``h`` is the local smoothing length. + +The effective dynamic viscosities are then computed as +```math +\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, +``` +and averaged as +```math +\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}. +``` + +This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. + +# Keywords +- `nu`: Standard kinematic viscosity. +- `C_S`: Smagorinsky constant. +- `epsilon=0.01`: Parameter to prevent singularities +""" +struct ViscosityMorrisSGS{ELTYPE} + nu::ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] + C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] + epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001] +end + +ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, epsilon) + +@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(particle_system, + neighbor_system, + v_particle_system, + v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, + m_b, rho_a, rho_b, grad_kernel) + epsilon = viscosity.epsilon + + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 + + nu_a = kinematic_viscosity(particle_system, + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle, sound_speed) + nu_b = kinematic_viscosity(neighbor_system, + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor, sound_speed) + + v_a = viscous_velocity(v_particle_system, particle_system, particle) + v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) + v_diff = v_a - v_b + + # SGS part: Compute the subgrid-scale eddy viscosity. + # See comments above for `ViscosityAdamiSGS`. + S_mag = norm(v_diff) / (distance + epsilon) + nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag + + # Effective viscosities include the SGS term. + nu_a_eff = nu_a + nu_SGS + nu_b_eff = nu_b + nu_SGS + + # For the Morris model, dynamic viscosities are: + mu_a = nu_a_eff * rho_a + mu_b = nu_b_eff * rho_b + + force_Morris = (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / + (distance^2 + epsilon * smoothing_length_average^2) * v_diff + return m_b * force_Morris +end + +function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_length, + sound_speed) + return viscosity.nu end diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index 4c557ee46..141d60278 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -25,18 +25,18 @@ clip_negative_pressure=true), "WCSPH with ViscosityAdami" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015),), + viscosity_fluid=ViscosityAdami(nu=0.0015),), "WCSPH with ViscosityMorris" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityMorris(nu=0.0015),), + viscosity_fluid=ViscosityMorris(nu=0.0015),), "WCSPH with ViscosityAdami and SummationDensity" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015), + viscosity_fluid=ViscosityAdami(nu=0.0015), fluid_density_calculator=SummationDensity(), clip_negative_pressure=true), "WCSPH with ViscosityMorris and SummationDensity" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityMorris(nu=0.0015), + viscosity_fluid=ViscosityMorris(nu=0.0015), fluid_density_calculator=SummationDensity(), clip_negative_pressure=true), "WCSPH with smoothing_length=1.3" => (smoothing_length=1.3,), @@ -55,7 +55,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity=viscosity_fluid, density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity))), @@ -63,7 +63,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity=viscosity_fluid, density_calculator=SummationDensity(), acceleration=(0.0, -gravity)),), diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index 9d551a527..85d657b67 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -178,13 +178,13 @@ end clip_negative_pressure=true), "WCSPH with ViscosityAdami" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015f0),), + viscosity_fluid=ViscosityAdami(nu=0.0015f0),), "WCSPH with ViscosityMorris" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityMorris(nu=0.0015f0),), + viscosity_fluid=ViscosityMorris(nu=0.0015f0),), "WCSPH with ViscosityAdami and SummationDensity" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015f0), + viscosity_fluid=ViscosityAdami(nu=0.0015f0), fluid_density_calculator=SummationDensity(), maxiters=38, # 38 time steps on CPU clip_negative_pressure=true), @@ -205,7 +205,7 @@ end smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity=viscosity_fluid, density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity))), @@ -213,7 +213,7 @@ end smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity=viscosity_fluid, density_calculator=SummationDensity(), acceleration=(0.0, -gravity)),) diff --git a/test/schemes/fluid/viscosity.jl b/test/schemes/fluid/viscosity.jl index f11f4b05b..852ebfbc3 100644 --- a/test/schemes/fluid/viscosity.jl +++ b/test/schemes/fluid/viscosity.jl @@ -9,8 +9,8 @@ fluid = rectangular_patch(particle_spacing, (3, 3), seed=1) - v_diff = [0.1, -0.75] - pos_diff = [-0.5 * smoothing_length, 0.75 * smoothing_length] + v_diff = [0.3, -1.0] + pos_diff = [-0.25 * smoothing_length, 0.375 * smoothing_length] distance = norm(pos_diff) rho_a = rho_b = rho_mean = 1000.0 @@ -29,8 +29,8 @@ rho_mean, rho_a, rho_b, smoothing_length, grad_kernel, 0.0, 0.0) - @test isapprox(dv[1], -0.007073849138494646, atol=6e-15) - @test isapprox(dv[2], 0.01061077370774197, atol=6e-15) + @test isapprox(dv[1], -0.02049217623299368, atol=6e-15) + @test isapprox(dv[2], 0.03073826434949052, atol=6e-15) end @testset verbose=true "`ViscosityMorris`" begin nu = 7e-3 @@ -41,13 +41,22 @@ grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, distance, 1) + v = fluid.velocity - dv = viscosity(sound_speed, v_diff, pos_diff, distance, - rho_mean, rho_a, rho_b, smoothing_length, - grad_kernel, nu, nu) + m_a = 0.01 + m_b = 0.01 - @test isapprox(dv[1], -0.00018421886647594437, atol=6e-15) - @test isapprox(dv[2], 0.0013816414985695826, atol=6e-15) + v[1, 1] = v_diff[1] + v[2, 1] = v_diff[2] + v[1, 2] = 0.0 + v[2, 2] = 0.0 + + dv = viscosity(system_wcsph, system_wcsph, + v, v, 1, 2, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + + @test isapprox(dv[1], -1.0895602048035404e-5, atol=6e-15) + @test isapprox(dv[2], 3.631867349345135e-5, atol=6e-15) end @testset verbose=true "`ViscosityAdami`" begin viscosity = ViscosityAdami(nu=7e-3) @@ -71,7 +80,59 @@ v, v, 1, 2, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - @test isapprox(dv[1], -1.8421886647594435e-6, atol=6e-15) - @test isapprox(dv[2], 1.3816414985695826e-5, atol=6e-15) + @test isapprox(dv[1], -1.089560204803541e-5, atol=6e-15) + @test isapprox(dv[2], 3.6318673493451364e-5, atol=6e-15) + end + @testset verbose=true "`ViscosityMorrisSGS`" begin + nu = 7e-3 + viscosity = ViscosityMorrisSGS(nu=nu) + system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(), + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity) + + grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, + distance, 1) + + v = fluid.velocity + + m_a = 0.01 + m_b = 0.01 + + v[1, 1] = v_diff[1] + v[2, 1] = v_diff[2] + v[1, 2] = 0.0 + v[2, 2] = 0.0 + + dv = viscosity(system_wcsph, system_wcsph, + v, v, 1, 2, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + + @test isapprox(dv[1], -2.032835697804103e-5, atol=6e-15) + @test isapprox(dv[2], 6.776118992680343e-5, atol=6e-15) + end + @testset verbose=true "`ViscosityAdamiSGS`" begin + viscosity = ViscosityAdamiSGS(nu=7e-3) + system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(), + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity) + + grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, + distance, 1) + v = fluid.velocity + + m_a = 0.01 + m_b = 0.01 + + v[1, 1] = v_diff[1] + v[2, 1] = v_diff[2] + v[1, 2] = 0.0 + v[2, 2] = 0.0 + + dv = viscosity(system_wcsph, system_wcsph, + v, v, 1, 2, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + + @test isapprox(dv[1], -2.0328356978041036e-5, atol=6e-15) + @test isapprox(dv[2], 6.776118992680346e-5, atol=6e-15) end end diff --git a/test/validation/validation.jl b/test/validation/validation.jl index d1c842084..70aceb4bf 100644 --- a/test/validation/validation.jl +++ b/test/validation/validation.jl @@ -50,7 +50,7 @@ if Sys.ARCH === :aarch64 # MacOS ARM produces slightly different pressure values than x86. # Note that pressure values are in the order of 1e5. - @test isapprox(error_edac_P1, 0, atol=4e-6) + @test isapprox(error_edac_P1, 0, atol=5e-6) @test isapprox(error_edac_P2, 0, atol=4e-11) @test isapprox(error_wcsph_P1, 0, atol=400.0) @test isapprox(error_wcsph_P2, 0, atol=0.03) diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index 73884e38c..454341230 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -2,6 +2,10 @@ using TrixiParticles # ========================================================================================== # ==== Resolution +# P. Ramachandran, K. Puri +# "Entropically damped artificial compressibility for SPH". +# In: Computers and Fluids, Volume 179 (2019), pages 579-594. +# https://doi.org/10.1016/j.compfluid.2018.11.023 # The paper provides reference data for particle spacings particle_spacings = [0.02, 0.01, 0.005] particle_spacing = 0.02