Errors in Numerical Solutions of Shock Physics Problems

1
Errors in Numerical Solutions of Shock Physics
Problems

• Presented by
• Yan Yu
• Advisor James Glimm

2
Presentation Outline

• Introduction
• Error Analysis for Planar Geometry
• Statistical Numerical Riemann Problem
• Composition Law for Errors
• Numerical Results
• Error Analysis for Spherical Geometry
• Statistical Numerical Riemann Problem
• Composition Law for Errors
• Error Decomposition
• Conclusions
• Future Objective

3

• Part ?
• Introduction

4
Introduction

• Our goal
• To formulate and validate a composition
  law to estimate errors in the solutions of
  complex problems in terms of the errors from the
  simpler ones. To understand the relative
  magnitude of the input uncertainty vs. the errors
  created within the numerical solution.
• References
• 1. Statistical riemann problems and a
  composition law for errors in numerical solutions
  of shock physics problems, J. Glimm, J. W. Grove,
  Y. Kang, T. Lee, X. Li, Y. Yu, K. Ye and M. Zhao,
  SISC(2003)
• 2. Errors in numerical solutions of
  spherically symmetric shock physics problems, J.
  Glimm, J. W. Grove, Y. Kang, T. Lee, X. Li, Y.
  Yu, K. Ye and M. Zhao, Contemporary Mathematics
  (2004)
• 3. Error Analysis of Composite Shock
  Interaction Problems, T. Lee, Y. Yu, M. Zhao, J.
  Glimm, X. Li and K. Ye, Conference Proceedings of
  PMC2004 (2004)

5
Graphic View
Introduction
From Study error models for individual
shock interactions To Estimate errors in
the complex multi-interaction problem
6
Uncertainty Quantification
Introduction

• Related approaches
• An early focus round off errors.
• Asymptotic analysis of numerical solution errors.
• The use of a posteriori error estimator.
• Mapping of input uncertainty to output random
  variables.
• Our approach
• We study the errors statistically, in a
  pre-asymptotic range.
• We allow for errors generated within the solution
  process.

7
Main Steps
Introduction

• Construct wave filters that decompose a complex
  flow into approximately independent elementary
  waves.
• Determine solution error models for a
  compre-hensive set of elementary wave
  interactions.
• Formulate a composition law that constructs the
  total solution error at any space-time point in
  terms of errors from repeated elementary
  interactions.

8
The Wave Filter
Wave profiles are reconstructed
For contact or shock waves
For rarefactions and compressions
9
Operation of A Wave Filter
contact or shock
rarefaction or compressions
10
SNRP

• Statistical Numerical Riemann Problem
• A statistical distribution of numerical incoming
  waves and starting states determines the SNRP.
• The waves in the SNRP have a finite width and the
  solution algorithm in the SNRP has only finite
  accuracy.
• The SNRP introduces errors in addition to
  propagating errors or uncertainty from input to
  output.

11

• Part ??
• Error Analysis for Planar Geometry

12
List of Interactions (SNRPs)
Planar Geometry

• Left frame shows the type and location of waves
  as determined by our
• wave filter analysis for the complex
  multi-interaction problem.
• Right frame is a schematic representation of the
  waves and the interactions. Ten Riemann problems
  are extracted from it.

13
Isolated Shock Waves

• Left frame shows the expected narrow and time
  independent shock width ( 2?x).
• Right frame shows an exponential approach of the
  numerical shock profile to its limiting values at
  x 8. The straight line is the asymptote to the
  exponential convergence rate.

