Fix dam break with corrections (#235)

* Fix dam break with corrections

* Add error check for `trixi_include`

* Fix tests

* Add test for `trixi_include`

* Fix typo

examples/fluid/dam_break_2d.jl

@@ -10,7 +10,7 @@ using OrdinaryDiffEq
 
 # Constants
 gravity = 9.81
-athmospheric_pressure = 100000.0
+atmospheric_pressure = 100000.0
 fluid_density = 1000.0
 
 # Simulation settings
@@ -51,8 +51,8 @@ move_wall(tank, tank.fluid_size[1])
 # ==== Fluid
 sound_speed = 20 * sqrt(gravity * initial_fluid_size[2])
 
-state_equation = StateEquationCole(sound_speed, 7, fluid_density, athmospheric_pressure,
-                                   background_pressure=athmospheric_pressure)
+state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure,
+                                   background_pressure=atmospheric_pressure)
 
 smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 

examples/fluid/dam_break_2d_corrections.jl

@@ -28,23 +28,33 @@ density_calculator_dict = Dict(
 )
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-              particle_spacing=particle_spacing, smoothing_length=smoothing_length,
+              fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length,
               boundary_density_calculator=ContinuityDensity(),
               fluid_density_calculator=ContinuityDensity(),
               correction=Nothing(), use_reinit=true,
               relaxation_step_file_prefix="relaxation_continuity_reinit",
               simulation_step_file_prefix="continuity_reinit",
               relaxation_tspan=relaxation_tspan, simulation_tspan=simulation_tspan)
 
+# Clip negative pressure to be able to use `SummationDensity`
+state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure,
+                                   background_pressure=atmospheric_pressure,
+                                   clip_negative_pressure=true)
+
 for correction_name in keys(correction_dict)
     local fluid_density_calculator = density_calculator_dict[correction_name]
     local correction = correction_dict[correction_name]
 
+    println("="^100)
+    println("fluid/dam_break_2d.jl with ", correction_name)
+
     trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-                  particle_spacing=particle_spacing, smoothing_length=smoothing_length,
+                  fluid_particle_spacing=particle_spacing,
+                  smoothing_length=smoothing_length,
                   boundary_density_calculator=boundary_density_calculator,
                   fluid_density_calculator=fluid_density_calculator,
                   correction=correction, use_reinit=false,
+                  state_equation=state_equation,
                   relaxation_step_file_prefix="relaxation_$(correction_name)",
                   simulation_step_file_prefix="$(correction_name)",
                   relaxation_tspan=relaxation_tspan, simulation_tspan=simulation_tspan)

src/util.jl

@@ -170,6 +170,16 @@ julia> redirect_stdout(devnull) do
 ```
 """
 function trixi_include(mod::Module, elixir::AbstractString; kwargs...)
+    # Check that all kwargs exist as assignments
+    code = read(elixir, String)
+    expr = Meta.parse("begin \n$code \nend")
+
+    for (key, val) in kwargs
+        # This will throw an error when `key` is not found
+        find_assignment(expr, key)
+    end
+
+    # Include `elixir` and replace assignments
     Base.include(ex -> replace_assignments(insert_maxiters(ex); kwargs...), mod, elixir)
 end
 
@@ -237,20 +247,26 @@ end
 function find_assignment(expr, destination)
     # declare result to be able to assign to it in the closure
     local result
+    found = false
 
     # find explicit and keyword assignments
     walkexpr(expr) do x
         if x isa Expr
             if (x.head === Symbol("=") || x.head === :kw) &&
                x.args[1] === Symbol(destination)
                 result = x.args[2]
+                found = true
                 # dump(x)
             end
         end
         return x
     end
 
-    result
+    if !found
+        throw(ArgumentError("assignment $destination not found in expression"))
+    end
+
+    return result
 end
 
 """

test/examples/examples.jl

@@ -102,7 +102,6 @@
             @test_nowarn trixi_include(@__MODULE__,
                                        joinpath(examples_dir(), "fsi",
                                                 "dam_break_gate_2d.jl"),
-                                       tspan_relaxing=(0.0, 2.0),
                                        tspan=(0.0, 0.4),
                                        dtmax=1e-3)
             @test sol.retcode == ReturnCode.Success

test/general/general.jl

@@ -1 +1,3 @@
+# `util.jl` is added here to avoid confusion with `test_util.jl`
+include("util.jl")
 include("initial_condition.jl")

test/general/util.jl

@@ -0,0 +1,23 @@
+@testset verbose=true "trixi_include" begin
+    example = """
+        x = 4
+        """
+
+    filename = tempname()
+    try
+        open(filename, "w") do file
+            write(file, example)
+        end
+
+        @test_nowarn trixi_include(filename)
+        @test @isdefined x
+        @test x == 4
+
+        @test_nowarn trixi_include(filename, x=7)
+        @test x == 7
+
+        @test_throws "assignment y not found in expression" trixi_include(filename, y=3)
+    finally
+        rm(filename, force=true)
+    end
+end