Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add tests for state equations #292

Merged
merged 1 commit into from
Nov 16, 2023
Merged

Add tests for state equations #292

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions src/schemes/fluid/weakly_compressible_sph/state_equations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,13 @@ function (state_equation::StateEquationIdealGas)(density)
return pressure
end

function inverse_state_equation(state_equation::StateEquationIdealGas, pressure)
(; sound_speed, reference_density, reference_pressure, background_pressure) = state_equation

return (pressure + background_pressure - reference_pressure) / sound_speed^2 +
reference_density
end

@doc raw"""
StateEquationCole(sound_speed, gamma, reference_density, reference_pressure;
background_pressure=0.0, clip_negative_pressure=false)
Expand Down
149 changes: 149 additions & 0 deletions test/schemes/fluid/weakly_compressible_sph/state_equation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,149 @@
@testset verbose=true "State Equations" begin
@testset verbose=true "StateEquationCole" begin
# The equation of state was designed by Cole to accurately describe the
# physical properties of water under pressures up to 25 kbar.
# We verify that it roughly coincides with online calculators (because I couldn't
# find the original data from 1942 that Cole used to design his equation).
@testset "Physical Properties of Water" begin
# Standard speed of sound for pure water at 20°C
sound_speed = 1484.0

# Density of pure water at 20°C
rest_density = 998.34

# Work with pressures in ATM
ATM = 101_325.0

# 7.15 is the value used by Cole (see p. 39 of Cole 1948)
state_equation = StateEquationCole(sound_speed, 7.15, rest_density, 1ATM)

# These densities differ from an online calculator by less than 0.05%.
# See https://www.omnicalculator.com/physics/water-density
@test TrixiParticles.inverse_state_equation(state_equation, 1ATM) == 998.34
@test TrixiParticles.inverse_state_equation(state_equation, 100ATM) ==
1002.8323123356663
@test TrixiParticles.inverse_state_equation(state_equation, 500ATM) ==
1019.8235062499685

# The online calculator says it's only accurate up to pressures of 1000 bar,
# while Cole designed his equation to be accurate up to 25 kbar.
# Therefore, these densities differ by up to 0.6%.
@test TrixiParticles.inverse_state_equation(state_equation, 1000ATM) ==
1038.8747989986027
@test TrixiParticles.inverse_state_equation(state_equation, 3000ATM) ==
1099.0607413267035
@test TrixiParticles.inverse_state_equation(state_equation, 9000ATM) ==
1210.4689472510186
end

@testset "Background Pressure and Clipping" begin
# Test background pressure
background_pressures = [0.0, 10_000.0, 100_000.0, 200_000.0]

for background_pressure in background_pressures
state_equation = StateEquationCole(10.0, 7, 1000.0, 100_000.0,
background_pressure=background_pressure)
@test state_equation(1000.0) == 100_000.0 - background_pressure
@test state_equation(1001.0) > 100_000.0 - background_pressure + 10
# No pressure clipping
@test state_equation(999.0) < 100_000.0 - background_pressure - 10
end

# Test pressure clipping
state_equation = StateEquationCole(10.0, 7, 1000.0, 100_000.0,
background_pressure=100_000.0,
clip_negative_pressure=true)
@test state_equation(999.0) == 0.0
@test state_equation(900.0) == 0.0
end
end

@testset verbose=true "StateEquationIdealGas" begin
# This equation of state does not accurately describe the physical properties of
# water under high pressures. However, it gives good results for low pressures,
# where the relation of pressure and density is almost linear.
# We verify again that it roughly coincides with online calculators.
@testset "Physical Properties of Water" begin
# Standard speed of sound for pure water at 20°C
sound_speed = 1484.0

# Density of pure water at 20°C
rest_density = 998.34

# Work with pressures in ATM
ATM = 101_325.0

state_equation = StateEquationIdealGas(sound_speed, rest_density, 1ATM)

# These densities differ from an online calculator by less than 0.1%.
# See https://www.omnicalculator.com/physics/water-density
@test TrixiParticles.inverse_state_equation(state_equation, 1ATM) == 998.34
@test TrixiParticles.inverse_state_equation(state_equation, 100ATM) ==
1002.8949541016121
@test TrixiParticles.inverse_state_equation(state_equation, 500ATM) ==
1021.2988090576209

# For higher pressures, this state equation fails. These densities differ
# by up to 16%.
@test TrixiParticles.inverse_state_equation(state_equation, 1000ATM) ==
1044.3036277526319
@test TrixiParticles.inverse_state_equation(state_equation, 3000ATM) ==
1136.3229025326757
@test TrixiParticles.inverse_state_equation(state_equation, 9000ATM) ==
1412.380726872807
end

@testset "Background Pressure and Clipping" begin
# Test background pressure
background_pressures = [0.0, 10_000.0, 100_000.0, 200_000.0]

for background_pressure in background_pressures
state_equation = StateEquationIdealGas(10.0, 1000.0, 100_000.0,
background_pressure=background_pressure)
@test state_equation(1000.0) == 100_000.0 - background_pressure
@test state_equation(1001.0) > 100_000.0 - background_pressure + 10
# No pressure clipping
@test state_equation(999.0) < 100_000.0 - background_pressure - 10
end

# Test pressure clipping
state_equation = StateEquationIdealGas(10.0, 1000.0, 100_000.0,
background_pressure=100_000.0,
clip_negative_pressure=true)
@test state_equation(999.0) == 0.0
@test state_equation(900.0) == 0.0
end
end

# Don't show all state equations in the final overview
@testset verbose=false "inverse_state_equation" begin
# Verify that the `inverse_state_equation` actually is the inverse
state_equations = [
StateEquationCole(1484.0, 7.15, 998.34, 101_325.0),
StateEquationCole(10.0, 7, 1000.0, 100_000.0, background_pressure=90_000.0),
StateEquationCole(10.0, 7, 1000.0, 100_000.0, background_pressure=100_000.0),
StateEquationCole(10.0, 7, 1000.0, 100_000.0, background_pressure=200_000.0),
StateEquationIdealGas(1484.0, 998.34, 101_325.0),
StateEquationIdealGas(10.0, 1000.0, 100_000.0),
StateEquationIdealGas(10.0, 1000.0, 100_000.0, background_pressure=10_000.0),
StateEquationIdealGas(10.0, 1000.0, 100_000.0, background_pressure=100_000.0),
StateEquationIdealGas(10.0, 1000.0, 100_000.0, background_pressure=200_000.0),
]

densities = [100.0, 500.0, 900.0, 990.0, 1000.0, 1005.0, 1100.0, 1600.0]
pressures = [-100.0, 0.0, 100.0, 10_000.0, 100_000.0, 100_000_000.0]

@testset "$state_equation" for state_equation in state_equations
for density in densities
pressure = state_equation(density)
@test TrixiParticles.inverse_state_equation(state_equation,
pressure) ≈ density
end

for pressure in pressures
density = TrixiParticles.inverse_state_equation(state_equation, pressure)
@test isapprox(state_equation(density), pressure, atol=2e-7, rtol=1e-10)
end
end
end
end
1 change: 1 addition & 0 deletions test/schemes/fluid/weakly_compressible_sph/weakly_compressible_sph.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
include("density_diffusion.jl")
include("state_equation.jl")