Paraview catalyst

docs/src/visualization.md

@@ -399,6 +399,48 @@ below:
   <style>.embed-container { position: relative; padding-bottom: 56.25%; height: 0; overflow: hidden; max-width: 100%; } .embed-container iframe, .embed-container object, .embed-container embed { position: absolute; top: 0; left: 0; width: 100%; height: 100%; }</style><div class='embed-container'><iframe src='https://www.youtube-nocookie.com/embed/UZtrqeDY1Fs' frameborder='0' allowfullscreen></iframe></div>
 ```
 
+### Visualizing results during a simulation using Paraview Catalyst
+With a Paraview Installation containing the Catalyst Library, you can use the ParaViewCatalystCallback to visualise your results in Paraview during the simulation.
+
+```julia
+    ParaViewCatalystCallback(; interval=0, nvisnodes = nothing, catalyst_pipeline = nothing)
+```
++ `interval` determines the amount of timesteps between calls of this callback
++ `nvisnodes` determines the number of visualization nodes. 
+  + Paraview can not handle the Polynomials for each grid cell, so an interpolation is necessary
+  + Visualization nodes are the nodes per dimension in each cell
+    + For example in a 3D Plot each cell has a nxnxn array of visualization nodes
++ `catalyst_pipeline` a path to the catalyst pipeline if the default should not be used
+  + The starting view shown by Paraview is determined by a python pipeline
+  + ParaViewCatalyst.jl includes a standard `catalyst_pipeline.py`
+  + In principle your current view in Paraview can be exported as a pipeline using File->Save Catalyst State...
+  but in practice the resulting pipelines did not work in our testing.
+  ```python
+  def catalyst_execute(info):
+    global input
+    global options
+    input.UpdatePipeline()
+    contour1.UpdatePipeline()
+    SaveExtractsUsingCatalystOptions(options)
+  ```
+  a block like this needed to replace the bottom of the exported pipeline, to get it working for us
+
+
+To be able to use the ParaViewCatalystCallback you need
++ A Paraview Installation containing the Catalyst Library
++ The ParaViewCatalyst.jl Julia package
++ The PARAVIEW_CATALYST_PATH Environment Variable set to the location of the Catalyst Library contained in Paraview
+
+To use the ParaViewCatalystCallback you have to
++ Include ParaViewCatalyst.jl in your julia Script via `using ParaViewCatalyst`
++ Create a ParaViewCatalystCallback and add it to your Callbackset
++ Before running the script:
+  + Open Paraview and press Catalyst->Connect... and use the standard Port of 22222 (if not specified otherwise in the pipeline)
+  + Press Catalyst->Pause Simulation
++ Now start running your Simulation
++ once the first set of values has been calculated and the Simulation pauses
+  + Adjust the view in Paraview to your liking
+  + Press Catalyst->Continue
 
 ## Trixi2Vtk
 Trixi2Vtk converts Trixi.jl's `.h5` output files to VTK files, which can be read

examples/p4est_2d_dgsem/elixir_advection_basic_catalyst.jl

@@ -0,0 +1,66 @@
+# The same setup as tree_2d_dgsem/elixir_advection_basic.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEq
+using Trixi
+using ParaViewCatalyst
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (1, 1)
+
+# Create P4estMesh with 8 x 8 trees and 16 x 16 elements
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 initial_refinement_level = 1)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 1.0));
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+catalyst_callback = ParaViewCatalystCallback(interval=100)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback, catalyst_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

examples/p4est_3d_dgsem/elixir_advection_cubed_sphere_catalyst.jl

@@ -0,0 +1,64 @@
+
+using OrdinaryDiffEq
+using Trixi
+using ParaViewCatalyst
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (20, -70, 50)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+initial_condition = initial_condition_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:inside => boundary_condition,
+                           :outside => boundary_condition)
+
+mesh = Trixi.P4estMeshCubedSphere(5, 3, 0.5, 0.5,
+                                  polydeg = 3, initial_refinement_level = 0)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+catalyst_callback = ParaViewCatalystCallback(interval=100, interpolation=true)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback, catalyst_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

examples/p4est_3d_dgsem/elixir_euler_circular_wind_nonconforming_catalyst.jl

@@ -0,0 +1,147 @@
+# A manufactured solution of a circular wind with constant angular velocity
+# on a planetary-sized non-conforming mesh
+#
+# Note that this elixir can take several hours to run.
+# Using 24 threads of an AMD Ryzen Threadripper 3990X (more threads don't speed it up further)
+# and `check-bounds=no`, this elixirs takes about 20 minutes to run.
+
+using OrdinaryDiffEq
+using Trixi
+using LinearAlgebra
+using ParaViewCatalyst
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations3D(gamma)
+
+function initial_condition_circular_wind(x, t, equations::CompressibleEulerEquations3D)
+    radius_earth = 6.371229e6
+    lambda, phi, r = cart_to_sphere(x)
+    rel_height = (r - radius_earth) / 30000
+    p = 1e5
+    rho = 1.0 + 0.1  * exp(-50 * ((lambda-1.0)^2/2.0 + (phi-0.4)^2)) +
+                0.08 * exp(-100 * ((lambda-0.8)^2/4.0 + (phi-0.5)^2))
+    v1 = -10 * x[2] / radius_earth
+    v2 = 10 * x[1] / radius_earth
+    v3 = 0.0
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+@inline function source_terms_circular_wind(u, x, t,
+                                            equations::CompressibleEulerEquations3D)
+    radius_earth = 6.371229e6
+    rho = u[1]
+
+    du1 = 0.0
+    du2 = -rho * (10 / radius_earth) * (10 * x[1] / radius_earth)
+    du3 = -rho * (10 / radius_earth) * (10 * x[2] / radius_earth)
+    du4 = 0.0
+    du5 = 0.0
+
+    return SVector(du1, du2, du3, du4, du5)
+end
+
+function cart_to_sphere(x)
+    r = norm(x)
+    lambda = atan(x[2], x[1])
+    if lambda < 0
+        lambda += 2 * pi
+    end
+    phi = asin(x[3] / r)
+
+    return lambda, phi, r
+end
+
+initial_condition = initial_condition_circular_wind
+
+boundary_conditions = Dict(:inside => boundary_condition_slip_wall,
+                           :outside => boundary_condition_slip_wall)
+
+# The speed of sound in this example is 374 m/s.
+surface_flux = FluxLMARS(374)
+# Note that a free stream is not preserved if N < 2 * N_geo, where N is the
+# polydeg of the solver and N_geo is the polydeg of the mesh.
+# However, the FSP error is negligible in this example.
+solver = DGSEM(polydeg = 4, surface_flux = surface_flux)
+
+# Other mesh configurations to run this simulation on.
+# The cylinder allows to use a structured mesh, the face of the cubed sphere
+# can be used for performance reasons.
+
+# # Cylinder
+# function mapping(xi, eta, zeta)
+#   radius_earth = 6.371229e6
+
+#   x = cos((xi + 1) * pi) * (radius_earth + (0.5 * zeta + 0.5) * 30000)
+#   y = sin((xi + 1) * pi) * (radius_earth + (0.5 * zeta + 0.5) * 30000)
+#   z = eta * 1000000
+
+#   return SVector(x, y, z)
+# end
+
+# # One face of the cubed sphere
+# mapping(xi, eta, zeta) = Trixi.cubed_sphere_mapping(xi, eta, zeta, 6.371229e6, 30000.0, 1)
+
+# trees_per_dimension = (8, 8, 4)
+# mesh = P4estMesh(trees_per_dimension, polydeg=3,
+#                  mapping=mapping,
+#                  initial_refinement_level=0,
+#                  periodicity=false)
+
+trees_per_cube_face = (2, 1)
+mesh = Trixi.P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000.0,
+                                  polydeg = 4, initial_refinement_level = 1)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_circular_wind,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 50 * 24 * 60 * 60.0) # time in seconds for 10 days
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 5000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 2000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 1.004,
+                                      max_level = 1, max_threshold = 1.11)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 2000,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+catalyst_callback = ParaViewCatalystCallback(interval=100)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        save_solution,
+                        catalyst_callback)
+
+###############################################################################
+# run the simulation
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+# Enable threading of the RK method for better performance on multiple threads
+sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1.0e-6,
+            reltol = 1.0e-6,
+            ode_default_options()..., callback = callbacks, maxiters=1e7);
+
+summary_callback() # print the timer summary