Custom Advection/diffusion-like equation with limiter, unstructured grid and Robin Boundary condition

[cgadal]

Hi all,

I'm moving here the discussion started on the Julia forum as it's easier to write here.

Problem Summary

I am trying to solve the following equation:

$$\frac{\partial \phi}{\partial t} + \nabla\cdot(\phi \boldsymbol{u}) - \frac{\partial}{\partial z} \left[\mathcal{F} \phi(1-\phi)\right] = \nabla \cdot (\mathcal{D} \nabla\phi),$$

where:

• $\phi$ is a scalar, $\boldsymbol{u}$ is a vector field advecting the scalar
• $\mathcal{D}$ and $\mathcal{F}$ are coefficients that can be taken constant, but that will probably become functions of $\phi$ but also $x$ and $z$.

The domain is bounded by two parabolas (red and brown) :

and the black lines represent streamlines of the advection field $\boldsymbol{u}$.

The boundary conditions are non-flux through the boundaries. As the vector field also respects this condition, it becomes:

$$ ( \mathcal{D} \nabla\phi ) \cdot \boldsymbol{n} + \mathcal{F} \phi(1-\phi) n_{z} = 0,$$

where $\boldsymbol{n} = (n_{x}, n_{z})$ is a vector normal to the boundary.

Workflow

As I am new to Trixi.jl, I will try to follow the current plan:

1. implement 1D equation: $\frac{\partial \phi}{\partial t} + \frac{\partial}{\partial z} (F \phi(1-\phi) = 0$
2. add the diffusion term $\nabla \cdot (\mathcal{D} \nabla\phi)$ and the Robin boundary condition
3. move to 2d and unstructured mesh
4. add the second convection term $\nabla\cdot(\phi \boldsymbol{u})$

Current state

I am currently trying to implement the 1D hyperbolic equation $\frac{\partial \phi}{\partial t} + \frac{\partial}{\partial z} (F \phi(1-\phi) = 0$.

Here is the code I've come up with based on tutorials and examples:

import Trixi as TR
import OrdinaryDiffEq as ODE
import Plots

# ### define equation
F = -1
phi_0 = 0.8

struct SegregationEquation <: TR.AbstractEquations{1, # number of spatial dimensions
    1} # number of primary variables, i.e. scalar
end

TR.flux(u, orientation, equation::SegregationEquation) = F .* u .* (TR.SVector(1) - u)
TR.varnames(_, ::SegregationEquation) = ("scalar",)
equation = SegregationEquation()

# characteristics
function TR.max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
    ::SegregationEquation)
    advection_velocity_ll = F .* (1 - 2 * u_ll[1])
    advection_velocity_rr = F .* (1 - 2 * u_rr[1])
    return max(TR.abs(advection_velocity_ll), TR.abs(advection_velocity_rr))
end

TR.max_abs_speeds(u, equations::SegregationEquation) = TR.abs(F .* (1 - 2 .* u[1]))

# ### set-up discretization
mesh = TR.TreeMesh(-1.0, 1.0, # min/max coordinates
    initial_refinement_level=5,
    n_cells_max=10^4,
    periodicity=false)

volume_flux = TR.flux_central
surface_flux = TR.flux_lax_friedrichs
solver = TR.DGSEM(polydeg=3, surface_flux=surface_flux) # set polynomial degree to 3

# #### initial coniditions
initial_condition_plain(x, t, equation::SegregationEquation) = phi_0

function boundary_condition_top(x, t, equation::SegregationEquation)
    return TR.SVector(0)
end

function boundary_condition_bot(x, t, equation::SegregationEquation)
    return TR.SVector(1)
end

# #### boundary condition
boundary_conditions = (x_neg=TR.BoundaryConditionDirichlet(boundary_condition_bot),
    x_pos=TR.BoundaryConditionDirichlet(boundary_condition_top))

semi = TR.SemidiscretizationHyperbolic(mesh, equation, initial_condition_plain, solver, boundary_conditions=boundary_conditions)


# ### Use temporal solver

# Create ODE problem with given time span
tspan = (0.0, 1)
ode = TR.semidiscretize(semi, tspan)

summary_callback = TR.SummaryCallback()
visualization_callback = TR.VisualizationCallback(solution_variables=TR.cons2cons, interval=5; clims=(0, 1), plot_data_creator=TR.PlotData1D)
stepsize_callback = TR.StepsizeCallback(cfl=0.5)
callbacks = TR.CallbackSet(summary_callback, visualization_callback, stepsize_callback)

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = ODE.solve(ode, ODE.SSPRK33(); dt=1.0, callback=callbacks);

Expected results: Two shocks should be generated at the top and bottom boundaries, travel and meet somewhere inside the domain.

Current issues:

1. volume flux: From this initial condition, instabilities quickly arise at the boundaries due to the discontinuity between the initial condition and the value imposed on the sides. I think I need to use a shock-capturing method, is that right? If so, this requires the definition of a two-point flux. Could you point me towards any literature that would help me derive this for my problem?

2. boundary conditions: Here, I imposed a constant $\phi$ value at the boundaries that correspond to the final state that I know a priori. However, this will not be possible in the 2D simulation with the additional convention term, for which I will need $\mathcal{F} \phi(1-\phi) n_{z} = 0$. Is there a way to use something like this instead?

Bonus question: In the past, I have successfully used Central upwind semi-discrete schemes to solve this kind of hyperbolic equations, such as described in Kurganov and Tadmor (2000), New High-Resolution Central Schemes
for Nonlinear Conservation Laws and Convection–Diffusion Equations. However, while I am far from being a specialist of numerical methods to solve PDE, I feel like their method could be used in Trixi by using the way they compute the flux as surface flux. Is that true ? More generally, could you briefly explain how their scheme relates to the methods used in Trixi?

[jlchan]

Hi @cgadal - forgive me if I'm misunderstanding something, but you're imposing boundary conditions of 0 and 1 but your initial condition is 0.8 here. Wouldn't you expect shocks/discontinuities to form since the BCs aren't consistent with the initial condition?

[jlchan]

In terms of instabilities - what do you observe specifically? Are you talking about the oscillations, or does the ODE solver crash? If it's oscillations, these are common for high order methods. Shock capturing will help suppress them, but they shouldn't crash a solver by themselves.

[cgadal]

Hi @jlchan,

Thanks for your quick answer, and sorry I should have been more explicit and uploaded a better image. Here is the same at longer time steps:

Hi @cgadal - forgive me if I'm misunderstanding something, but you're imposing boundary conditions of 0 and 1 but your initial condition is 0.8 here. Wouldn't you expect shocks/discontinuities to form since the BCs aren't consistent with the initial condition?

Yes, two shocks should form at the boundaries (at the left, a shock between 1 and 0.8, and at the right a shock between 0 and 0.8) and meet somewhere in the centre of the domain.

In terms of instabilities - what do you observe specifically? Are you talking about the oscillations, or does the ODE solver crash? If it's oscillations, these are common for high order methods. Shock capturing will help suppress them, but they shouldn't crash a solver by themselves.

The right-moving shock forms at the left and moves to the right with reasonable shape and amplitude. However, oscillations develop and grow on the right side until the solver returns nans.
The fluxes should go to 0 as soon as phi gits 0 or 1, so $\phi$ should not go outside of $[0,1]$.

I'm also confused by the non-symmetry of the result, because the equations are invariant under $\phi \to 1 - \phi$ and $x \to -x$.

But the initial condition is not symmetric, so oscillations grow and get crazy when the shock gets too big. If I use 0.2 as a starting condition, then oscillations develop in at the other end. So I guess it is linked to the size of the shock, which would point towards the use of a limiter?

Making the grid coarser indeed help as it makes the shock less sharp.

[ranocha]

The first thing to try is to set polydeg = 0 (instead of polydeg = 3) when constructing the solver. High-order methods will generate oscillations because of the Gibbs phenomenon. If that removes the oscillations as you like, you can change the solver to use a higher polynomial degree and the shock-capturing volume integral, see https://trixi-framework.github.io/Trixi.jl/stable/tutorials/shock_capturing/

[cgadal]

Setting polydeg = 0 and increasing the number of grid points indeed give the right solution:

• Is the point of using high order methods to be able to getter better approximation of solutions with coarser meshes, hence improving performances?

• The documentation specifies that shock capturing methods require the use of "two-point" fluxes, and lists "flux_ranocha, flux_shima_etal, flux_chandrashekar or flux_kennedy_gruber".

This is getting me confused, as I thought that central fluxes ( F(u_L​, u_R​) = 0.5 * ​ (f(u_L​) + f(u_R​) ) ) or Lax-Friedrichs fluxes were also two-point fluxes. Could you clarify this for me?

[ranocha]

Is the point of using high order methods to be able to getter better approximation of solutions with coarser meshes, hence improving performances?

Yes. When you have smooth (parts of the) solutions, high-order methods will perform better since they give you much better results on coarse meshes.

The documentation specifies that shock capturing methods require the use of "two-point" fluxes, and lists "flux_ranocha, flux_shima_etal, flux_chandrashekar or flux_kennedy_gruber".

You can use flux_central for now. The shock-capturing volume integral requires a non-dissipative, symmetric flux (like flux_central) and a dissipative, non-symmetric flux (like flux_lax_friedrichs). flux_shima_etal etc. replace flux_central and can be more efficient/robust/...

[cgadal]

Thanks, it works amazingly (code below). I'll think about a more efficient and robust non-dissipative flux later in the process if needed. I'm now moving to adding the parabolic part :)

