implementing an equation

[rveltz]

Dear developpers,

I have my own FVM method for solving the eqaution below but it is slow (~10us per time step) because a lot of space is wasted as the solution is quite localized. Hence, I was hoping to use mesh adaptation in Trixi. BUT:

I am intimidated by the framework. Basically, I want to simulate

$$\partial_t \mu +\partial_x(F^1(x,y)\mu)+ \partial_x(F^2(x,y)\mu) + \partial_x^2\mu=0$$

Looking at the docs, I understand that I could use https://trixi-framework.github.io/Trixi.jl/stable/tutorials/adding_nonconservative_equation/

but this is written in non-conservative form. I guess my question is "Is it worth it?" Would it be faster than my in-house code (explicit euler with upwind fvm)?

[jlchan]

I don't think that's a non-conservative from? That looks like conservative form with spatially varying coefficients and a parabolic (e.g., 2nd order derivative) term.

Just checking, is the sign on $\partial_x\mu$ correct? Usually I see a negative sign on that term.

[rveltz]

yes it is conservative, I mentioned the link because the vector field is not constant There is no sign - because mu is a measure but we can put the minus if this measure has a density.

Is this an equation handled by trixi?

[jlchan]

Thanks for the clarification - yes, I think you'd need to follow the approach in the docs for spatially varying coefficients, which is to just include those as an extra set of unknowns. However, this isn't too difficult - you'd just create a set of equations with both mu and F as unknowns, and initialize them both in the initial_condition function. When you define the flux, you would just set the flux for F to zero so that it's not evolved.

In terms of speed, it probably depends on what your solution looks like. If it's smooth and differentiable, yes, I think Trixi should be faster than a 1st or 2nd order FVM for higher levels of accuracy. If it has a lot of shocks or discontinuities, I'm not as sure.

[rveltz]

no it is smooth

[jlchan]

Great, then I think it should do fairly well. Both using a higher order timestepper and using a higher order spatial discretization should improve accuracy.

If you plan on implementing this using Trixi, please don't hesitate to contact us if we can help clarify anything about the interface.

[rveltz]

ok thanks, I will give it a try

[rveltz]

This is tough...

[rveltz]

I dont know how to write the function flux_conservative. If you can give me a push, it'd be lovely

using Trixi
using Trixi: AbstractEquations, get_node_vars
import Trixi: varnames, default_analysis_integrals, flux, max_abs_speed_naive,
              have_nonconservative_terms

# Since there is no native support for variable coefficients, we use 3
# variables: one for the basic unknown `u` and 2 for the coefficient `a`
struct IsolatedNeuron <: AbstractEquations{2 #= spatial dimension =#,
                                           3 #= 3 variables (u,F1,F2) =#}
end

varnames(::typeof(cons2cons), ::IsolatedNeuron) = ("scalar", "F1", "F2")

default_analysis_integrals(::IsolatedNeuron) = ()

# The conservative part of the flux is zero
flux(u, orientation, equation::IsolatedNeuron) = zero(u)

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, ::IsolatedNeuron)
    _, F1_ll, F2_ll = u_ll
    _, F1_rr, F2_rr = u_rr
    return max(abs(F1_ll), abs(F2_ll), abs(F1_rr), abs(F1_rr))
end

# We use conservative terms
have_nonconservative_terms(::IsolatedNeuron) = Trixi.False()

function flux_conservative(u_mine, u_other, orientation,
                              equations::IsolatedNeuron)
    _, F1, F2 = u_mine
    scalar, _ = u_other
    return SVector(F1 * scalar + F2 * scalar, zero(scalar), zero(scalar))
end

[jlchan]

I don't think this is a non-conservative equation, so you should only need to construct flux, which I believe should look like

function flux(u, orientation, equations::IsolatedNeuron)
    scalar, F1, F2 = u
    if orientation == 1
        return SVector(F1 * scalar, zero(scalar), zero(scalar))
    else # if orientation == 2
        return SVector(F2 * scalar, zero(scalar), zero(scalar))
end

[jlchan]

Also, have_nonconservative_terms defaults to False() so you don't need to include that.

[jlchan]

I should note that this is enough for a standard nodal DG discretization. This is usually enough for a scalar transport equation like yours, but if you observe any instability, there are "flux differencing" discretizations which address this.

[rveltz]

