Write metadata to a dedicated JSON file

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).

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},

examples/fluid/dam_break_2d.jl

@@ -47,18 +47,29 @@ 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)
 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)
 

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,

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)
 

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)
 

examples/fluid/hydrostatic_water_column_2d.jl

@@ -32,15 +32,15 @@ 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()
 
 # This is to set acceleration with `trixi_include`
 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)
 

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

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)

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)
 

examples/n_body/n_body_system.jl

@@ -105,7 +105,7 @@ end
 
 TrixiParticles.vtkname(system::NBodySystem) = "n-body"
 
-function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem; write_meta_data=true)
+function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem)
     (; mass) = system
 
     vtk["velocity"] = v

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

src/callbacks/post_process.jl

@@ -216,6 +216,8 @@ function initialize_postprocess_callback!(cb::PostprocessCallback, u, t, integra
     # Apply the callback
     cb(integrator)
 
+    write_meta_data(cb, integrator)
+
     return cb
 end
 

src/callbacks/solution_saving.jl

@@ -2,8 +2,7 @@
     SolutionSavingCallback(; interval::Integer=0, dt=0.0, save_times=Array{Float64, 1}([]),
                            save_initial_solution=true, save_final_solution=true,
                            output_directory="out", append_timestamp=false, prefix="",
-                           verbose=false, write_meta_data=true, max_coordinates=2^15,
-                           custom_quantities...)
+                           verbose=false, max_coordinates=2^15, custom_quantities...)
 
 
 Callback to save the current numerical solution in VTK format in regular intervals.
@@ -28,7 +27,6 @@ To ignore a custom quantity for a specific system, return `nothing`.
 - `append_timestamp=false`:     Append current timestamp to the output directory.
 - 'prefix=""':                  Prefix added to the filename.
 - `custom_quantities...`:       Additional user-defined quantities.
-- `write_meta_data=true`:       Write meta data.
 - `verbose=false`:              Print to standard IO when a file is written.
 - `max_coordinates=2^15`:       The coordinates of particles will be clipped if their
                                 absolute values exceed this threshold.
@@ -70,7 +68,6 @@ mutable struct SolutionSavingCallback{I, CQ}
     save_times            :: Vector{Float64}
     save_initial_solution :: Bool
     save_final_solution   :: Bool
-    write_meta_data       :: Bool
     verbose               :: Bool
     output_directory      :: String
     prefix                :: String
@@ -84,8 +81,8 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
                                 save_times=Float64[],
                                 save_initial_solution=true, save_final_solution=true,
                                 output_directory="out", append_timestamp=false,
-                                prefix="", verbose=false, write_meta_data=true,
-                                max_coordinates=Float64(2^15), custom_quantities...)
+                                prefix="", verbose=false, max_coordinates=Float64(2^15),
+                                custom_quantities...)
     if (dt > 0 && interval > 0) || (length(save_times) > 0 && (dt > 0 || interval > 0))
         throw(ArgumentError("Setting multiple save times for the same solution " *
                             "callback is not possible. Use either `dt`, `interval` or `save_times`."))
@@ -101,8 +98,8 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
 
     solution_callback = SolutionSavingCallback(interval, Float64.(save_times),
                                                save_initial_solution, save_final_solution,
-                                               write_meta_data, verbose, output_directory,
-                                               prefix, max_coordinates, custom_quantities,
+                                               verbose, output_directory, prefix,
+                                               max_coordinates, custom_quantities,
                                                -1, Ref("UnknownVersion"))
 
     if length(save_times) > 0
@@ -133,6 +130,8 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in
     solution_callback.latest_saved_iter = -1
     solution_callback.git_hash[] = compute_git_hash()
 
+    write_meta_data(solution_callback, integrator)
+
     # Save initial solution
     if solution_callback.save_initial_solution
         solution_callback(integrator)
@@ -151,8 +150,8 @@ end
 
 # `affect!`
 function (solution_callback::SolutionSavingCallback)(integrator)
-    (; interval, output_directory, custom_quantities, write_meta_data, git_hash,
-     verbose, prefix, latest_saved_iter, max_coordinates) = solution_callback
+    (; interval, output_directory, custom_quantities, git_hash, verbose,
+     prefix, latest_saved_iter, max_coordinates) = solution_callback
 
     vu_ode = integrator.u
     semi = integrator.p
@@ -175,8 +174,8 @@ function (solution_callback::SolutionSavingCallback)(integrator)
 
     @trixi_timeit timer() "save solution" trixi2vtk(vu_ode, semi, integrator.t;
                                                     iter, output_directory, prefix,
-                                                    write_meta_data, git_hash=git_hash[],
-                                                    max_coordinates, custom_quantities...)
+                                                    git_hash=git_hash[], max_coordinates,
+                                                    custom_quantities...)
 
     # Tell OrdinaryDiffEq that `u` has not been modified
     u_modified!(integrator, false)