Currently, the positivity limiter works in prevent phi from going negative. I suspect I could also use it to prevent phi from going above 1. Can I define another variable to pass to the limiter as phi_in = 1- phi ? Should I expand my system of equations to two variables?

import Trixi as TR
import OrdinaryDiffEq as ODE
import Plots

# #### parameter values
F = -1
phi_0 = 0.8
D = 0.01

#############################
# ##### Set-up equation ##### 
#############################

# ###### hyperbolic part ######

struct SegregationEquation <: TR.AbstractEquations{1, # number of spatial dimensions
    1} # number of primary variables, i.e. scalar
end

TR.flux(u, orientation, equation::SegregationEquation) = F .* u .* (TR.SVector(1) - u)
TR.varnames(_, ::SegregationEquation) = ("phi",)
equation_hyperbolic = SegregationEquation()

# characteristics
function TR.max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
    ::SegregationEquation)
    advection_velocity_ll = F .* (1 - 2 * u_ll[1])
    advection_velocity_rr = F .* (1 - 2 * u_rr[1])
    return max(TR.abs(advection_velocity_ll), TR.abs(advection_velocity_rr))
end

TR.max_abs_speeds(u, equations::SegregationEquation) = TR.abs(F .* (1 - 2 .* u[1]))

#################################
# ##### Boundary conditions ##### 
#################################

boundary_conditions_hyperbolic = (x_neg=TR.BoundaryConditionDirichlet((x, t, equation) -> TR.SVector(1.0)),
    x_pos=TR.BoundaryConditionDirichlet((x, t, equation) -> TR.SVector(0.0)))

###################################
# ##### Create discretization ##### 
###################################

mesh = TR.TreeMesh(-1.0, 1.0, # min/max coordinates
    initial_refinement_level=4,
    n_cells_max=10^4,
    periodicity=false)

volume_flux = TR.flux_central
surface_flux = TR.flux_lax_friedrichs
basis = TR.LobattoLegendreBasis(3)

#_______ part specific to shock capturing method
indicator_sc = TR.IndicatorHennemannGassner(equation, basis,
    alpha_max=0.5,
    alpha_min=0.001,
    alpha_smooth=true,
    variable=TR.first)

volume_integral = TR.VolumeIntegralShockCapturingHG(indicator_sc;
    volume_flux_dg=volume_flux,
    volume_flux_fv=surface_flux)