Thank you!

I now have the following but it complains LoadError: UndefVarError: solve

using Trixi
using Trixi: AbstractEquations, get_node_vars
import Trixi: varnames, default_analysis_integrals, flux, max_abs_speed_naive,
              have_nonconservative_terms

# Since there is no native support for variable coefficients, we use 3
# variables: one for the basic unknown `u` and 2 for the coefficient `a`
struct IsolatedNeuron <: AbstractEquations{2 #= spatial dimension =#,
                                           3 #= 3 variables (u,F1,F2) =#}
end

varnames(::typeof(cons2cons), ::IsolatedNeuron) = ("scalar", "F1", "F2")

default_analysis_integrals(::IsolatedNeuron) = ()

# The conservative part of the flux is zero
flux(u, orientation, equation::IsolatedNeuron) = zero(u)

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, ::IsolatedNeuron)
    _, F1_ll, F2_ll = u_ll
    _, F1_rr, F2_rr = u_rr
    return max(abs(F1_ll), abs(F2_ll), abs(F1_rr), abs(F1_rr))
end

# We use conservative terms
have_nonconservative_terms(::IsolatedNeuron) = Trixi.False()

function flux(u, orientation, equations::IsolatedNeuron)
    scalar, F1, F2 = u
    if orientation == 1
        return SVector(F1 * scalar, zero(scalar), zero(scalar))
    else # if orientation == 2
        return SVector(F2 * scalar, zero(scalar), zero(scalar))
    end
end

using OrdinaryDiffEq

equation = IsolatedNeuron()

function initial_condition(x, t, equation::IsolatedNeuron)
    scalar = exp( -((x[1]-1)^2+x[2]^2) / (2*0.3^2))
    F1 = -x[1]
    F2 = -x[2]
    SVector(scalar, F1, F2)
end

# Create a uniform mesh in 1D in the interval [-π, π] with periodic boundaries
mesh = TreeMesh((-2.,-2.), (2., 2.), # min/max coordinates
                initial_refinement_level=4, n_cells_max=10^4)

