Move SWEs from Trixi.jl to TrixiShallowWater.jl

[patrickersing]

Accompanying PR to trixi-framework/Trixi.jl#2379.

This PR moves the following equations and related functionality from Trixi.jl to TrixiShallowWater.jl:

• ShallowWaterEquations1D
• ShallowWaterEquations2D
• ShallowWaterEquationsQuasi1D

The ShallowWaterEquationsWetDry are then renamed to ShallowWaterEquations. This should only be merged together with the accompanying PR trixi-framework/Trixi.jl#2379 in Trixi.jl.

For testing purposes the ci.yml has been modified to use the Trixi.jl version from the corresponding PR. Before merging, this needs to be cleaned up and the compat bound for Trixi.jl in the Project.toml needs to be updated.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.30%. Comparing base (dd51a30) to head (8048a78).
Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #96      +/-   ##
==========================================
+ Coverage   99.17%   99.30%   +0.12%     
==========================================
  Files          77       82       +5     
  Lines        3635     4263     +628     
==========================================
+ Hits         3605     4233     +628     
  Misses         30       30              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[patrickersing]

After moving the SWE code to TrixiSW.jl, elixir_shallowwater_perturbation_wet_dry_amr.jl is now crashing after the third time step, due to negative water height in max_abs_speeds, which is called in the StepsizeCallback. This seems to be caused by the AMR callback, which is first activated in the third time step.
The reason why this did not happen before is that max_abs_speeds was called from Trixi.jl and therefore used Trixi.sqrt(), which returns NaN and this somehow allowed the run to progress. Now max_abs_speeds in TrixiSW.jl uses Base.sqrt() which throws an error if negative water heights are encountered.
I have played around with a smaller CFL and different thresholds, but the elixir always crashes the first time the AMR callback is activated, which indicates that there might be a positivity issue with the AMR.

This is probably related to similar problems mentioned in trixi-framework/Trixi.jl#1766 and trixi-framework/Trixi.jl#2064.

[andrewwinters5000]

Below is my (experimental) attempt to implement a function that limits the mortars on P4estMesh. This redoes some of the work from prolong2mortars! so the call chain might need modified in the assembly of the right hand side. Also, this strategy to "fix" the positivity issues on AMR relies on the fact that the run has a stage_limiter!. Thus, this strategy would not fix the instabilities encountered in trixi-framework/Trixi.jl#1766 as the MHD rotor does no use the limiter.

limit_shallow_water_mortars.jl

# Could be placed in the limiter call before the "safety" application of the velocity desingularization
# or at least somewhere after the refinement/coarsening occurs
if hasproperty(cache, :mortars)
    limit_shallow_water_mortars!(u, threshold, variable, mesh, equations, dg, cache)
end

