Write meta data to vtk files (#199)

* write some meta data

* add type2string function additional data

* write the pvd file

* format

* merge error

* add additional output for viscosity models

* make meta data optional

* fix

* format

* Ignore markdown file changes

* [skip ci] review

* [skip ci] unpack all models

* [skip ci] add additional meta data

---------

Co-authored-by: Erik Faulhaber <[email protected]>

.github/workflows/ci.yml

@@ -4,7 +4,25 @@ on:
   push:
     branches:
       - main
+    paths-ignore:
+      - 'AUTHORS.md'
+      - 'CITATION.bib'
+      - 'CONTRIBUTING.md'
+      - 'LICENSE.md'
+      - 'NEWS.md'
+      - 'README.md'
+      - '.github/workflows/CompatHelper.yml'
+      - 'docs/**'
   pull_request:
+      paths-ignore:
+      - 'AUTHORS.md'
+      - 'CITATION.bib'
+      - 'CONTRIBUTING.md'
+      - 'LICENSE.md'
+      - 'NEWS.md'
+      - 'README.md'
+      - '.github/workflows/CompatHelper.yml'
+      - 'docs/**'
   workflow_dispatch:
 
 # Cancel redundant CI tests automatically

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)
+function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem; write_meta_data=true)
     (; mass) = system
 
     vtk["velocity"] = v

src/TrixiParticles.jl

@@ -15,7 +15,7 @@ using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
 using ThreadingUtilities
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
-using WriteVTK: vtk_grid, MeshCell, VTKCellTypes
+using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
 using ForwardDiff
 
 # util needs to be first because of macro @trixi_timeit

src/util.jl

@@ -316,3 +316,7 @@ macro autoinfiltrate(condition=true)
                 lnn,
                 esc(condition))
 end
+
+function type2string(type)
+    return string(nameof(typeof(type)))
+end

src/visualization/write2vtk.jl

@@ -47,7 +47,8 @@ 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), custom_quantities...)
+                   iter=nothing, system_name=vtkname(system), write_meta_data=true,
+                   custom_quantities...)
     mkpath(output_directory)
 
     # handle "_" on optional pre/postfix strings
@@ -58,14 +59,20 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix
                     add_opt_str_pre(prefix) * "$system_name"
                     * add_opt_str_post(iter))
 
+    collection_file = joinpath(output_directory,
+                               add_opt_str_pre(prefix) * "$system_name")
+
+    pvd = paraview_collection(collection_file; append=true)
+
     points = periodic_coords(current_coordinates(u, system), periodic_box)
     cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)]
 
     vtk_grid(file, points, cells) do vtk
-        write2vtk!(vtk, v, u, t, system)
+        write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data)
 
         # Store particle index
         vtk["index"] = eachparticle(system)
+        vtk["time"] = t
 
         # Extract custom quantities for this system
         for (key, quantity) in custom_quantities
@@ -74,7 +81,11 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix
                 vtk[string(key)] = value
             end
         end
+
+        # Add to collection
+        pvd[t] = vtk
     end
+    vtk_save(pvd)
 end
 
 function custom_quantity(quantity::AbstractArray, v, u, t, system)
@@ -128,52 +139,124 @@ vtkname(system::FluidSystem) = "fluid"
 vtkname(system::TotalLagrangianSPHSystem) = "solid"
 vtkname(system::BoundarySPHSystem) = "boundary"
 
-function write2vtk!(vtk, v, u, t, system::FluidSystem)
+function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true)
     vtk["velocity"] = view(v, 1:ndims(system), :)
     vtk["density"] = [particle_density(v, system, particle)
                       for particle in eachparticle(system)]
     vtk["pressure"] = [particle_pressure(v, system, particle)
                        for particle in eachparticle(system)]
 
+    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
+            vtk["state_equation"] = type2string(system.state_equation)
+            vtk["state_equation_rho0"] = system.state_equation.reference_density
+            vtk["state_equation_p0"] = system.state_equation.reference_pressure
+            vtk["state_equation_pa"] = system.state_equation.background_pressure
+            vtk["state_equation_c"] = system.state_equation.sound_speed
+            if system.state_equation isa StateEquationCole
+                vtk["state_equation_gamma"] = system.state_equation.gamma
+            end
+
+            vtk["solver"] = "WCSPH"
+        else
+            vtk["solver"] = "EDAC"
+            vtk["sound_speed"] = system.sound_speed
+        end
+    end
+
     return vtk
 end
 
-function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem)
+function write2vtk!(vtk, viscosity::ViscosityAdami)
+    vtk["viscosity_nu"] = viscosity.nu
+    vtk["viscosity_epsilon"] = viscosity.epsilon
+end
+
+function write2vtk!(vtk, viscosity::ArtificialViscosityMonaghan)
+    vtk["viscosity_alpha"] = viscosity.alpha
+    vtk["viscosity_beta"] = viscosity.beta
+    vtk["viscosity_epsilon"] = viscosity.epsilon
+end
+
+function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_data=true)
     n_fixed_particles = nparticles(system) - n_moving_particles(system)
 
     vtk["velocity"] = hcat(view(v, 1:ndims(system), :),
                            zeros(ndims(system), n_fixed_particles))
     vtk["material_density"] = system.material_density
+    vtk["young_modulus"] = system.young_modulus
+    vtk["poisson_ratio"] = system.poisson_ratio
+    vtk["lame_lambda"] = system.lame_lambda
+    vtk["lame_mu"] = system.lame_mu
 
-    write2vtk!(vtk, v, u, t, system.boundary_model, system)
+    write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data)
 end
 
-function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem)
-    write2vtk!(vtk, v, u, t, system.boundary_model, system)
+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)
 end
 
-function write2vtk!(vtk, v, u, t, model, system)
+function write2vtk!(vtk, v, u, t, model, system; write_meta_data=true)
     return vtk
 end
 
-function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system)
-    write2vtk!(vtk, v, u, t, model, model.viscosity, system)
+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
 end
 
-function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, viscosity, system)
+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(system.boundary_model.smoothing_kernel)
+        vtk["smoothing_length"] = system.boundary_model.smoothing_length
+        vtk["density_calculator"] = type2string(system.boundary_model.density_calculator)
+        vtk["state_equation"] = type2string(system.boundary_model.state_equation)
+    end
+
+    write2vtk!(vtk, v, u, t, model, model.viscosity, system,
+               write_meta_data=write_meta_data)
+end
+
+function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, viscosity, system;
+                    write_meta_data=true)
     vtk["hydrodynamic_density"] = [particle_density(v, system, particle)
                                    for particle in eachparticle(system)]
     vtk["pressure"] = model.pressure
 
+    if write_meta_data
+        vtk["viscosity_model"] = type2string(viscosity)
+    end
+
     return vtk
 end
 
 function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles,
-                    viscosity::ViscosityAdami, system)
+                    viscosity::ViscosityAdami, system; write_meta_data=true)
     vtk["hydrodynamic_density"] = [particle_density(v, system, particle)
                                    for particle in eachparticle(system)]
     vtk["pressure"] = model.pressure
     vtk["wall_velocity"] = view(model.cache.wall_velocity, 1:ndims(system), :)
 
+    if write_meta_data
+        vtk["viscosity_model"] = "ViscosityAdami"
+    end
+
     return vtk
 end