14
Isolated Contact Waves
Top step down contact (flow from high density to
low) Contact width grows with a rate
asymptotically proportional to t1/3.
i.e. wc cct1/3, and we find cc
1. Bottom step up contact (the reverse)
We find wc min5, cct1/3. Difference
the spreading is associated with the up stream
side of the contact, and the continued spreading
depends on the up stream flow being subsonic. We
used the MUSCL scheme here.
15
Multinomial Expansion for SNRP Output
First, we represent the wave properties as
Then, the multinomial expansion for the output is
defined as (for wave strength)
Given a statistical ensemble of input and output
values ?i and ?o, we use a least square algorithm
to determine the best fitting model parameters a
K,J.
16
Expansion Coefficients
Shock contact interaction
coef variable constant ?1i ?2i model error model error
coef variable constant ?1i ?2i L8 STD
?1o -0.208 0.454 0.251 0.47 0.001
?2o -0.042 0.000 0.912 0.03 0.0001
?3o -0.286 1.004 0.346 0.30 0.001
?1o 2.184 -0.563 0.000 122 0.240
?2o 4.725 0.110 -1.466 0.67 0.010
?3o 2.197 0.068 0.106 5.35 0.057
p1o 0.221 -0.014 0.023 27.1 0.022
p2o 0.426 0.001 -0.092 1.78 0.002
p3o 0.332 -0.004 -0.099 3.47 0.005
17
Expansion Coefficients
Shock reflection
coef variable constant ?1i model error model error
coef variable constant ?1i L8 STD
?1o -0.002 0.716 0.014 0.00003
?1o 2.291 -0.422 7.923 0.062
p1o 0.060 -0.039 5065 0.009
?2o 0.057 0.0003 24.4 0.005
?2o 5.9 50 0.7
wall errors
18
Expansion Coefficients
Contact reshock interaction
coef variable constant ?1i ?2i model error model error
coef variable constant ?1i ?2i L8 STD
?1o 0.282 -0.314 0.645 0.57 0.0008
?2o 0.013 0.819 0.118 0.20 0.0003
?3o -0.128 0.143 0.458 0.41 0.0004
?1o 2.383 0.754 -1.307 5.47 0.038
?2o 0.909 0.011 0.216 1.00 0.005
?3o 3.619 0.151 -0.974 14.8 0.138
p1o 0.242 0.043 0.042 10.4 0.014
p2o -0.036 0.045 0.066 75.5 0.008
p3o -0.447 0.078 -0.036 15.7 0.029
19
The Composition Law
The composition law A formula for
combining the wave interaction errors defined for
single Riemann problems to yield the error for
arbitrary points in the complex wave interaction
problem. Which pieces to add up? The initial
waves and errors located inside the domain of
dependence.
Domain of dependence
20
Multi-path Integral
The composition law
The total error at the final time can be
formulated as
G a planar graph with all the Riemann
problems and traveling waves. V the set of
vertices of G, VV(G) . B a subset of V,
the vertices of which are inside the domain of
dependence. Iv the interaction
coefficient, taken from the table we obtained
before. d?B the multi-path propagator, a
product of the individual propagator for
each single path.
21
Transmission of Position Errors
The composition law
For interaction with a wall


For interaction of two incoming waves

22
Errors in Fully Resolved Calculations
The composition law
Simulation Prediction
mean wave strengths mean wave strengths mean wave strengths
?1o 0.451 0.452
?2o 0.704 0.703
?3o 0.999 0.998
wave strength errors wave strength errors wave strength errors
Var(?1o) 0.0008 0.0008
Var(?2o) 0.0019 0.0018
Var(?3o) 0.0035 0.0036
Simulation Prediction
wave width errors wave width errors wave width errors
?1o 1.630 1.622
?2o 3.636 3.635
?3o 2.346 2.352
wave position errors wave position errors wave position errors
P1o 0.220 0.226
P2o 0.313 0.312
p3o 0.200 0.202
for interaction 1
23
Errors in Fully Resolved Calculations
The composition law
Simulation Prediction
mean wave strengths mean wave strengths mean wave strengths
?1o 0.713 0.714
wave strength errors wave strength errors wave strength errors
Var(?1o ) 0.0018 0.0018
wave width errors wave width errors wave width errors
?1o 1.868 1.859
wave position errors wave position errors wave position errors
P1o -0.118 -0.092
for interaction 2
24
Errors in Fully Resolved Calculations
The composition law
Simulation Prediction
mean wave strengths mean wave strengths mean wave strengths
?1o 0.520 0.519
?2o 0.674 0.674
?3o 0.306 0.305
wave strength errors wave strength errors wave strength errors
Var(?1o) 0.0009 0.0010
Var(?2o) 0.0012 0.0013
Var(?3o) 0.0004 0.0004
Simulation Prediction
wave width errors wave width errors wave width errors
?1o 2.097 1.982
?2o 5.027 4.918
?3o 2.875 3.033
wave position errors wave position errors wave position errors
P1o -0.097 -0.105
P2o -0.003 0.013
p3o -0.151 -0.134
for interaction 3
25

