A-posteriori estimators for conservation laws A NDREAS D EDNER AND J AN G ISSELMANN A.S.Dedner@warwick.ac.uk, www2.warwick.ac.uk/fac/sci/maths/people/staff/andreas_dedner Mathematics Institute, Birmingham, Januar 5, 2016
Structure of solution Scalar non linear conservation law: Find u : R d × R + → R solution of ∂ t u ( x , t ) + ∇ · f ( u ( x , t )) = 0 u ( x , 0 ) = u 0 ( x ) with e.g. f ( u ) = 1 2 u 2 Viscosity limit: let u ε be a classical solution of the regularized problem: ∂ t u ε ( x , t ) + ∇ · f ( u ε ( x , t )) = ε △ , u ε ( x , t ) , u ε ( x , 0 ) = u 0 ( x ) There exists u = lim ε → 0 u ε (a.e.) and u is weak solution. u is physically relevant weak solution Equivalent: Entropy Solution � � � − ( S ( u ) ∂ t φ + F S ( u ) · ∇ φ ) dt dx − S ( u 0 ) φ ( x , 0 ) dx ≤ 0 R d R + R d for all entropy pairs ( S , F S ) , i.e., S convex and F ′ S = S ′ f ′
First order finite-volume scheme System of conservation law (e.g. Euler Equations): Find U : R d × R + → R m entropy solution ∂ t U ( x , t ) + ∇ · F ( U ( x , t )) = 0 U ( x , 0 ) = U 0 ( x ) Integrate over T ∈ T with tessellation T of Ω : � � � ∂ t U ( · , t ) = − ∇ · F ( U ( · , t )) = − F ( U ( · , t )) · n ∂ T T T � 1 Piecewise constant approximation U T ( t ) ≈ T U ( · , t ) : | T | � dtU T ( t ) = − 1 d F h ( t ) | T | ∂ T with numerical flux F h ( t ) = F T , T ′ ( U T ( t ) , U T ′ ( t )) on intersection between neighboring elements T , T ′ : � dtU T ( t ) = − 1 d F T , T ′ ( U T ( t ) , U T ′ ( t )) | T | T ′ ∈ N ( T ) N ( T ) is set of all neighbors of T .
First order finite-volume scheme System of conservation law (e.g. Euler Equations): Find U : R d × R + → R m entropy solution ∂ t U ( x , t ) + ∇ · F ( U ( x , t )) = 0 U ( x , 0 ) = U 0 ( x ) Semi discrete scheme: � dtU T ( t ) = − 1 d F T , T ′ ( U T ( t ) , U T ′ ( t )) | T | T ′ ∈ N ( T ) N ( T ) is set of all neighbors of T . Forward Euler in time for time steps t n and ∆ t n = t n + 1 − t n : � T − ∆ t n U n + 1 = U n F T , T ′ ( U n T , U n T ′ ) T | T | T ′ ∈ N ( T )
First order finite-volume scheme System of conservation law (e.g. Euler Equations): Find U : R d × R + → R m entropy solution ∂ t U ( x , t ) + ∇ · F ( U ( x , t )) = 0 U ( x , 0 ) = U 0 ( x ) Fully discrete scheme: � T − ∆ t n U n + 1 = U n F T , T ′ ( U n T ( t ) , U n T ′ ( t )) T | T | T ′ ∈ N ( T ) Define U h ( x , t ) := U n T for x ∈ T , t ∈ [ t n , t n + 1 ) . A-priori error estimate for scalar case Let u h be a first order finite-volume approximation then under suitable con- ditions on the numerical flux: 1 � u h ( · , t ) − u ( · , t ) � L 1 ( R d ) ≤ C ( u ) h max 4 t 1 2 , only proven for structured grids. Should be h
Euler System of Hydrodynamics ∂ t ρ + ∇ · ( ρ� u ) = 0 (conservation of mass) , u T + P ) = 0 ∂ t ( ρ� u ) + ∇ · ( ρ� u � (conservation of momentum) , ∂ t ( ρ e ) + ∇ · ( ρ e � u + P � u ) = 0 (conservation of energy) , e = ε + 1 u | 2 2 | � (total energy) , P = p ( ρ, ε ) I (equation of state for pressure) , ( ρ, ρ� u , ρ e )( · , 0 ) = ( ρ 0 , ρ 0 � u 0 , ρ 0 e 0 ) (initial conditions) ρ : density , u : momentum , ρ e : total energy density ρ�
Numerical results Riemann problem (Euler) u 0 = u L ( x < 0 ) , u 0 = u R ( x > 0 ) . Solution: rarefaction, contact, shock 1 0.001 0.8 1 contact shock 0.9 0.7 0.8 0.1 L 1 convergence rate 0.7 0.6 density L 1 -error 0.6 error 0.01 0.0001 0.5 contact 0.5 0.4 0.001 shock 0.3 0.4 0.2 solution error approx eoc 0.1 0.0001 1e-05 0.3 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.0001 0.001 0.01 x x h Single contact (Euler) u 0 = u L ( x < 0 ) , u 0 = u R ( x > 0 ) . Solution: u ( x , t ) = u 0 ( x − at ) 1 0.01 0.6 1 0.9 0.55 0.8 0.1 L 1 convergence rate 0.5 0.7 density L 1 -error 0.6 error 0.01 0.001 0.45 0.5 0.4 0.4 0.001 0.3 0.35 0.2 solution error approx eoc 0.1 0.0001 0.0001 0.3 0.2 0.25 0.3 0.35 0.4 0.45 0.2 0.25 0.3 0.35 0.4 0.45 0.0001 0.001 0.01 x x h resolution in 3d requires 130.000.000 elements for results shown
Numerical results Forward facing step (Euler) Right moving Mach 3 flow structure: 15.000 elements 230.000 element
High order methods • highly efficient for smooth solution • Loss of efficiency for non-smooth solution ... • ... and unstable for non linear discontinuities (shocks) No stabilization 1.4 exact solution unlimited DG solution With simple stabilization 1.2 1 0.1 1 0.9 0.8 0.8 L 1 -error u(x) L 1 -eoc 0.01 0.7 0.6 0.6 p=3 p=3 0.4 0.5 p=2 p=2 p=1 p=1 0.001 0.4 100 1000 10000 100000 100 1000 10000 100000 0.2 grid size grid size 0.1 1.04 1.02 0 1 0.98 0.01 0.96 L 1 -error L 1 -eoc -1 -0.5 0 0.5 1 0.94 0.92 x 0.9 0.001 0.88 p=3 p=3 0.86 p=2 p=2 0.84 p=1 p=1 0.0001 0.82 100 1000 10000 100000 100 1000 10000 100000 grid size grid size
High order methods • highly efficient for smooth solution • Loss of efficiency for non-smooth solution ... • ... and unstable for non linear discontinuities (shocks) Basis Algorithm: 1 determine troubled cells where the error is high or the scheme is unstable 2 for each troubled cell either increase the order (if solution is smooth) or reduce the order und refine the grid Determine troubled cells heuristically or by error estimate.
Higher order (Discontinuous Galerkin) Riemann problem (Euler) u 0 = u L ( x < 0 ) , u 0 = u R ( x > 0 ) . Solution: rarefaction, contact, shock 1 1 p=0 error p=0 1.2 contact shock p=2 eoc p=0 0.9 error p=2 1.1 eoc p=2 0.8 0.001 0.1 L 1 convergence rate 1 0.7 density L 1 -error 0.9 0.6 error 0.01 0.5 0.0001 0.8 contact 0.4 0.7 shock 0.001 0.3 0.6 solution 0.2 p=0 1e-05 0.5 p=2 0.1 0.0001 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.0001 0.001 0.01 x x h Single contact (Euler) u 0 = u L ( x < 0 ) , u 0 = u R ( x > 0 ) . Solution: u ( x , t ) = u 0 ( x − at ) 1 1 error p=0 1.1 eoc p=0 0.9 error p=2 1 0.001 eoc p=2 0.8 0.1 L 1 convergence rate 0.7 0.9 density L 1 -error 0.6 error 0.01 0.8 0.5 0.0001 0.7 0.4 0.001 0.3 0.6 solution 0.2 p=0 p=0 1e-05 0.5 p=2 p=2 0.1 0.0001 0.2 0.25 0.3 0.35 0.4 0.45 0.2 0.25 0.3 0.35 0.4 0.45 0.0001 0.001 0.01 x x h
Heuristic strategy for local mesh adaption • Find interface between cells where solution has a large jump • Mark the two elements at that intersection • Mark a neighborhood of marked elements • Mark elements for coarsening where the jump is very small Possibly look at curvature of solution (i.e. jumps between gradients) Idea: Error is where the discontinueties are Problems? • Can not distinguish between contacts and shocks • Could coarsen wrongly (kinks at ends of rarefactions) • Indicator does not get smaller with reduction of grid size • No mathematical proof that it works only many people using successfuly...
Higher order schemes on adaptive grids Estimator Nead indictor η K for "smoothness" of solution � O ( h q K ) smooth region η K = O ( h − 1 K ) troubled region
A posteriori error estimate Kruzkov framework (semi implicit) D, Makridakis, Ohlberger ’06 Structure of a posteriori error estimate: � � � � � � || ( u − u h )( T ) || 2 R n T , j + R n S , jl + R n u n u n L 1 ( B R ( x 0 )) ≤ K h j + || � j − � j || L ∞ L , j n j ∈ J n • R n T , j : Element residual • R n S , jl : Jump residual (numerical viscosity) • R n L , j : Coarsening error • || � u n u n j − � j || L ∞ : difference between aerage and higher order polynomial Numerical test show: � O ( h p − 1 � � ) h j solution is smooth ¯ R n R n T , j + R n S , jl + R n j j := = L , j | T j | ∆ t n O ( 1 ) ... discontinuous
Strategy for local mesh adaption First step: Define the set of grid cells I s with a significant contribution to the overall error indicator η h . Second step: Use an equal distribution strategy to refine or coarsen the elements in I s according to the error estimate. Third step: For elements that are marked for coarsening, check if the projection error is small enough.
Forward facing step Approximately 14 . 000 elements at t = T
We now come to a short commercial break... at University of Warwick 4th to 8th of July, 2016
We now come to a short commercial break... at University of Warwick 4th to 8th of July, 2016
A very different approach Issues with Kruzkov: only works with scalars Issue with proof: only done for semi discrete scheme Giesselmann, Makridakis, Pryer: use of relative entropy (RE) framework Given one convex entropy η then η ( U | V ) := η ( U ) − η ( V ) − η ′ ( V )( U − V ) ≈ � U − V � L 2 and if V is a Lipschitz solution and U a weak solution then � � d η ( U | V ) ≤ C ( V ) η ( U | V ) dt This can also be used to study perturbed solutions U : Advantage: works with system of equations (need only one entropy) Issue with RE: only works in the case that a Lipschitz solution extists Issue with proof: semi discrete scheme in 1d, very restrictive on flux function
Recommend
More recommend