#_______

solver = TR.DGSEM(basis, surface_flux, volume_integral)

# initial conditions
initial_condition_plain(x, t, equation::SegregationEquation) = phi_0

# Create discretization
semi = TR.SemidiscretizationHyperbolic(mesh,
    equation_hyperbolic,
    initial_condition_plain,
    solver,
    boundary_conditions=boundary_conditions_hyperbolic)


#################################
# ##### Use temporal solver ##### 
#################################

# Create ODE problem with given time span
tspan = (0.0, 2)
ode = TR.semidiscretize(semi, tspan)

summary_callback = TR.SummaryCallback()
visualization_callback = TR.VisualizationCallback(solution_variables=TR.cons2cons, interval=5; clims=(0, 1), plot_data_creator=TR.PlotData1D)
stepsize_callback = TR.StepsizeCallback(cfl=0.5)
callbacks = TR.CallbackSet(summary_callback, visualization_callback, stepsize_callback)

stage_limiter! = TR.PositivityPreservingLimiterZhangShu(thresholds=(5e-6,), variables=(TR.first,))

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = ODE.solve(ode, ODE.SSPRK33(stage_limiter!); dt=1.0, callback=callbacks);

Solution for polydeg=3 and initial_refinement_level=4:

[cgadal]

I have an issue with the definition of boundary conditions when adding the parabolic term. I tried to inspire myself from the
14: Adding new parabolic terms tutorial and the 1d_laplace equation definition to come up with this definition for a 0 diffusive flux at the boundary:

struct BoundaryConditionConstantNeumann{T<:Real}
    boundary_value::T
end


@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::TR.AbstractVector,
    x, t,
    operator_type::TR.Divergence,
    equations_parabolic::Diffusion1D)
    return boundary_condition.boundary_value
end

@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::TR.AbstractVector,
    x, t,
    operator_type::TR.Gradient,
    equations_parabolic::Diffusion1D)
    return flux_inner
end

boundary_conditions_parabolic = (x_neg=BoundaryConditionConstantNeumann(0.0),
    x_pos=BoundaryConditionConstantNeumann(0.0))

but this returns the following error:

ERROR: MethodError: no method matching (::BoundaryConditionConstantNeumann{…})(::StaticArraysCore.SVector{…}, ::StaticArraysCore.SVector{…}, ::Int64, ::Int64, ::StaticArraysCore.SVector{…}, ::Float64, ::Trixi.Gradient, ::Diffusion1D{…})
The object of type `BoundaryConditionConstantNeumann{Float64}` exists, but no method is defined for this combination of argument types when trying to treat it as a callable object.

Closest candidates are:
  (::BoundaryConditionConstantNeumann)(::Any, ::Any, ::AbstractVector, ::Any, ::Any, ::Trixi.Divergence, ::Diffusion1D)
   @ Main ~/Documents/Work/Research/Projects/Granular_manchester/seggregation/stripes/numerical/test_trixi.jl/test_1D_segregation_hyperbolic_parabolic.jl:66
  (::BoundaryConditionConstantNeumann)(::Any, ::Any, ::AbstractVector, ::Any, ::Any, ::Trixi.Gradient, ::Diffusion1D)
   @ Main ~/Documents/Work/Research/Projects/Granular_manchester/seggregation/stripes/numerical/test_trixi.jl/test_1D_segregation_hyperbolic_parabolic.jl:74

Full code (click to expand)

import Trixi as TR
import OrdinaryDiffEq as ODE
import Plots

# #### parameter values
F = -1
phi_0 = 0.8
D = 0.01

#############################
# ##### Set-up equation ##### 
#############################

# ###### hyperbolic part ######

struct SegregationEquation <: TR.AbstractEquations{1, # number of spatial dimensions
    1} # number of primary variables, i.e. scalar
end

TR.flux(u, orientation, equation::SegregationEquation) = F .* u .* (TR.SVector(1) - u)
TR.varnames(_, ::SegregationEquation) = ("phi",)
equation_hyperbolic = SegregationEquation()

# characteristics
function TR.max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
    ::SegregationEquation)
    advection_velocity_ll = F .* (1 - 2 * u_ll[1])
    advection_velocity_rr = F .* (1 - 2 * u_rr[1])
    return max(TR.abs(advection_velocity_ll), TR.abs(advection_velocity_rr))
end

TR.max_abs_speeds(u, equations::SegregationEquation) = TR.abs(F .* (1 - 2 .* u[1]))

# ###### parabolic part
struct Diffusion1D{E,T} <: TR.AbstractEquationsParabolic{1,1,TR.GradientVariablesConservative}
    diffusivity::T
    equation_hyperbolic::E
end

function TR.varnames(variable_mapping, equation_parabolic::Diffusion1D)
    TR.varnames(variable_mapping, equation_parabolic.equation_hyperbolic)
end

function TR.flux(u, gradients, orientation::Integer,
    equation_parabolic::Diffusion1D)
    TR.@unpack diffusivity = equation_parabolic
    return TR.SVector(diffusivity * gradients[1])
end

equation_parabolic = Diffusion1D(D, equation_hyperbolic)

#################################
# ##### Boundary conditions ##### 
#################################

# hyperbolic
boundary_conditions_hyperbolic = (x_neg=TR.BoundaryConditionDirichlet((x, t, equation) -> TR.SVector(1.0)),
    x_pos=TR.BoundaryConditionDirichlet((x, t, equation) -> TR.SVector(0.0)))

# parabolic
struct BoundaryConditionConstantNeumann{T<:Real}
    boundary_value::T
end


@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::TR.AbstractVector,
    x, t,
    operator_type::TR.Divergence,
    equations_parabolic::Diffusion1D)
    return boundary_condition.boundary_value
end

@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::TR.AbstractVector,
    x, t,
    operator_type::TR.Gradient,
    equations_parabolic::Diffusion1D)
    return flux_inner
end

