StructuredMesh parabolic equations supported Trixi

Specific Domains
1

Hi, i was using Trixi.jl for the Brown and Minion test case and i was trying to use the StructuredMesh function to refine the mesh cells by cells (20x20 to 25x25 and so on) but it seems that the create_cache_parabolic function is not defined for a StructuredMesh type

ERROR: LoadError: MethodError: no method matching create_cache_parabolic(::StructuredMesh{2, Float64}, ::CompressibleEulerEquations2D{Float64}, ::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive, Float64, Float64, CompressibleEulerEquations2D{Float64}}, ::DGSEM{LobattoLegendreBasis{Float64, 5, SVector{5, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 5, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{typeof(flux_hllc)}, VolumeIntegralWeakForm}, ::ViscousFormulationBassiRebay1, ::Type{Float64}, ::Type{Float64})

Closest candidates are:
  create_cache_parabolic(!Matched::P4estMesh{3}, ::Trixi.AbstractEquations, ::Trixi.AbstractEquationsParabolic, ::DG, ::Any, ::Any, ::Any)
   @ Trixi /scratch/stg-cfds/renoux/Trixi.jl/src/solvers/dgsem_p4est/dg_3d_parabolic.jl:11
  create_cache_parabolic(!Matched::P4estMesh{2}, ::Trixi.AbstractEquations, ::Trixi.AbstractEquationsParabolic, ::DG, ::Any, ::Any, ::Any)
   @ Trixi /scratch/stg-cfds/renoux/Trixi.jl/src/solvers/dgsem_p4est/dg_2d_parabolic.jl:11
  create_cache_parabolic(!Matched::TreeMesh{3, TreeType} where TreeType<:Trixi.AbstractTree{3}, ::Trixi.AbstractEquations, ::Trixi.AbstractEquationsParabolic, ::DG, ::Any, ::Any, ::Any)
   @ Trixi /scratch/stg-cfds/renoux/Trixi.jl/src/solvers/dgsem_tree/dg_3d_parabolic.jl:986

On the Trixi.jl documentation, i saw that CompressibleNavierStokesEquations() should be supported for StructuredMesh

Capture d’écran 2024-07-19 à 13.57.28
Capture d’écran 2024-07-19 à 13.57.28425×699 35.8 KB

So i just wanted to know if it was a mistake on my part or just a future implementation ?

2

Here is the following julia code

using Trixi, OrdinaryDiffEq

prandtl_number() = 0.72
mu = 0.00005

dg = DGSEM(polydeg = 4, surface_flux = flux_hllc)

equations = CompressibleEulerEquations2D(1.4)
equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number(), gradient_variables = GradientVariablesPrimitive())
"""
A compressible version of the double shear layer initial condition. Adapted from
Brown and Minion (1995). See Section 3, equations (27)-(28) for the original
incompressible version.

- David L. Brown and Michael L. Minion (1995)
  Performance of Under-resolved Two-Dimensional Incompressible Flow Simulations.
  [DOI: 10.1006/jcph.1995.1205](https://doi.org/10.1006/jcph.1995.1205)
"""
function initial_condition_BM_vortex(x, t, equations::CompressibleEulerEquations2D)

  L = 1.0
  x0 = 0.0
  Ma0 = 0.2        # Maximum Mach number
  rho0 = 1.2         # density outside the shear layer
  
  n_vortex=1
  l_ratio = 100.0      # length scale ratio (30.0 default)
  d_ratio = 1.0       # density ratio (4.0 default)
  v_ratio = 0.05       # velocity ratio (0.2 default)
  shear = 200.0        # shear value (15.0 default)

  y0 = L/2
  ya = y0 - L/4
  yb = y0 + L/4
  sigma = L / l_ratio
  u0 = shear*sigma/2
  v0 = v_ratio * u0
  rho1 = d_ratio * rho0

  B = tanh((x[2]-ya)/sigma)/2 - tanh((x[2]-yb)/sigma)/2

  rho = rho0 * (1.0 + B*(d_ratio-1.0))
  u = u0 * (2*B - 1.0)
  v = v0 * cos(2.0 * pi * n_vortex * (x[1]-x0)/L)# + u0

  c0 = u0 / Ma0
  temp0 = c0^2 / (equations.gamma*8.314)
  p0 = rho1 * 8.314 * temp0

    return prim2cons(SVector(rho, u, v, p0), equations)
end

initial_condition = initial_condition_BM_vortex

coordinates_min = (0.0,0.0) # minimum coordinates (min(x), min(y))
coordinates_max = (1.0,1.0) # maximum coordinates (max(x), max(y))

trees_per_dimension = (1,1)
nb_cells = 67
#mesh = P4estMesh(trees_per_dimension,
                 #polydeg = 1, initial_refinement_level = 6,
                 #coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                 #periodicity = true)

mesh = StructuredMesh((nb_cells, nb_cells), coordinates_min, coordinates_max, periodicity = true)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations,equations_parabolic), initial_condition, dg)

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

alive_callback = AliveCallback(alive_interval = 1000)

save_solution = SaveSolutionCallback(interval = 200,
                                     save_initial_solution = true,
                                     save_final_solution = true,
                                     solution_variables = cons2prim, output_directory="out_shearlayer")

analysis_interval = 100

analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

stepsize_callback = StepsizeCallback(cfl = 0.1)

callbacks = CallbackSet(summary_callback, save_solution, analysis_callback, stepsize_callback)

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

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

summary_callback() # print the timer summary
3

Yes, you are right. The StructuredMesh does not support parabolic terms yet. This was an error in our documentation, which will be fixed soon. Thanks for pointing that out!

1 Like
4

Thanks for the report. @JoshuaLampert has already prepared a hotfix for the docs. We should update the code to work with StructuredMeshes, too. In the meantime, you should be able to use the P4estMesh.

3 Likes
5

Thanks ! Glad i could help. I also had another question about the fonctionnalities for the DGmulti solver, I have tried to get an output in the format .h5 but the SaveSolutionCallback returns an error. Is it implemented in Trixi.jl for the DGmulti solver ?

6

No, the SaveSolutionCallback is not available for DGMulti. But you can plot numerical solutions using Plots.jl or Makie.jl (in an interactive session).

1 Like
7

Thanks, I am trying to get the data from Julia to Python for learning purposes (Through CSV.jl), I was wondering how you manage the discontinuities in the plot of the solution ? Because we have (P+1)^(NDIM) points in a cell and the faces that are not at the boundaries share points with theirs neighbors.
Thanks again for the quick response !

8

It depends on the setup. With Trixi2Vtk.jl and similar tools, we can interpolate solutions to a uniform grid of nodes strictly inside the elements for plotting. For other purposes (in particular in 1D), it can be very helpful to move the nodes at the element interfaces a tiny bit inside the element (like nextfloat and prevfloat) to really show the discontinuous solutions.

9

Oh I see, to visualise the riemann problem at the faces of elements, thanks !