Merge branch 'main' into validate-open-boundaries

Author: LasNikas

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

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"]
-version = "0.3.1-dev"
+version = "0.3.2-dev"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -47,7 +47,7 @@ DelimitedFiles = "1"
 DiffEqCallbacks = "4"
 FastPow = "0.1"
 FileIO = "1"
-ForwardDiff = "0.10"
+ForwardDiff = "1"
 GPUArraysCore = "0.2"
 JSON = "0.21"
 KernelAbstractions = "0.9"

docs/src/preprocessing/preprocessing.md

@@ -250,6 +250,51 @@ Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_hormann.jl")
 Modules = [TrixiParticles]
 Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_jacobson.jl")]
 ```
+# [Read geometries from file](@id read_geometries_from_file)
+Geometries can be imported using the [`load_geometry`](@ref) function.
+- For 3D geometries, we support the binary (`.stl`) format.
+- For 2D geometries, the recommended format is DXF (`.dxf`), with optional support for a simple ASCII (`.asc`) format.
+
+## ASCII Format (.asc)
+An .asc file contains a list of 2D coordinates, space-delimited, one point per line,
+where the first column are the x values and the second the y values.
+For example:
+```plaintext
+# ASCII
+0.0 0.0
+1.0 0.0
+1.0 1.0
+0.0 1.0
+```
+It is the user’s responsibility to ensure the points are ordered correctly.
+This format is easy to generate and inspect manually.
+
+## DXF Format (.dxf) – recommended
+The DXF (Drawing Exchange Format) is a widely-used CAD format for 2D and 3D data.
+We recommend this format for defining 2D polygonal.
+Only DXF entities of type `LINE` or `POLYLINE` are supported.
+To create DXF files from scratch, we recommend using the open-source tool [FreeCAD](https://www.freecad.org/).
+For a less technical approach, we recommend using [Inkscape](https://inkscape.org/) to create SVG files and then export them as DXF.
+
+### Creating DXF Files with FreeCAD
+
+- Open FreeCAD and create a new document.
+- Switch to the Sketcher workbench and create a new sketch.
+- Choose the XY-plane orientation and draw your geometry.
+- Select the object to be exported.
+- Go to "File > Export..." and choose "AutoDesk DXF (*.dxf)" as the format.
+- Ensure the following Import-Export options are enabled:
+    - "Use legacy Python exporter".
+    - "Project exported objects along current view direction".
+
+### Creating DXF Files with Inkscape
+SVG vector graphics can also be used as a basis for geometry.
+Use Inkscape to open or create the SVG.
+You can simply draw a Bezier curve by pressing "B" on your keyboard.
+We reommend the mode "Create spiro paths" for a smooth curve.
+Select the relevant path and export it as DXF:
+- Go to "File > Save As...".
+- Choose "Desktop Cutting Plotter (AutoCAD DXF R12)(*.dxf)" format.
 
 ```@autodocs
 Modules = [TrixiParticles]

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/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/pipe_flow_2d.jl

@@ -155,7 +155,8 @@ saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
 
 extra_callback = nothing
 
-callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback)
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(),
+                        ParticleShiftingCallback(), extra_callback)
 
 sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)

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