boundary_conditions_parabolic = (x_neg=BoundaryConditionConstantNeumann(0.0),
    x_pos=BoundaryConditionConstantNeumann(0.0))

###################################
# ##### Create discretization ##### 
###################################

mesh = TR.TreeMesh(-1.0, 1.0, # min/max coordinates
    initial_refinement_level=6,
    n_cells_max=10^4,
    periodicity=false)

volume_flux = TR.flux_central
surface_flux = TR.flux_lax_friedrichs
basis = TR.LobattoLegendreBasis(3)

#_______ part specific to shock capturing method
indicator_sc = TR.IndicatorHennemannGassner(equation, basis,
    alpha_max=0.5,
    alpha_min=0.001,
    alpha_smooth=true,
    variable=TR.first)

volume_integral = TR.VolumeIntegralShockCapturingHG(indicator_sc;
    volume_flux_dg=volume_flux,
    volume_flux_fv=surface_flux)
#_______

solver = TR.DGSEM(basis, surface_flux, volume_integral)

# #### initial coniditions
initial_condition_plain(x, t, equation) = phi_0

# #### boundary condition

semi = TR.SemidiscretizationHyperbolicParabolic(mesh,
    (equation_hyperbolic, equation_parabolic),
    initial_condition_plain,
    solver;
    boundary_conditions=(boundary_conditions_hyperbolic, boundary_conditions_parabolic))


# ### Use temporal solver

# Create ODE problem with given time span
tspan = (0.0, 2)
ode = TR.semidiscretize(semi, tspan)

summary_callback = TR.SummaryCallback()
visualization_callback = TR.VisualizationCallback(solution_variables=TR.cons2cons, interval=5; clims=(0, 1), plot_data_creator=TR.PlotData1D)
stepsize_callback = TR.StepsizeCallback(cfl=0.5)
callbacks = TR.CallbackSet(summary_callback, visualization_callback, stepsize_callback)

stage_limiter! = TR.PositivityPreservingLimiterZhangShu(thresholds=(5e-6,), variables=(TR.first,))

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = ODE.solve(ode, ODE.SSPRK33(stage_limiter!); dt=1.0, callback=callbacks);

[jlchan]

I think you have TR.AbstractVector instead of just AbstractVector in your Neumann BCs:

@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::TR.AbstractVector,
    x, t,
    operator_type::TR.Divergence,
    equations_parabolic::Diffusion1D)
    return boundary_condition.boundary_value
end

@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::TR.AbstractVector,
    x, t,
    operator_type::TR.Gradient,
    equations_parabolic::Diffusion1D)
    return flux_inner
end

[cgadal]

Changing to

@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::AbstractVector,
    x, t,
    operator_type::TR.Divergence,
    equations_parabolic::Diffusion1D)
    return boundary_condition.boundary_value
end

@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    normal::AbstractVector,
    x, t,
    operator_type::TR.Gradient,
    equations_parabolic::Diffusion1D)
    return flux_inner
end

results in the same error

[ranocha]

You're using a Cartesian TreeMesh, I guess? If so, you need to accept an orientation::Integer and a direction::Integer instead of a normal::AbstractVector, see, e.g.,

• https://trixi-framework.github.io/Trixi.jl/stable/conventions/#Spatial-dimensions-and-directions
• https://github.com/trixi-framework/Trixi.jl/blob/f586b0dc837cea99c1eff805065c3a1c5c39f550/src/solvers/dgsem_tree/dg_1d_parabolic.jl#L322C56-L323C45

[cgadal]

The tutorial https://trixi-framework.github.io/Trixi.jl/stable/tutorials/adding_new_parabolic_terms/ also uses a TreeMesh and defines the boundary condition with an normal::AbstractVector, which works.

What is the difference in my case?

Edit: Oh, I noticed by looking at https://github.com/trixi-framework/Trixi.jl/blob/f586b0dc837cea99c1eff805065c3a1c5c39f550/src/solvers/dgsem_tree/dg_2d_parabolic.jl that

you need to accept an orientation::Integer and a direction::Integer instead of a normal::AbstractVector

is valid just for the 1D case of cartesian TreeMesh it seems

I'm leaving the working 1D code here in case it can help someone in the future, and moving towards the 2D part :)

1D working Code (click to expand)

import Trixi as TR
import OrdinaryDiffEq as ODE
import Plots

# #### parameter values
F = -1
phi_0 = 0.8
D = 0.1

#############################
# ##### Set-up equation ##### 
#############################

# ###### hyperbolic part ######

struct SegregationEquation <: TR.AbstractEquations{1, # number of spatial dimensions
    1} # number of primary variables, i.e. scalar
end

TR.flux(u, orientation, equation::SegregationEquation) = F .* u .* (TR.SVector(1) - u)
TR.varnames(_, ::SegregationEquation) = ("phi",)
equation_hyperbolic = SegregationEquation()

# characteristics
function TR.max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
    ::SegregationEquation)
    advection_velocity_ll = F .* (1 - 2 * u_ll[1])
    advection_velocity_rr = F .* (1 - 2 * u_rr[1])
    return max(TR.abs(advection_velocity_ll), TR.abs(advection_velocity_rr))
end

TR.max_abs_speeds(u, equations::SegregationEquation) = TR.abs(F .* (1 - 2 .* u[1]))

# ###### parabolic part
struct Diffusion1D{E,T} <: TR.AbstractEquationsParabolic{1,1,TR.GradientVariablesConservative}
    diffusivity::T
    equation_hyperbolic::E
end

function TR.varnames(variable_mapping, equation_parabolic::Diffusion1D)
    TR.varnames(variable_mapping, equation_parabolic.equation_hyperbolic)
end

function TR.flux(u, gradients, orientation::Integer,
    equation_parabolic::Diffusion1D)
    dudx = gradients
    return equation_parabolic.diffusivity * dudx
end

equation_parabolic = Diffusion1D(D, equation_hyperbolic)

#################################
# ##### Boundary conditions ##### 
#################################