• Part ???
• Error Analysis for Spherical Geometry

26
Main Steps
Spherical Geometry
Composition
Wave
Law
Filter
Data
Analysis
27
List of Interactions (SNRPs)
Spherical Geometry
This plot shows the type and location
of waves as determined by our wave filter
analysis for the complex multi-interaction
problem. It is also a schematic representation of
the waves and the interactions. Five Riemann
problems are extracted from it.
28
1D Spherical Geometry
New issues The solution is not constant
between two waves. Waves do not have constant
strength between interactions. Causes As the
shock wave strengthens while moving inward, it
generates a moving out rarefaction of the
opposite family. Challenge Need a simple
model for the growth of a wave as it moves
radially inward or the decrease as it moves out.
29
Single Propagating Shock Waves
Inward Shock
The Guderley Solution

• Left frame shows the Mach number vs. radius for a
  single inward propagating shock.
• Right frame shows the same data plotted on a
  log-log scale.

30
Single Propagating Shock Waves
Outward Shock

• Left frame shows the Mach number vs. radius for a
  single outward propagating shock wave starting at
  different radii r0.
• Right frame shows the same data plotted on a
  log-log scale, the dashed lines represent the
  power law model.

31
Multinomial Expansion for SNRP Output
First, we represent the wave properties as
Then, the multinomial expansion for the output is
defined as (for wave strength)
Given a statistical ensemble of input and output
values ?i and ?o, we use a least square algorithm
to determine the best fitting model parameters a
K,J.
32
Expansion Coefficients
Shock contact interaction
coef variable constant ?1i ?2i model error model error
coef variable constant ?1i ?2i STD STD/ ?o
?1o -33.353 19.521 2.501 0.860 0.954
?2o 0.374 0.200 0.0003 0.042 7.650
?3o 3.568 0.402 -0.045 0.009 0.463
d10 -2.039 -3.200 -0.01 0.157 0.174
d20 0.236 0.016 -0.002 0.021 3.825
d30 0.053 0.003 -0.001 0.0008 0.041
?1o 1.675 0.305 0.017 0.085
?2o 7.093 0.482 -0.146 0.239
?3o 2.829 0.302 -0.024 0.107
p1o -0.247 0.242 0.005 0.009
p2o 0.643 0.065 -0.011 0.192
p3o -0.042 0.062 0.004 0.009
33
Expansion Coefficients
Shock reflection
coef variable constant ?1i model error model error
coef variable constant ?1i STD STD/ ?o
?1o -242.394 5.606 1.137 0.468
d10 -3.27 0.031 0.112 0.045
?1o 1.221 0.018 0.099
p1o 0.474 0.001 0.012
34
Expansion Coefficients
Contact reshock interaction
coef variable constant ?1i ?2i model error model error
coef variable constant ?1i ?2i STD STD/ ?o
?1o 0.097 -0.108 0.436 0.031 13.305
?2o 0.103 -0.192 1.168 0.007 1.116
?3o 0.988 0.195 -0.225 0.003 0.262
d10 -0.291 0.161 -0.468 0.017 7.296
d20 -0.067 0.142 -0.125 0.006 0.957
d30 -0.030 0.107 -0.0004 0.001 0.087
?1o 9.776 -6.372 5.091 0.484
?2o 1.903 0.156 -0.677 0.534
?3o 4.088 -1.401 1.549 0.168
p1o 4.782 -3.602 2.372 0.379
p2o -0.453 0.409 -0.054 0.177
p3o -0.199 -0.685 3.213 0.052
35
Composite Shock Interaction Problems
The composition law
The total error at the final time can be
formulated as
G a planar graph with all the Riemann
problems and traveling waves. V the set of
vertices of G, VV(G) . B a subset of V,
the vertices of which are inside the domain of
dependence. Iv the interaction
coefficient, taken from the table we obtained
before. d?B the multi-path propagator, a
product of the individual propagator for
each single path.
36
Errors in Fully Resolved Calculations
The composition law
Simulation Prediction Simulation Prediction
100 vs. 2000 mesh 100 vs. 2000 mesh 500 vs. 2000 mesh 500 vs. 2000 mesh
wave strength errors and propagated initial uncertainties wave strength errors and propagated initial uncertainties wave strength errors and propagated initial uncertainties wave strength errors and propagated initial uncertainties
d1o 0.042(0.03) 0.032(0.02) 0.012(0.02) 0.0092(0.01)
d2o 0.142(0.05) 0.122(0.02) 0.032(0.01) 0.032(0.008)
d3o -0.022(0.02) -0.022(0.01) -0.0062(0.005) -0.0072(0.004)
mean wave width errors mean wave width errors mean wave width errors mean wave width errors
?1o 3.04 2.83 2.63 2.72
?2o 5.36 6.11 5.56 6.08
?3o 2.71 3.04 2.92 2.98
mean wave position errors mean wave position errors mean wave position errors mean wave position errors
P1o 1.25 0.23 0.12 0.18
P2o 0.43 0.06 0.05 0.04
p3o -0.73 -0.15 -0.08 -0.11
37
Decomposition of the Error
Left frame the type and location
of waves as determined by our wave filter
analysis for the complex multi-interaction
problem. Right frame Schematic diagram shows
different sources contribute to the error at the
output of interaction 3
38
Main Results from Data Analysis
Pie charts show the relative
importance of transmitted input uncertainty and
created error (100 and 500 mesh units,
respectively). The plots are based on the
variance of error or uncertainty associated with
each of the six wave integral paths contributing
to the post reshock wave strength for the contact
wave. The two paths associated with input
uncertainty are shown in hatched gray and the
others are solid gray scales.
39

• Part ?V
• Conclusions

40
Conclusions

• A simple linear model of solution error is
  sufficient for the study of a highly nonlinear
  problem.
• A composition law for combining errors and
  predicting errors for composite interactions on
  the basis of an error model of the simple
  constituent interactions has been formulated and
  validated.
• For spherically symmetric shock physics problems,
  the main difficulties encountered were the
  non-constancy of waves and errors between
  interactions. The errors grow by a power law in
  the radius for a inward moving shock. Similarly
  outward moving shocks and their errors weaken by
  a power law.

41
Conclusions
Continued

• For planar case, although our formulism allows
  for statistical errors in the ensemble, in fact
  the dominant part of all errors studied were
  deterministic.
• For spherical geometry, we observed that for a
  500 cell grid, the dominant error comes from the
  initial uncertainty, while for the 100 grid over
  75 of the error arises within the numerical
  simulations.
• The wave filter performs well as a diagnostic
  tool, but its limitation lies in assuming well
  separated waves.

42

• Part V
• Future Objective

43
Future Objective

• Error analysis for 2D perturbed
  simulations
• Start with the solution convergence of direct
  numerical simulations (DNS) through mesh
  refinement, regardless of different numerical
  algorithms.
• Determine the solution sensitivity to mesh size,
  algorithm, mass diffusion and etc.
• Conduct error analysis in 2D perturbed
  simulations (both single mode perturbed interface
  and multi-mode perturbed interface).
• Code comparison. For instance, FronTier code
  (tracked and untracked), AMR code, and AMR code
  with curvilinear coordinate, etc.

44
Perturbed Interface
Future work
Single-mode perturbation
Multi-mode perturbation
45
Preliminary Simulations
Future work
Without offset
With offset
46
Code Comparison
Future work
FronTier simulation
AMR simulation
47
Thank You