From db890674f44c307639f2db82f3aed14174782789 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 24 Mar 2025 08:39:28 +0100 Subject: [PATCH 01/18] Integrate write2json for initial metadata export --- src/TrixiParticles.jl | 3 +- src/callbacks/solution_saving.jl | 5 + src/visualization/write2json.jl | 152 +++++++++++++++++++++++++++++++ 3 files changed, 159 insertions(+), 1 deletion(-) create mode 100644 src/visualization/write2json.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 07fa175744..165b29a48b 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -52,6 +52,7 @@ include("schemes/schemes.jl") include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") +include("visualization/write2json.jl") include("visualization/recipes_plots.jl") include("preprocessing/preprocessing.jl") @@ -77,7 +78,7 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx export BoundaryMovement export examples_dir, validation_dir -export trixi2vtk +export trixi2vtk, trixi2json export RectangularTank, RectangularShape, SphereShape, ComplexShape export ParticlePackingSystem, SignedDistanceField export WindingNumberHormann, WindingNumberJacobson diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index fa06e8b6d3..742fd8a460 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -133,6 +133,11 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in solution_callback.latest_saved_iter = -1 solution_callback.git_hash[] = compute_git_hash() + # Write metadata to JSON-file + if integrator.iter == 0 + trixi2json(solution_callback, integrator) + end + # Save initial solution if solution_callback.save_initial_solution # Update systems to compute quantities like density and pressure diff --git a/src/visualization/write2json.jl b/src/visualization/write2json.jl new file mode 100644 index 0000000000..225138bfd3 --- /dev/null +++ b/src/visualization/write2json.jl @@ -0,0 +1,152 @@ +""" + trixi2json(solution_callback, integrator) + +Write simulation metadata to a JSON-file. + +# Arguments +- `solution_callback`: Callback storing metadata and output settings. +- `integrator`: ODE integrator containing simulation data. + +# Example +```jldoctest; output = false, saving_callback = SolutionSavingCallback(dt = 0.1, output_directory = "output", prefix = "solution"), +setup = :(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), tspan = (0.0, 0.01), callbacks = saving_callback)) + +# output + +``` +""" + +function trixi2json(solution_callback, integrator) + + semi = integrator.p + (; systems) = semi + + output_directory = solution_callback.output_directory + prefix = solution_callback.prefix + git_hash = solution_callback.git_hash + + filenames = system_names(systems) + + foreach_system(semi) do system + system_index = system_indices(system, semi) + + trixi2json(system; system_name = filenames[system_index], output_directory, prefix, git_hash) + end +end + +function trixi2json(system; system_name, output_directory, prefix, git_hash) + mkpath(output_directory) + + meta_data = Dict{String, Any}( + "solver_version" => git_hash, + "julia_version" => string(VERSION), + ) + + get_meta_data!(meta_data, system) + + # handle "_" on optional prefix strings + add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") + + # Write metadata to JSON-file + json_file = joinpath(output_directory, + add_opt_str_pre(prefix) * "$(system_name)_metadata.json") + + open(json_file, "w") do file + JSON.print(file, meta_data, 2) + end +end + +function get_meta_data!(meta_data, system::FluidSystem) + meta_data["acceleration"] = system.acceleration + meta_data["viscosity"] = type2string(system.viscosity) + get_meta_data!(meta_data, system.viscosity) + meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) + meta_data["smoothing_length"] = system.smoothing_length + meta_data["density_calculator"] = type2string(system.density_calculator) + + if system isa WeaklyCompressibleSPHSystem + 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["solver"] = "WCSPH" + + 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 + else + 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 + + return meta_data +end + +get_meta_data!(meta_data, viscosity::Nothing) = meta_data + +function get_meta_data!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) + meta_data["viscosity_nu"] = viscosity.nu + meta_data["viscosity_epsilon"] = viscosity.epsilon +end + +function get_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 get_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"] = system.smoothing_length + + get_meta_data!(meta_data, system.boundary_model, system) +end + +function get_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["flow_direction"] = system.flow_direction + 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 get_meta_data!(meta_data, system::BoundarySPHSystem) + get_meta_data!(meta_data, system.boundary_model, system) +end + +function get_meta_data!(meta_data, model, system) + return meta_data +end + +function get_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 get_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 get_meta_data!(meta_data, system::BoundaryDEMSystem) + return meta_data +end From d8906ee1c1377ad114f14045e7fd4c16a015cd21 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 24 Mar 2025 10:58:24 +0100 Subject: [PATCH 02/18] Remove old `write_meta_data` parameter from VTK writing functions and callbacks --- examples/fluid/falling_water_spheres_2d.jl | 3 +- examples/fsi/falling_sphere_3d.jl | 1 - examples/fsi/falling_spheres_2d.jl | 3 +- examples/n_body/n_body_system.jl | 2 +- src/callbacks/solution_saving.jl | 21 ++-- src/preprocessing/particle_packing/system.jl | 4 +- src/visualization/write2vtk.jl | 112 +++---------------- 7 files changed, 31 insertions(+), 115 deletions(-) diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index a59d25e259..0eaa5d4e70 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/fsi/falling_sphere_3d.jl b/examples/fsi/falling_sphere_3d.jl index d44911409e..990a8bdf5b 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 a06056f8f7..4444df5b72 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 2add374d26..9c40b415d4 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -110,7 +110,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/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index 742fd8a460..871ea2ef4b 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 @@ -162,8 +159,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 @@ -186,8 +183,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/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index a5b9db6007..6ac8713c08 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -142,10 +142,8 @@ end update_callback_used!(system::ParticlePackingSystem) = system.update_callback_used[] = true -function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true) - if write_meta_data +function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem) vtk["signed_distances"] = system.signed_distances - end end write_v0!(v0, system::ParticlePackingSystem) = v0 .= zero(eltype(system)) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index dade585ede..fa2a09e4e6 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.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 @@ -70,15 +68,14 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix trixi2vtk(v, u, t, system, 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(v_, u_, t, system_, periodic_box; output_directory="out", prefix="", - iter=nothing, system_name=vtkname(system_), write_meta_data=true, - max_coordinates=Inf, git_hash=compute_git_hash(), - custom_quantities...) + iter=nothing, system_name=vtkname(system_), max_coordinates=Inf, + git_hash=compute_git_hash(), custom_quantities...) mkpath(output_directory) # Skip empty systems @@ -118,17 +115,12 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre @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) vtk["time"] = t - 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, v, u, t, system) @@ -233,13 +225,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::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"] = [particle_density(v, system, particle) @@ -253,42 +245,6 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["neighbor_count"] = system.cache.neighbor_count 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"] = system.smoothing_length - vtk["density_calculator"] = type2string(system.density_calculator) - - if system isa WeaklyCompressibleSPHSystem - 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 @@ -305,7 +261,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), :), @@ -324,19 +280,10 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d vtk["material_density"] = system.material_density - if write_meta_data - vtk["young_modulus"] = system.young_modulus - vtk["poisson_ratio"] = system.poisson_ratio - vtk["lame_lambda"] = system.lame_lambda - vtk["lame_mu"] = system.lame_mu - vtk["smoothing_kernel"] = type2string(system.smoothing_kernel) - vtk["smoothing_length"] = system.smoothing_length - 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"] = [particle_density(v, system, particle) @@ -344,45 +291,22 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data vtk["pressure"] = [particle_pressure(v, system, particle) for particle in active_particles(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, system; write_meta_data=true) +function write2vtk!(vtk, v, u, t, model, 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"] = [particle_density(v, system, particle) for particle in eachparticle(system)] vtk["pressure"] = model.pressure @@ -398,6 +322,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 From 8558b13f4fde576d922745a727ff69d4c93fa6b2 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 24 Mar 2025 11:00:57 +0100 Subject: [PATCH 03/18] Add `solver_name` to metadata and format files --- src/preprocessing/particle_packing/system.jl | 2 +- src/visualization/write2json.jl | 177 ++++++++++--------- 2 files changed, 90 insertions(+), 89 deletions(-) diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index 6ac8713c08..a89a5c63a3 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -143,7 +143,7 @@ end update_callback_used!(system::ParticlePackingSystem) = system.update_callback_used[] = true function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem) - vtk["signed_distances"] = system.signed_distances + vtk["signed_distances"] = system.signed_distances end write_v0!(v0, system::ParticlePackingSystem) = v0 .= zero(eltype(system)) diff --git a/src/visualization/write2json.jl b/src/visualization/write2json.jl index 225138bfd3..dc7fb344ac 100644 --- a/src/visualization/write2json.jl +++ b/src/visualization/write2json.jl @@ -17,136 +17,137 @@ setup = :(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrosta """ function trixi2json(solution_callback, integrator) + semi = integrator.p + (; systems) = semi - semi = integrator.p - (; systems) = semi + output_directory = solution_callback.output_directory + prefix = solution_callback.prefix + git_hash = solution_callback.git_hash - output_directory = solution_callback.output_directory - prefix = solution_callback.prefix - git_hash = solution_callback.git_hash + filenames = system_names(systems) - filenames = system_names(systems) + foreach_system(semi) do system + system_index = system_indices(system, semi) - foreach_system(semi) do system - system_index = system_indices(system, semi) - - trixi2json(system; system_name = filenames[system_index], output_directory, prefix, git_hash) - end + trixi2json(system; system_name=filenames[system_index], output_directory, prefix, + git_hash) + end end function trixi2json(system; system_name, output_directory, prefix, git_hash) - mkpath(output_directory) + mkpath(output_directory) - meta_data = Dict{String, Any}( - "solver_version" => git_hash, - "julia_version" => string(VERSION), - ) + meta_data = Dict{String, Any}( + "solver_name" => "TrixiParticles.jl", + "solver_version" => git_hash, + "julia_version" => string(VERSION) + ) - get_meta_data!(meta_data, system) + get_meta_data!(meta_data, system) - # handle "_" on optional prefix strings - add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") + # handle "_" on optional prefix strings + add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") - # Write metadata to JSON-file - json_file = joinpath(output_directory, - add_opt_str_pre(prefix) * "$(system_name)_metadata.json") + # Write metadata to JSON-file + json_file = joinpath(output_directory, + add_opt_str_pre(prefix) * "$(system_name)_metadata.json") - open(json_file, "w") do file - JSON.print(file, meta_data, 2) - end + open(json_file, "w") do file + JSON.print(file, meta_data, 2) + end end function get_meta_data!(meta_data, system::FluidSystem) - meta_data["acceleration"] = system.acceleration - meta_data["viscosity"] = type2string(system.viscosity) - get_meta_data!(meta_data, system.viscosity) - meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) - meta_data["smoothing_length"] = system.smoothing_length - meta_data["density_calculator"] = type2string(system.density_calculator) - - if system isa WeaklyCompressibleSPHSystem - 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["solver"] = "WCSPH" - - 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 - else - 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 - - return meta_data + meta_data["acceleration"] = system.acceleration + meta_data["viscosity"] = type2string(system.viscosity) + get_meta_data!(meta_data, system.viscosity) + meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) + meta_data["smoothing_length"] = system.smoothing_length + meta_data["density_calculator"] = type2string(system.density_calculator) + + if system isa WeaklyCompressibleSPHSystem + 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["solver"] = "WCSPH" + + 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 + else + 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 + + return meta_data end get_meta_data!(meta_data, viscosity::Nothing) = meta_data function get_meta_data!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) - meta_data["viscosity_nu"] = viscosity.nu - meta_data["viscosity_epsilon"] = viscosity.epsilon + meta_data["viscosity_nu"] = viscosity.nu + meta_data["viscosity_epsilon"] = viscosity.epsilon end function get_meta_data!(meta_data, viscosity::ArtificialViscosityMonaghan) - meta_data["viscosity_alpha"] = viscosity.alpha - meta_data["viscosity_beta"] = viscosity.beta - meta_data["viscosity_epsilon"] = viscosity.epsilon + meta_data["viscosity_alpha"] = viscosity.alpha + meta_data["viscosity_beta"] = viscosity.beta + meta_data["viscosity_epsilon"] = viscosity.epsilon end function get_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"] = system.smoothing_length - - get_meta_data!(meta_data, system.boundary_model, system) + 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"] = system.smoothing_length + + get_meta_data!(meta_data, system.boundary_model, system) end function get_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["flow_direction"] = system.flow_direction - meta_data["velocity_function"] = type2string(system.reference_velocity) - meta_data["pressure_function"] = type2string(system.reference_pressure) - meta_data["density_function"] = type2string(system.reference_density) + meta_data["boundary_zone"] = type2string(system.boundary_zone) + meta_data["width"] = round(system.boundary_zone.zone_width, digits=3) + meta_data["flow_direction"] = system.flow_direction + 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 get_meta_data!(meta_data, system::BoundarySPHSystem) - get_meta_data!(meta_data, system.boundary_model, system) + get_meta_data!(meta_data, system.boundary_model, system) end function get_meta_data!(meta_data, model, system) - return meta_data + return meta_data end function get_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 + meta_data["boundary_model"] = "BoundaryModelMonaghanKajtar" + meta_data["boundary_spacing_ratio"] = model.beta + meta_data["boundary_K"] = model.K end function get_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) + 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 get_meta_data!(meta_data, system::BoundaryDEMSystem) - return meta_data + return meta_data end From 8803ef900375a12680102e8962c94de76e252663 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 14 Apr 2025 11:46:22 +0200 Subject: [PATCH 04/18] update `write2json.jl` to match merge --- src/visualization/write2json.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/visualization/write2json.jl b/src/visualization/write2json.jl index dc7fb344ac..8dd6f07e5f 100644 --- a/src/visualization/write2json.jl +++ b/src/visualization/write2json.jl @@ -62,7 +62,7 @@ function get_meta_data!(meta_data, system::FluidSystem) meta_data["viscosity"] = type2string(system.viscosity) get_meta_data!(meta_data, system.viscosity) meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) - meta_data["smoothing_length"] = system.smoothing_length + meta_data["smoothing_length_factor"] = system.cache.smoothing_length_factor meta_data["density_calculator"] = type2string(system.density_calculator) if system isa WeaklyCompressibleSPHSystem @@ -111,7 +111,7 @@ function get_meta_data!(meta_data, system::TotalLagrangianSPHSystem) 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"] = system.smoothing_length + meta_data["smoothing_length_factor"] = initial_smoothing_length(system) / particle_spacing(system, 1) get_meta_data!(meta_data, system.boundary_model, system) end @@ -129,7 +129,7 @@ function get_meta_data!(meta_data, system::BoundarySPHSystem) get_meta_data!(meta_data, system.boundary_model, system) end -function get_meta_data!(meta_data, model, system) +function get_meta_data!(meta_data, model::Nothing, system) return meta_data end From b86e2d64f3e115966e744c0e773fb9c18b33dede Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 19 May 2025 11:13:13 +0200 Subject: [PATCH 05/18] rename file --- src/{visualization/write2json.jl => io/write_json.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{visualization/write2json.jl => io/write_json.jl} (100%) diff --git a/src/visualization/write2json.jl b/src/io/write_json.jl similarity index 100% rename from src/visualization/write2json.jl rename to src/io/write_json.jl From 164b56e4cb1c83f3bdba4fde9102d7133a6f322b Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 19 May 2025 11:22:49 +0200 Subject: [PATCH 06/18] update references --- src/io/io.jl | 1 + src/io/write_vtk.jl | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io/io.jl b/src/io/io.jl index 380151bcff..85c484c8e4 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -1,2 +1,3 @@ include("write_vtk.jl") +include("write_json.jl") include("read_vtk.jl") diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index ad9a9bfcbf..d0886c8616 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -365,8 +365,7 @@ 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) +function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system) vtk["hydrodynamic_density"] = current_density(v, system) vtk["pressure"] = model.pressure From 2165d3747a9683a7cac96e3092c627ea10869330 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 19 May 2025 11:29:48 +0200 Subject: [PATCH 07/18] rename function --- src/io/write_json.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/io/write_json.jl b/src/io/write_json.jl index 8dd6f07e5f..199bb65cd0 100644 --- a/src/io/write_json.jl +++ b/src/io/write_json.jl @@ -43,7 +43,7 @@ function trixi2json(system; system_name, output_directory, prefix, git_hash) "julia_version" => string(VERSION) ) - get_meta_data!(meta_data, system) + write2json!(meta_data, system) # handle "_" on optional prefix strings add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") @@ -57,10 +57,10 @@ function trixi2json(system; system_name, output_directory, prefix, git_hash) end end -function get_meta_data!(meta_data, system::FluidSystem) +function write2json!(meta_data, system::FluidSystem) meta_data["acceleration"] = system.acceleration meta_data["viscosity"] = type2string(system.viscosity) - get_meta_data!(meta_data, system.viscosity) + write2json!(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) @@ -92,20 +92,20 @@ function get_meta_data!(meta_data, system::FluidSystem) return meta_data end -get_meta_data!(meta_data, viscosity::Nothing) = meta_data +write2json!(meta_data, viscosity::Nothing) = meta_data -function get_meta_data!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) +function write2json!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) meta_data["viscosity_nu"] = viscosity.nu meta_data["viscosity_epsilon"] = viscosity.epsilon end -function get_meta_data!(meta_data, viscosity::ArtificialViscosityMonaghan) +function write2json!(meta_data, viscosity::ArtificialViscosityMonaghan) meta_data["viscosity_alpha"] = viscosity.alpha meta_data["viscosity_beta"] = viscosity.beta meta_data["viscosity_epsilon"] = viscosity.epsilon end -function get_meta_data!(meta_data, system::TotalLagrangianSPHSystem) +function write2json!(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 @@ -113,10 +113,10 @@ function get_meta_data!(meta_data, system::TotalLagrangianSPHSystem) meta_data["smoothing_kernel"] = type2string(system.smoothing_kernel) meta_data["smoothing_length_factor"] = initial_smoothing_length(system) / particle_spacing(system, 1) - get_meta_data!(meta_data, system.boundary_model, system) + write2json!(meta_data, system.boundary_model, system) end -function get_meta_data!(meta_data, system::OpenBoundarySPHSystem) +function write2json!(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["flow_direction"] = system.flow_direction @@ -125,21 +125,21 @@ function get_meta_data!(meta_data, system::OpenBoundarySPHSystem) meta_data["density_function"] = type2string(system.reference_density) end -function get_meta_data!(meta_data, system::BoundarySPHSystem) - get_meta_data!(meta_data, system.boundary_model, system) +function write2json!(meta_data, system::BoundarySPHSystem) + write2json!(meta_data, system.boundary_model, system) end -function get_meta_data!(meta_data, model::Nothing, system) +function write2json!(meta_data, model::Nothing, system) return meta_data end -function get_meta_data!(meta_data, model::BoundaryModelMonaghanKajtar, system) +function write2json!(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 get_meta_data!(meta_data, model::BoundaryModelDummyParticles, system) +function write2json!(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 @@ -148,6 +148,6 @@ function get_meta_data!(meta_data, model::BoundaryModelDummyParticles, system) meta_data["viscosity_model"] = type2string(model.viscosity) end -function get_meta_data!(meta_data, system::BoundaryDEMSystem) +function write2json!(meta_data, system::BoundaryDEMSystem) return meta_data end From 51317e3e31d2637483258c68f81b9a248d6ff00a Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 19 May 2025 11:30:24 +0200 Subject: [PATCH 08/18] Refactor write2vtk! --- src/preprocessing/particle_packing/system.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index c9e8c5c296..cc63142c66 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -217,12 +217,13 @@ 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 +end + +function write2json!(vtk, v, u, t, system::ParticlePackingSystem) + meta_data["signed_distances"] = system.signed_distances end # Skip for fixed systems From 8f6eda96fc6f0a043ac47ab1b1ba6660189f6a5a Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 19 May 2025 11:35:39 +0200 Subject: [PATCH 09/18] formatting --- src/callbacks/solution_saving.jl | 2 +- src/io/write_json.jl | 3 ++- src/io/write_vtk.jl | 5 +++-- src/preprocessing/particle_packing/system.jl | 2 +- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index 4a6af94c69..8a978be518 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -154,7 +154,7 @@ end # `affect!` function (solution_callback::SolutionSavingCallback)(integrator) (; interval, output_directory, custom_quantities, git_hash, verbose, - prefix, latest_saved_iter, max_coordinates) = solution_callback + prefix, latest_saved_iter, max_coordinates) = solution_callback vu_ode = integrator.u semi = integrator.p diff --git a/src/io/write_json.jl b/src/io/write_json.jl index 199bb65cd0..4ddab8894b 100644 --- a/src/io/write_json.jl +++ b/src/io/write_json.jl @@ -111,7 +111,8 @@ function write2json!(meta_data, system::TotalLagrangianSPHSystem) 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) + meta_data["smoothing_length_factor"] = initial_smoothing_length(system) / + particle_spacing(system, 1) write2json!(meta_data, system.boundary_model, system) end diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index d0886c8616..0799184dde 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -71,7 +71,8 @@ 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_), max_coordinates=Inf, git_hash=compute_git_hash(), + prefix="", iter=nothing, system_name=vtkname(system_), + max_coordinates=Inf, git_hash=compute_git_hash(), custom_quantities...) mkpath(output_directory) @@ -123,7 +124,7 @@ function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_direc vtk["ndims"] = ndims(system) vtk["particle_spacing"] = [particle_spacing(system, particle) - for particle in active_particles(system)] + for particle in active_particles(system)] # Extract custom quantities for this system for (key, quantity) in custom_quantities diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index cc63142c66..88a1eaafdf 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -223,7 +223,7 @@ function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem) end function write2json!(vtk, v, u, t, system::ParticlePackingSystem) - meta_data["signed_distances"] = system.signed_distances + meta_data["signed_distances"] = system.signed_distances end # Skip for fixed systems From 7825eb9688fc4621b27dc7a3e4c4e809ef1803c9 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Tue, 20 May 2025 09:49:05 +0200 Subject: [PATCH 10/18] refactor functions --- src/TrixiParticles.jl | 3 +- src/callbacks/post_process.jl | 1 + src/callbacks/solution_saving.jl | 5 +- src/io/io.jl | 140 +++++++++++++++++++++++++++- src/io/write_json.jl | 154 ------------------------------- 5 files changed, 143 insertions(+), 160 deletions(-) delete mode 100644 src/io/write_json.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 410dbd076b..e0d3268979 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -83,7 +83,8 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx export HertzContactModel, LinearContactModel export BoundaryMovement export examples_dir, validation_dir -export trixi2vtk, trixi2json, vtk2trixi +export trixi2vtk, vtk2trixi +export write_meta_data export RectangularTank, RectangularShape, SphereShape, ComplexShape export ParticlePackingSystem, SignedDistanceField export WindingNumberHormann, WindingNumberJacobson diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 9370bcdadc..3d2404c0e2 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -216,6 +216,7 @@ 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 8a978be518..deab0a0d13 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -130,10 +130,7 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in solution_callback.latest_saved_iter = -1 solution_callback.git_hash[] = compute_git_hash() - # Write metadata to JSON-file - if integrator.iter == 0 - trixi2json(solution_callback, integrator) - end + write_meta_data(solution_callback, integrator) # Save initial solution if solution_callback.save_initial_solution diff --git a/src/io/io.jl b/src/io/io.jl index 85c484c8e4..57f7daf2c9 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -1,3 +1,141 @@ include("write_vtk.jl") -include("write_json.jl") include("read_vtk.jl") + +function write_meta_data(callback::Union{SolutionSavingCallback, PostprocessCallback}, + integrator) + semi = integrator.p + (; systems) = semi + + output_directory = callback.output_directory + prefix = callback.prefix + git_hash = callback.git_hash + + filenames = system_names(systems) + + foreach_system(semi) do system + system_index = system_indices(system, semi) + + write_meta_data(system; system_name=filenames[system_index], output_directory, + prefix, + git_hash) + end +end + +function write_meta_data(system; system_name, output_directory, prefix, git_hash) + mkpath(output_directory) + + meta_data = Dict{String, Any}( + "solver_name" => "TrixiParticles.jl", + "solver_version" => git_hash, + "julia_version" => string(VERSION) + ) + + get_meta_data!(meta_data, system) + + # handle "_" on optional prefix strings + add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") + + # Write metadata to JSON-file + json_file = joinpath(output_directory, + add_opt_str_pre(prefix) * "$(system_name)_metadata.json") + + open(json_file, "w") do file + JSON.print(file, meta_data, 2) + end +end + +function get_meta_data!(meta_data, system::FluidSystem) + meta_data["acceleration"] = system.acceleration + meta_data["viscosity"] = type2string(system.viscosity) + get_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) + + if system isa WeaklyCompressibleSPHSystem + 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["solver"] = "WCSPH" + + 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 + else + 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 + + return meta_data +end + +get_meta_data!(meta_data, viscosity::Nothing) = meta_data + +function get_meta_data!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) + meta_data["viscosity_nu"] = viscosity.nu + meta_data["viscosity_epsilon"] = viscosity.epsilon +end + +function get_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 get_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) + + get_meta_data!(meta_data, system.boundary_model, system) +end + +function get_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["flow_direction"] = system.flow_direction + 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 get_meta_data!(meta_data, system::BoundarySPHSystem) + get_meta_data!(meta_data, system.boundary_model, system) +end + +function get_meta_data!(meta_data, model::Nothing, system) + return meta_data +end + +function get_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 get_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 get_meta_data!(meta_data, system::BoundaryDEMSystem) + return meta_data +end diff --git a/src/io/write_json.jl b/src/io/write_json.jl deleted file mode 100644 index 4ddab8894b..0000000000 --- a/src/io/write_json.jl +++ /dev/null @@ -1,154 +0,0 @@ -""" - trixi2json(solution_callback, integrator) - -Write simulation metadata to a JSON-file. - -# Arguments -- `solution_callback`: Callback storing metadata and output settings. -- `integrator`: ODE integrator containing simulation data. - -# Example -```jldoctest; output = false, saving_callback = SolutionSavingCallback(dt = 0.1, output_directory = "output", prefix = "solution"), -setup = :(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), tspan = (0.0, 0.01), callbacks = saving_callback)) - -# output - -``` -""" - -function trixi2json(solution_callback, integrator) - semi = integrator.p - (; systems) = semi - - output_directory = solution_callback.output_directory - prefix = solution_callback.prefix - git_hash = solution_callback.git_hash - - filenames = system_names(systems) - - foreach_system(semi) do system - system_index = system_indices(system, semi) - - trixi2json(system; system_name=filenames[system_index], output_directory, prefix, - git_hash) - end -end - -function trixi2json(system; system_name, output_directory, prefix, git_hash) - mkpath(output_directory) - - meta_data = Dict{String, Any}( - "solver_name" => "TrixiParticles.jl", - "solver_version" => git_hash, - "julia_version" => string(VERSION) - ) - - write2json!(meta_data, system) - - # handle "_" on optional prefix strings - add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") - - # Write metadata to JSON-file - json_file = joinpath(output_directory, - add_opt_str_pre(prefix) * "$(system_name)_metadata.json") - - open(json_file, "w") do file - JSON.print(file, meta_data, 2) - end -end - -function write2json!(meta_data, system::FluidSystem) - meta_data["acceleration"] = system.acceleration - meta_data["viscosity"] = type2string(system.viscosity) - write2json!(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) - - if system isa WeaklyCompressibleSPHSystem - 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["solver"] = "WCSPH" - - 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 - else - 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 - - return meta_data -end - -write2json!(meta_data, viscosity::Nothing) = meta_data - -function write2json!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) - meta_data["viscosity_nu"] = viscosity.nu - meta_data["viscosity_epsilon"] = viscosity.epsilon -end - -function write2json!(meta_data, viscosity::ArtificialViscosityMonaghan) - meta_data["viscosity_alpha"] = viscosity.alpha - meta_data["viscosity_beta"] = viscosity.beta - meta_data["viscosity_epsilon"] = viscosity.epsilon -end - -function write2json!(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) - - write2json!(meta_data, system.boundary_model, system) -end - -function write2json!(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["flow_direction"] = system.flow_direction - 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 write2json!(meta_data, system::BoundarySPHSystem) - write2json!(meta_data, system.boundary_model, system) -end - -function write2json!(meta_data, model::Nothing, system) - return meta_data -end - -function write2json!(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 write2json!(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 write2json!(meta_data, system::BoundaryDEMSystem) - return meta_data -end From 251eaf290c5f4f61d7c903d1bed94988fc1b6a6c Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 2 Jun 2025 13:46:43 +0200 Subject: [PATCH 11/18] implement suggestions --- src/TrixiParticles.jl | 1 - src/callbacks/post_process.jl | 1 + src/io/io.jl | 30 ++++++++++---------- src/io/write_vtk.jl | 2 +- src/preprocessing/particle_packing/system.jl | 5 +--- 5 files changed, 18 insertions(+), 21 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index e0d3268979..bab7e7a813 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -84,7 +84,6 @@ export HertzContactModel, LinearContactModel export BoundaryMovement export examples_dir, validation_dir export trixi2vtk, vtk2trixi -export write_meta_data export RectangularTank, RectangularShape, SphereShape, ComplexShape export ParticlePackingSystem, SignedDistanceField export WindingNumberHormann, WindingNumberJacobson diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 3d2404c0e2..7fdce08bfe 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -217,6 +217,7 @@ function initialize_postprocess_callback!(cb::PostprocessCallback, u, t, integra cb(integrator) write_meta_data(cb, integrator) + return cb end diff --git a/src/io/io.jl b/src/io/io.jl index 57f7daf2c9..71210e72fc 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -30,7 +30,7 @@ function write_meta_data(system; system_name, output_directory, prefix, git_hash "julia_version" => string(VERSION) ) - get_meta_data!(meta_data, system) + add_meta_data!(meta_data, system) # handle "_" on optional prefix strings add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") @@ -44,10 +44,10 @@ function write_meta_data(system; system_name, output_directory, prefix, git_hash end end -function get_meta_data!(meta_data, system::FluidSystem) +function add_meta_data!(meta_data, system::FluidSystem) meta_data["acceleration"] = system.acceleration meta_data["viscosity"] = type2string(system.viscosity) - get_meta_data!(meta_data, 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) @@ -79,20 +79,20 @@ function get_meta_data!(meta_data, system::FluidSystem) return meta_data end -get_meta_data!(meta_data, viscosity::Nothing) = meta_data +add_meta_data!(meta_data, viscosity::Nothing) = meta_data -function get_meta_data!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) +function add_meta_data!(meta_data, viscosity::Union{ViscosityAdami, ViscosityMorris}) meta_data["viscosity_nu"] = viscosity.nu meta_data["viscosity_epsilon"] = viscosity.epsilon end -function get_meta_data!(meta_data, viscosity::ArtificialViscosityMonaghan) +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 get_meta_data!(meta_data, system::TotalLagrangianSPHSystem) +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 @@ -101,10 +101,10 @@ function get_meta_data!(meta_data, system::TotalLagrangianSPHSystem) meta_data["smoothing_length_factor"] = initial_smoothing_length(system) / particle_spacing(system, 1) - get_meta_data!(meta_data, system.boundary_model, system) + add_meta_data!(meta_data, system.boundary_model, system) end -function get_meta_data!(meta_data, system::OpenBoundarySPHSystem) +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["flow_direction"] = system.flow_direction @@ -113,21 +113,21 @@ function get_meta_data!(meta_data, system::OpenBoundarySPHSystem) meta_data["density_function"] = type2string(system.reference_density) end -function get_meta_data!(meta_data, system::BoundarySPHSystem) - get_meta_data!(meta_data, system.boundary_model, system) +function add_meta_data!(meta_data, system::BoundarySPHSystem) + add_meta_data!(meta_data, system.boundary_model, system) end -function get_meta_data!(meta_data, model::Nothing, system) +function add_meta_data!(meta_data, model::Nothing, system) return meta_data end -function get_meta_data!(meta_data, model::BoundaryModelMonaghanKajtar, system) +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 get_meta_data!(meta_data, model::BoundaryModelDummyParticles, system) +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 @@ -136,6 +136,6 @@ function get_meta_data!(meta_data, model::BoundaryModelDummyParticles, system) meta_data["viscosity_model"] = type2string(model.viscosity) end -function get_meta_data!(meta_data, system::BoundaryDEMSystem) +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 0799184dde..f48b46db15 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -115,7 +115,7 @@ function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_direc end @trixi_timeit timer() "write to vtk" vtk_grid(file, points, cells) do vtk - # dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem + # Dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem write2vtk!(vtk, v, u, t, system) # Store particle index diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index bf52a78427..311ed6554e 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -226,10 +226,7 @@ update_callback_used!(system::ParticlePackingSystem) = system.update_callback_us function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem) vtk["velocity"] = [advection_velocity(v, system, particle) for particle in active_particles(system)] -end - -function write2json!(vtk, v, u, t, system::ParticlePackingSystem) - meta_data["signed_distances"] = system.signed_distances + vtk["signed_distances"] = system.signed_distances end # Skip for fixed systems From fecc01fc9fe76125e52e5c17dcfed6b5075cd8d8 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 2 Jun 2025 15:16:47 +0200 Subject: [PATCH 12/18] fix bugs --- src/io/io.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/io/io.jl b/src/io/io.jl index 71210e72fc..051c0d641e 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -7,7 +7,7 @@ function write_meta_data(callback::Union{SolutionSavingCallback, PostprocessCall (; systems) = semi output_directory = callback.output_directory - prefix = callback.prefix + prefix = hasproperty(callback, :prefix) ? callback.prefix : "" git_hash = callback.git_hash filenames = system_names(systems) @@ -44,6 +44,14 @@ function write_meta_data(system; system_name, output_directory, prefix, git_hash 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::FluidSystem) meta_data["acceleration"] = system.acceleration meta_data["viscosity"] = type2string(system.viscosity) @@ -107,7 +115,6 @@ 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["flow_direction"] = system.flow_direction meta_data["velocity_function"] = type2string(system.reference_velocity) meta_data["pressure_function"] = type2string(system.reference_pressure) meta_data["density_function"] = type2string(system.reference_density) From 3a399fff6315b63fb9b276d52b8df5a94a56124b Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Tue, 3 Jun 2025 11:01:33 +0200 Subject: [PATCH 13/18] add add_meta_data! function for `system::ParticlePackingSystem` --- src/preprocessing/particle_packing/particle_packing.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/preprocessing/particle_packing/particle_packing.jl b/src/preprocessing/particle_packing/particle_packing.jl index c0ffd2b8d0..1ea1856d87 100644 --- a/src/preprocessing/particle_packing/particle_packing.jl +++ b/src/preprocessing/particle_packing/particle_packing.jl @@ -2,3 +2,7 @@ include("nhs_faces.jl") include("signed_distance.jl") include("system.jl") include("rhs.jl") + +function add_meta_data!(meta_data, system::ParticlePackingSystem) + return meta_data +end From 9868c4aff5556397d133d223df0ff7865371f7ca Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Tue, 3 Jun 2025 11:17:10 +0200 Subject: [PATCH 14/18] relocate `add_meta_data!` to `system.jl` --- src/preprocessing/particle_packing/particle_packing.jl | 4 ---- src/preprocessing/particle_packing/system.jl | 4 ++++ 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/preprocessing/particle_packing/particle_packing.jl b/src/preprocessing/particle_packing/particle_packing.jl index 1ea1856d87..c0ffd2b8d0 100644 --- a/src/preprocessing/particle_packing/particle_packing.jl +++ b/src/preprocessing/particle_packing/particle_packing.jl @@ -2,7 +2,3 @@ include("nhs_faces.jl") include("signed_distance.jl") include("system.jl") include("rhs.jl") - -function add_meta_data!(meta_data, system::ParticlePackingSystem) - return meta_data -end diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index 311ed6554e..8b713bd407 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -229,6 +229,10 @@ function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem) vtk["signed_distances"] = system.signed_distances end +function add_meta_data!(meta_data, system::ParticlePackingSystem) + return meta_data +end + # Skip for fixed systems write_u0!(u0, system::ParticlePackingSystem{<:Any, true}) = u0 From 45940d97c0a75a29812939ff559828cc62002f22 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Wed, 4 Jun 2025 09:37:22 +0200 Subject: [PATCH 15/18] formatting --- src/io/io.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/io/io.jl b/src/io/io.jl index 051c0d641e..2511c8e825 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -16,8 +16,7 @@ function write_meta_data(callback::Union{SolutionSavingCallback, PostprocessCall system_index = system_indices(system, semi) write_meta_data(system; system_name=filenames[system_index], output_directory, - prefix, - git_hash) + prefix, git_hash) end end From 3ed23c5851e3d6b53a295ca6a4c1f0b7f25a5e40 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 13 Jun 2025 15:26:48 +0200 Subject: [PATCH 16/18] Add simple SGS to Adami and Morris Viscosity (#753) * add sgs * update * fix * remove * fix doc * fix * fixes * typo * tpo * format * add explanation * adjust tests * add to validation setup for wcsph * review fixes * Update taylor_green_vortex_2d.jl * Update validation_taylor_green_vortex_2d.jl * Update dam_break_2d.jl * fix and add citations * add news * implement suggestions * format * rename viscosity * too much renaming * implement suggestions * Update validation_taylor_green_vortex_2d.jl * fix tests again --------- Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- NEWS.md | 2 + docs/src/refs.bib | 24 +- examples/fluid/dam_break_2d.jl | 28 +- examples/fluid/dam_break_2phase_2d.jl | 2 +- examples/fluid/dam_break_oil_film_2d.jl | 3 +- examples/fluid/hydrostatic_water_column_2d.jl | 4 +- examples/fluid/taylor_green_vortex_2d.jl | 2 +- src/TrixiParticles.jl | 3 +- src/io/write_vtk.jl | 4 +- src/schemes/fluid/viscosity.jl | 254 ++++++++++++++++-- test/examples/examples_fluid.jl | 12 +- test/examples/gpu.jl | 10 +- test/schemes/fluid/viscosity.jl | 83 +++++- test/validation/validation.jl | 2 +- .../validation_taylor_green_vortex_2d.jl | 4 + 15 files changed, 381 insertions(+), 56 deletions(-) diff --git a/NEWS.md b/NEWS.md index f004f1410e..fc9b832e30 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 10c1b2ef8c..b857c7a128 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 35493a8c2b..950f8c8060 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 de99a5f397..247462b52b 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 320f185092..f301f99644 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/hydrostatic_water_column_2d.jl b/examples/fluid/hydrostatic_water_column_2d.jl index b492290756..f60a76f965 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 c1f0f72fd4..75e6387184 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/src/TrixiParticles.jl b/src/TrixiParticles.jl index bab7e7a813..ff23f94902 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/io/write_vtk.jl b/src/io/write_vtk.jl index 0e5aaf4c15..6a9663c68c 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -357,7 +357,9 @@ 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 diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index e39e7205d0..ddd74a7481 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 4c557ee465..141d602787 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 9d551a5273..85d657b678 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 f11f4b05b7..852ebfbc34 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 d1c8420842..70aceb4bf4 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 73884e38c2..4543412306 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 From 7d34fe52aa0564d864275a0a65d60268b2d86524 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 16 Jun 2025 14:46:55 +0200 Subject: [PATCH 17/18] dispatch `FluidSystem` --- src/io/io.jl | 56 ++++++++++++++++++++++++++++++---------------------- 1 file changed, 32 insertions(+), 24 deletions(-) diff --git a/src/io/io.jl b/src/io/io.jl index 2511c8e825..a84ba1b4be 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -51,7 +51,8 @@ function add_meta_data!(meta_data, system::DEMSystem) return meta_data end -function add_meta_data!(meta_data, system::FluidSystem) +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) @@ -59,31 +60,38 @@ function add_meta_data!(meta_data, system::FluidSystem) meta_data["smoothing_length_factor"] = system.cache.smoothing_length_factor meta_data["density_calculator"] = type2string(system.density_calculator) - if system isa WeaklyCompressibleSPHSystem - 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["solver"] = "WCSPH" - - 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 - else - 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 + # 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 - return meta_data +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 From ff5651fb533afe57d3e6edf3c20615e2696151c0 Mon Sep 17 00:00:00 2001 From: Marcel Schurer Date: Mon, 16 Jun 2025 14:48:18 +0200 Subject: [PATCH 18/18] Write one unified metadata file for all systems --- src/io/io.jl | 50 ++++++++++++++++++++++++++------------------------ 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/src/io/io.jl b/src/io/io.jl index a84ba1b4be..4a0f7e51c6 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -3,43 +3,45 @@ include("read_vtk.jl") function write_meta_data(callback::Union{SolutionSavingCallback, PostprocessCallback}, integrator) - semi = integrator.p - (; systems) = semi + # handle "_" on optional prefix strings + add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") - output_directory = callback.output_directory - prefix = hasproperty(callback, :prefix) ? callback.prefix : "" git_hash = callback.git_hash + prefix = hasproperty(callback, :prefix) ? callback.prefix : "" - filenames = system_names(systems) + semi = integrator.p + names = system_names(semi.systems) + # fill `systems` with metadata for each system + systems = Dict{String, Any}() foreach_system(semi) do system - system_index = system_indices(system, semi) + idx = system_indices(system, semi) + name = add_opt_str_pre(prefix) * "$(names[idx])" - write_meta_data(system; system_name=filenames[system_index], output_directory, - prefix, git_hash) - end -end + meta_data = Dict{String, Any}() + add_meta_data!(meta_data, system) -function write_meta_data(system; system_name, output_directory, prefix, git_hash) - mkpath(output_directory) + systems[name] = meta_data + end - meta_data = Dict{String, Any}( - "solver_name" => "TrixiParticles.jl", - "solver_version" => git_hash, - "julia_version" => string(VERSION) + # 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 ) - add_meta_data!(meta_data, system) - - # handle "_" on optional prefix strings - add_opt_str_pre(str) = (str === "" ? "" : "$(str)_") + # write JSON-file + output_directory = callback.output_directory + mkpath(output_directory) - # Write metadata to JSON-file - json_file = joinpath(output_directory, - add_opt_str_pre(prefix) * "$(system_name)_metadata.json") + json_file = joinpath(output_directory, "simulation_metadata.json") open(json_file, "w") do file - JSON.print(file, meta_data, 2) + JSON.print(file, simulation_meta_data, 2) end end