# hyperbolic
boundary_conditions_hyperbolic = (x_neg=TR.BoundaryConditionDirichlet((x, t, equation) -> TR.SVector(1.0)),
    x_pos=TR.BoundaryConditionDirichlet((x, t, equation) -> TR.SVector(0.0)))

# parabolic
struct BoundaryConditionConstantNeumann{T<:Real}
    boundary_value::T
end


@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    orientation::Integer,
    direction::Integer,
    x, t,
    operator_type::TR.Divergence,
    equation_parabolic::Diffusion1D)
    return boundary_condition.boundary_value
end

@inline function (boundary_condition::BoundaryConditionConstantNeumann)(flux_inner, u_inner,
    orientation::Integer,
    direction::Integer,
    x, t,
    operator_type::TR.Gradient,
    equation_parabolic::Diffusion1D)
    return flux_inner
end

boundary_conditions_parabolic = (x_neg=BoundaryConditionConstantNeumann(0.0),
    x_pos=BoundaryConditionConstantNeumann(0.0))

###################################
# ##### Create discretization ##### 
###################################

mesh = TR.TreeMesh(-1.0, 1.0, # min/max coordinates
    initial_refinement_level=6,
    n_cells_max=10^4,
    periodicity=false)

volume_flux = TR.flux_central
surface_flux = TR.flux_lax_friedrichs
basis = TR.LobattoLegendreBasis(3)

#_______ part specific to shock capturing method
indicator_sc = TR.IndicatorHennemannGassner(equation_hyperbolic, basis,
    alpha_max=0.5,
    alpha_min=0.001,
    alpha_smooth=true,
    variable=TR.first)

volume_integral = TR.VolumeIntegralShockCapturingHG(indicator_sc;
    volume_flux_dg=volume_flux,
    volume_flux_fv=surface_flux)
#_______

solver = TR.DGSEM(basis, surface_flux, volume_integral)

# #### initial coniditions
initial_condition_plain(x, t, equation) = phi_0

# #### boundary condition

semi = TR.SemidiscretizationHyperbolicParabolic(mesh,
    (equation_hyperbolic, equation_parabolic),
    initial_condition_plain,
    solver;
    boundary_conditions=(boundary_conditions_hyperbolic, boundary_conditions_parabolic))


# ### Use temporal solver

# Create ODE problem with given time span
tspan = (0.0, 2)
ode = TR.semidiscretize(semi, tspan)

# summary_callback = TR.SummaryCallback()
# visualization_callback = TR.VisualizationCallback(solution_variables=TR.cons2cons, interval=5; clims=(0, 1), plot_data_creator=TR.PlotData1D)
# stepsize_callback = TR.StepsizeCallback(cfl=0.1)
# callbacks = TR.CallbackSet(summary_callback, visualization_callback, stepsize_callback)

# stage_limiter! = TR.PositivityPreservingLimiterZhangShu(thresholds=(5e-6,), variables=(TR.first,))

# # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
# sol = ODE.solve(ode, ODE.SSPRK33(stage_limiter!); dt=1.0, callback=callbacks);


summary_callback = TR.SummaryCallback()
visualization_callback = TR.VisualizationCallback(solution_variables=TR.cons2cons, interval=5; clims=(0, 1), plot_data_creator=TR.PlotData1D)
stepsize_callback = TR.StepsizeCallback(cfl=0.5)
callbacks = TR.CallbackSet(summary_callback, visualization_callback, stepsize_callback)

stage_limiter! = TR.PositivityPreservingLimiterZhangShu(thresholds=(5e-6,), variables=(TR.first,))

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
time_int_tol = 1.0e-6
sol = ODE.solve(ode, ODE.RDPK3SpFSAL49(stage_limiter!); callback=callbacks, abstol=time_int_tol, TR.ode_default_options()..., reltol=time_int_tol);

[ranocha]

Yeah... We should probably fix/improve this in the internals 👍

[cgadal]

Okay, so I have moved to 2D and everything is going very well on cartesian meshes, created from TreeMesh or P4estMesh. I have tried to use then meshes created from gmsh but, while the solver runs, I have trouble understanding what I observe, and hence going further.

For simplicity, I have removed the parabolic part and kept only the hyperbolic part.

2D hyperbolic code on cartesian P4estMesh (click to expand)

import Trixi as TR
import OrdinaryDiffEq as ODE
import Plots

# #### parameter values
S = TR.SVector{2}([0.0, -1.0])
phi_0 = 0.8

#############################
# ##### Set-up equation ##### 
#############################

# ###### hyperbolic part ######
struct SegregationEquation{S} <: TR.AbstractEquations{2,1}
    segregation::S
end

function TR.flux(u, normal_direction::AbstractVector, equation::SegregationEquation)
    TR.@unpack segregation = equation
    a = TR.dot(segregation, normal_direction) # segregation component in normal direction
    # a = segregation[1] * normal_direction[1] + segregation[2] * normal_direction[2]
    return a .* u .* (TR.SVector(1.0) - u)
end

# characteristics
function characteristics(u, S)
    return S .* (1 - 2 .* u)
end

function TR.max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
    equation::SegregationEquation)
    TR.@unpack segregation = equation
    advection_velocity_ll = TR.dot(characteristics(u_ll[1], segregation), normal_direction)
    advection_velocity_rr = TR.dot(characteristics(u_rr[1], segregation), normal_direction)
    return max(TR.abs(advection_velocity_ll), TR.abs(advection_velocity_rr))
end

function TR.max_abs_speeds(u, equation::SegregationEquation)
    TR.@unpack segregation = equation
    return TR.abs.(characteristics(u[1], segregation))
end

TR.varnames(_, ::SegregationEquation) = ("phi",)
equation_hyperbolic = SegregationEquation(S)

#################################
# ##### Boundary conditions ##### 
#################################

# hyperbolic
function zeroflux_boundarycondition(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equation::SegregationEquation)
    return 0.0
end