# Create a DGSEM solver with polynomials of degree `polydeg`
# Remember to pass a tuple of the form `(conservative_flux, nonconservative_flux)`
# as `surface_flux` and `volume_flux` when working with nonconservative terms
volume_flux  = (flux_central, flux)
surface_flux = (flux_lax_friedrichs, flux)
solver = DGSEM(polydeg=3, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# Setup the spatial semidiscretization containing all ingredients
semi = SemidiscretizationHyperbolic(mesh, equation, initial_condition, solver)

# Create an ODE problem with given time span
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

# Set up some standard callbacks summarizing the simulation setup and computing
# errors of the numerical solution
summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi, interval=50)
callbacks = CallbackSet(summary_callback, analysis_callback)

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes
# the passed callbacks
sol = solve(ode, Tsit5(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks)

[jlchan]

Gotcha - you will need to load OrdinaryDiffEq.jl to run this.

Since you aren't changing Trixi internals, you could just create a new run project directory with Trixi.jl and OrdinaryDiffEq.jl, define your new equations in this driver file, and run from within this run directory.

[rveltz]

but OrdinaryDiffEq is in the code above....

[jlchan]

Oh, I see it now, I missed that at first.

That is strange. Can you post the full stacktrace?

[rveltz]

I fixed it by importing it before Trixi.

[rveltz]

The solve errors

1-element ExceptionStack:
LoadError: MethodError: objects of type Tuple{typeof(flux_central), typeof(flux)} are not callable
Stacktrace:
  [1] flux_differencing_kernel!
    @ ~/.julia/packages/Trixi/wpKZk/src/solvers/dgsem_tree/dg_2d.jl:262 [inlined]
  [2] flux_differencing_kernel!
    @ ~/.julia/packages/Trixi/wpKZk/src/solvers/dgsem_tree/dg_2d.jl:248 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/Trixi/wpKZk/src/solvers/dgsem_tree/dg_2d.jl:236 [inlined]
  [4] #565
    @ ~/.julia/packages/Polyester/W70Gq/src/closure.jl:274 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/Polyester/W70Gq/src/batch.jl:159 [inlined]
  [6] _batch_no_reserve
    @ ~/.julia/packages/Polyester/W70Gq/src/batch.jl:112 [inlined]
  [7] batch
    @ ~/.julia/packages/Polyester/W70Gq/src/batch.jl:231 [inlined]
  [8] macro expansion
    @ ~/.julia/packages/Polyester/W70Gq/src/closure.jl:387 [inlined]
  [9] macro expansion
    @ ~/.julia/packages/Trixi/wpKZk/src/auxiliary/auxiliary.jl:246 [inlined]
 [10] calc_volume_integral!(du::StrideArraysCore.PtrArray{Float64, 4, (1, 2, 3,

[JoshuaLampert]

Try not to pass a tuple, but only

volume_flux  = flux_central
surface_flux = flux_lax_friedrichs

You need a Tuple for nonconservative equations. Since your equation is conservative only one flux is needed.

[rveltz]

Great! it returns something. I need to plot and add the laplacian wrt x Ill try the plot first

[jlchan]

Great. The Laplacian can probably just use LaplaceDiffusion2D though maybe the sign is wrong.

[rveltz]

OK, I think I got it. But the mass is not conserved and the density becomes negative :(

using Trixi: AbstractEquations, get_node_vars
import Trixi: varnames, default_analysis_integrals, flux, max_abs_speed_naive,
              have_nonconservative_terms

# Since there is no native support for variable coefficients, we use 3
# variables: one for the basic unknown `u` and 2 for the coefficient `a`
struct IsolatedNeuron <: AbstractEquations{2 #= spatial dimension =#,
                                           3 #= 3 variables (u,F1,F2) =#}
end

varnames(::typeof(cons2cons), ::IsolatedNeuron) = ("scalar", "F1", "F2")

default_analysis_integrals(::IsolatedNeuron) = ()

# The conservative part of the flux is zero
flux(u, orientation, equation::IsolatedNeuron) = zero(u)

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, ::IsolatedNeuron)
    _, F1_ll, F2_ll = u_ll
    _, F1_rr, F2_rr = u_rr
    return max(abs(F1_ll), abs(F2_ll), abs(F1_rr), abs(F1_rr))
end

# We use conservative terms
have_nonconservative_terms(::IsolatedNeuron) = Trixi.False()

function flux(u, orientation, equations::IsolatedNeuron)
    scalar, F1, F2 = u
    if orientation == 1
        return SVector(F1 * scalar, zero(scalar), zero(scalar))
    else # if orientation == 2
        return SVector(F2 * scalar, zero(scalar), zero(scalar))
    end
end

function initial_condition(x, t, equation::IsolatedNeuron)
    scalar = exp( -((x[1]-1)^2+x[2]^2) / (2*.2^2))
    F1 = -x[1]
    F2 = -x[2]
    SVector(scalar, F1, F2)
end

equation_transport = IsolatedNeuron()
equations = LaplaceDiffusion2D((1.,1.), equation_transport);

# Create a uniform mesh in 1D in the interval [-π, π] with periodic boundaries
mesh = TreeMesh((-4.,-4.), (4., 4.), # min/max coordinates
                initial_refinement_level=4, n_cells_max=10^4)

# Create a DGSEM solver with polynomials of degree `polydeg`
# Remember to pass a tuple of the form `(conservative_flux, nonconservative_flux)`
# as `surface_flux` and `volume_flux` when working with nonconservative terms
volume_flux  = flux_central
surface_flux = flux_lax_friedrichs
solver = DGSEM(polydeg=3, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# Setup the spatial semidiscretization containing all ingredients
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

# Create an ODE problem with given time span
tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

# Set up some standard callbacks summarizing the simulation setup and computing
# errors of the numerical solution
summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi, interval=50)
callbacks = CallbackSet(summary_callback, analysis_callback)

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes
# the passed callbacks
sol = solve(ode, Tsit5(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks)

plot(sol)

pd = PlotData2D(sol)

pd = PlotData2D(sol(0.), semi);plot(pd["scalar"], color = :viridis, title = "t=0")

pd = PlotData2D(sol(2.), semi);plot(pd["scalar"], color = :viridis, title = "t=2")

[rveltz]

[rveltz]

[rveltz]

Given the vector field, the density should converge to a dirac at (x,y)=0 without the Laplacian. The Laplacian prevents this so I am expecting an invariant distribution around 0

[rveltz]

I have an issue with the above code. When adding the laplacian, what do I declare for the initial condition?

equation_transport = IsolatedNeuron()
equations = LaplaceDiffusion2D((1.,1.), equation_transport);

Should I just change into ?

function initial_condition(x, t, equation::LaplaceDiffusion2D{IsolatedNeuron, 3, Tuple{Float64, Float64}})

[jlchan]

Good question, at the moment we require that you do not add a type annotation for equations in the initial condition when using parabolic terms

[rveltz]

OK thank you. solve now complain that

ERROR: MethodError: no method matching flux(::SVector{3, Float64}, ::Int64, ::LaplaceDiffusion2D{IsolatedNeuron, 3, Tuple{Float64, Float64}})

Closest candidates are:
  flux(::Any, ::Any, ::Integer, ::LaplaceDiffusion1D)
   @ Trixi ~/.julia/packages/Trixi/wpKZk/src/equations/laplace_diffusion_1d.jl:21
  flux(::Any, ::Any, ::Integer, ::LaplaceDiffusion2D)
   @ Trixi ~/.julia/packages/Trixi/wpKZk/src/equations/laplace_diffusion_2d.jl:22
  flux(::Any, ::Any, ::Integer, ::CompressibleNavierStokesDiffusion1D)
   @ Trixi ~/.julia/packages/Trixi/wpKZk/src/equations/compressible_navier_stokes_1d.jl:152
  ...

[rveltz]

The current version of the code being:

using OrdinaryDiffEq, Trixi, Plots

using Trixi: AbstractEquations, get_node_vars
import Trixi: varnames, default_analysis_integrals, flux, max_abs_speed_naive,
              have_nonconservative_terms

# Since there is no native support for variable coefficients, we use 3
# variables: one for the basic unknown `u` and 2 for the coefficient `a`
struct IsolatedNeuron <: AbstractEquations{2 #= spatial dimension =#,
                                           3 #= 3 variables (u,F1,F2) =#}
end

varnames(::typeof(cons2cons), ::IsolatedNeuron) = ("scalar", "F1", "F2")

default_analysis_integrals(::IsolatedNeuron) = ()

# The conservative part of the flux is zero
flux(u, orientation, equation::IsolatedNeuron) = zero(u)

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, ::IsolatedNeuron)
    _, F1_ll, F2_ll = u_ll
    _, F1_rr, F2_rr = u_rr
    return max(abs(F1_ll), abs(F2_ll), abs(F1_rr), abs(F1_rr))
end

# We use conservative terms
have_nonconservative_terms(::IsolatedNeuron) = Trixi.False()

function flux(u, orientation, equations::IsolatedNeuron)
    scalar, F1, F2 = u
    if orientation == 1
        return SVector(F1 * scalar, zero(scalar), zero(scalar))
    else # if orientation == 2
        return SVector(F2 * scalar, zero(scalar), zero(scalar))
    end
end

function initial_condition(x, t, equation)
    scalar = exp( -((x[1]-1.2)^2+(x[2]-4)^2) / (2*.2^2))
    # scalar = exp(-(x[1]-v0b)^2/(2 * σv0^2) -(x[1]-w0b)^2/(2 * σw0^2) )
    F1 = -x[1]
    F2 = -x[2]
    SVector(scalar, F1, F2)
end

equation_transport = IsolatedNeuron()
equations = LaplaceDiffusion2D((1.,1.), equation_transport);

# Create a uniform mesh in 1D in the interval [-π, π] with periodic boundaries
mesh = TreeMesh((-2.5,-2.), (5., 12.), # min/max coordinates
                initial_refinement_level=4, n_cells_max=10^5)

# Create a DGSEM solver with polynomials of degree `polydeg`
# Remember to pass a tuple of the form `(conservative_flux, nonconservative_flux)`
# as `surface_flux` and `volume_flux` when working with nonconservative terms
volume_flux  = flux_central
surface_flux = flux_lax_friedrichs
solver = DGSEM(polydeg=3, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# Setup the spatial semidiscretization containing all ingredients
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

# Create an ODE problem with given time span
tspan = (0.0, 4.0)
ode = semidiscretize(semi, tspan)

# Set up some standard callbacks summarizing the simulation setup and computing
# errors of the numerical solution
summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi, interval=50)
callbacks = CallbackSet(summary_callback, analysis_callback)

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes
# the passed callbacks
sol = solve(ode, Tsit5(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks)

[jlchan]

Taking a look...

[JoshuaLampert]

Have you had a closer look at this elixir or a similar one? Maybe this helps as an inspiration.

[jlchan]

Yeah, the elixir @JoshuaLampert linked would help. The main issue is that you are still using SemidiscretizationHyperbolic instead of SemidiscretizationHyperbolicParabolic https://github.com/trixi-framework/Trixi.jl/blob/dfdaef5adc243fcd6c55262a69ede59fe3a32a8f/examples/tree_2d_dgsem/elixir_advection_diffusion.jl#L45-L49

[JoshuaLampert]

Exactly. Then you can also just keep function initial_condition(x, t, equation::IsolatedNeuron), I guess.

[jlchan]

I think you still need to remove the annotation from initial_condition because it will be called with both equation and equation_parabolic in different parts of the code (like callbacks)

[jlchan]

Wait, I'm wrong - sorry @JoshuaLampert, you're correct, you can keep the type annotation in that case then.

[rveltz]

Thank you for the elixr. It helped me with the boundary condition. I still have an error from solve. These are not easy to debug

using OrdinaryDiffEq, Trixi, Plots
# using AbbreviatedStackTraces

# using FHNMF
# neuron = FHNMF.NeuronParam2d(a = 0.01, b = 0.1, Iext = 4.0, J = 0.05)
# mesh_n = MeshPDE2d(Nw=300, param_neuron = neuron)

using Trixi: AbstractEquations, get_node_vars
import Trixi: varnames, default_analysis_integrals, flux, max_abs_speed_naive,
              have_nonconservative_terms

# Since there is no native support for variable coefficients, we use 3
# variables: one for the basic unknown `u` and 2 for the coefficient `a`
struct IsolatedNeuron <: AbstractEquations{2 #= spatial dimension =#,
                                           3 #= 3 variables (u,F1,F2) =#}
end

varnames(::typeof(cons2cons), ::IsolatedNeuron) = ("scalar", "F1", "F2")

default_analysis_integrals(::IsolatedNeuron) = ()

# The conservative part of the flux is zero
flux(u, orientation, equation::IsolatedNeuron) = zero(u)

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, ::IsolatedNeuron)
    _, F1_ll, F2_ll = u_ll
    _, F1_rr, F2_rr = u_rr
    return max(abs(F1_ll), abs(F2_ll), abs(F1_rr), abs(F1_rr))
end

# We use conservative terms
have_nonconservative_terms(::IsolatedNeuron) = Trixi.False()

function flux(u, orientation, equations::IsolatedNeuron)
    scalar, F1, F2 = u
    if orientation == 1
        return SVector(F1 * scalar, zero(scalar), zero(scalar))
    else # if orientation == 2
        return SVector(F2 * scalar, zero(scalar), zero(scalar))
    end
end

function initial_condition(x, t, equation::IsolatedNeuron)
    scalar = exp( -((x[1]-1.2)^2+(x[2]-4)^2) / (2*.2^2))
    F1 = -x[1]
    F2 = -x[2]
    # F1 = FHNMF.f(x[1], x[2], neuron)
    # F2 = FHNMF.f(x[1], x[2], neuron)
    SVector(scalar, F1, F2)
end

equations_transport = IsolatedNeuron()
equations_parabolic = LaplaceDiffusion2D((1.,1.), equations_transport);

# Create a uniform mesh in 1D in the interval [-π, π] with periodic boundaries
mesh = TreeMesh((-2.5,-2.), (5., 12.), # min/max coordinates
                initial_refinement_level=4, n_cells_max=10^5)

# Create a DGSEM solver with polynomials of degree `polydeg`
# Remember to pass a tuple of the form `(conservative_flux, nonconservative_flux)`
# as `surface_flux` and `volume_flux` when working with nonconservative terms
volume_flux  = flux_central
surface_flux = flux_lax_friedrichs
solver = DGSEM(polydeg=3, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# Setup the spatial semidiscretization containing all ingredients
# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))

boundary_conditions_hyperbolic = (; x_neg = boundary_condition_do_nothing,
                                    y_neg = boundary_condition_do_nothing,
                                    y_pos = boundary_condition_do_nothing,
                                    x_pos = boundary_condition_do_nothing)

boundary_conditions_parabolic = (; x_neg = boundary_condition_zero_dirichlet,
                                   y_neg = boundary_condition_zero_dirichlet,
                                   y_pos = boundary_condition_zero_dirichlet,
                                   x_pos = boundary_condition_zero_dirichlet);

semi = SemidiscretizationHyperbolicParabolic(mesh,
                                             (equations_transport, equations_parabolic),
                                             initial_condition, solver;
                                             boundary_conditions=(boundary_conditions_hyperbolic,
                                                                  boundary_conditions_parabolic))

# Create an ODE problem with given time span
tspan = (0.0, 4.0)
ode = semidiscretize(semi, tspan)

# Set up some standard callbacks summarizing the simulation setup and computing
# errors of the numerical solution
summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi, interval=50)
callbacks = CallbackSet(summary_callback, analysis_callback)

sol = solve(ode, Tsit5(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks)

The error being:

ERROR: MethodError: no method matching *(::Tuple{Float64, Float64}, ::SVector{3, Float64})

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:578
  *(::StaticArraysCore.StaticArray{Tuple{N, M}, T, 2} where {N, M, T}, ::StaticArraysCore.StaticArray{Tuple{N}, T, 1} where {N, T})
   @ StaticArrays ~/.julia/packages/StaticArrays/cZ1ET/src/matrix_multiply.jl:10
  *(::Union{StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}, LinearAlgebra.Adjoint{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}, LinearAlgebra.Diagonal{T, <:StaticArraysCore.StaticArray{Tuple{s1}, T, 1}}, LinearAlgebra.Hermitian{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}, LinearAlgebra.LowerTriangular{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}, LinearAlgebra.Symmetric{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}, LinearAlgebra.Transpose{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}, LinearAlgebra.UnitLowerTriangular{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}, LinearAlgebra.UnitUpperTriangular{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}, LinearAlgebra.UpperTriangular{T, <:StaticArraysCore.StaticArray{Tuple{s1, s2}, T, 2}}} where {s1, s2, T}, ::StaticArraysCore.StaticArray{Tuple{N}, T, 1} where {N, T})
   @ StaticArrays ~/

[jlchan]

Would you mind posting the unabbreviated stacktrace?

[JoshuaLampert]

I am not sure if LaplaceDiffusion2D accepts a Tuple for the diffusivity. Looks like that could be the problem at first glance.

[jlchan]

One possible source of error:

boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))

should probably be

boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0, 0.0, 0.0))

since you have 3 "variables" in your equation

[JoshuaLampert]

Looks like it does not...

[rveltz]

Is it easy to add just a diffusion on the x variable?

[jlchan]

Yes, if you don't mind a hacky way, you could just redefine https://github.com/trixi-framework/Trixi.jl/blob/dfdaef5adc243fcd6c55262a69ede59fe3a32a8f/src/equations/laplace_diffusion_2d.jl#L22-L29

so that it returns zero if orientation == 2

[rveltz]

Judging from the opened issues, it seems I cannot use AMR with diffusion, is it true or not?

[jlchan]

Unfortunately, yes, though we are very close to merging the PR for that https://github.com/trixi-framework/Trixi.jl/pull/1629

[rveltz]

Arg, I was not sure I was doing something wrong.

[rveltz]

I have one more question if I may. That the equation is in fact

$$\partial_t \mu +\partial_x([F^1(x,y) + \int g\mu]\ \mu)+ \partial_x(F^2(x,y)\mu) + \partial_x^2\mu=0$$

you see that the vector field is modified by the integral of $g\mu$ for some known g. Would it be possible to tackle this with a callback for example? Or maybe adding a nonlinearity

[ranocha]

Trixi.jl is not really designed to include integrals in the fluxes. You could try to use something similar to our tutorial https://trixi-framework.github.io/Trixi.jl/stable/tutorials/custom_semidiscretization/ to compute the integral first and use it inside the flux afterwards. You can use integrate to compute the integral numerically (see https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#PolynomialBases.integrate-Union{Tuple{Func},%20Tuple{Func,%20Any,%20Trixi.AbstractSemidiscretization}}%20where%20Func).

[JoshuaLampert]

Can this issue be closed or are there any open questions, @rveltz?

[rveltz]

Yes we can close. It was really helpful.