Choice of `ldg_switch` for `ViscousFormulationLocalDG`

[jlchan]

Currently we use

Trixi.jl/src/solvers/solvers_parabolic.jl

Lines 121 to 126 in 1f32f50

	function flux_parabolic(u_ll, u_rr, normal_direction,
	::Gradient, equations_parabolic,
	parabolic_scheme::ViscousFormulationLocalDG)
	ldg_switch = sign(sum(normal_direction)) # equivalent to sign(dot(normal_direction, ones))
	return 0.5f0 * (u_ll + u_rr - ldg_switch * (u_rr - u_ll))
	end

This can be numerically unstable on certain meshes if sum(normal_direction) is approximately zero. I think one can avoid this just by taking max(normal_direction[1], normal_direction[2]) instead, which should preserve that the value of ldg_switch is equal and opposite across interfaces.

[ranocha]

Would you like to make a PR, ideally with a mesh having some interfaces close to the problematic case?

[jlchan]

Yes, I'll probably do so after I finish my talk for this week 😅

[DanielDoehring]

I briefly experimented with this setup

using OrdinaryDiffEqLowStorageRK
using Trixi

###############################################################################
# semidiscretization of the linear advection-diffusion equation

diffusivity() = 1.0e-2
advection_velocity = (-1.0, 1.0)
equations = LinearScalarAdvectionEquation2D(advection_velocity)
equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)

function initial_condition_gauss_damped(x, t, equations)
    damping_factor = 1 + 4 * diffusivity() * t
    return SVector(exp(-(x[1]^2 + x[2]^2)/damping_factor) / damping_factor)
end
initial_condition = initial_condition_gauss_damped

solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)

# This maps the domain [-1, 1]^2 to a 45-degree rotated increased square
square_size() = 5.0
function mapping(xi, eta)
    x = square_size() * xi
    y = square_size() * eta
    return SVector((x - y) / sqrt(2), (x + y) / sqrt(2))
end

trees_per_dimension = (23, 23)
mesh = P4estMesh(trees_per_dimension,
                 polydeg = 3, initial_refinement_level = 0,
                 mapping = mapping, periodicity = true)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), 
                                             initial_condition, solver;
                                             solver_parabolic = ViscousFormulationLocalDG(),
                                             #solver_parabolic = ViscousFormulationBassiRebay1(),
                                             boundary_conditions = (boundary_condition_periodic,
                                                                    boundary_condition_periodic))

###############################################################################
# ODE solvers, callbacks etc.

n_passes = 2
tspan = (0.0, n_passes * square_size() * sqrt(2))
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)

###############################################################################
# run the simulation

time_int_tol = 1.0e-6
sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
            ode_default_options()..., callback = callbacks)

which gives you the rotated square

and the normals are indeed fluctuating around 0 (I printed them sum from the parabolic solver):

-1.1102230246251565e-16
-8.326672684688674e-17
5.551115123125783e-17
5.828670879282072e-16
-3.885780586188048e-16
8.326672684688674e-17
0.0
-3.3306690738754696e-16

Current implementation LDG:

 L2 error:       4.84514880e-05
 Linf error:     6.28771998e-04

Proposed implementation:

 L2 error:       5.81842989e-05
 Linf error:     5.92882222e-04

So L2 seems to grow and LInf seems to be smaller. I also tested ldg_switch = sign(max(normal_direction[1], normal_direction[2])) which gives

 L2 error:       4.85992565e-05
 Linf error:     6.44842303e-04

[jlchan]

Sorry, I meant to say sign(max(normal_direction...)) otherwise it's not a switch. Thanks for catching that.

Seems like the new approach behaves similarly to the old one (maybe a little less accurate)?