boundary_conditions_hyperbolic = Dict(
    :x_neg => zeroflux_boundarycondition,
    :y_neg => zeroflux_boundarycondition,
    :x_pos => zeroflux_boundarycondition,
    :y_pos => zeroflux_boundarycondition,)

###################################
# ##### Create discretization ##### 
###################################
# #### create mesh
coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
trees_per_dimension = (4, 4)
mesh = TR.P4estMesh(trees_per_dimension,
    polydeg=3, initial_refinement_level=5,
    coordinates_min=coordinates_min, coordinates_max=coordinates_max,
    periodicity=false)

# ### set-up discretization
volume_flux = TR.flux_central
surface_flux = TR.flux_lax_friedrichs
basis = TR.LobattoLegendreBasis(3)

#_______ part specific to shock capturing method
indicator_sc = TR.IndicatorHennemannGassner(equation_hyperbolic, basis,
    alpha_max=0.5,
    alpha_min=0.001,
    alpha_smooth=true,
    variable=TR.first)

volume_integral = TR.VolumeIntegralShockCapturingHG(indicator_sc;
    volume_flux_dg=volume_flux,
    volume_flux_fv=surface_flux)
#_______

solver = TR.DGSEM(basis, surface_flux, volume_integral)

# #### initial conditions
initial_condition_plain(x, t, equation) = phi_0

semi = TR.SemidiscretizationHyperbolic(mesh,
    equation_hyperbolic,
    initial_condition_plain,
    solver;
    boundary_conditions=boundary_conditions_hyperbolic)


# ### Use temporal solver
# Create ODE problem with given time span
tspan = (0.0, 2)
ode = TR.semidiscretize(semi, tspan)

summary_callback = TR.SummaryCallback()
visualization_callback = TR.VisualizationCallback(solution_variables=TR.cons2cons, interval=5; clims=(0, 1))
alive_callback = TR.AliveCallback(alive_interval=20)
stepsize_callback = TR.StepsizeCallback(cfl=0.5)
callbacks = TR.CallbackSet(summary_callback, visualization_callback, alive_callback, stepsize_callback)

stage_limiter! = TR.PositivityPreservingLimiterZhangShu(thresholds=(5e-6,), variables=(TR.first,))

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
time_int_tol = 1.0e-6
sol = ODE.solve(ode, ODE.RDPK3SpFSAL49(stage_limiter!); callback=callbacks, abstol=time_int_tol, TR.ode_default_options()..., reltol=time_int_tol)

# sol = ODE.solve(ode, ODE.SSPRK83(stage_limiter!); dt=1.0, callback=callbacks, abstol=time_int_tol, TR.ode_default_options()..., reltol=time_int_tol);

Two shocks are generated at the bottom and at the top and they travel to meet somewhere in the center.

Code on P4estMesh imported from .inp file coming from gmsh (click to expand)

import Trixi as TR
import OrdinaryDiffEq as ODE
import Plots

# #### parameter values
S = TR.SVector{2}([0.0, -1.0])
phi_0 = 0.8

#############################
# ##### Set-up equation ##### 
#############################

# ###### hyperbolic part ######
struct SegregationEquation{S} <: TR.AbstractEquations{2,1}
    segregation::S
end

function TR.flux(u, normal_direction::AbstractVector, equation::SegregationEquation)
    TR.@unpack segregation = equation
    a = TR.dot(segregation, normal_direction) # segregation component in normal direction
    # a = segregation[1] * normal_direction[1] + segregation[2] * normal_direction[2]
    return a .* u .* (TR.SVector(1.0) - u)
end

# characteristics
function characteristics(u, S)
    return S .* (1 - 2 .* u)
end

function TR.max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
    equation::SegregationEquation)
    TR.@unpack segregation = equation
    advection_velocity_ll = TR.dot(characteristics(u_ll[1], segregation), normal_direction)
    advection_velocity_rr = TR.dot(characteristics(u_rr[1], segregation), normal_direction)
    return max(TR.abs(advection_velocity_ll), TR.abs(advection_velocity_rr))
end

function TR.max_abs_speeds(u, equation::SegregationEquation)
    TR.@unpack segregation = equation
    return TR.abs.(characteristics(u[1], segregation))
end

TR.varnames(_, ::SegregationEquation) = ("phi",)
equation_hyperbolic = SegregationEquation(S)

#################################
# ##### Boundary conditions ##### 
#################################

# hyperbolic
function zeroflux_boundarycondition(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equation::SegregationEquation)
    return 0.0
end

boundary_conditions_hyperbolic = Dict(
    :ExtBoundary => zeroflux_boundarycondition,
)

###################################
# ##### Create discretization ##### 
###################################
# #### load mesh
mesh_file = "numerical/test_trixi.jl/mesh_test.inp"
boundary_symbols = [:ExtBoundary]
mesh = TR.P4estMesh{2}(mesh_file, boundary_symbols=boundary_symbols, polydeg=3)

# ### set-up discretization
volume_flux = TR.flux_central
surface_flux = TR.flux_lax_friedrichs
basis = TR.LobattoLegendreBasis(3)

#_______ part specific to shock capturing method
indicator_sc = TR.IndicatorHennemannGassner(equation_hyperbolic, basis,
    alpha_max=0.5,
    alpha_min=0.001,
    alpha_smooth=true,
    variable=TR.first)

volume_integral = TR.VolumeIntegralShockCapturingHG(indicator_sc;
    volume_flux_dg=volume_flux,
    volume_flux_fv=surface_flux)
#_______

solver = TR.DGSEM(basis, surface_flux, volume_integral)

# #### initial coniditions
initial_condition_plain(x, t, equation) = phi_0

semi = TR.SemidiscretizationHyperbolic(mesh,
    equation_hyperbolic,
    initial_condition_plain,
    solver;
    boundary_conditions=boundary_conditions_hyperbolic)


# ### Use temporal solver
# Create ODE problem with given time span
tspan = (0.0, 2)
ode = TR.semidiscretize(semi, tspan)

