Skip to content

DGMulti CGSEM example with nonperiodic boundaries #2599

@ranocha

Description

@ranocha

It would be nice to have such an example. Since I was asked, here is a setup that should work:

julia> using SummationByPartsOperators, Trixi, OrdinaryDiffEqLowStorageRK, Plots

julia> equations = LinearScalarAdvectionEquation1D(1.0)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation1D                                                                  │
│ ═══════════════════════════════                                                                  │
│ #variables: ………………………………………………… 1                                                                │
│ │ variable 1: …………………………………………… scalar                                                           │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> D = couple_continuously(                                # CGSEM
               legendre_derivative_operator(-1.0, 1.0, 3 + 1), # polynomial degree 3
               UniformMesh1D(xmin = -1.0, xmax = 1.0, Nx = 8)) # 8 elements
First derivative operator {T=Float64} on 4 Lobatto Legendre nodes in [-1.0, 1.0]
coupled continuously on UniformMesh1D{Float64} with 8 cells in (-1.0, 1.0)

julia> dg = DGMulti(element_type = Line(),
                    approximation_type = D,
                    surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
                    volume_integral = VolumeIntegralWeakForm())
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64}                                                                                      │
│ ═══════════                                                                                      │
│ basis: ……………………………………………………………… RefElemData{UniformNonperiodicCoivative:1, accuracy:3), Line()} │
│ mortar: …………………………………………………………… nothing                                                          │
│ surface integral: ………………………………… SurfaceIntegralWeakForm                                          │
│ │ surface flux: ……………………………………… FluxLaxFriedrichs(max_abs_speed)                                 │
│ volume integral: …………………………………… VolumeIntegralWeakForm                                           │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> mesh = DGMultiMesh(dg, (1,), # only 1 element - otherwise, CG/DG would be mixed 
                          coordinates_min = (-1.0,), coordinates_max = (1.0,), 
                          periodicity = false)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DGMultiMesh{1, Trixi.Affine},                                                                    │
│ ══════════════════════════════                                                                   │
│ number of elements: …………………………… 1                                                                │
│ number of boundaries: ……………………… 1                                                                │
│ │ nfaces on entire_boundary: …… 2                                                                │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, dg;
           boundary_conditions = (; :entire_boundary => 
               BoundaryConditionDirichlet(initial_condition_convergence_test)))
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic                                                                     │
│ ════════════════════════════                                                                     │
│ #spatial dimensions: ………………………… 1                                                                │
│ mesh: ………………………………………………………………… Trixi.Affine DGMultiMesh with NDIMS = 1.                         │
│ equations: …………………………………………………… LinearScalarAdvectionEquation1D                                  │
│ initial condition: ……………………………… initial_condition_convergence_test                               │
│ boundary conditions: ………………………… 1                                                                │
│ │ entire_boundary: ……………………………… BoundaryConditionDirichlet{typeoal_condition_convergence_test)} │
│ source terms: …………………………………………… nothing                                                          │
│ solver: …………………………………………………………… DG                                                               │
│ total #DOFs per field: …………………… 25                                                               │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> ode = semidiscretize(semi, (0.0, 1.7))
ODEProblem with uType RecursiveArrayTools.VectorOfArray{Float64, 2, StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 1.7)
u0: VectorOfArray{Float64,2}:
25×1 StructArray(::Matrix{Float64}) with eltype SVector{1, Float64}:
 [0.9999999999999999]
 [0.8923110900798575]
 [0.7308917846662063]
 [0.6464466094067264]
 [0.5785966677399663]
 [0.5117346022087871]
 [0.5]
 [0.511734602208787]
 [0.5785966677399662]
 [0.6464466094067262]
 [0.7308917846662062]
 [0.8923110900798575]
 [1.0]
 [1.1076889099201424]
 [1.269108215333794]
 [1.3535533905932737]
 [1.421403332260034]
 [1.488265397791213]
 [1.5]
 [1.488265397791213]
 [1.421403332260034]
 [1.3535533905932737]
 [1.269108215333794]
 [1.1076889099201426]
 [1.0]

julia> sol = solve(ode, RDPK3SpFSAL35(); ode_default_options()...);

julia> plot(sol)

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions