removing top wall in RectangularTank

[Nour-Touil]

Hi,

Is there a recommended way to simulate an open-top tank using RectangularTank in TrixiParticles.jl—specifically by removing the top boundary ? thanks in advance .

[LasNikas]

Yes, you can set the faces by passing a tuple to the faces kwarg. For an example, take a look at the hydrostatic_water_column_2d.jl example and modify the following code:

tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
                       n_layers=boundary_layers, acceleration=(0.0, -gravity),
                       state_equation=state_equation, 
                       faces=(true, true, true, false))

This configuration allows you to specify which faces of the tank are active. For instance, (true, true, true, false) activates the left, right, and bottom faces, while deactivating the top face.

Additionally, you can visualize the initial condition using the Plots.jl package:

using Plots

plot(tank.fluid, tank.boundary)

[Nour-Touil]

it's working now , just a last question , is it possible to turn a 2D simulation like the moving_wall to 3D ?

[LasNikas]

You're welcome! Of course. You can take a look at the example files dam_break_2d.jl and moving_wall_2d.jl, for guidance. Use the function reset_wall! for the tank in the dam break example and then define the movement function. Feel free to reach out if you encounter any difficulties or have further questions.

[LasNikas]

Here is the modified moving wall example extended to 3D.

using TrixiParticles
using OrdinaryDiffEq

# ==========================================================================================
# ==== Resolution
fluid_particle_spacing = 0.05

# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
boundary_layers = 3

# ==========================================================================================
# ==== Experiment Setup
gravity = 9.81
tspan = (0.0, 2.0)

# Boundary geometry and initial fluid particle positions
initial_fluid_size = (1.0, 0.8, 0.8)
tank_size = (4.0, 1.0, 0.8)

fluid_density = 1000.0
sound_speed = 10 * sqrt(gravity * initial_fluid_size[2])
state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
                                   exponent=7)

tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
                       n_layers=boundary_layers, spacing_ratio=1.0,
                       acceleration=(0.0, -gravity, 0.0), state_equation=state_equation)

reset_wall!(tank, (false, true, false, false, false, false), (0.0, tank.fluid_size[1], 0.0, 0.0, 0.0,0.0))

# Movement function
movement_function(t) = SVector(0.5t^2, 0.0, 0.0)

is_moving(t) = t < 1.5

boundary_movement = BoundaryMovement(movement_function, is_moving,
                                     moving_particles=tank.face_indices[2])

# ==========================================================================================
# ==== Fluid
smoothing_length = 1.2 * fluid_particle_spacing
smoothing_kernel = SchoenbergCubicSplineKernel{ndims(tank.boundary)}()

fluid_density_calculator = ContinuityDensity()
viscosity = ArtificialViscosityMonaghan(alpha=0.1, beta=0.0)

fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                           state_equation, smoothing_kernel,
                                           smoothing_length, viscosity=viscosity,
                                           acceleration=(0.0, -gravity, 0.0))

# ==========================================================================================
# ==== Boundary
boundary_density_calculator = AdamiPressureExtrapolation()
boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                             state_equation=state_equation,
                                             boundary_density_calculator,
                                             smoothing_kernel, smoothing_length)

boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
                                    movement=boundary_movement)

# ==========================================================================================
# ==== Simulation
semi = Semidiscretization(fluid_system, boundary_system)
ode = semidiscretize(semi, tspan)

info_callback = InfoCallback(interval=100)
saving_callback = SolutionSavingCallback(dt=0.02, prefix="")

callbacks = CallbackSet(info_callback, saving_callback)

# Use a Runge-Kutta method with automatic (error based) time step size control.
# Limiting of the maximum stepsize is necessary to prevent crashing.
# When particles are approaching a wall in a uniform way, they can be advanced
# with large time steps. Close to the wall, the stepsize has to be reduced drastically.
# Sometimes, the method fails to do so because forces become extremely large when
# fluid particles are very close to boundary particles, and the time integration method
# interprets this as an instability.
sol = solve(ode, RDPK3SpFSAL35(),
            abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
            dtmax=1e-2, # Limit stepsize to prevent crashing
            save_everystep=false, callback=callbacks);

[Nour-Touil]

Hi again @LasNikas and i'm sorry if i'm asking a lot , i visualized a simulation based on TrixiParticle and what i noticed that the parametres i can analyse are like density , velocity and pressure but i need also the water high , i tried to use the density to determine the water high but i didn't get it , is there any other solution ? thank you in advance.

[LasNikas]

Do you need the maximum or the average water height? For the maximum height, you can simply check the particle coordinates to determine the highest point. However, if you want to detect the free surface, this would require a free-surface detection method, which is not (yet) implemented in TrixiParticles.jl.

I hope this helps! Let me know if you have any further questions

[Nour-Touil]

Thank you for replying what i really need is plot the high of fluid particles over time in a precise position , like for x=5 i need to plot a graph of all the high in that point over the simulation time

[LasNikas]

You can extract the position using the post-processing callback (see PostprocessCallback documentation) and then visualize it using a suitable tool, such as Plots.jl or ParaView.

[svchb]

Alternatively you can use the tools provided in ParaView. For example a combination of calculators, line extraction and than plot it in a plot view. Or you can use python in ParaView. https://docs.paraview.org/en/latest/index.html

Or you can use python/julia to read in the vtk files and extract the particle positions. https://github.com/JuliaVTK/ReadVTK.jl https://community.geodynamics.org/t/loading-vtk-data-into-python/1661

[Nour-Touil]

Hi @svchb , can you explain more pls because i didn't get to make it work

[svchb]

If you are not familiar with ParaView it might be simpler to do what @LasNikas described by using the PostprocessCallback.

For example write every 0.02s the maximum height assuming height is the z-axis:

function max_z_coord(system, data, t)
    return maximum(j -> data.coordinates[3, j], axes(data.coordinates, 2))
end

postprocessing_cb = PostprocessCallback(; dt=0.02, output_directory="out",
                                        filename="max_height", max_z_coord)

[Nour-Touil]

yeah I got it , but it's not used for a specific point !

[svchb]

You can use in this function any code:

function test_func(system, data, t)

  for coord_z in data.coordinates[3, :]
   if coord_z> 0.5 && coord_z < 0.4
     ....
   end
end

    return value
end

postprocessing_cb = PostprocessCallback(; dt=0.02, output_directory="out",
                                        filename="max_height", test_func)

[Nour-Touil]

Hello again,

To extract the fluid height at a specific point, I previously tried using ParaView and Python, but I couldn't obtain valid results. So now, I would like to use a PostprocessCallback instead. However, I wasn’t able to get it working correctly.

Could you please help me integrate a working PostprocessCallback into my script?

https://drive.google.com/file/d/1-C2CdLc53ny82wXTvd2dtVEj4nxd4RTc/view?usp=sharing

Thank you in advance!

[svchb]

using TrixiParticles
using OrdinaryDiffEq
using Plots
using StaticArrays

# =======================
# Paramètres du cas Madsen (1970)
# =======================
h = 0.46
A = 0.1006
h_over_L = 0.132
L = h / h_over_L
k = 2π / L
ω = sqrt(9.81 * k)
T = 2π / ω
n_cycles = 11

t_stop = n_cycles * T
tspan = (0.0, t_stop + 5.0)

# =======================
# Domaine et fluide
# =======================
fluid_particle_spacing = 0.05
boundary_layers = 3
gravity = 9.81

initial_fluid_size = (22.0, h)
tank_size = (22.0, 2.0)

fluid_density = 1000.0
sound_speed = 10 * sqrt(gravity * h)
state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=7)

tank = RectangularTank(
    fluid_particle_spacing,
    initial_fluid_size,
    tank_size,
    fluid_density,
    n_layers=boundary_layers,
    spacing_ratio=1.0,
    acceleration=(0.0, -gravity),
    state_equation=state_equation,
    faces=(true, true, true, false)
)

reset_wall!(tank, (false, true, false, false), (0.0, tank_size[1], 0.0, 0.0))

# =======================
# Mouvement du piston
# =======================
movement_function(t) = SVector(A * cos(ω * t / 2)^2, 0.0)
is_moving(t) = t ≤ t_stop

boundary_movement = BoundaryMovement(
    movement_function,
    is_moving,
    moving_particles = tank.face_indices[1]
)

# =======================
# Système SPH
# =======================
smoothing_length = 1.2 * fluid_particle_spacing
smoothing_kernel = SchoenbergCubicSplineKernel{2}()

fluid_density_calculator = ContinuityDensity()
viscosity = ArtificialViscosityMonaghan(alpha=0.01, beta=0.0)

fluid_system = WeaklyCompressibleSPHSystem(
    tank.fluid,
    fluid_density_calculator,
    state_equation,
    smoothing_kernel,
    smoothing_length,
    viscosity=viscosity,
    acceleration=(0.0, -gravity)
)

boundary_density_calculator = AdamiPressureExtrapolation()
boundary_model = BoundaryModelDummyParticles(
    tank.boundary.density,
    tank.boundary.mass,
    state_equation=state_equation,
    boundary_density_calculator,
    smoothing_kernel,
    smoothing_length
)

boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, movement=boundary_movement)

# =======================
# Résolution de l’ODE
# =======================
semi = Semidiscretization(fluid_system, boundary_system)
ode = TrixiParticles.semidiscretize(semi, tspan)

info_callback = InfoCallback(interval=100)
saving_callback = SolutionSavingCallback(dt=0.01, prefix="")

function test_func(system, data, t)

  for coord_z in data.coordinates[3, :]
   if coord_z> 0.5 && coord_z < 0.4
     ....
   end
end

    return value
end

postprocessing_cb = PostprocessCallback(; dt=0.02, output_directory="out",
                                        filename="max_height", test_func)

callbacks = CallbackSet(info_callback, saving_callback, postprocessing_cb)

sol = solve(ode, RDPK3SpFSAL35();
    abstol=1e-6,
    reltol=1e-4,
    dtmax=1e-2,
    save_everystep=false,
    callback=callbacks
)

plot(sol)