summary_callback = TR.SummaryCallback()
visualization_callback = TR.VisualizationCallback(solution_variables=TR.cons2cons, interval=5; clims=(0, 1))
alive_callback = TR.AliveCallback(alive_interval=20)
stepsize_callback = TR.StepsizeCallback(cfl=0.5)
callbacks = TR.CallbackSet(summary_callback, visualization_callback, alive_callback, stepsize_callback)

stage_limiter! = TR.PositivityPreservingLimiterZhangShu(thresholds=(5e-6,), variables=(TR.first,))

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
time_int_tol = 1.0e-6
sol = ODE.solve(ode, ODE.RDPK3SpFSAL49(stage_limiter!); callback=callbacks, abstol=time_int_tol, TR.ode_default_options()..., reltol=time_int_tol)

# sol = ODE.solve(ode, ODE.SSPRK83(stage_limiter!); dt=1.0, callback=callbacks, abstol=time_int_tol, TR.ode_default_options()..., reltol=time_int_tol);

The result that I do not know how to interpret is:

The mesh is created using julia API from gmsh:

Code to create .inp file from gmsh julia API (click to expand)

import Gmsh: gmsh

xmin, xmax = -1, 1
ymin, ymax = -1, 1
n_width = 60
typical_meshsize = (xmax - xmin) / n_width

gmsh.initialize()
gmsh.model.add("domain")

bot_left = gmsh.model.geo.addPoint(xmin, ymin, 0, typical_meshsize)
bot_right = gmsh.model.geo.addPoint(xmax, ymin, 0, typical_meshsize)
top_right = gmsh.model.geo.addPoint(xmax, ymax, 0, typical_meshsize)
top_left = gmsh.model.geo.addPoint(xmin, ymax, 0, typical_meshsize)

left_line = gmsh.model.geo.addSpline([bot_left, top_left])
top_line = gmsh.model.geo.addSpline([top_left, top_right])
right_line = gmsh.model.geo.addSpline([top_right, bot_right])
bot_line = gmsh.model.geo.addSpline([bot_right, bot_left])

gmsh.model.geo.removeAllDuplicates()
gmsh.model.geo.addCurveLoop([left_line, top_line, right_line, bot_line], 1)
surface = gmsh.model.geo.addPlaneSurface([1], 1)

gmsh.model.geo.synchronize()
gmsh.model.mesh.setRecombine(2, 1) # Recombine triangles into quadrilaterals

gmsh.model.addPhysicalGroup(1, [left_line, right_line, bot_line, top_line], 1, "ExtBoundary")

gmsh.option.setNumber("Mesh.SaveAll", 1)
gmsh.option.setNumber("Mesh.SaveGroupsOfElements", 1)
gmsh.option.setNumber("Mesh.SaveGroupsOfNodes", 1) # Important for p4est and Trixi

gmsh.model.mesh.generate(2)

gmsh.fltk.run()
gmsh.write("mesh_test.inp")
gmsh.finalize()

Any suggestions on where to start looking?

[ranocha]

Can you just evaluate an initial condition and check what it looks like?

[cgadal]

I was not sure how to check the initial condition so I ran the code with tspan = (0.0, 0.0), and then all(sol.u[1] .== phi_0) returns true. So all cells are properly set to phi_0.

Things seem to be a bit better when I decrease polydeg in both P4estMesh and LobattoLegendreBasis to 1.

It suggests a wavy behaviour similar to what I had at the beginning in 1D without the shock-capturing method developing where the shape shock is (at the top, 0 to 1).

This is confirmed by changing phi_0 to 0.2, hence putting the sharp shock at the bottom and leading to:

where the oscillating behaviour is now at the bottom.

[jlchan]

Since this behavior is near the boundary, what happens if you use the Gmsh BCs with the Cartesian P4estMesh?

[cgadal]

Since this behaviour is near the boundary, what happens if you use the Gmsh BCs with the Cartesian P4estMesh?

Sorry, I do not understand what you mean by this. The boundary condition function is the same for both, and also the same everywhere around the domain. The only things that changes from this point is how this function is passed to the discretization:

Cartesian P4estMesh:

boundary_conditions_hyperbolic = Dict(
    :x_neg => zeroflux_boundarycondition,
    :y_neg => zeroflux_boundarycondition,
    :x_pos => zeroflux_boundarycondition,
    :y_pos => zeroflux_boundarycondition)

GMSH:

  boundary_conditions_hyperbolic = Dict(
      :ExtBoundary => zeroflux_boundarycondition,
  )

[jlchan]

Sorry, I read your code incorrectly. I am trying to see if this is a mathematical issue or a bug. Is it possible to generate a uniform mesh using Gmsh to see if it recovers the Cartesian P4estMesh case?

[cgadal]

No worries.

It is indeed possible ! Computing the structured cartesian grid using gmsh still gets wavy (although a bit less quickly and strongly than with the unstructured grid).

Code for generating structured grid with gmsh (click to expand)

import Gmsh: gmsh

xmin, xmax = -1, 1
ymin, ymax = -1, 1
n_width = 60
typical_meshsize = (xmax - xmin) / n_width

gmsh.initialize()
gmsh.model.add("domain")

bot_left = gmsh.model.geo.addPoint(xmin, ymin, 0, typical_meshsize)
bot_right = gmsh.model.geo.addPoint(xmax, ymin, 0, typical_meshsize)
top_right = gmsh.model.geo.addPoint(xmax, ymax, 0, typical_meshsize)
top_left = gmsh.model.geo.addPoint(xmin, ymax, 0, typical_meshsize)

left_line = gmsh.model.geo.addLine(bot_left, top_left)
top_line = gmsh.model.geo.addLine(top_left, top_right)
right_line = gmsh.model.geo.addLine(top_right, bot_right)
bot_line = gmsh.model.geo.addLine(bot_right, bot_left)

gmsh.model.geo.removeAllDuplicates()
gmsh.model.geo.addCurveLoop([left_line, top_line, right_line, bot_line], 1)
surface = gmsh.model.geo.addPlaneSurface([1], 1)

