Nonconforming / AMR improvements

[andrewwinters5000]

Functionality for well-balanced mortars and nonconforming AMR approximations for the shallow water equations with wet/dry interfaces on P4estMesh was prototyped in #45. Everything works well for fully wet flow configurations and well-balancedness is maintained for fully wet or wet/dry lake-at-rest steady-states. However, the combination of well-balanced mortars + AMR + wet/dry leads to an appreciable loss in conservation of the total water height. The conservation errors are between 1e-4 to 1e-5, compared to ≈1e-13 for the fully wet case with AMR. Also, robustness of the wet/dry + AMR is very sensitive to velocity desingularization.

So, of principle interest is how one can make the well-balanced mortars + AMR + wet/dry configuration conservative as well.

Other concerns to improve the AMR and nonconforming functionality are:

• Add modified Audusse et al. hydrostatic reconstruction to work for wet/dry transitions and add tests.
  • Done in Move SWEs from Trixi.jl to TrixiShallowWater.jl #96
• Adapt the implementation to work on TreeMesh. Dropping to this simpler mesh type might also help track down where the conservation loss occurs
• Sometimes including an extra velocity desingularization in the limiter "safety" step (at line 81 in src/callback_stage/positivity_shallow_water_dg2d.jl) helps with robustness; why?
• Apply limiting of the solution on the mortars in a similar fashion used in TrixiLW.jl. A version of this limiting is given as source code below, although in testing it leads to a loss of well-balancedness
• Apply limiting and/or desingularization after the refine or coarsen in the AMR process to ensure admissible solution states.
  • Done in Add call to the positivity limiter after coarsen / refine #97
• Specialize the coarsen and refine in the AMR steps to avoid introducing discontinuities in the bottom topography. One idea would be to re-evaluate the bottom topography function (via the initial_condition) after either process is performed. However, these functions do not have "knowledge" of the initial_condition function as it is stored in the Semidiscretization struct.
• Adapt the strategy to test it with the MultilayerEquations with a single layer. It carries the pressure contribution directly in the nonconservative term and it should be possible to simplify the back projection of the fluxes from the mortar back to the parent.
  • Done in Well-balanced mortars for ShallowWaterMultiLayerEquations2D #102
• Investigate the limiter and what occurs when the u_mean drops below zero. This is likely directly related to the conservation loss.

limit_mortars.jl

# call signature at the end of the limiter call at line 70

    if hasproperty(cache, :mortars)
        limit_shallow_water_mortars!(u, threshold, variable, mesh, equations, dg, cache)
    end

# !!! 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::ShallowWaterEquationsWetDry2D, 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
</details>