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

Hi all,

I'm moving here [the discussion](https://discourse.julialang.org/t/solving-2d-pde-similar-to-advection-diffusion-on-unstructured-mesh-with-trixi/125462/5) 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:

```julia
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_](https://www.sciencedirect.com/science/article/pii/S0021999100964593). 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?

Start with the custom SegregationEquation entry points, including TR.flux, TR.max_abs_speed_naive, TR.max_abs_speeds, boundary conditions, and SemidiscretizationHyperbolic, then compare them with the tutorials and examples referenced in the issue. Done means determining how the 1D equation, shock treatment, boundary condition, and proposed Kurganov–Tadmor flux map onto Trixi before extending the workflow to diffusion, Robin boundaries, and an unstructured 2D mesh.

Written by the indexing model from the issue text.