# #### for structured mesh
gmsh.model.geo.mesh.setTransfiniteCurve(left_line, n_width)
gmsh.model.geo.mesh.setTransfiniteCurve(right_line, n_width)
gmsh.model.geo.mesh.setTransfiniteCurve(top_line, n_width)
gmsh.model.geo.mesh.setTransfiniteCurve(bot_line, n_width)

# gmsh.model.geo.mesh.setTransfiniteSurface(surface, "Left", [1, 2, 3, 4])
gmsh.model.geo.mesh.setTransfiniteSurface(surface, "Left")

gmsh.model.geo.synchronize()
gmsh.model.mesh.setRecombine(2, 1)  # Recombine triangles into quadrilaterals

gmsh.model.addPhysicalGroup(1, [left_line, right_line, bot_line, top_line], 1, "ExtBoundary")

gmsh.option.setNumber("Mesh.SaveAll", 1)
gmsh.option.setNumber("Mesh.SaveGroupsOfElements", 1)
gmsh.option.setNumber("Mesh.SaveGroupsOfNodes", 1)  # Important for p4est and Trixi

gmsh.model.mesh.generate(2)

gmsh.fltk.run()
gmsh.write("mesh_test.inp")
gmsh.finalize()

For comparison, the result when building the cartesian mesh directly from p4est is the following, as expected:

EDIT:
I'm however suspicious about the properties of the two grids being similar.

• cartesian grid from P4estMesh

julia> mesh.boundary_names
4×1 Matrix{Symbol}:
 :x_neg
 :x_pos
 :y_neg
 :y_pos

julia> mesh.nodes
4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 -1.0
 -0.4472135954999579
  0.4472135954999579
  1.0

• cartesian mesh from gmsh:

julia> mesh.boundary_names
4×2401 Matrix{Symbol}:
 :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   …  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")
 Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")     :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary   :ExtBoundary
 :ExtBoundary   Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")     Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")
 Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")     Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  :ExtBoundary

julia> mesh.nodes
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 -1.0
  1.0

[jlchan]

That's quite odd; you'd expect them to be the same if things were processed correctly. I'm not sure why this would be the case, however.

Naive question - is the BC actually being used in both cases? E.g., if you put a print statement out inside is it being called?

[cgadal]

That's quite odd; you'd expect them to be the same if things were processed correctly. I'm not sure why this would be the case, however.

Agreed ... See my edit on my previous answer, I might have done something wrong with the meshes ...

Naive question - is the BC actually being used in both cases? E.g., if you put a print statement out inside is it being called?

Yes, I just checked adding a print statement and it's called in both cases at each timestep.

[jlchan]

Your edits on your previous answer do look a bit odd. I can't say for sure, but it almost looks like :ExtBoundary is being set in the interior of the mesh...which might explain the odd wave behavior.

[cgadal]

I just loaded the foil mesh from https://github.com/trixi-framework/Trixi.jl/blob/main/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl

and it gives pretty similar results to the mesh I generated using gmsh:

julia> mesh.boundary_names
4×876 Matrix{Symbol}:
 Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  …  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")
 Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")     Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")
 Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")     Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")
 Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")     Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")  Symbol("---")

julia> mesh.nodes
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 -1.0
  1.0

The only difference is that the actual boundary names :PhysicalLine1, :PhysicalLine2, :PhysicalLine3, :PhysicalLine4 do not appear anymore in mesh.boundary_names.

[ranocha]

@DanielDoehring You have worked with this type of setup, haven't you?

[DanielDoehring]

Are you supplying the boundary symbols to the mesh constructor, i.e.,

Trixi.jl/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl

Lines 74 to 76 in 1e1f643

	boundary_symbols = [:PhysicalLine1, :PhysicalLine2, :PhysicalLine3, :PhysicalLine4]
	mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols)

?

[cgadal]

Yes, as

boundary_symbols = [:ExtBoundary]
mesh = TR.P4estMesh{2}(mesh_file, boundary_symbols=boundary_symbols, polydeg=polydeg)

I just tested to actually run my code using the .inp file from the foil from the example you linked, and the code runs as expected:

So the problem comes from my .inp files or how I load them. I compared the .inp files I generated myself to the one of the foil in the example, and they look very similar ...

[DanielDoehring]

So the mesh you generate should be fine. One thing I can imagine is that your BCs assume something on the mesh which is true for the p4est generated one but not for the imported one. You could check this by using one of the tested, known to work BCs.

[cgadal]

Okay, so I figured it out. I was basically defining the outer domain by following the corners clockwise instead of anticlockwise/trigo.

It looks like if you start from

 // outer bounding box
 Point(1) = {-1, -1, 0, 1.0};
 Point(2) = {1, -1, 0, 1.0};
 Point(3) = {1, 1, 0, 1.0};
 Point(4) = {-1, 1, 0, 1.0};

Then this would be right:

 // lines of the bounding box
 Line(1) = {1, 2};
 Line(2) = {2, 3};
 Line(3) = {3, 4};
 Line(4) = {4, 1};

but this would be wrong:

 // lines of the bounding box
 Line(1) = {1, 4};
 Line(2) = {4, 3};
 Line(3) = {3, 2};
 Line(4) = {2, 1};

It does not prevent the generation of the mesh by gmsh, and the nodes are simply ordered differently in the file. Trixi or `p4est may rely on the ordering of the nodes to determine what is the "inward" direction of the boundary ?

[jlchan]

Ah, that makes sense. I'm not sure of the P4estMesh specifics but the mesh orientation typically matters. If the nodes were flipped then it might be more than just the "inward" direction that got reversed, some geometric scaling terms might also have been negated.

[DanielDoehring]

Okay, so I figured it out. I was basically defining the outer domain by following the corners clockwise instead of anticlockwise/trigo.

Yeah that makes sense, p4est assumes that things are ordered right-handed /anti clock-wise.

[ranocha]

Can we detect this when reading in the mesh?