# This also performs the projection, so this actually redoes the work
# of the `prolong2mortars!` routine and should be called after it to ensure the positivity
#
# !!! warning "Experimental code"
#     This is an experimental feature and may change in future releases.
function limit_shallow_water_mortars!(u, threshold::Real, variable, mesh::P4estMesh{2},
                                      equations::ShallowWaterEquations2D, dg::DGSEM, cache)
    mortar_l2 = dg.mortar

    @unpack neighbor_ids, node_indices = cache.mortars
    @unpack mortars = cache
    @unpack weights = dg.basis
    @unpack inverse_jacobian = cache.elements
    index_range = eachnode(dg)

    Trixi.@threaded for mortar in Trixi.eachmortar(dg, cache)
        # Copy solution data from the small elements using "delayed indexing" with
        # a start value and a step size to get the correct face and orientation.

        small_indices = node_indices[1, mortar]
        small_direction = Trixi.indices2direction(small_indices)

        # Buffer to copy solution values of the large element in the correct orientation
        # before interpolating
        u_buffer = cache.u_threaded[Threads.threadid()]

        # Copy solution of large element face to buffer in the
        # correct orientation
        large_indices = node_indices[2, mortar]
        large_direction = Trixi.indices2direction(large_indices)

        i_large_start, i_large_step = Trixi.index_to_start_step_2d(large_indices[1],
                                                                   index_range)
        j_large_start, j_large_step = Trixi.index_to_start_step_2d(large_indices[2],
                                                                   index_range)

        i_large = i_large_start
        j_large = j_large_start
        element = neighbor_ids[3, mortar] # large element
        for i in eachnode(dg)
            # for v in eachvariable(equations)
            #     u_buffer[v, i] = u[v, i_large, j_large, element]
            # end
            if u[1, i_large, j_large, element] >= 2 * (threshold + eps())
                u_buffer[1, i] = u[1, i_large, j_large, element] +
                                 u[4, i_large, j_large, element]
            else
                u_buffer[1, i] = equations.H0 # from Benov et al.
            end
            u_buffer[2:4, i] = u[2:4, i_large, j_large, element]

            i_large += i_large_step
            j_large += j_large_step
        end

        # compute mean value
        u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
        total_volume = zero(eltype(u))
        for j in eachnode(dg), i in eachnode(dg)
            volume_jacobian = abs(inv(Trixi.get_inverse_jacobian(inverse_jacobian, mesh,
                                                                 i, j, element)))
            u_node = get_node_vars(u, equations, dg, i, j, element)
            u_mean += u_node * weights[i] * weights[j] * volume_jacobian
            total_volume += weights[i] * weights[j] * volume_jacobian
        end
        # normalize with the total volume
        u_mean = u_mean / total_volume

        # We compute the value directly with the mean values, as we assume that
        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
        value_mean = variable(u_mean, equations)

        # Interpolate large element face data from buffer to small face locations
        Trixi.multiply_dimensionwise!(view(mortars.u, 2, :, 1, :, mortar),
                                      mortar_l2.forward_lower,
                                      u_buffer)
        Trixi.multiply_dimensionwise!(view(mortars.u, 2, :, 2, :, mortar),
                                      mortar_l2.forward_upper,
                                      u_buffer)

        for i in eachnode(dg)
            mortars.u[2, 1, 1, i, mortar] = max(mortars.u[2, 1, 1, i, mortar] -
                                                mortars.u[2, 4, 1, i, mortar], threshold)
            mortars.u[2, 1, 2, i, mortar] = max(mortars.u[2, 1, 2, i, mortar] -
                                                mortars.u[2, 4, 2, i, mortar], threshold)
        end

        value_min = typemax(eltype(u))
        for i in eachnode(dg)
            u_node1 = @view mortars.u[2, :, 1, i, mortar]
            u_node2 = @view mortars.u[2, :, 2, i, mortar]
            value_min = min(value_min, variable(u_node1, equations),
                            variable(u_node2, equations))
        end

        # detect if limiting is necessary
        value_min < threshold || continue

        # Correct u so that mortar values will be positive
        # theta = min(abs((value_mean - threshold) / (value_mean - value_min)), 1)
        theta = (value_mean - threshold) / (value_mean - value_min)
        for j in eachnode(dg), i in eachnode(dg)
            u_node = get_node_vars(u, equations, dg, i, j, element)

            # Cut off velocity in case that the water height is smaller than the threshold

            h_node, h_v1_node, h_v2_node, b_node = u_node
            h_mean, h_v1_mean, h_v2_mean, _ = u_mean # b_mean is not used as it must not be overwritten

            if h_node <= threshold
                h_v1_node = zero(eltype(u))
                h_v2_node = zero(eltype(u))
                h_v1_mean = zero(eltype(u))
                h_v2_mean = zero(eltype(u))
            end

            u_node = SVector(h_node, h_v1_node, h_v2_node, b_node)
            u_mean = SVector(h_mean, h_v1_mean, h_v2_mean, b_node)

            # When velocities are cut off, the only averaged value is the water height,
            # because the velocities are set to zero and this value is passed.
            # Otherwise, the velocities are averaged, as well.
            # Note that the auxiliary bottom topography variable `b` is never limited.
            set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
                           equations, dg, i, j, element)
        end
    end
end

[patrickersing]

The results for elixir_shallowwater_ec_float32.jl look reasonable, but vary slightly from the expected test values which causes the tests to fail. @ranocha do you have an idea what could cause these values to differ?

[andrewwinters5000]

Just one last small suggestion, but overall it looks good to me!

[ranocha]

Currently, there are a significant number of CI errors in Trixi.jl. We need to resolve them to ensure that everything is functioning properly before making the breaking release. I apologize for the continued wait...

[andrewwinters5000]

With the CI fixed due to the Trixi PR #2415 this is good to go once the next breaking release of Trixi is ready.

[ranocha]

Just a few reminders

[ranocha]

Let's re-start CI: Trixi.jl v0.12 has been registered.

[ranocha]

This should be fine when CI passes. Thanks!

[patrickersing]

The CI is looking good now, so this should be ready to merge.

[ranocha]

Thanks for all your work!