Numerical Estimation of DiscretizationError on Unstructured MeshesbyKai Kin Gary YanBASc., University of Toronto, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)March 2018c© Kai Kin Gary Yan 2018AbstractA numerical estimation of discretization error is performed for solutions tosteady and unsteady models for compressible flow. An accurate and reliableestimate of discretization error is useful in obtaining a more accurate defectcorrected solution, as well as a tight uncertainty bound as error bars. Theerror estimation procedure is performed by solving an auxiliary problem,known as the error transport equation (ETE), solved on the same mesh asthe original model equations. Unlike unsteady adjoint methods for error infunctionals, the ETE requires only one other set of equations to be solved,agnostic to the choice and number of output functionals, including commonaerodynamic quantities such as lift or drag. Furthermore, co-advancing theETE in time only requires the storage of local solutions in time and notthe entire history, reducing memory requirements. This method of errorestimation is performed in the context of higher order finite-volume meth-ods on unstructured meshes. Approaches based on solving the ETE canbe found in the literature for uniform or smooth meshes, but this has notbeen well studied for unstructured meshes. Such meshes necessarily havenonsmooth geometric features, which create many difficulties in accurateerror estimation. These difficulties in accurate discretizations of the ETE areinvestigated, including the discretization of the ETE source term, which iscritical to error estimate accuracy. It was found that the proposed schemesby the ETE approach can be more efficient and robust compared to solvingthe higher order problem. The choices of discretization schemes need tobe made carefully, and these results demonstrate how it is possible, alongwith justification, to obtain asymptotically accurate, efficient, and robust er-ror estimates that can be used with vast possibilities of model equations inpractice.iiLay SummaryThe use of computer simulations for scientific predictions and engineeringdesign has become commonplace nowadays. However, the quality of a par-ticular solution, measured in terms of computational error, is not well stud-ied, but is nevertheless a critical measure of how reliable the computed so-lution is in agreeing with observation. The current work is concerned withobtaining efficient and robust error estimates through solving an auxiliaryequation for the solution error, which can be used either as a correctionor as a bound on the error magnitude, so that effective validation can beperformed.iiiPrefaceThe developments of the methods and results of this dissertation are theoutcomes of the research conducted by the author, Kai Kin Gary Yan, un-der the supervision of Dr. Carl Ollivier-Gooch, as part of the PhD programrequirements in the department of Mechanical Engineering. The researchcode used for 2D problems is an ongoing development to which members ofthe research group contribute. The parts directly relevant to discretizationerror estimation, as well as a 1D code in its entirety, have been written bythe author, Kai Kin Gary Yan.Excerpts from the dissertation have been published. A version of Chapter5 for unsteady manufactured test problems has been published in Comput-ers & Fluids [121]. Parts of Chapters 3 and 4 have been submitted to theJournal of Verification, Validation, and Uncertainty Quantification and thearticle has been accepted for publication [120]. An article containing the re-sults for applications to unsteady problems with periodic functionals in thelatter part of Chapter 5 has been submitted for publication. Results for relin-earization and applications to turbulent models in the latter part of Chapter4 are currently in preparation for submission for publication. Conferenceversions of the papers can be found in [115–119].ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Mesh Properties . . . . . . . . . . . . . . . . . . . . . . . . . 51.3 Classification of Error . . . . . . . . . . . . . . . . . . . . . . 81.4 Approaches to Error Estimation . . . . . . . . . . . . . . . . . 111.4.1 Method of Manufactured Solutions . . . . . . . . . . . 111.4.2 Richardson Extrapolation . . . . . . . . . . . . . . . . 111.4.3 Adjoint Methods . . . . . . . . . . . . . . . . . . . . . 121.4.4 Error Transport Equation . . . . . . . . . . . . . . . . 141.4.5 Higher Order Methods . . . . . . . . . . . . . . . . . 151.4.6 Quality of Error Estimates . . . . . . . . . . . . . . . . 151.5 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161.6 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.1 Discretization in Space . . . . . . . . . . . . . . . . . . . . . 202.2 Finite-Volume Discretization . . . . . . . . . . . . . . . . . . 22vTable of Contents2.3 Order of Accuracy . . . . . . . . . . . . . . . . . . . . . . . . 262.3.1 Error Accuracy . . . . . . . . . . . . . . . . . . . . . . 272.3.2 Higher Order Accuracy via k-Exact Reconstruction . . 282.3.3 Residual Evaluation . . . . . . . . . . . . . . . . . . . 322.4 Discretization in Time . . . . . . . . . . . . . . . . . . . . . . 342.4.1 Residual Evaluation . . . . . . . . . . . . . . . . . . . 372.5 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 372.6 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . 383 ETE Accuracy on Smooth and Nonsmooth Meshes . . . . . . . . 393.1 Spectral Properties . . . . . . . . . . . . . . . . . . . . . . . . 403.2 Uniform Meshes . . . . . . . . . . . . . . . . . . . . . . . . . 413.3 Smooth Meshes . . . . . . . . . . . . . . . . . . . . . . . . . 443.4 Perturbed Meshes . . . . . . . . . . . . . . . . . . . . . . . . 493.4.1 Nested Grids . . . . . . . . . . . . . . . . . . . . . . . 533.5 Nonlinear Model Problems . . . . . . . . . . . . . . . . . . . 533.6 Other Attempts at Truncation Error Estimation . . . . . . . . 573.6.1 Residual Filtering . . . . . . . . . . . . . . . . . . . . 573.6.2 h-Truncation Error Estimate . . . . . . . . . . . . . . . 583.6.3 Projection onto Uniform Mesh . . . . . . . . . . . . . 603.7 Error Estimation for Other Discretization Schemes . . . . . . 623.7.1 Other Discretization Orders . . . . . . . . . . . . . . . 623.7.2 Finite-Difference Method . . . . . . . . . . . . . . . . 643.7.3 Finite-Element Methods . . . . . . . . . . . . . . . . . 674 ETE Methodology Applied to Steady Problems . . . . . . . . . . 694.1 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.1.1 Relinearization of the ETE . . . . . . . . . . . . . . . 744.2 Properties of ETE Discretization . . . . . . . . . . . . . . . . 754.2.1 Exact Discrete Residual Source Term . . . . . . . . . . 754.2.2 Expected Order of Accuracy for the Error Estimate . . 784.2.3 Equivalence of Defect Correction and Higher OrderScheme when q = r . . . . . . . . . . . . . . . . . . . 794.2.4 Equivalence of Nonlinear ETE and Order Ramping whenq = r . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.3 Compressible Inviscid Flow . . . . . . . . . . . . . . . . . . . 804.3.1 Supersonic Vortex . . . . . . . . . . . . . . . . . . . . 814.3.2 Gaussian Bump . . . . . . . . . . . . . . . . . . . . . 844.4 Compressible Viscous Flow . . . . . . . . . . . . . . . . . . . 844.4.1 Manufactured Solution . . . . . . . . . . . . . . . . . 87viTable of Contents4.4.2 NACA 0012 Airfoil . . . . . . . . . . . . . . . . . . . . 874.4.2.1 Efficiency Considerations . . . . . . . . . . . 894.4.2.2 Robustness Considerations . . . . . . . . . . 914.5 Simulation of Turbulent Flows . . . . . . . . . . . . . . . . . 934.6 Zero Pressure Gradient Flat Plate . . . . . . . . . . . . . . . . 954.7 Flow over a NACA 0012 Airfoil . . . . . . . . . . . . . . . . . 975 ETE Methodology Applied to Unsteady Problems . . . . . . . . 1055.1 Time-Dependent Residual Evaluation . . . . . . . . . . . . . 1055.1.1 Source Term for Exact Error Estimate . . . . . . . . . 1055.1.2 T = 0 Gives the Trivial Zero Estimate . . . . . . . . . 1065.1.3 Higher Order Error Estimate in Space Only . . . . . . 1065.1.4 MEA Source Term . . . . . . . . . . . . . . . . . . . . 1075.1.5 Computing the Time Derivative in −(∂tu˜+R(u˜)) Di-rectly by Finite-Difference . . . . . . . . . . . . . . . . 1115.2 Expected Error Estimate Accuracy Using FD ETE Source . . . 1125.3 Co-Advancing the ETE . . . . . . . . . . . . . . . . . . . . . . 1145.4 Accuracy Comparisons for Error in Solution . . . . . . . . . . 1155.4.1 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . 1165.4.2 Advection . . . . . . . . . . . . . . . . . . . . . . . . . 1175.4.3 Navier-Stokes Equations . . . . . . . . . . . . . . . . . 1205.5 Accuracy Comparisons for Error in Functionals . . . . . . . . 1215.5.1 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . 1255.5.2 Advection . . . . . . . . . . . . . . . . . . . . . . . . . 1265.5.3 Navier-Stokes Equations . . . . . . . . . . . . . . . . . 1325.6 Comparisons to the Higher Order Primal Problem . . . . . . . 1375.7 Smoothly Varying Time Steps . . . . . . . . . . . . . . . . . . 1386 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . 1416.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1416.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 1426.2.1 Extensions to 3D . . . . . . . . . . . . . . . . . . . . . 1426.2.2 PDEs of Higher Order . . . . . . . . . . . . . . . . . . 1456.2.3 Turbulence Modeling Extensions . . . . . . . . . . . . 1466.2.4 Weak Solutions . . . . . . . . . . . . . . . . . . . . . 147Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150viiTable of ContentsAppendicesA Richardson Extrapolation . . . . . . . . . . . . . . . . . . . . . . 162A.1 Known p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162A.2 Unknown p . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163B CFD Model Equations . . . . . . . . . . . . . . . . . . . . . . . . 164C ANSLib Command Line Options . . . . . . . . . . . . . . . . . . 166C.1 Steady Problems . . . . . . . . . . . . . . . . . . . . . . . . . 167C.1.1 Inviscid Supersonic Vortex . . . . . . . . . . . . . . . 167C.1.2 Inviscid Gaussian Bump . . . . . . . . . . . . . . . . . 167C.1.3 NACA 0012 Laminar Navier-Stokes . . . . . . . . . . 167C.1.4 Flat plate with SA-neg Turbulence Model . . . . . . . 168C.1.5 NACA 0012 airfoil with SA-neg Turbulence Model . . 168C.2 Unsteady Problems . . . . . . . . . . . . . . . . . . . . . . . . 169C.2.1 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . 169C.2.2 Advection . . . . . . . . . . . . . . . . . . . . . . . . . 169C.2.3 Translating Vortex . . . . . . . . . . . . . . . . . . . . 170C.2.4 Vortex shedding . . . . . . . . . . . . . . . . . . . . . 170viiiList of Tables3.1 Empirical asymptotic sizes of the maximum eigenvalues andsingular values of the discretization Jacobian and its inverse. . 403.2 Asymptotic sizes for discretization error and truncation error. 414.1 Summary of values used for the supersonic vortex test case. . 814.2 Accuracy of entropy (error) for the different schemes for thelaminar Euler equations test case on the Gaussian bump ge-ometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.3 Accuracy of CD,Viscous for the different schemes for the lami-nar Navier-Stokes test case on the NACA 0012 geometry. . . . 894.4 Empirical values for the Spalart-Allmaras turbulence model. . 954.5 Accuracy of CD for the different schemes for the flat plate testcase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 964.6 Accuracy of CL for the NACA 0012 test case, along with ref-erence values. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045.1 Summary of constants used for the translating vortex test case.1215.2 Orders of accuracy for measures of functionals in the diffu-sion test case, using primal solution only and ETE correctedsolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1285.3 Orders of accuracy for measures of the functional in the ad-vection test case using perturbed ETE initial conditions. . . . 1315.4 Flow parameters in vortex shedding test case. . . . . . . . . . 1325.5 Orders of accuracy for the angular frequency of the functionalin the vortex shedding test case. . . . . . . . . . . . . . . . . . 134ixList of Figures1.1 The role of computational methods in validation of models. . 21.2 Classifications of meshes with quadrilateral cells. . . . . . . . 71.3 Typical unstructured mesh used on a square domain for thisstudy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.1 Conservation law in a domain. . . . . . . . . . . . . . . . . . 192.2 Discretization Jacobian fill comparison of the current finite-volume method and a DG scheme using FEniCS on the samemesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.3 Required reconstruction stencil for various orders of accuracy. 312.4 Schematic of second order interface flux evaluation based ona numerical flux function of contributions from either side. . . 312.5 Exact truncation error and estimate using a higher order dis-crete operator. . . . . . . . . . . . . . . . . . . . . . . . . . . 353.1 Accuracy of ETE estimate for uniform (2, 2, 4), N = 40. . . . . 423.2 Discretization error and error in error estimate for (2, 2, 4)using a uniform mesh. . . . . . . . . . . . . . . . . . . . . . . 433.3 Schematic of multigrid scheme for the ETE. . . . . . . . . . . 453.4 Fine and coarse grid numbering. . . . . . . . . . . . . . . . . 453.5 Discretization error and error in error estimate for (2, 2, 4)multigrid scheme using uniform meshes. . . . . . . . . . . . . 463.6 Mesh that is smoothly expanded by a factor 1 + s. . . . . . . . 473.7 Mesh points for a smoothly mapped mesh with N = 10, andits underlying grid mapping function. . . . . . . . . . . . . . . 483.8 Discretization error and error in error for (2, 2, 4) using smoothlystretched meshes. . . . . . . . . . . . . . . . . . . . . . . . . . 503.9 Exact discretization error and estimate for (2, 2, 4) and (2, 4, 4)schemes using a randomly perturbed N = 40 mesh. . . . . . . 513.10 Discretization error and error in error for (2, 2, 4) and (2, 4, 4)schemes using randomly perturbed meshes. . . . . . . . . . . 52xList of Figures3.11 Initial coarse perturbed mesh (N = 10) and its nested refine-ment until N = 40. . . . . . . . . . . . . . . . . . . . . . . . . 543.12 Error accuracy for the (2, 2, 4) scheme using a nested per-turbed grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.13 Error accuracy for the (2, 4, 4) scheme using a nested per-turbed grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.14 Discretization error and error in error using nested perturbedgrids. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.15 Discretization error and error in error using uniform and ran-dom meshes for Burgers’ equation. . . . . . . . . . . . . . . . 583.16 Discretization error and error estimate resulting from the ap-plication of a Gaussian filter to the residual source term. . . . 593.17 Discretization error and error estimate for a (2, 2, 2∗) schemeusing h-truncation error estimate. . . . . . . . . . . . . . . . . 603.18 Projection of primal solution from a perturbed mesh to a uni-form mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.19 Projection of residual from uniform mesh to original perturbedmesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.20 Discretization error and (2, 2, 4) error estimate using a uni-form mesh to compute residual, then solving ETE on per-turbed mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.21 Discretization error and (2, 2, 4) error estimate using a uni-form mesh to compute residual, then solving ETE on the sameuniform mesh, and projecting back to the perturbed mesh. . . 633.22 Discretization error and estimate for uniform mesh for the 1DPoisson equation. . . . . . . . . . . . . . . . . . . . . . . . . . 653.23 Discretization error and estimate for randomly perturbed meshfor the 1D Poisson equation. . . . . . . . . . . . . . . . . . . . 663.24 Discretization error and error in error for finite-differenceschemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.1 Results for the Poisson problem using the ETE with (2, 2, 4). . 704.2 Error estimate and error in error using the ETE with (2, 4, 4). . 714.3 Accuracy of the second order solution, and error estimatesusing the (2, 2, 4) and (2, 4, 4) schemes for the Poisson problem. 714.4 Convergence in error estimate with mesh length scale usingthe fully nonlinear error flux, the linearized flux, and the pri-mal flux for the 1D viscous Burgers’ equation. . . . . . . . . . 744.5 Accuracy of relinearization for 1D Burgers’ equation. . . . . . 76xiList of Figures4.6 Exact error, error estimate, and difference for the energy vari-able using the (2, 4, 4) discretization for the supersonic vortextest case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.7 Convergence of the error estimate for the supersonic vortexproblem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.8 Entropy error from primal solutions and ETE corrected solu-tion for the Gaussian bump test case. . . . . . . . . . . . . . . 854.9 Accuracy of entropy norm for the Gaussian bump test case. . . 864.10 Convergence of the error estimate in nondimensionalized en-ergy for the 2D Navier-Stokes equations using manufacturedsolutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.11 Comparison of solution error plotted against overall compu-tational time required for the 2D Navier-Stokes equations us-ing manufactured solutions. . . . . . . . . . . . . . . . . . . . 904.12 Comparison of an output functional (CD,Viscous) plotted againstmesh length scale and overall computational time requiredfor 2D Navier-Stokes equations for the NACA 0012 test case. 904.13 Particularly unsuitable meshes which converge for the lin-earized (2, 4, 4) scheme but not (4, 0, 0) or (2, 4, 4). . . . . . . 924.14 Geometry and boundary conditions for the flat plate test case. 974.15 Plot of the turbulence variable µ′ν˜ for N = 8160 and compar-ison with reference solutions. . . . . . . . . . . . . . . . . . . 984.16 Accuracy of CD for the flat plate test case. . . . . . . . . . . . 994.17 Functionals for the flat plate test case. . . . . . . . . . . . . . 1004.18 Turbulence variable ν˜, error estimates, and corrected solutions.1014.19 Estimate of error in the turbulence variable ν˜ for the secondorder primal solution and the linearized ETE corrected solution.1024.20 Functional accuracy for the NACA 0012 test case. The refer-ences use N = 230529. . . . . . . . . . . . . . . . . . . . . . . 1035.1 The diffusion problem using p = 2, q = 4, pt = qt = 2 andfq(U)− fp(U) as the ETE source. . . . . . . . . . . . . . . . . 1085.2 Accuracy of error estimate using the MEA source for a simpleODE system. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1105.3 The diffusion problem using p = 2, q = 4, pt = qt = 2 andMEA ETE source. . . . . . . . . . . . . . . . . . . . . . . . . . 1115.4 Norms of the MEA source term ||k2f3(U)/12|| at the initialtime, keeping k/h approximately fixed. . . . . . . . . . . . . . 112xiiList of Figures5.5 Plot of the evolution in time of the L2 norm of the exact error,and error estimates using the MEA and FD source terms forthe ETE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1135.6 The diffusion problem using (2, 4, 4; 2, 2) FD time derivativeETE source, with h ≈ 0.0059, k = 0.025. . . . . . . . . . . . . 1175.7 A comparison of lower order primal error, error in error us-ing ETE, and higher order primal error for diffusion using(2, 4, 4; 2, 2), with refinement in both k and h, k/h ≈ 4.2 fixed. 1185.8 The advection problem using (2, 4, 4; 2, 2) FD time derivativeETE source, with h ≈ 0.0059, k = 0.0025. . . . . . . . . . . . 1195.9 A comparison of lower order primal error, error in error us-ing ETE, and higher order primal error for advection using(2, 4, 4; 2, 2), with refinement in both k and h, k/h ≈ 0.42 fixed.1195.10 Plots of the norms of exact primal error and error in error esti-mate for the advection problem using different discretizationorders. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1205.11 Energy variable for the Navier-Stokes test case using (2, 4, 4; 2, 2)FD time derivative ETE source, with h ≈ 0.0059, k = 0.0025. . 1225.12 A comparison of lower order primal error, error in error usingETE, and higher order primal error for Navier-Stokes using(2, 4, 4; 2, 2), with refinement in both k and h, k/h ≈ 0.42fixed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1225.13 Solution and errors at the final time for the diffusion test case. 1265.14 Functionals J1, J2 for the diffusion test case. . . . . . . . . . 1275.15 Accuracy of measures of primal and corrected diffusion func-tional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1285.16 Solution and errors at the final time for the advection test case.1295.17 Functional for the advection test case using perturbed initialconditions for the ETE. . . . . . . . . . . . . . . . . . . . . . 1305.18 Accuracy of measures of primal and corrected advection func-tional using perturbed initial conditions for the ETE. . . . . . 1315.19 Solution and errors at the final time for the vortex sheddingtest case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1335.20 CL for the vortex shedding case. . . . . . . . . . . . . . . . . 1355.21 CL for the vortex shedding case using non-reflecting bound-ary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . 1365.22 Accuracy of Strouhal number computed from the frequencyof the lift coefficient functional using primal and ETE cor-rected solutions. . . . . . . . . . . . . . . . . . . . . . . . . . 136xiiiList of Figures5.23 Comparison of the efficiency of the ETE approach (combinedtime) versus solving only the primal problem. . . . . . . . . . 1375.24 Stability region of the BDF family of time discretization schemeswith different orders of accuracy. . . . . . . . . . . . . . . . . 1395.25 Error accuracy with uniform in space, time step varied smoothlyand randomly. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1406.1 Results for a 3D Poisson problem using a uniform mesh. . . . 1446.2 Accuracy of the primal solution and error estimate for a 3DPoisson problem, on uniform meshes. . . . . . . . . . . . . . . 1456.3 Lorenz solution using initial conditions that differ by a 0.1%perturbation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1486.4 Discretization error for stable and chaotic Lorenz system. . . . 149xivList of Symbols1 Identity matrixA Jacobian of discretizationCD Drag coefficientCf Skin friction coefficientCL Lift coefficientCP Pressure Coefficientcp Specific heat at constant pressurecv Specific heat at constant volumeD Finite-difference operatorD DiameterD Destruction term in negative Spalart-Allmaras turbulence modeld Number of derivativesd Number of dimensionsE Energye Exact continuous discretization errorF FluxF Interface fluxf Right hand side of the primal equation after discretization in spaceg Right hand side of the ETE after discretization in spaceh Length scale in spaceIh2h1 Mesh transfer operatorJ Functional in time domainJ Functional in frequency domaink Length scale in timeL Differential operatorL∗ Adjoint differential operatorM Mach numberN Number of control volumesNT Number of time stepsnˆ Outward unit normalxvList of Symbolsn Time levelP PressureP Production term in SA-neg modelp Primal order of accuracypt Primal order of accuracy in timePr Prandtl numberPrT Turbulent Prandtl numberq Heat fluxq Error transport equation order of accuracyqt Error transport equation order of accuracy in timeR Spatial residualR Discrete primal spatial residualRˆ Discrete ETE spatial residualr Reconstruction order for p-truncation errorrt Time accuracy for ETE source termRe Reynolds numberS Discrete source termS Vorticitys Continuous source termSt Strouhal numberT Time-dependent discrete residualT Temperaturet TimeU Discrete solutionU Exact solution on discrete spaceUp Discrete solution to order pUn Discrete solution at time level nu Exact continuous exact solutionu˜ Approximate solution interpolated onto exact solution spacev Velocity vectorw Reconstruction weightsx Collection of spatial coordinatesx, y Spatial coordinatesα Angle of attackα Jump termβ Reconstruction weight exponentγ Ratio of specific heats Discrete discretization errorp Exact discrete error estimate to order pxviList of Symbolspqr Discrete error estimate using (p,q,r) schemen Discrete error estimate at time level nκ Thermal conductivityµ ViscosityµT Turbulent viscosityν˜ Spalart Allmaras turbulence variableξ Grid mapping functionρ Densityτ Stress tensorτ Truncation errorΩ Domain in spaceΩT Domain in timeω FrequencyxviiAcknowledgmentsIt is with great humility that I express gratitude for my research supervisor,Dr. Carl Ollivier-Gooch, who is always looking out for the best interestsof his students. His unfaltering support and confidence in my academic,research, and personal development is demonstrably more than what I couldhave expected.I also wish to acknowledge the individuals who have taken the time togive impartial evaluations and validations of this work, including membersof the supervisory committee, the examining committee, and the scientificcommunity.I am also grateful for the financial support that has made this researchpossible, which was provided in part by the Natural Sciences and Engineer-ing Research Council of Canada, the Killam Trusts, the Four Year Fellowshipfrom the University of British Columbia, and the Air Force Office of ScientificResearch.Finally, I thank my colleagues at the Advanced Numerical SimulationLaboratory, with whom professional and personal discussions have alwaysproven constructive.xviiiChapter 1Introduction1.1 MotivationThe aspiration to predict, with confidence, observed phenomena by com-putation has existed since ancient times, even before the development offormal computational methods. In this endeavor, the first step is to devise,either from first principles or some more heuristic basis, a model whichmathematically represents or approximates the situation. Next comes anattempt to solve these mathematical equations. In computational fluid dy-namics (CFD), the current field of interest, while examples of analytical solu-tions for compressible flow can be found in the literature (see, for instance,[58]), they are for the most part only valid under restrictive assumptionson the physics or geometry. For virtually all realistic applications, the mod-els are complex enough that solutions cannot be obtained analytically. Analternative in making insightful predictions is to use physical experiments,which isolates some testing parameters, or a scaled version of the actualphenomenon. Although experimental observation is the gold standard forevaluating predictions, relying solely on experiment is not feasible in manyinstances. Experiments are often costly, unsafe, and come with their ownuncertainties such as in measurements or the question of repeatability in acontrolled environment. In light of these alternatives, computational meth-ods can be seen as a bridge between mathematical models and phenomeno-logical observation, providing a means for validation, supplementing ana-lytical and experimental approaches. The role of computational methods inthe context of validation and prediction is summarized in Figure 1.1.Computational advances have been made through history, from the smallbeginnings of CFD, born of Lewis Fry Richardson in the early twentiethcentury with his efforts in numerical weather prediction. His concept ofsolving the simplified flow equations by a "forecast factory" [86], consist-ing of an arrangement of humans performing calculations, has the sameblueprints of a modern parallel computer. Notably, the many inadequaciesin this approach are concerns mostly still familiar in today’s state-of-the-art.The weather forecast simulations for the next few hours would have taken11.1. MotivationFigure 1.1: The role of computational methods in validation of models.days, highlighting the computational time considerations still prevalent to-day. Concerns of computational resources are also common: Richardson’scomputational thought experiment would have required tens of thousandsof people; the similar problem exists today with computer resources man-ifested in physical memory limitations. His particular choice of numericalscheme turned out to be unstable, and the lack of resolution in the inher-ently chaotic model equations were also issues that plagued his efforts. Nev-ertheless, Richardson’s work was of pioneering nature to CFD, even in theshortcomings which are still of important relevance. Through an under-standing in many of these issues as well as developments on techniquesin how to approach them, computational methods have become an indis-pensable tool nowadays in CFD, and also in many domains of science andengineering, marking many triumphs in recent decades.Different mathematical models have been developed and used in CFD.The current work focuses on the compressible Navier-Stokes equations, oneof the most widely used models, but there is also much investigation in theliterature on using other model equations such as computations based onthe Lattice Boltzmann models [106], which consider the collision of fluid21.1. Motivationparticles, with varying degrees of success.The areas to which this general methodology of modeling and computa-tion can be applied are vast and diverse, and not limited to the context ofCFD. In finance, the Black-Scholes equation [19] stochastically models thepricing of contracts on futures of an underlying asset, known as options. Inthis case, there are strong assumptions embedded into the model that fail tocapture observed financial markets in their entirety, such as government pol-icy and human psychological factors. While the model has been validatedwith benchmark data having consistent assumptions, there are some limita-tions in its use for real time or historical observations [40]. The Schrödingerequation [97] is a model in modern physics governing the wavefunction,commonly interpreted as the probability amplitude, whose numerical so-lutions have been found to be consistent with experimental observationon a fundamental level of particles [33]. If there were unlimited compu-tational resources, it would be possible to solve the Schrödinger equationfor every particle, and validation of the model would be possible with anyphenomenon of the observable universe. For macroscopic applications, thisis infeasible even with algorithms and hardware in the foreseeable future.These two examples highlight the role and characteristics of computationalmethods with modeling in that although there can be limitations of mod-els, computational results for the model equations can be useful in someinstances.Having motivated the important role of computational methods, thequantification of errors is an important but sometimes overlooked consid-eration. This is important since CFD practitioners are often interested inthe quality and fidelity of a particular computational solution [75], as wellas estimating if the current mesh is fine enough for a desired threshold ac-curacy. The accuracy level required in aerodynamic applications are quitestringent; for wind tunnel tests the error is often on the order of one percent[14]. In another instance, Vassberg et al. [110] found that even an error ofone percent in the prediction of drag is not only discernible, but can have asignificant impact on flight economics. Therefore, it is imperative to quan-tify and understand the errors that arise from CFD simulations, which is thefocus of this work.In each stage of the procedures in Figure 1.1, there are associated errors.The overall deviations of computed solutions from physical observations canbe attributed to two main sources of error: modeling error and numericalerror, with the former arising from the inadequacy of the derived mathe-matical equations to describe the underlying physics, and the latter arisingfrom the inability to solve these equations exactly. A more precise classi-31.1. Motivationfication of errors and the current scope is described in Section 1.3, but abrief discussion is given here. In many fluid flow applications, numericalerrors are often comparable in magnitude to modeling errors. Therefore, toaccount for the overall discrepancies between predictions and observations,one must first be able to accurately estimate and control the numerical errorin solving the modeled equations. Reliable estimates of the numerical errorscan give an assessment of the quality of the associated computed solution,much akin to error bars in an experimental study; a computed quantity with-out an estimate of error is weak at best. In addition to appraising solutionaccuracy, an accurate estimate of error is useful in improving the solutionthrough defect correction [104], as well as guiding mesh adaptation. Thisexploration is important because reliable and robust error estimators in thecontext of higher order finite-volume methods for unstructured methods areessential for real applications but lacking in the literature. The current workoffers a robust and efficient method to accurately estimate these errors thatmeasure how the computed solutions fail to satisfy the model equations.The approach is agnostic to the choice of model equation, hence versatilefor different applications in science and engineering.Since the advent of computational machines, development in hardwarehas seen an incredible rate of increase of processing power, doubling ap-proximately every 18 months in recent decades. This observation is oftencited as Moore’s law, attributed to his original work on trends of packingcomponents in integrated circuits [71]. As more computational resourcesbecome available, scientists and engineers are interested in solving largerproblems, rather than problems of the same size in a shorter time. Althougha formal analysis of the run-time complexities of our algorithms will notbe performed, the importance of algorithmic efficiency is motivated andemphasized. In the CFD community, the expectation is that problems withrealistic geometry and complex physics such as combustion will be common-place by the year 2030 with grid points on the order 10 to 100 billion [99],implying unknowns being in the trillions. In fact, problem sizes of this mag-nitude are being reached today, requiring an incredible amount of compu-tational resources [16]. It is evident that software efficiency will grow evermore important as problem size increases. To illustrate this point, for theconservative estimate of n = 1012 inputs, an improvement of a particular al-gorithm that takes computational time proportional to n2, to one that takescomputational time proportional to n log n, would far dominate advances inhardware components, requiring approximately 53 years of equivalent hard-ware improvement. This difference in magnitude only increases as problemsize increases. Furthermore, with evidence of a flattening trend in Moore’s41.2. Mesh Propertieslaw, and certainly the exponential increase becoming inevitably unsustain-able due to physical limitations, improvement in algorithms is extremelyimportant. Therefore, simply waiting for more powerful hardware to solvelarge problems will prove insufficient. Throughout the current work on er-ror estimation, emphasis on computational efficiency is hence warranted.If more computational resources, including processing units and mem-ory, are available at a given time, simulation times can be reduced throughparallel computing. This has been a focus in modern scientific computing,and has given rise to supercomputing clusters with many cores. Recently,some parallel applications have shown good scalability in the context ofCFD problems to a number of cores on the order of hundreds of thousandsto millions [9, 16]. The current work on error estimation is parallel-safein the sense that it preserves all of the parallel capabilities of the currentsolver, which uses the message passing interface (MPI) protocol, the ori-gins of which can be found in the literature [36]. As we will see, the cur-rent error estimation approach builds on the methodology used for the dis-cretization of the original partial differential equation (PDE) problem, whichmeans that any parallel efficiency improvements made on a fundamentallevel would be readily transferable.1.2 Mesh PropertiesMost discretization schemes approximate solutions on a computational mesh.For a domain Ω ⊂ Rd, where d is the number of spatial dimensions, amesh, or synonymously a grid, is a tessellation of its closure Ω¯ into regionsΩi, i = 1, ..N called cells or control volumes, whose interiors are pairwisedisjoint, and whose union is Ω¯. Throughout this thesis, closed domains areused so that Ω¯ = Ω.One of the main objectives of this work is to perform error estimation inthe context of general unstructured meshes. It is instructive, then, to elab-orate on the classification of meshes based on their topology and geometry.Examples of each mesh type are shown in Figure 1.2, which are categorizedas:• Structured meshesThis type of mesh can be identified by their regular connectivity. Thisproperty makes it natural to refer to individual cells by indices (i, j, k)in each of the independent directions. Further classification is possiblebased on geometrical features, as follows.51.2. Mesh Properties– Uniform meshesThe simplest and most restrictive subclassification of structuredmeshes are uniform meshes, where all cells have the same shapeand volume. Discretization schemes are simple to implement forthese meshes, with correspondingly smooth profiles of error.– Smooth meshesRelaxing the uniformity requirement, meshes with cells whosegeometrical features vary smoothly can be constructed. In Figure1.2b, the cell widths expand by constant multiplicative factor ineach direction. This implies that there exists a mapping functionbetween a reference uniform mesh and the actual smooth mesh.– Nonsmooth meshesThe most general type of structured mesh is classified as nons-mooth meshes, where geometric features do not vary smoothlyacross cells, but the regular connectivity is retained.• Unstructured MeshesHaving topologically regular and smoothly varying features is very re-strictive in meshing, especially in complex 3D geometries. Mesheswith general connectivity, known as unstructured meshes, relax therequirements for mesh generation. This is the most general notion ofa mesh that is considered, where connectivity as well as element typecan be arbitrary. As a result, cell connectivity information must beexplicitly stored to be used during computation. For such a generalmesh, a mapping to a reference mesh cannot be defined. Therefore,meshes with smoothly varying features are not possible with generaltopology.Our general approach for error estimation targets the nonsmooth class ofmeshes, a category in which unstructured meshes naturally fall, which alsocontains structured nonsmooth meshes. While solver implementation forstructured meshes is simpler, a larger amount of time for mesh generationis required for realistic geometries. It is often the case that mesh genera-tion time is a significant proportion of the overall time to obtain a solution,and should be an important consideration when evaluating discretizationmethods. Furthermore, without the geometric flexibility, the automationof generating even block structured meshes is extremely difficult, and thepossibility of generating a valid mesh in 3D for realistic geometries withoutexplicit user intervention becomes questionable.61.2. Mesh Properties(a) Uniform mesh (b) Smoothly varying mesh(c) Perturbed mesh (nonsmooth) (d) Unstructured mesh with general topol-ogyFigure 1.2: Classifications of meshes with quadrilateral cells. Triangular andmixed meshes can be analogously classified.71.3. Classification of ErrorEven though unstructured meshes alleviate some difficulties by remov-ing some of the constraints introduced by structured meshes, automaticmesh generation is a difficult problem in general, especially in 3D. Manyother classes of techniques are studied in the literature, including Octree,Advancing Front, and Delaunay. An overview of such mesh generation tech-niques can be found in [80]. Currently, meshes used for model problemsare generated by GRUMMP [77], our in-house Delaunay-based mesh gener-ation code. The unstructured framework also allows mixed element meshesnaturally. There are instances where quadrilateral cells are suitable, such assteep gradients in boundary layers, and others where simplicial (triangularin 2D) elements are suitable. For most model problems, triangular meshesare used, but the conclusions are unchanged when mixed element meshesare used in Chapter 4. Hence, we wish to establish an approach for effi-cient error estimation for general meshes, applicable broadly in practice. Tostudy the accuracy of solutions and error estimates, a sequence of mesheswith decreasing length scale is used. In the current context, we note that themeshes are generated with no natural nesting from one mesh to the next,meaning that cells are not simply subdivided for the next mesh resolution.A typical unstructured mesh used in these studies is shown in Figure 1.3.1.3 Classification of ErrorThere are many notions of error that will be relevant to our study. A briefdescription is given as follows.• Modeling Error.This is defined as the difference between the observation and the ex-act solution to the derived mathematical model equations. This canbe further subdivided into errors arising from the failure of the equa-tions to capture the observable physics, and errors arising from thediscrepancy of the model to capture the geometrical domain of theobservation with sufficient accuracy. This is still an open problem inmany areas of application. Specific to CFD applications, modeling tur-bulence as chaotic behavior on different scales of the Navier-Stokesequations poses one of the greatest challenges in the formulation ofmodel equations [69].• Discretization Error.This belongs to a category of error known as numerical error, whichwill be the focus of this work. Discretization error is defined as the81.3. Classification of Error0 101Figure 1.3: Typical unstructured mesh used on a square domain for thisstudy.91.3. Classification of Errordifference between the exact solution to the model equations and thesolution of the discretized approximate algebraic equations. This no-tion of error is what is meant typically when referring to error in thecomputed solution. The focus of the current work is to solve an aux-iliary equation, the error transport equation (ETE), which serves as amodel for the discretization error.• Truncation Error.The truncation error is defined as the difference between the exact op-erator applied to the exact solution and the discrete operator appliedto the exact solution. As derived later, the truncation error entersthe ETE as a source term, and needs to be accurately estimated. Forsmooth meshes the truncation error is asymptotically the same orderas the discretization error, but this is not true for nonsmooth meshes.This turns out to be a major concern, as the ETE models an asymptot-ically small discretization error, but the error made in evaluating theETE source term can be large. Note that the truncation error is thediscrete residual evaluated at the exact solution. The technical dis-tinction between continuous and discrete residual is seldom made inthe current context.• Iteration Error.Iterative methods are essential in solving the large systems of algebraicequations after discretization. The question of whether the iterationscheme has fully converged to a fixed point, within some specifiedtolerance, arises naturally. Iteration error is the lack of convergencein this sense of reaching the fixed point. There are many strategies toquantify this and it is a source of error that is relatively well studied[38, 39].• Round-off Error.Another relatively benign and well studied [46] source of error isround-off error, which is concerned with the preservation of preci-sion when arithmetic operations are performed on numbers stored ona computer. This propagation of error will not be of great concernthroughout this work.With this taxonomy of errors, some general remarks can be made. Modelingerrors are difficult to estimate, due to dependencies on measurement uncer-tainties, along with other forms of experimental error as well. Iteration and101.4. Approaches to Error Estimationround-off error are relatively well controlled and understood. We proceedunder the premise that a deeper understanding of discretization error anda framework to obtain a robust, accurate estimate is a necessary step toperform reliable validation to evaluate modeling errors.1.4 Approaches to Error EstimationHere we present an overview of existing methods for discretization errorestimation, highlighting the features and limitations of each approach, andprovide justification of the current ETE approach in the context of higherorder finite-volume methods.1.4.1 Method of Manufactured SolutionsTo measure the error, the easiest method is to subtract the computed solu-tion from the exact solution. In many cases the exact solution is not readilyavailable. However, different techniques are available to test the accuracy ofthe numerical scheme in such cases. The method of manufactured solutionsis one such method. This is done by creating a desired exact solution. Theanalytical source term that is required when the desired exact solution issubstituted in the model equations is computed, and used in the test case.This method is well known, and examples and descriptions can be found inthe literature [96].1.4.2 Richardson ExtrapolationThis classical method of Richardson extrapolation [85, 87] uses a sequenceof approximate solutions with varying resolutions to obtain an asymptot-ically better solution estimate and order properties. There are abundantexamples of its use in the literature [26, 76]. Romberg integration [90] isa classical application of this technique applied to evaluating a definite in-tegral, where repeated applications of the quadrature approximation usingdifferent resolutions can be combined for a more accurate estimate. Thiscan also be used to compute an error estimate on the coarsest mesh. If theorder of accuracy is known, then at least two solutions are needed for anextrapolated solution and error estimate. In the general case of an unknownorder of accuracy, at least three solutions are required. The details can befound in Appendix A, but the main result is that the order of accuracy p canbe obtained from solutions U1(h), U1(h2 ), U1(h4 ) on three meshes with the111.4. Approaches to Error Estimationlength scale parameter reduced each time. A convenient, but not necessaryapproach is to divide this by a factor of two each time, for which p can beobtained asp ≈ log[(U1(h2)− U1(h)) / (U1 (h4 )− U1 (h2 ))]log 2(1.1)which gives the extrapolated solutionU2(h) :=2p2p − 1U1(h2)− 12p − 1U1(h) (1.2)from which the error estimate on the coarsest mesh is ≈ 2p2p − 1(U1(h2)− U1(h)). (1.3)An advantage to this approach is its simplicity and robustness; only solu-tions on successively finer meshes need to be obtained, so that no intrusivemodification to the solver code is necessary. However, several assumptionsneed to be made when using such extrapolation methods, and this warrantsfurther examination of their validity for different applications. Perhaps themost severe drawback lies in the assumption that the asymptotic range hasbeen reached for the coarsest mesh; in practice there are many instanceswhere it is evident that the higher order terms are larger than the leadingorder terms for the given mesh resolution. In many realistic problems, thesuccessively finer meshes are not feasible to solve on due to computationalresource and time constraints. Furthermore, extra work is needed to applythis analysis for problems that exhibit oscillatory convergence behavior.1.4.3 Adjoint MethodsIn this subsection, a popular alternative approach for error estimation isdiscussed. However, a distinction first needs to be made about error quan-tities. The ETE approach is able to obtain an estimate of the discretizationerror, which is the error on the solution variables. Adjoint methods, onthe other hand, are formulated to yield error estimates on functionals ofthe solution. Studies in the literature on adjoint-based approaches for suchoutput-based error estimation [111], as well as aerodynamic optimization[61], have been available for some time. The basic ideas and key features ofusing adjoints for output error estimation is described here. For a more de-tailed description the reader is directed to [44]. To illustrate, consider the121.4. Approaches to Error EstimationCauchy problem for a diffusion partial differential equation (PDE), posedfor u(x, t), which vanishes on the boundary of a domain Ω× ΩT , as∂tu− ∂2xu = s. (1.4)Using L for the differential operator, this isLu = s. (1.5)Using the standard L2 inner product notation, a specific target functional Jcan be written asJ = (u, s∗), (1.6)and s∗ can be determined. The corresponding adjoint operator is L∗, definedas(w1,Lw2) = (w2,L∗w1) (1.7)for all elements w1 and w2 in the solution space. In particular(v,Lu) = (u,L∗v) (1.8)where the adjoint variable v satisfiesL∗v = s∗. (1.9)Applying integration by parts on the left side of Equation 1.8, we get(v,Lu) : =∫ΩT∫Ωv(∂tu− ∂2xu)dxdt=∫ΩT∫Ωu(−∂tv − ∂2xv) dxdt= (u,L∗v) (1.10)and the adjoint operator isL∗ = −∂t − ∂2x. (1.11)Therefore, the adjoint problem for this case, Equation 1.9, is not well posedif initial conditions are given and the solution is advanced to the final time.This is similar in nature to backwards diffusion models such as Black-Scholesequation. A standard von Neumann analysis [28] would show that smallperturbations to initial data do not stay small for modes with large wavenum-bers. Instead, it is only well posed if final conditions are imposed and the131.4. Approaches to Error Estimationequations are advanced backwards in time, which is the approach taken inconventional unsteady adjoint applications.From this overview, we see that there are some issues in using the ad-joint method for output error estimation. First, one adjoint problem needsto be formulated and solved for each functional of interest. For design prob-lems where many functionals are of interest, this would be computationallycostly. Furthermore, for problems with a large time domain, many inte-gration steps would be needed, and the entire solution history needs tobe stored for backwards integration of the adjoint problem, which can bememory prohibitive. There are more sophisticated treatments to reducestorage requirements [47], but a versatile and efficient framework for un-steady adjoints is still an open problem. There has been some recent work[107] which performs comparisons and presents results that reinforce thecompetitiveness of ETE-based methods against traditional adjoint methods.Throughout the present work, where the ETE approach is the focus, thereare some open questions which relate to implications for analogous accuracyorder properties for adjoint methods, but these will not be further exploredin this thesis.1.4.4 Error Transport EquationThe proposed approach addresses the objective of an accurate, robust dis-cretization error estimate, by means of solving the ETE, an auxiliary equa-tion, which can be solved together with the primal on the same grid, over-coming some of the drawbacks of alternatives such as Richardson extrapola-tion and unsteady adjoint methods. One of the key features of this methodis that an error estimate is computed for a single grid, removing the needfor additional cost in mesh generation, which is often a large proportion ofoverall simulation time. For unsteady problems, the ETE can be co-evolvedwith the primal equation to the final simulation time, requiring only localquantities in time, so that the solution history does not need to be stored.One desirable property of the ETE methodology, as demonstrated later, isthat it is not model-specific; in these studies we apply the ETE to compress-ible flow equations, but in principle this can be applied to other differentmathematical models as well.Another feature of the ETE approach is that it produces error estimatesof the solution variables, agnostic to the choice of functionals of the solu-tion. This is unlike the formulation of adjoint methods, where only errorsin functionals are obtained. The ETE approach is more flexible in the sensethat corrected functionals can be a postprocessing step using the solution141.4. Approaches to Error Estimationcorrected by the error estimate. Furthermore, if many functionals are ofinterest, equations for a new adjoint variable need to be formulated andsolved for every functional examined. In the ETE approach, there is onlyone auxiliary problem for the discretization error.Comparisons of different error estimation methods can be found in theliterature. In one instance [57], it was found that the ETE approach is astrong candidate in accurate error estimation against Richardson extrapo-lation based methods and another approach using solution reconstruction,over a range of subsonic and supersonic flow test cases.The development of ETE is quite mature for smooth meshes, charac-terized by the smoothness in the right hand side forcing term of the dis-cretized equations. Many studies have been conducted with success in ob-taining accurate error estimates in the literature [10, 25]. Recently, moresophisticated techniques for smooth right hand sides have brought signifi-cant advances such as higher order accurate corrections [29, 37]. Despitethis progress, the development of the ETE approach to obtain accurate errorestimates for nonsmooth unstructured meshes is lacking, and we attempt toaddress this research gap.1.4.5 Higher Order MethodsMethods that have a formal order of accuracy higher than two are referredto as higher order methods. With higher order discretizations, the error for-mally decreases at a faster rate as the mesh is refined. However, the pricefor the added accuracy needs to be paid with more computational work insolving the associated equations, as well as a possible compromise in stabil-ity properties. Nevertheless, there is convincing evidence in the literature[113] which suggests that in many applications, higher order methods areasymptotically more efficient in computational work and memory in obtain-ing accurate answers than their lower order counterparts.1.4.6 Quality of Error EstimatesWe measure the quality of an error estimate by comparing it with the exactdiscretization error. The difference between the exact discretization errorand the error estimate is referred to as the "error in error" throughout thisthesis. The objective is to devise discretization schemes that will result inan error in error that is asymptotically small. A common approach in theliterature to evaluate the quality of the measure uses what is known as the151.5. Objectiveseffectivity index [57]. This quantity is the ratio of the norm of the discretiza-tion error estimate to the norm of the exact discretization error, and usuallyonly its convergence to unity as the mesh is refined is examined. If such anapproach is taken, this only shows that the error in the error estimate goesto zero faster than the error itself. That is, if the exact discretization hasformal order of accuracy p, the effectivity index approaching one only saysthat the error in error has formal order of accuracy larger than p. We focus,rather, on obtaining higher order accurate error estimates, which gives moreasymptotic information. When used for defect correction, we show that ahigher order accurate error estimate is equivalent to a higher order accuratecorrected solution. This makes it asymptotically more efficient than lowerorder schemes for a given level of error [113], suggesting the existence of acrossover point beyond which the higher order method is more efficient.1.5 ObjectivesThe novelty of the current work can be summarized as follows:• Establishing a framework for solving the ETE on unstructured meshesfor steady and unsteady models of compressible fluid flow, and eval-uating to what extent strategies that exist in the literature for smoothmeshes can be used to obtain asymptotically higher order accuratediscretization error estimates.• Evaluating methods of computing the residual, which is needed as thesource term for ETE, that gives accurate results.• Exploring expected theoretical orders of accuracy for the error esti-mate, and confirming these results from computation.• Discerning what choices should be made in the discretization of theprimal problem, ETE, and residual, to get asymptotically higher orderaccurate and robust error estimates. A common comparison for accu-racy and robustness will be made between the corrected solution andthe solution to the higher order primal problem.• Investigate if it is possible to efficiently compute these error estimatesso that they can be used practically, preserving the parallel-safe prop-erties of the current discretizations and linear algebra operations.161.6. Outline1.6 OutlineThe thesis is arranged as follows. Chapter 2 discusses the methodologyfor both solving the primal problem and ETE that is followed throughout,highlighting details such as the ETE source term evaluation, and generaltreatments of initial and boundary conditions. Chapter 3 explores the dif-ficulties in 1D of using the ETE on nonsmooth meshes, which will persistas we move to multiple dimensions. Discretization requirements in obtain-ing higher order accurate error estimates will be demonstrated, which alsocarries over in the later chapters. Chapter 4 reinforces these concepts byextending to steady test cases in 2D, as well as some theoretical results ofthe expected accuracy of the error estimates using the ETE. More difficultmodel problems, including the Reynolds Averaged Navier-Stokes (RANS)equations with the Spalart-Allmaras turbulence model, will be considered.Chapter 5 extends the ETE methodology further to unsteady problems, andinvestigates the possibility of obtaining higher order error estimates in bothspace and time. Directions of future work and perspectives are given inChapter 6.17Chapter 2MethodologyWriting mathematical model equations as conservation laws encompassesa very broad classification of models that can be applied to many fields ofstudy, an overview of which is discussed in Appendix B. Conservation lawsare based on the observation that the rate of change of the total amount ofa certain quantity u(x, t) in a stationary spatial domain Ω is equal to the netamount crossing the boundary, combined with any sources or sinks s withinthe interior. This bookkeeping can be stated asddt∫Ωudx = −∮∂ΩF(u) · nˆdσ +∫Ωs(u)dx, (2.1)where F is the flux in the coordinate directions, ∂Ω is the boundary of thedomain, and nˆ is the outward unit normal on the boundary. This is some-times referred to as the integral form of the model equations. A schematicis shown in Figure 2.1. Applying the divergence theorem and rearranginggives ∫Ω(∂tu+∇ · F(u)− s(u)) dx = 0. (2.2)Since the domain in space is arbitrary, the integrand must be identicallyzero, which gives the differential form of the governing mathematical model∂tu+∇ · F(u) = s(u), (2.3)where the space-time domain is {x, t} ∈ Ω×ΩT ⊂ Rd×[0,∞), and the num-ber of spatial dimensions is d = 1, 2, or 3. A particular numerical methodaims to solve the model equations approximately, which involves discretelyenforcing one of these forms of the equations approximately over the com-putational domain. The equations for u will be referred to as the primalequations, since this is the quantity that is of interest. For a chosen numeri-cal method, there will be an associated approximate solution. Let u˜(x, t) bethe projection of this approximate solution onto the space to which the ex-act solution u(x, t) belongs. The corresponding discretization error is thendefined as18Chapter 2. MethodologyFigure 2.1: Conservation law in a domain.e(x, t) := u(x, t)− u˜(x, t), (2.4)and the objective is to derive an equation that governs the behavior of thisquantity. Solving for u in Equation 2.4 and substituting that into Equation2.3, we get∂t(e+ u˜) +∇ · F(e+ u˜) = s(e+ u˜). (2.5)Subtracting the quantity ∂tu˜ +∇ · F(u˜) − s(u˜) from both sides, then rear-ranging, this yields∂te+∇ · (F(e+ u˜)− F(u˜))− (s(e+ u˜)− s(u˜)) = −(∂tu˜+R(u˜)), (2.6)whereR(u˜) := ∇ · F(u˜)− s(u˜) (2.7)is defined as the spatial residual. Equation 2.6 is the continuous ETE, wherethe discretization error is governed by a differential operator that is in gen-eral different than that for the solution; these are only the same in the caseof a linear primal operator. Using the primal flux erroneously for the ETE ofnonlinear problems severely degrades the error estimate.The ETE is driven by a source term that is the time-dependent residualof the approximate primal solution. It is possible to show that this sourceterm corresponds to the truncation error if an analogous derivation for thediscretized primal equations is performed. Physically, Equation 2.6 says thatthe error depends on a local source term which is the residual, and possibly192.1. Discretization in Spacetransported from other regions of the domain as governed by the differen-tial operators on the left side of the ETE. The strategy is to solve the ETE,Equation 2.6, in conjunction with the primal equation, Equation 2.3, in or-der to get an approximate solution, along with an accompanying accurateerror estimate by appropriate choices of discretizations in space and time.2.1 Discretization in SpaceThere are several possible choices to be made when discretizing models inthe form of Equation 2.3 and Equation 2.6 in space. Several common ap-proaches are considered, and the current choice of the finite-volume methodwill be justified.The finite-difference method, where the derivatives in the differentialform are directly approximated at grid points, is a popular choice of dis-cretization due to its simplicity on uniform grids. However, there are defi-ciencies which make these methods quite problematic in practice for modelsencountered in CFD applications, even in obtaining higher order accuracythrough a reconstruction procedure. First, much effort is required to extendthe methodology to unstructured meshes, especially in multiple dimensions,thus diminishing the benefits of its simple discretization. Perhaps a moresevere problem is that the finite-difference approach solves for pointwisesolution values using the differential form of the equations, which tends tofail for finding weak solutions to nonlinear problems, because the deriva-tives are ill-defined when there is a discontinuity in the solution, knownas a shock. Furthermore, even if the shock can somehow be captured atan instant, due the lack of conservation built in to the scheme by default,the computation of shock speeds can be incorrect in many instances [65],leading to unusable results.The finite-volume discretization method mimics the global conservationargument in Equation 2.1 for each cell of the mesh, also known as a controlvolume. Instead of pointwise values, the discrete unknowns are the controlvolume averages. Unlike some other schemes, conservation is naturally en-forced, which is an important property in many physical applications. Thefinite-volume approach taken throughout this work reconstructs polynomi-als piecewise for each cell, without the constraint of enforcing continuityat interfaces, which has proved to be quite successful in CFD applications.Finite-volume methods have been traditionally successful in solving CFDmodel equations, encompassing developments in solving the Riemann prob-lems that naturally arise across the boundary of cells. In these situations,202.1. Discretization in Spacemethods based on the Godunov scheme [45] for satisfying entropy condi-tions for shock solutions, or the widely used Roe’s approximate Riemannsolver [89] used currently, have been thoroughly studied and applied to awide range of problems in the literature [4]. However, from an error esti-mation point of view, the availability of accurate and robust error estimatesfor unstructured meshes, either a priori or a posteriori, is lacking. The goalof the current work is to address the gaps in error estimation in the contextof unstructured finite-volume methods, which is well justified in light of therelevance of this class of discretization methods.There are also some strong competitors to finite-volume methods. Inrecent decades, there has been a proliferation of finite-element methodsin CFD applications, especially the discontinuous Galerkin (DG) family ofschemes. In finite-element methods, the integral form of the equations isenforced on a smaller space of functions, where enforcing the conditions onthe basis functions of that space is sufficient by the linearity properties. Thediscrete unknowns are the coordinates, in the linear algebra sense, of eachbasis function for the chosen solution space to approximate the actual solu-tion space. Finite-element methods in general have seen an abundance ofdevelopments in theoretical results, including error estimates, due to theirmathematical nature. The major features are highlighted here; more de-tailed comparisons of finite-volume and DG schemes can be found in theliterature [7, 74]. Similar to the current finite-volume method, one of thecharacteristics of DG methods is the removal of a continuity constraint be-tween elements. There are tradeoffs between DG and our finite-volumescheme. For instance, using a higher order discretization, the DG schemesrequires many degrees of freedom per cell, due to increased dimension ofthe basis for each element, whereas finite-volume only requires one un-known per cell. However, the number of cell neighbors required for a higherorder stencil in finite-volume methods is large in comparison to DG, whichmeans that the bandwidth of the discretization Jacobian can be wider, whichaffects the linear convergence of some iterative schemes. Some optimiza-tion can be performed in either case such as reordering, but certainly finite-volume methods have the advantage of requiring fewer degrees of freedomfor higher order discretizations on the same mesh. Since DG has more de-grees of freedom per cell, the number of nonzeros per degree of freedomfor the two schemes are similar. The fill of the Jacobian of the discretiza-tion for a typical Poisson problem in 2D is compared in Figure 2.2 betweenthe current finite-volume scheme and a DG scheme using FEniCS [6], a freesoftware package for finite-elements, both for second order and fourth orderaccuracy. The mesh that is used has 1,806 control volumes. In these cases,212.2. Finite-Volume Discretizationthe mesh unknowns were reordered using the reverse Cuthill-McKee (RCM)method [31, 43]. It can be seen that the degrees of freedom for higher orderDG is much larger than for the current higher order finite-volume method.Another reason to pursue error estimation techniques for finite-volumemethods is their ubiquity in commercial CFD software. Even though DG hasgained a vast amount of traction in CFD research codes, due to the lack ofconclusive and convincing evidence that it is consistently more efficient thanfinite-volume methods, commercial codes have not shown much enthusiasmto abandon traditional finite-volume methods. Recent to the time of writing,some estimates peg the proportion of commercial codes that use a finite-volume discretization as opposed to other schemes at about 80% [112].However, robust, accurate error estimation techniques for these schemes areoften not well studied, which makes the fidelity of many simulation resultsquestionable.These are the reasons that justify the current choice of error estimationin the context of higher order finite-volume methods. The ETE approach,however, is not specific to the finite-volume method, and results shown laterare convincing in showing that the ETE approach for nonsmooth meshes canbe equally applied to other discretization schemes.2.2 Finite-Volume DiscretizationWe now turn our attention to an overview of the finite-volume discretiza-tion. To proceed, the integral form of the equations is enforced for eachcontrol volume. This can be obtained by integrating the PDE Equation 2.3over a control volume Ωi, i = 1, ..., N , and then dividing by its measure |Ωi|,giving1|Ωi|∫Ωi∂tu dx+1|Ωi|∫Ωi∇ · F dx = 1|Ωi|∫Ωis dx. (2.8)The finite-volume formulation is obtained by applying the divergence theo-rem to Equation 2.8, from which the formdUidt+1|Ωi|∮∂ΩiF · nˆ dσ = Si (2.9)can be obtained, with Ui and Si denoting averages of the solution and sourceover the control volume i, asUi :=1|Ωi|∫Ωiudx, (2.10)222.2. Finite-Volume Discretization0 1000 2000 3000 4000 50000500100015002000250030003500400045005000nz = 58574(a) DG second order0 5000 10000 15000020004000600080001000012000140001600018000nz = 519158(b) DG fourth order0 500 1000 1500020040060080010001200140016001800nz = 17872(c) finite-volume second order0 500 1000 1500020040060080010001200140016001800nz = 60367(d) finite-volume fourth orderFigure 2.2: Discretization Jacobian fill comparison of the current finite-volume method and a DG scheme using FEniCS on the same mesh. Notethe difference in matrix size.232.2. Finite-Volume DiscretizationandSi :=1|Ωi|∫Ωisdx. (2.11)The discrete spatial residual is defined asRi(U) :=1|Ωi|∮∂ΩiF · nˆ dσ − Si, (2.12)where U is the collection of averages for all the control volumes. Theseare the unknown quantities to be solved for in finite-volume methods. Thequantity of Equation 2.12 is also known as the flux integral, including sourceterms, for finite-volume methods. Writing R as the collection of flux inte-grals for each control volume, the following system of ordinary differentialequations (ODEs) can then written for the averages:dUdt= −R(U). (2.13)Note that no numerical approximation has yet taken place. The manner inwhich the fluxes are computed from cell averages and integrated to evaluatethe right hand side of Equation 2.13 will dictate the spatial accuracy of thescheme. For steady problems, the discrete solution satisfies the nonlinearalgebraic systemR(U) = 0. (2.14)Practically, a common approach is to use pseudo time-stepping, also knownas pseudo-transient continuation [62], where Equation 2.13 is solved inpseudo time instead of Equation 2.14. This can be considered as a regu-larized root finding problem. For steady problems, this pseudo time vari-able is physically meaningless, so the time accuracy is irrelevant as a designrequirement. In the absence of such a time accuracy constraint, a popularchoice is to use the backward Euler method, which is written asUn+1 − Unk= −R(Un+1), (2.15)where k is the time step. Expanding the right side about time level n andneglecting terms second order and higher in time, this becomes(1k+∂R∂U)∣∣∣∣n (Un+1 − Un) = −R(Un) (2.16)where 1 is the identity matrix. Combining this with stable choices of timestep k that tend to infinity will reach convergence to the steady-state solu-tion. In the limit of large k, Equation 2.16 coincides with Newton’s method242.2. Finite-Volume Discretizationfor a nonlinear algebraic system. The interpretation for using pseudo-transientcontinuation can be justified as starting with an initial guess that can benavigated towards the solution, taking quite conservative steps far awayfrom the solution. Near convergence, larger time steps are taken, and theNewton-like iterations accelerate the trajectory to steady-state with a near-quadratic rate. In unsteady problems, we advance Equation 2.13 by se-lecting a suitable time integration scheme that is appropriate for both timeaccuracy and stability, as discussed later. This manner of discretizing first inspace, then in time is known as a method-of-lines approach.For practical applications, especially in multiple dimensions, these linearsystems that arise from the spatial discretization are sparse, but very largein size. Due to the arbitrary topology of unstructured meshes, the matrix isin general not symmetric. A Krylov subspace method known as generalizedminimal residual (GMRES) can be used to efficiently solve the linear sys-tems [95]. Throughout, the linear algebra routines are handled externallythrough PETSc [101]. The discretization Jacobian for unstructured meshesdoes not have much structure to exploit in general, such as symmetry ordefiniteness, so that GMRES is used instead of methods such as conjugategradient [52], BiCGStab [108], or MINRES [81].The ETE, Equation 2.6, can be discretized in space in the same way asthe primal equation, asddt= −(R(U + )−R(U)) + τ(U) (2.17)where the terms on the right make up the residual for the ETE, wherenow U is the approximate solution, the space-discrete version of the dis-cretization error e, and τ the space-discrete version of the ETE source term−(∂tu˜ + R(u˜)). This can be solved numerically in the same way as theprimal equation, whether by pseudo-transient continuation for steady-statesolutions or an appropriate time accurate scheme for unsteady problems.For steady problems, the ETE is solved after the primal problem, and onlythe spatial part remains in the source term. In the unsteady case, since theETE at every time step depends on the time-dependent residual of the solu-tion, co-advancing both sets of equations removes the memory requirementof storing the full solution history in time, which would be needed if theETE were solved after the primal problem. The pseudo-transient approachfor the ETE source term can of course be taken for steady problems, but thiswould require more computational resources, and any stability benefits overusing only the steady state solution would remain to be seen.252.3. Order of Accuracy2.3 Order of AccuracySince we are interested in characterizing and reducing discretization errors,it is instructive to discuss what metrics can be used as an assessment. Forexample, one such measure would be to predict how the error in some quan-tity will behave asymptotically as finer meshes are used. The spatial lengthscale for a mesh is defined in terms of the number of cells N in the mesh, ash := N−1/d, (2.18)where d is the number of spatial dimensions. Here it is assumed that somemeasure on cell variation remains fixed, so that Equation 2.18 is valid. Weopt to use the standard big-oh notation, where|||| = O(hp) (2.19)means that there exist positive numbers C and δ such that || ≤ Chp when-ever 0 < h < δ. Loosely, this gives an asymptotic bound on the behavior ofthe error.1 This means that the leading order term is the one that dominatesas the length scale approaches zero. The norm in Equation 2.19 is taken tobe the typical L1, L2, and L∞ norms defined on grid vectors, as||||Lp =(N∑i=1|i|phdi)1/p, (2.20)which are analogous to the continuous form||e(x)||Lp =(∫Ω|e(x)|pdx)1/p. (2.21)In the order of accuracy studies, we opt to use the L2 norm as a measureof the size of the error to strike a balance between being able to detectasymptotically larger areas on small regions of the mesh, and being overlystringent. Using the other norms to measure error was found not to alterthe conclusions materially. A discussion of this can be found in [66]. Typicalof asymptotic analyses, we are mindful of the fact that although this givesan indication of how the error behaves as the small parameter is reduced,little information is known for a given magnitude. In general, it can wellbe the case that the mesh under testing is not yet in the asymptotic region.1Note that the equality notation does not constitute an equivalence relation, since h2 =O(h) but h 6= O(h2) shows that it is not symmetric.262.3. Order of AccuracyNevertheless, this analysis is still useful in practice, because as more com-putational resources are available, larger problems, or equivalently smallerlength scales, can be solved, hence the order of accuracy gives a good senseof the how much return can be expected for the allocation of more compu-tational resources.In these studies, the order of accuracy is determined by a sequence ofmeshes with successively smaller length scales and time steps, which givesa good indication of the error trend. It is noted that the sequence of mesheswith different spatial resolutions are not simply subdivided in a nested man-ner, which is investigated in Section 3.4.1. Rather, the length scale used forthe Delaunay-based mesh generation process is reduced by a constant factorwith each mesh resolution. This is more general than the nested refinementapproach, in that the accuracy assessment can just as well apply to a se-quence of meshes generated by another algorithm. Nevertheless, the resultsdemonstrate that the current methodology gives the desired orders of accu-racy with the sequence of unnested meshes. More sophisticated statisticalstudies can be done as well in general for unstructured meshes, consideringthe possible variation in mesh quality of the different resolutions. However,evidence suggests that the sampled approach is in fact quite representative[116].2.3.1 Error AccuracyIf an exact solution is available, its projection onto the discrete space, U ,can be used to measure the error in the discrete solution. For a particulardiscrete solution Up to the primal problem that is accurate to order p, let theexact discretization error on the mesh be p, and we wish to estimate thisby ˜p. If we apply the error estimate as a defect correction, then ensuringhigher order accuracy in the corrected solution is equivalent to higher orderaccuracy of the error estimate, becauseU − (Up + ˜p) = (U − Up)− ˜p = p − ˜p. (2.22)In the results to follow, we measure the order of accuracy of the error es-timate by either computing p − ˜p or the corrected solution U − (Up + ˜p)at the final time as the spatial length scale h and time step k are refined,keeping k/h roughly fixed for each case. The asymptotic accuracy is easilyvisualized when the error and the small parameter are plotted on a log-logplot, where the power law relationship can be readily seen as the slopes.If an exact solution is not available from the model equations, an exact272.3. Order of Accuracysolution can be constructed for a particular test problem using the methodof manufactured solutions, where the resulting source term from imposingthe desired exact solution is added to the model equations. Examples of thiscan be found in the literature [96].If there is no exact solution and one does not wish to use the methodof manufactured solutions, solutions on at least three different mesh resolu-tions that are asymptotic can be used for error estimation and determinationof the order of accuracy using Richardson extrapolation [24]. These mea-surements of accuracy can also be performed in the same way on functionalsof the solution, such as the coefficient of lift or drag, which are often quan-tities that are reported in reference data.2.3.2 Higher Order Accuracy via k-Exact ReconstructionThe approach that is taken currently is part of a class of finite-volume meth-ods based on reconstructing a piecewise high degree polynomial for thesolution in order to compute the fluxes at the control volume boundariesaccurately. A method known as k-exact least-squares reconstruction is used.A detailed discussion is presented in the literature [11, 78], and some the-oretical aspects of these finite-volume schemes are described by [12]. Thefeature of this reconstruction procedure is that polynomials of degree p− 1can be reconstructed exactly from control volume average data, and thereconstruction is order p accurate for general smooth functions. The proce-dure for reconstructing one solution component in the 1D case is outlinedhere. The application to multidimensional systems is straightforward. For adesired order of accuracy p, a piecewise polynomial of degree p−1 is soughtwithin control volume i. A Taylor series expansion of the reconstruction incontrol volume i, taken about its reference point, usually at its centroid, iswritten asui(x) = ui + (x− xi) · ∇ui + 12!(x− xi)2 : ∇∇ui + ...+1(p− 1)!(x− xi)p−1[·]p−1∇p−1ui=p−1∑m=01m!(x− xi)m[·]m∇mui (2.23)with unknown coefficients{1m! ∇mu i}. The notation is chosen such that(x−xi)m is a tensor of rank m and [·]m is an extension of the inner productcontraction operator that returns a scalar from two tensors of rank m. The282.3. Order of Accuracycoefficients are determined by enforcing the mean as a constraint in controlvolume i exactly, and predicting the mean in n neighboring control volumesin the stencil of i in a least-squares sense, with n ≥ p. Practically the neigh-bors are added one layer at a time, until a sufficient number is reached forthe overdetermined least-squares system.The mean in control volume i should be preserved by the reconstructionpolynomial in control volume i. This is expressed asui(x) :=1|Ωi|∫Ωiui(x)dx (2.24)= ui + xi · ∇ui + ...+ 1(p− 1)!xp−1i [·]p−1∇p−1ui (2.25)withxmi :=1|Ωi|∫Ωi(x− xi)mdx (2.26)defined as the control volume moments. If the cell centroids are chosen asthe reference points, then xi = 0.The reconstruction in control volume i should be able to predict theaverages in another control volume j in its stencil reasonably accurately.This quantity isui,j(x) :=1|Ωj |∫Ωjui(x)dx, (2.27)and is computed asui,j(x) = ui +1|Ωj |∫Ωj(x− xi) · ∇uidx+ ...+1(p− 1)!|Ωj |∫Ωj(x− xi)p−1[·]p−1∇p−1uidx= ui +1|Ωj |∫Ωj(x− xj + xj − xi) · ∇uidx+ ...+1(p− 1)!|Ωj |∫Ωj(x− xj + xj − xi)p−1[·]p−1∇p−1uidx= ui + (xj + xj − xi) · ∇ui + 1(p− 1)!|Ωj | · ...∫Ωjp−1∑m=0(p− 1m)(x− xj)p−1−m(xj − xi)m[·]p−1∇p−1uidx292.3. Order of Accuracy= ui + x̂ij · ∇ui + ...+ 1(p− 1)! x̂p−1ij [·]p−1∇p−1ui (2.28)withx̂p−1ij :=p−1∑m=0(p− 1)!m!(p− 1−m)!xp−1−mj (xj − xi)m (2.29)defined as generalized moments. These geometric terms can be precom-puted for a particular mesh.For a topologically structured triangular mesh having three face neigh-bors, the required reconstruction stencil for different orders of accuracy isshown in Figure 2.3. The resulting linear system symbolically takes the form1 xi x2i ...wi1 wi1x̂i1 wi1x̂2i1 ...wi2 wi2x̂i2 wi2x̂2i2 ............win winx̂in winx̂2in ...u∇u12!∇∇u...i=uiwi1uiwi2ui...winui (2.30)withwij :=1||xj − xi||β , β ≥ 0 (2.31)as weights to emphasize nearby control volumes. The first row of Equation2.30 is the mean constraint. This is in the form of a standard constrainedweighted least-squares problem for the unknown coefficients.After reconstructing a piecewise high degree polynomial in each controlvolume by this procedure, the fluxes can be evaluated at control volumeboundaries using values of the solution and derivatives, and integrated overthe boundary using Gauss quadrature. Since continuity of the reconstructionacross faces is not enforced, for a particular point on a face, there are dif-ferent values of the solution (UL and UR) and derivatives (∇UL and ∇UR)from the reconstructions on either side. A unique flux is evaluated as afunction of these two sets of values through what is known as the numericalflux function, a schematic of which is shown in Figure 2.4. In general, theflux integral will depend on the fluxes on both sides of a face, which meansthat the flux integral depends on the reconstruction stencil of the cell’s firstneighbors. For advective fluxes, the numerical flux function is computed byupwinding via Roe’s flux difference splitting scheme [89],302.3. Order of AccuracyXY0.4 0.45 0.5 0.55 0.6 0.650.30.350.40.450.522 23 33333444444444Figure 2.3: Required reconstruction stencil for various orders of accuracy.Neighbors are added one layer at a time until a sufficient number is obtainedfor the least-squares system.Figure 2.4: Schematic of second order interface flux evaluation based on anumerical flux function of contributions from either side. For second orderaccuracy the Gauss quadrature point xq is at the midpoint of the face.312.3. Order of AccuracyFadv =F (UL) + F (UR)2+12|A˜|(UL − UR) (2.32)with |A˜| = P−1|Λ|P , P the matrix of eigenvectors, and |Λ| = diag{|λj |} theeigenvalues of the flux Jacobian A. In a similar way, for diffusive fluxes,an arithmetic average of the values on either side is used, with additionaldamping known as the jump term [73]. The total flux at the face is thenFdiff =F (UL) + F (UR)2+ αUL − UR|rij · nˆ| (2.33)where α is the constant jump term, rij is the vector connecting the referencepoints of the control volumes i and j sharing the face, nˆ is the outward unitnormal of that face. The jump term couples adjacent control volumes andprovides damping to stabilize odd-even modes, which is necessary for thestability of second order schemes for pure diffusion problems. For purediffusion problems, α = 1 is chosen. With this choice, for a second orderscheme on a uniform grid, the discretization coincides with the standardsecond order stencil. In 1D, this is a tridiagonal discretization matrix insteadof pentadiagonal for other choices of α.2.3.3 Residual EvaluationThere are several choices at hand to discretize the primal equation and ETE.We can choose separately the discretization order of the primal equation p,the discretization order of the ETE q, and the order r to which the primalsolution is reconstructed for residual evaluation. For steady problems, wewill refer to the discretization orders by the ordered triple (p, q, r). We alsodenote schemes solving only the primal problem by (p, 0, 0).The continuous residual R, which is the source term for the ETE, isdefined as the primal differential operator applied to the discrete solutionmapped onto the continuous space, subtracting physical source terms. Thisquantity is well defined, but obtaining an accurate estimate of its discreteanalog, the truncation error, is not trivial for unstructured meshes. Hay andVisonneau [50] discussed several approaches for mapped Cartesian mesheswith embedded smoothness, along with applications in CFD such as mod-eling turbulent flows [49]. However, obtaining accurate estimates to usein the ETE for general unstructured meshes has not been well studied inthe literature. An obvious method that works well for smooth meshes isto compute the residual by first projecting the converged primal solutiononto a higher order discrete space, a method known as p-truncation error322.3. Order of Accuracyestimation. In the current context, this translates to taking the convergedsolution of control volume averages of the primal problem, and performinga least-squares reconstruction to order r. This yields a piecewise polynomialof degree r−1 within each control volume, from which the flux integral canbe evaluated. It was found in the present study that if the ETE source termestimate is within O(hr) of the exact ETE source term, this is sufficient forthe discretization error estimate to be within O(hr) of the exact error also.However, it can be shown that using a reconstruction to order r of the lowerorder solution does not give an ETE source term that is accurate toO(hr) fornonsmooth meshes, and no other methods are known to be able to estimatethis accurately.For finite-volume methods on unstructured meshes whose geometry doesnot vary smoothly, the truncation error is asymptotically larger than the dis-cretization error and is inherently not smooth. Furthermore, estimation ofthis by means of higher order reconstruction of the primal solution doesnot produce an accurate estimate, whereas this is not a problem for uni-form or smooth meshes. To illustrate this, the exact truncation error usinga second order scheme and its estimates for a periodic Poisson problem in1D are shown in Figure 2.5. Estimates are obtained by a fourth order re-construction of the second order primal solution. This is performed usinga uniform and a nonsmooth mesh, which is generated by random pertur-bations of grid nodes proportional to a fraction of the nominal cell size, sothat no cells overlap. The truncation error on the smooth mesh is smoothand asymptotically small, whereas the truncation error on the nonsmoothmesh is O(1) and noisy. Furthermore, increasing the order of reconstructionof the primal solution in the p-truncation error estimate does not asymptot-ically give better results; the high amplitude noise components cannot bewell represented as the mesh is refined.One interpretation of the success of solving for the error estimate as acorrection for smooth meshes is that if the primal solution is accurate, thecorresponding residual will be small, and the correction should in turn besmall, which leads to a corrected solution that is not contaminated whenthe correction is performed. For nonsmooth meshes, the observation is thataccuracy in the primal solution, characterized by small discretization error,does not necessarily translate to a small residual using the p-truncation er-ror estimate. Nevertheless, we use r as the reconstruction order for residualevaluation and recognize that this is not the order of accuracy of the trunca-tion error estimate as in the literature for structured meshes. This method ofp-truncation error estimation is a common approach for uniform or smoothmeshes, as is an analogous h-truncation error in which projection of the pri-332.4. Discretization in Timemal solution onto a finer mesh is used for the estimation of the truncationerror on the coarse mesh. The technique of [41], known as τ -estimation,cannot be readily used in this context since the nonsmooth nature of thetruncation error violates many of the smoothness assumptions there. In-deed, the authors of that work recognize that for the randomly perturbed"highly distorted" meshes considered here, it is difficult to estimate trunca-tion error accurately.Proceeding in this way, values of p, q, and r are left free to comparethe various combinations. For smooth meshes considered in [50, 122], anycombination with the constraint r > p seems to give good error estimates.It can be shown also that r > p is necessary for accurate error estimatesusing unstructured meshes, for the same reason shown in the literature foruniform meshes: the residual source term will not be asymptotically correctotherwise. However, as we have seen for unstructured meshes, this condi-tion is not sufficient, because the ETE source term cannot be accurately es-timated using a higher order reconstruction of the primal solution. It turnsout that most combinations of p, q, and r give poor error estimates, but wewill investigate which particular combinations lead to good error estimates.2.4 Discretization in TimeThe time derivative is discretized separately from the spatial discretizationas a method-of-lines approach. Other methods which are possible such asspace-time formulations, where time is discretized as an extra space dimen-sion, can be found in the literature [56]. Under some circumstances, space-time formulations are found to be equivalent to a method-of-lines treatment.The current work focuses on using uniform or smoothly varying time steps,in hopes of exploiting the superior accuracy properties of the error estimatedue to smooth discretizations. This is justified in the time discretization,since there is usually little reason to use nonsmooth time steps, unlike inspace where there are often stringent geometrical constraints in multiple di-mensions when discretizing the domain. By discretizing in this way, smooth-ness in the error attributed to the time discretization can be preserved.In Equation 2.13 and Equation 2.17, we write the ODE systems asdUdt= f(U) (2.34)ddt= g() + T (2.35)342.4. Discretization in Time0 0.2 0.4 0.6 0.8 1−0.5−0.4−0.3−0.2−0.100.10.20.30.40.5xτ,Uniform ExactR4(U2)(a) Uniform mesh truncation error estimate0 0.2 0.4 0.6 0.8 1−8−6−4−202468 x 10−3xτ2−R4(U2)(b) Difference for uniform mesh0 0.2 0.4 0.6 0.8 1−15−10−5051015xτ,Perturbed ExactR4(U2)(c) Perturbed mesh truncation error estimate0 0.2 0.4 0.6 0.8 1−15−10−5051015xτ2−R4(U2)(d) Difference for perturbed meshFigure 2.5: Exact truncation error and estimate using a higher order discreteoperator.352.4. Discretization in Timewhere f and g represent the flux integral, and T is an estimate of the ETEsource term. The time discretization order is denoted as pt and qt for theprimal problem and ETE, respectively. We demonstrate the methodologyfor two classical implicit Runge-Kutta type schemes for time integration:backward Euler (pt = 1) and Crank-Nicolson (pt = 2). Using uniform timesteps k, the backward Euler discretization isUn+1 − Unk= f(Un+1) (2.36)n+1 − nk= g(n+1) + Tn+1 (2.37)and the Crank-Nicolson discretization isUn+1 − Unk=12(f(Un+1) + f(Un)) (2.38)n+1 − nk=12(g(n+1) + g(n)) +12(Tn+1 + Tn) (2.39)where the superscript indicates that the quantity is evaluated at that timelevel.Our objective is to pick Tn, n = 0, 1, . . . NT , to estimate τ in Equation2.17 so that the corrected solution V n = Un + n is closer than Un is to Un,the exact continuous solution projected onto the discrete space. With f = fp,and g = gq = fq(U + ·) − fq(U) as the right hand sides corresponding toorder p and q, adding Equation 2.36 and Equation 2.37 gives the correctedsolution asV n+1 − V nk= fq(Vn+1) + fp(Un+1)− fq(Un+1) + Tn+1 (2.40)for backward Euler, and adding Equation 2.38 and Equation 2.39 givesV n+1 − V nk=12(fq(Vn+1) + fq(Vn) + fp(Un+1) + fp(Un)−fq(Un+1)− fq(Un) + Tn+1 + Tn) (2.41)for Crank-Nicolson.For unsteady problem, the accuracy of solution and error estimates hasboth time and space contributions, depending on the spatial discretization362.5. Initial Conditionsorder p and temporal discretization order pt. For the primal problem it willbe in the form = O(hp) +O(kpt). (2.42)It will typically be the case that these are chosen so that the terms are asymp-totically of equal size. For example, if p = pt = 2, then the accuracy studywould usually proceed by taking k/h fixed as the meshes are refined. Thereis also the concern that some of the terms dominate over others, hence sucha refinement study may not measure the behavior of all the terms. Thiswas found not to be the case in the results; the space and time terms werequite comparable in order of magnitude. In these asymptotic expansions,note that more exotic forms, such as terms involving log(h), are not con-sidered, the justification being that our discretization schemes are based onTaylor series expansions, and our small parameter does not appear in themodel equations which would commonly be the source of non-polynomialexpansions that arise such as in singular perturbation problems.2.4.1 Residual EvaluationThe continuous time-dependent residual is −(∂tu˜ + R(u˜)), which is thesource term for the ETE. Having discussed the evaluation of the spatial partof the residual using a higher order reconstruction to order r of the primalsolution, we now examine the techniques for computing the time-dependentpart of the residual. The order of accuracy in time of the residual will bedenoted by rt. Then the discretization orders can be altogether specified bythe ordered 6-tuple (p, q, r; pt, qt, rt). The combinations of these values thatgives accurate error estimates will be investigated. For example, we maywish to examine if (2, 4, 4; 2, 2, 4) gives fourth order accuracy in both spaceand time using uniform time steps. This is solving the primal problem tosecond order in both space and time, and solving the ETE fourth order inspace and second order in time, combined with evaluating the ETE sourceterm using a fourth order reconstruction in space and fourth order in time.2.5 Initial ConditionsFor given initial conditions, a cell-weighted average of the specified valuesis used at the initial discrete step. The initial condition for the ETE is takento be 0 = U0 − U0 = 0, which is reasonable considering that the primalproblem uses the averages of the exact initial condition, or U0 = U0. For372.6. Boundary Conditionssteady problems, if the converged solution is unique, it will be independentof the initial condition.2.6 Boundary ConditionsFor boundary conditions specifying fluxes, those values are used at quadra-ture points when integrating over the boundary faces. When applying Dirich-let conditions, values are enforced as constraints in the reconstruction pro-cedure. For linear advection, at the inflow boundary the velocity is specified,and at the outflow boundary the reconstructed velocity from the interior isused. In addition, wall conditions are enforced by imposing zero velocity(hence zero flux). For the Navier-Stokes equations, characteristic boundaryconditions are used for inflow and outflow, similar to approaches found inthe literature [53]. This is enforced by using the value reconstructed fromthe interior if the characteristic of that physical variable comes from the in-terior, and specifying its value if the characteristic comes from the exterior.For unsteady problems, the boundary condition enforcement is performedin the same way. There are cases in the compressible flow models wherenon-reflecting boundary conditions need to be used to suppress oscillationsthat are not physical. Conventional approaches are used whose detailed de-scription can be found in the literature [83]. For the boundary conditionsin the ETE, they are enforced by the same type of boundary condition asin the primal problem, with the nonhomogeneous value replaced by 0. Forexample, if a nonhomogeneous Dirichlet condition U = UB is enforced inthe primal solution, then = 0 is enforced in the same way in the errorvariable.38Chapter 3ETE Accuracy on Smooth andNonsmooth MeshesThe purpose of this chapter is first to demonstrate the ETE methodology ondifferent types of meshes in 1D. The comparison to be made is the asymp-totic accuracy of the primal solution, U − Up, and the accuracy of the errorestimate, p− ˜p. The ETE methodology for smooth meshes is studied in theliterature [25, 50], as well as a range of other error estimation techniques[91]. The current focus will be on using the ETE with a p-truncation errorsource term, computed as a higher order operator applied to the primal so-lution. Other attempts at truncation error estimation, such as h-truncationerror estimation, where a finer mesh discrete operator is used, will be exam-ined later in this chapter. The asymptotic order of the error in error estimateusing the ETE on uniform meshes has been studied [10] and found to be||p − ˜p|| = O(hmin(p+q,r)). (3.1)A typical choice, then, is p = q = 2, r = 4, so that the error accuracy,hence the corrected solution accuracy, will be fourth order accurate fromtwo second order problems. This will be first verified using uniform meshes,and the extent to which this can be carried over from smooth meshes togeneral perturbed meshes will be examined. This will emulate the natureof general nonsmooth unstructured meshes one expects to see in multipledimensions by making apparent the difficulties of using nonsmooth meshes.Errors are dependent, in general, on the solution and its derivatives on themesh as well as the underlying function that maps the mesh to a uniformone. On nonsmooth perturbed meshes, some of these quantities may not bewell defined, highlighting the difficulties of using such meshes.393.1. Spectral PropertiesUniform Perturbedmax |λ| of A 1/h2 1/h2max |σ| of A 1/h2 1/h2max |λ| of A−1 1 1max |σ| of A−1 1 1Table 3.1: Empirical asymptotic sizes of the maximum eigenvalues and sin-gular values of the discretization Jacobian and its inverse.3.1 Spectral PropertiesThe Poisson equation is considered first for verification, asd2udx2= −4pi2 sin 2pix, Ω = {x ∈ [0, 1]} (3.2)with boundary conditionsu(0) = u(1) = 0. (3.3)The exact solution isu(x) = sin 2pix. (3.4)As discussed, the truncation error on nonsmooth meshes are asymptoticallylarge even for schemes where the discretization error is second order. Theapparent contradiction between the observed asymptotically large trunca-tion error and second order discretization error is discussed here in moredetail using some insight from the discretization Jacobian and its spectrum.For this linear problem the primal solution satisfies the discretized equationAU = S, (3.5)where A := ∂R∂U is the discretization matrix. Using = U − U , this becomesA = −τ. (3.6)Note that the classical stability analysis of the primal discretization schemeonly requires ||A−1|| to be bounded by a constant independent of h. This isinvestigated further, using information about how the maximum and mini-mum magnitude eigenvalues and singular values of A scale asymptoticallywith h empirically, which are tabulated in Table 3.1. Next, consider the403.2. Uniform MeshesUniform Perturbed||||2 h2 h2||τ ||2 h2 1Table 3.2: Asymptotic sizes for discretization error and truncation error.induced matrix 2-norm||A||2 = sup ||Az||2||z||2 . (3.7)From this we deduce||τ ||2||||2 =||A||2||||2 ≤ sup||Az||2||z||2 = ||A||2. (3.8)Through a singular value decomposition, for the matrix norm we have also||A||2 = ||UΣVT ||2 = ||Σ||2 = max |σ| = O(1h2). (3.9)Combining Equation 3.8 and Equation 3.9, the condition||τ ||2||||2 = O(1h2)(3.10)follows. By the same argument for A−1, the result is that||||2||τ ||2 = O(1). (3.11)The asymptotic behavior of the norms of discretization error and truncationerror for uniform and perturbed meshes is shown in Table 3.2. It is evidentthat neither of the conditions Equation 3.10 or Equation 3.11 are violatedfor the uniform as well as the perturbed meshes; the observed truncationerror for uniform meshes is smaller than the guaranteed asymptotic boundsdue to error cancellation.3.2 Uniform MeshesFor N = 40 and a (2, 2, 4) scheme, the exact discretization error, estimateusing the ETE, and the difference is shown in Figure 3.1. It can be seen413.2. Uniform Meshes0 0.2 0.4 0.6 0.8 1−2.5−2−1.5−1−0.500.511.522.5 x 10−3xDiscretization error Exact(2,2,4) estimate(a) Computed and exact error estimate0 0.2 0.4 0.6 0.8 1−3−2−10123 x 10−5Error in errorx(b) Error in error estimateFigure 3.1: Accuracy of ETE estimate for uniform (2, 2, 4), N = 40.that the error estimate is extremely close to the exact discretization error,and the difference is about two orders of magnitude smaller than the erroritself on this mesh resolution. The order of accuracy as N is varied is shownin Figure 3.2, which demonstrates that the error in error is close to fourthorder as expected. For these grid studies on the accuracy of error and er-ror estimates throughout this thesis, the primal schemes (p, 0, 0) will showthe exact discretization error along the vertical axis, while the primal andETE scheme (p, q, r) will show the error in error. The two can be simultane-ously labelled as the corrected discretization error, due to the relationship ofEquation 2.22. In this way, the vertical distance between these two sets ofdata, which widens as finer meshes are taken, can be interpreted as the im-provement in accuracy that would be obtained from the second order primalsolution if the error estimate is applied as a correction by solving anotherproblem that is the ETE. This is consistent with the relation in Equation 3.1for the accuracy of the error estimate using discretization order p, q, and r.Owing to the smoothness in the mesh, the truncation error is also smooth,and more sophisticated improvements can be made in reducing the compu-tational work even further. One such approach is shown here, which takesinspiration from what are known as multigrid schemes [21], which involvesprojecting the problem to nested coarser meshes, which has the effect ofdramatically reducing the computational work per iteration, as well as re-quiring fewer iterations due to the altered spectrum on the coarser meshes.A schematic of this is shown in Figure 3.3, and the fine and coarse grid num-423.2. Uniform MeshesCorrecteddiscretizationerror10−3 10−2 10−110−910−810−710−610−510−410−310−210−1h (2,0,0)(2,2,4)second orderfourth orderFigure 3.2: Discretization error and error in error estimate for (2, 2, 4) usinga uniform mesh.433.3. Smooth Meshesbering is shown in Figure 3.4. In the current context, the computational costcan be reduced by projecting the residual computed on the fine mesh ontothe coarse mesh, and solving the ETE on the coarse mesh to obtain the er-ror estimate, and finally projecting the error estimate back to the fine mesh.The fine to coarse mesh operator is a simple injection, and the coarse to fineoperator is a second order accurate interpolation, given byUf,2i−1 =34Uc,i +14Uc,i−1Uf,2i =34Uc,i +14Uc,i+1 (3.12)As in the traditional multigrid approach, multiple coarse mesh levels canbe combined together. Results for one and two levels are shown in Figure3.5. For a linear problem, the fine mesh (2, 2, 4) scheme obtains the errorestimate which can be represented as224 = A−12 (−R)= (−A−12 A4 + 1)U2 (3.13)where pqr is the error estimate resulting from the (p, q, r) scheme, Ap is theJacobian of the pth order discretization, and Up is the pth order discrete so-lution. Using one coarse level with multigrid, the error estimate is computedas224 = Ih2hA−12c (−I2hh R)= (−Ih2hA−12c I2hh A4 + Ih2hA−12c I2hh A2)U2 (3.14)where Ih2h1 is the transfer operator from a mesh with length scale h1 to a meshwith length scale h2, chosen to preserve the nominal order of the scheme.More accurate and more complex transfer operators could be used, but forthe uniform mesh these were not needed.3.3 Smooth MeshesNext, nonuniform smooth meshes are examined. This is taken to mean thatthere exists an underlying smooth grid mapping ξ(x), in the limit of zerolength scale, from the given mesh to a uniform mesh, usually referred toas physical space and computational space, respectively. If such a mappingexists, the differential equation can be rewritten through a coordinate trans-formation. The first derivative can be written as ∂x = ξx∂ξ by the chain rule.443.3. Smooth MeshesFigure 3.3: Schematic of multigrid scheme for the ETE.Figure 3.4: Fine and coarse grid numbering.453.3. Smooth MeshesCorrecteddiscretizationerror10−3 10−2 10−110−910−810−710−610−510−410−310−210−1h (2,0,0)(2,2,4)(2,2,4) multigrid(2,2,4) multigrid two levelsSecond orderFourth orderFigure 3.5: Discretization error and error in error estimate for (2, 2, 4) multi-grid scheme using uniform meshes.463.3. Smooth MeshesFigure 3.6: Mesh that is smoothly expanded by a factor 1 + s.The second derivative can then be written as∂2xu = ∂x(∂xu)= ∂x(ξxuξ)= ξx∂xuξ + ξxxuξ= (ξx)2uξξ + ξxxuξ. (3.15)This means that the equation can be discretized in either physical or com-putational space, and if the latter is chosen then this is equivalent to solvinga different linear PDE on a uniform mesh.One simple example in 1D is a sequence of meshes where each cell isexpanded by a factor (1 + s) compared to the preceding cell, with the firstcell having size h0, and s = O(1/N). This is shown in Figure 3.6. The gridmapping is found as follows. First, the domain constraint isN−1∑j=0h0(1 + s)j = 1 (3.16)which can be expanded ash01− (1 + s)N1− (1 + s) = 1, (3.17)from which the conditionh0 = − s1− (1 + s)N (3.18)473.3. Smooth Meshes0 0.2 0.4 0.6 0.8 100.20.40.60.81xξ(x) UniformSmoothUnderlying grid functionFigure 3.7: Mesh points for a smoothly mapped mesh with N = 10, and itsunderlying grid mapping function.483.4. Perturbed Meshesis required. For a given grid, we haveξ(xi) =i−1∑j=0h0(1 + s)j= h0(1− (1 + s)i1− (1 + s))= −h0s(1− (1 + s)i)=1− (1 + s)i1− (1 + s)N . (3.19)Using xi = ih, s = h = 1/N , this becomesξ(x) =1− (1 + s)x/h1− (1 + s)1/h=1− (1 + h)x/h1− (1 + h)1/h . (3.20)In the continuous limit, h→ 0, this isξ(x) = limh→01− (1 + h)x/h1− (1 + h)1/h= limh→011− (1 + h)1/h − limh→0(1 + h)x/h1− (1 + h)1/h=11− e −ex1− e=1− ex1− e (3.21)An illustration of the N = 10 mesh and the underlying mapping is shown inFigure 3.7. For such smooth meshes, the asymptotic accuracy of the errorestimate is similar to the uniform case. This can be seen in Figure 3.8, whereagain a (2, 2, 4) scheme produces a fourth order accurate error estimate.3.4 Perturbed MeshesTo simulate the nonsmooth characteristics of an unstructured mesh geome-try, we can randomly perturb the nodes of a uniform grid. A uniform random493.4. Perturbed Meshes10−3 10−2 10−110−910−810−710−610−510−410−310−210−1hCorrected discretization error (2,0,0)(2,2,4)Second orderFourth orderFigure 3.8: Discretization error and error in error for (2, 2, 4) using smoothlystretched meshes.503.4. Perturbed Meshes0 0.2 0.4 0.6 0.8 1−1012345 x 10−3xDiscretization Error Exact(2,2,4)(2,4,4)Figure 3.9: Exact discretization error and estimate for (2, 2, 4) and (2, 4, 4)schemes using a randomly perturbed N = 40 mesh.distribution is chosen between −h/3 and h/3, so that the possibility of meshentanglement is eliminated. The same procedure as before is applied to amesh with N = 40, the exact error and estimates of which are shown inFigure 3.9.It can be seen that even though (2, 2, 4) captures much of the trends ofthe exact error, it is not quite accurate, as quantified by asymptotic behaviorof the norm of the error in error as N is increased, shown in Figure 3.10.Also shown in these figures is an alternative scheme of (2, 4, 4), which turnsout to give accurate error estimates.It seems necessary to have p < q = r for an accurate error estimate,which is our conjecture for the time being; this is discussed in more detail inSection 3.7. This requirement does not seem more economical than usinga higher order solution, a priori, since a second order problem and a fourthorder problem need to be solved in either case. Indeed this is true if theprimal problem is linear, but the purpose in this chapter is to first discernviable combinations of schemes, and thereafter improvements in efficiencyand robustness are explored for nonlinear and unsteady problems.513.4. Perturbed Meshes10−3 10−2 10−110−910−810−710−610−510−410−310−210−1hCorrected discretization error (2,0,0)(2,2,4)(2,4,4)Second orderFourth orderFigure 3.10: Discretization error and error in error for (2, 2, 4) and (2, 4, 4)schemes using randomly perturbed meshes.523.5. Nonlinear Model Problems3.4.1 Nested GridsA subclass of nonsmooth meshes are ones that result from nested refine-ments of a coarse nonsmooth mesh. These meshes have only C0 continu-ity in the mapping, which might cause problems if the existence of secondderivatives in the mapping is required. For simplicial meshes, this involvessplitting each cell in half in 1D, and each triangle into four by adding ver-tices at the midpoint of each edge in 2D. This nested splitting operation, ineffect, preserves the initial nonsmooth characteristics. In 1D this convergesto a piecewise linear mapping ξ(x), which is insufficient for the Poissonequation, since two derivatives of the grid function would be required.Some mesh operations such as adaptation will generate cells in this man-ner. However, in general, global nested refinements would be restrictive inmultiple dimensions; in some instances there is no underlying coarse meshfrom which to start. In others, the geometry imposes constraints whichwould prohibit these operations. The underlying grid function for a N = 40mesh obtained from nested refinement of aN = 10 perturbed mesh is shownin Figure 3.11. The discretization error, error estimate, and error in error forthe (2, 2, 4) and (2, 4, 4) schemes are shown in Figure 3.12 and Figure 3.13,respectively. The accuracy of the solution and error estimates are summa-rized in Figure 3.14. Although the (2, 2, 4) error estimate seems qualitativelygood, a grid study reveals that it is not higher order accurate, unlike the(2, 4, 4) scheme.3.5 Nonlinear Model ProblemsThis example illustrates that the results of this chapter also extend to non-linear problems. Consider the steady viscous Burgers’ equation,ududx= νd2udx2+ s, Ω = {x ∈ [0, 1]}. (3.22)A manufactured solution ofu(x) = sin(2pix) (3.23)is used, which satisfies the boundary conditionsu(0) = u(1) = 0. (3.24)Figure 3.15 shows the accuracy of the error estimates for uniform and per-turbed meshes using the (2, 2, 4) and (2, 4, 4) schemes. The same conclusion533.5. Nonlinear Model Problems0 0.2 0.4 0.6 0.8 100.20.40.60.81xξ(x) Initial perturbed gridNested refinementFigure 3.11: Initial coarse perturbed mesh (N = 10) and its nested refine-ment until N = 40.543.5. Nonlinear Model Problems0 0.2 0.4 0.6 0.8 1−3−2−10123 x 10−3xDiscretization error ExactEstimate(a) Discretization error and error estimate0 0.2 0.4 0.6 0.8 1−8−6−4−20246 x 10−5xError in error(b) Error in error estimateFigure 3.12: Error accuracy for the (2, 2, 4) scheme using a nested perturbedgrid.553.5. Nonlinear Model Problems0 0.2 0.4 0.6 0.8 1−3−2−10123 x 10−3xDiscretization error ExactEstimate(a) Discretization error and error estimate0 0.2 0.4 0.6 0.8 1−2−1.5−1−0.500.511.52 x 10−5xError in error(b) Error in error estimateFigure 3.13: Error accuracy for the (2, 4, 4) scheme using a nested perturbedgrid.563.6. Other Attempts at Truncation Error Estimation10−3 10−2 10−110−910−810−710−610−510−410−310−210−1hCorrected discretization error (2,0,0)(2,2,4)(2,4,4)Second orderFourth orderFigure 3.14: Discretization error and error in error using nested perturbedgrids.is reached as with the linear case: (2, 2, 4) is not sufficient to give accu-rate error estimates for nonsmooth meshes, and (2, 4, 4) needs to be usedinstead.3.6 Other Attempts at Truncation Error EstimationIn this section several other plausible methods of computing the residualand solving the ETE are explored for nonsmooth meshes.3.6.1 Residual FilteringFrom the observation that the residual on nonsmooth meshes has large highfrequency content, a natural resolution would be to apply a filter as a wayto smooth the p-truncation error estimate. Many approaches exist in otherfields of study, such as low pass filters in signal processing, as well Laplaciansmoothing [51], Gaussian filter [20], and many other variants, which have573.6. Other Attempts at Truncation Error Estimation10−3 10−2 10−110−910−810−710−610−510−410−310−210−1hCorrected discretization error (2,0,0) uniform(2,0,0) perturbed(2,2,4) uniform(2,2,4) perturbed(2,4,4) perturbedSecond orderFourth orderFigure 3.15: Discretization error and error in error using uniform and ran-dom meshes for Burgers’ equation.corresponding operations in the frequency domain. An attempt at smooth-ing using a Gaussian filter is examined, with the error estimate resultingfrom filtering the residual source term shown in Figure 3.16. It is evidentthat this is not an effective approach to obtain accurate error estimates; itseems that the high frequency components are important in the sense of aGreen’s function for the error estimate. Other more sophisticated filteringtechniques give qualitatively similar results.3.6.2 h-Truncation Error EstimateAnother plausible method of evaluating the residual source term is throughwhat is known as the h-truncation error estimate. This draws from theinsight that the continuous ETE has first been derived, and while it is conve-nient to use the same mesh (which is in fact the focus of the current work),it is also possible to use a different mesh to evaluate residuals. This is infact quite common in finite-element methods, where the residual is orthog-onal to the space of test functions, and an enriched space must be used to583.6. Other Attempts at Truncation Error Estimation0 0.2 0.4 0.6 0.8 1−0.025−0.02−0.015−0.01−0.00500.005xDiscretization error ExactEstimateFigure 3.16: Discretization error and error estimate resulting from the ap-plication of a Gaussian filter to the residual source term.evaluate the residual. As in the uniform case, it is possible to derive theinterpolation operator for a nonuniform mesh, by finding coefficients a andb for the fine and coarse meshes related byUf,2i−1 = aUc,i + bUc,i−1 (3.25)Uf,2i = aUc,i + bUc,i+1. (3.26)By imposing the control volume average conditions which should be pre-served after applying the transfer operators, the respective coefficients canbe found to bea =hi + 2hi−12(hi + hi−1), b =hi2(hi + hi−1), (3.27)a =hi + 2hi+12(hi + hi+1), b =hi2(hi + hi+1), (3.28)which of course reduces to a = 34 , b =14 for a uniform mesh. This schemeis denoted by (2, 2, 2∗), since the residual is computed using r = 2, buton a finer mesh. The resulting discretization error and error estimate areshown in Figure 3.17. These results are analogous to the p-truncation error593.6. Other Attempts at Truncation Error Estimation0 0.2 0.4 0.6 0.8 1−20246810 x 10−3xDiscretization error ExactEstimateFigure 3.17: Discretization error and error estimate for a (2, 2, 2∗) schemeusing h-truncation error estimate.estimates in the same way that (2, 2, 4) schemes do not produce higher orderaccurate discretization error estimates for nonsmooth meshes.3.6.3 Projection onto Uniform MeshDrawing once again from the insight that the continuous ETE has first beenderived, another approach to is to evaluate the residual on a uniform mesh,motivated by the desirable properties that were observed for such meshes.Thereafter, the ETE can be solved on either the original perturbed mesh orthe uniform mesh and the error estimate projected back to the perturbedmesh.The projection of the primal solution from the perturbed mesh to a uni-form mesh is shown in Figure 3.18, and the projection of the residual com-puted using this solution is projected back to the perturbed mesh is shownin Figure 3.19. The ETE can then be solved on either the original perturbedmesh, shown in Figure 3.20, or solved on the uniform mesh with the errorestimate projected back to the original perturbed mesh, shown in Figure3.21. It is evident that the error estimates are unacceptable, suggesting thatthis mesh to mesh transfer fails to capture the important nonsmooth features603.6. Other Attempts at Truncation Error Estimation0 0.2 0.4 0.6 0.8 1−1−0.8−0.6−0.4−0.200.20.40.60.81xU2 Perturbed meshProjection on uniform meshFigure 3.18: Projection of primal solution from a perturbed mesh to a uni-form mesh.0 0.2 0.4 0.6 0.8 1−100−50050100150xR4(U2) Uniform meshProjection on perturbed meshFigure 3.19: Projection of residual from uniform mesh to original perturbedmesh.613.7. Error Estimation for Other Discretization Schemes0 0.2 0.4 0.6 0.8 1−0.04−0.03−0.02−0.0100.010.020.030.04xDiscretization error ExactEstimateFigure 3.20: Discretization error and (2, 2, 4) error estimate using a uniformmesh to compute residual, then solving ETE on perturbed mesh.in the truncation error, and leads to a discretization error estimate that isnot accurate. In these results the order of accuracy of the error estimate isno higher than second order, the accuracy of the primal solution.3.7 Error Estimation for Other DiscretizationSchemes3.7.1 Other Discretization OrdersHere, we explore the other choices of discretization orders p, q, r and showempirically that p < q = r is the only combination that gives an asymptoti-cally accurate error estimate. For the Poisson problem with p = 2, the choiceof q is varied in {2, 4, 6} and r varied in {4, 6}, keeping r > p, as we showto be a requirement in Chapter 4. The resulting error estimate and error inerror plots are shown in Figure 3.22. It can be seen that Equation 3.1 holdsfor these combinations. The corresponding plots for perturbed meshes areshown in Figure 3.23. It can be seen that the error accuracy behaves verydifferently here, and it is only those schemes where p < q = r that give error623.7. Error Estimation for Other Discretization Schemes0 0.2 0.4 0.6 0.8 1−0.04−0.03−0.02−0.0100.010.020.030.04xDiscretization error ExactEstimateFigure 3.21: Discretization error and (2, 2, 4) error estimate using a uniformmesh to compute residual, then solving ETE on the same uniform mesh, andprojecting back to the perturbed mesh.633.7. Error Estimation for Other Discretization Schemesestimates that are asymptotically higher order, suggesting the form||p − ˜p|| = O(hp+(q−p)δqr). (3.29)It can be also observed when comparing the (2, 2, 4) and (2, 2, 6) schemes,that for the nonsmooth meshes, using p-truncation error with p = q, theerror estimate is only accurate to order p and not higher order, regardless ofthe value of r.3.7.2 Finite-Difference MethodThe ETE methodology can be applied to other discretization schemes aswell, and similar conclusions result. This suggests that difficulties are dueto mesh properties rather than choice of discretization. The Poisson testcase is used with a finite-difference discretization. The perturbed grids aregenerated in the same way, with the only difference being the way in whichEquation 3.2 is discretized. The approach is to approximate derivatives atgrid points to the desired accuracy using Taylor series expansions. For uni-form meshes,d2Udx2∣∣∣∣i=Ui+1 − 2Ui + Ui−1h2+O(h2) (3.30)is used for a second order accurate approximate to the second derivative.However, for perturbed meshes, Equation 3.30 is not a second order accu-rate approximation. Instead, in seeking the coefficients ford2Udx2∣∣∣∣i= c−2Ui−2 + c−1Ui−1 + c0Ui + c1Ui+1 + c2Ui+2 +O(h2), (3.31)Taylor series expansions give the coefficients as follows:c−2 =−2(2hi−1hi − hihi+1 + hi−1hi+1 − h2i )hi−2(hi−2 + hi−1)(hi−2 + hi−1 + hi)(hi−2 + hi−1 + hi + hi+1)c−1 =2(2hi−1hi + 2hi−2hi − hihi+1 + hi−1hi+1 + hi−2hi+1 − h2i )hi−2hi−1(hi−1 + hi)(hi−1 + hi + hi+1)c0 =−2hi−1hi(hi + hi+1)(hi−2 + hi−1)(4hi−1hi + 2hi−2hi − hihi+1 − ...hi−2hi−1 + 2hi−1hi+1 + hi−2hi+1 − h2i−1 − h2i )c1 =2(2hi−1hi + hi−2hi − hi−2hi−1 + 2hi−1hi+1hi−2hi+1 − h2i−1)hihi+1(hi−1 + hi)(hi−2 + hi−1 + hi)643.7. Error Estimation for Other Discretization Schemes0 1−202x 10−3xExact error Computed error(a) Exact and computed discretization error for (2, q, r) schemes0 1−202x 10−3xExact error Computed error0 1−33−5x−33x 10−5−33x 10−5−55x 10−7−33x 10−5−55x 10−7O(h4) O(h6)(b) Error in error estimate for (2, q, r) schemes colored by asymptotic behaviorsFigure 3.22: Discretization error and estimate for uniform mesh for the 1DPoisson equation. 653.7. Error Estimation for Other Discretization Schemes0 1024x 10−3xExact error Computed error t rr r t rr r(a) Exact and computed discretization error for (2, q, r) schemes0 1−202x 10−3xExact error Computed error0 1−33−5x−33x 10−5−33x 10−5−55x 10−7−33x 10−5−55x 10−7O(h4) O(h6)55 0355 3 355 3 (h2) (h4) O(h6)(b) Error in error estimate for (2, q, r) schemes colored by asymptotic behaviorsFigure 3.23: Discretization error and estimate for randomly perturbed meshfor the 1D Poisson equation. 663.7. Error Estimation for Other Discretization Schemes10−3 10−2 10−110−1010−910−810−710−610−510−410−310−2hCorrected discretization error (2,0,0) uniform(2,0,0) perturbed(2,2,4) uniform(2,2,4) perturbed(2,4,4) perturbedSecond orderFourth orderFigure 3.24: Discretization error and error in error for finite-differenceschemes.c2 =−2(2hi−1hi + hi−2hi − hi−2hi−1 − h2i−1).hi+1(hi + hi+1)(hi−1 + hi + hi+1)(hi−2 + hi−1 + hi + hi+1)(3.32)It can be seen that very complicated expressions result for nonuniform meshes.If the perturbation goes to zero, then the expression becomes the stan-dard fourth order accurate approximation of the second derivative using fivepoints. Using analogous expressions for the other derivatives, the same ETEapproach is applied as before for the finite-volume method. The resultingerror estimates have accuracy behavior which is shown in Figure 3.24. Asbefore, for uniform meshes (2, 2, 4) is sufficient for a fourth order accurateerror estimate, but for perturbed meshes (2, 4, 4) is needed instead.3.7.3 Finite-Element MethodsFinite-element methods are somewhat different in formulation from finite-difference and finite-volume, in the sense that the solution is optimal insome norm such as the energy norm or L2 norm in the chosen functionspaces. For instance, using continuous linear elements, no perturbation673.7. Error Estimation for Other Discretization Schemeschosen from that space can be added to the primal solution to make animprovement. The notion of truncation error is therefore not defined in astraightforward way, possibly requiring enriching the primal solution spaceby means of a higher order operator, or using a different mesh. However, er-ror estimation techniques are abundant using some theoretical techniques,and such estimates can give robust error bounds, as found in the literature[3, 8], and more recently [1].68Chapter 4ETE Methodology Applied toSteady ProblemsThe difficulties with accurate error estimation encountered in one dimen-sion exist equally in multiple dimensions. For a second order PDE, thediscretization error estimate computed from the ETE using a p-truncationerror source term from the previous chapter is not asymptotically higher or-der accurate when (2, 2, 4) is used with unstructured meshes. This is firstconfirmed using the 2D Poisson equation, posed as−∆u = s, on Ω = {(x, y) ∈ [0, 1]2}, (4.1)with Dirichlet boundary conditions u|∂Ω = 0. The exact solution isu(x, y) = sin(pix) sin(piy), (4.2)and the source term is s(x, y) = 2pi2 sin(pix) sin(piy). Results for the (2, 2, 4)scheme are shown in Figure 4.1. Even though the error estimate and theexact discretization error are qualitatively similar, their peak difference iscomparable to the error itself at this resolution, and is still second orderasymptotically, which only provides, at best, an improvement by a constantfactor. In Figure 4.2, it can be seen that using a (2, 4, 4) scheme instead givesan error in error that is approximately two orders of magnitude smaller atthis resolution, and is shown to be fourth order asymptotically in Figure4.3. This difference may be hard to discern qualitatively by the plots of theerror estimate, hence it is important to compute the error in error. In thischapter, the generalization of this requirement will be explored for schemeswith other combinations of (p, q, r).4.1 LinearizationIn the previous section, there was evidence that when using the p-truncationerror as a source term in the ETE for nonsmooth meshes, a fourth order694.1. Linearization0 101 0.005 0.105 0.205 0.305(a) solution0 101 -0.00018 -7E-05 4E-05(b) exact error0 101 -0.00018 -7E-05 4E-05(c) (2, 2, 4) error estimate0 101 -1.9E-05 -4.5E-06 1E-05(d) error in errorFigure 4.1: Results for the Poisson problem using the ETE with (2, 2, 4).704.1. Linearization0 101 -0.00018 -7E-05 4E-05(a) (2, 4, 4) error estimate0 101 -1E-07 1E-07 3E-07(b) error in errorFigure 4.2: Error estimate and error in error using the ETE with (2, 4, 4).10−3 10−2 10−110−910−810−710−610−510−410−310−2hCorrected discretization error (2,0,0)(2,2,4)(2,4,4)Second orderFourth orderFigure 4.3: Accuracy of the second order solution, and error estimates usingthe (2, 2, 4) and (2, 4, 4) schemes for the Poisson problem.714.1. Linearizationaccurate error estimate for the second order primal solution can only beobtained with (2, 4, 4). This suggests that accurate error estimates will beasymptotically expensive, requiring the solution of an auxiliary fourth or-der problem. For a nonlinear primal problem, the asymptotic ETE cost canbe reduced by means of solving a linearized ETE problem. Starting fromEquation 2.6, which is the proper, non-approximated form of the ETE, weconsider the most common linearization approach, as in [25, 50, 84, 122],sometimes referred to as Newton linearization. This method works well forsmooth solutions without discontinuities.Assuming the discretization error is small in comparison to the solution,||e|| ||u˜||, the series expansion of the flux and source term in the contin-uous ETE of Equation 2.6, taken about u˜, can be written asF(u˜+ e) = F(u˜) + ∂uF(u˜)e+O(||e||2) (4.3)s(u˜+ e) = s(u˜) + ∂us(u˜)e+O(||e||2). (4.4)The O(||e||2) terms are neglected, and the result is substituted back intoEquation 2.6, which becomes∂te+∇ · (∂uF(u˜)e)− ∂us(u˜)e ≈ −(∂tu˜+R(u˜)). (4.5)Since we neglected terms that are on the order of the square of the dis-cretization error, we expect this linearization to be sufficiently accurate ifthe ETE is discretized to an order that is not more than twice the discretiza-tion order of the primal problem, which we confirm in the results. Severalfurther simplifications of Equation 4.5, and the failures of some of theseapproaches can be found in the literature, as seen in [10] and [50]. Thediscretized equations of Equation 2.16 for the primal equation, arising fromthe pseudo-transient continuation formulation, can be written for the ETE,as (1k+∂Rˆ∂)∣∣∣∣∣n(n+1 − n) = −Rˆ(n), (4.6)whereRˆ() = R(Up + )−R(Up)− T (Up) (4.7)is the spatial residual for the ETE, defined in terms of T , which is the ETEsource term, and Up which is the converged primal solution to order p. Thefirst two terms on the right side of Equation 4.7 are the primal residuals tothe discretization order of the ETE, which individually include the physicalsource terms, as in Equation 2.12. The analytic Jacobian is724.1. Linearization∂Rˆ∂=∂R(Up + )∂U, (4.8)where the derivatives on R and Rˆ are taken simply with respect to their ownarguments. For viscous flows, the term involving the gradient of the errorshould also be included, as∂Rˆ∂∇ =∂R(Up + ,∇(Up + ))∂∇U . (4.9)When linearization is performed, we haveRˆ() = R(Up) +∂R(Up)∂U+O(||||2)−R(Up)− T (Up)=∂R(Up)∂U− T (Up) +O(||||2). (4.10)This ETE residual involves the primal Jacobian evaluated at the convergedprimal solution, which is a constant matrix. The Jacobian for the linearizedETE can be directly computed as∂Rˆ∂=∂R(Up)∂U, (4.11)and also if there is a gradient contribution∂Rˆ∂∇ =∂R(Up,∇Up)∂∇U . (4.12)If the source term also depends on solution or its gradient, as we will see inthe models for turbulent flows, there will also be a source Jacobian contri-bution. For the nonlinear ETE, we perform pseudo time-stepping to reachsteady-state as in the primal problem. For the linearized ETE, the linearityof the problem can be exploited, requiring only the solution of one linear al-gebraic system. This can be seen from Equation 4.10, which at steady-stategives∂R(Up)∂U = T (Up). (4.13)The matrix on the left is simply the primal Jacobian (to the order of theETE discretization), and the ETE residual source term is on the right. Theletter "L" will be appended to the (p, q, r) notation to indicate that the lin-earized ETE is solved. Reducing the ETE problem to a single linear algebraic734.1. Linearization10−2 10−110−1210−1010−810−610−410−2100length scale||ǫ2−ǫ244||4.114.012.21 Primal FluxNewton LinearizationFully Nonlinear(a) The (2, 4, 4) scheme10−2 10−110−1210−1010−810−610−410−2100length scale||ǫ2−ǫ266||5.814.512.21 Primal FluxNewton LinearizationFully Nonlinear(b) The (2, 6, 6) schemeFigure 4.4: Convergence in error estimate with mesh length scale using thefully nonlinear error flux, the linearized flux, and the primal flux for the 1Dviscous Burgers’ equation.system has two important implications. First, only a single linear algebrasystem needs to be solved, resulting in efficiency advantages over the fullynonlinear problem. This is done by the same GMRES iterative method. Fur-thermore, linearization removes any stability requirements in solving theETE by pseudo time-stepping a higher order problem. Comparisons of bothefficiency and robustness with the higher order primal benchmark will beexamined in the results.To observe the effect of linearization error for a 1D viscous Burgers’equation of Equation 3.22, Figure 4.4 shows that (2, 4, 4)L gives an errorestimate whose accuracy is fourth order, whereas (2, 6, 6)L gives an errorestimate that is also only fourth order, since it is limited by the order 2plinearization error. It can also be seen that using the primal flux for the ETEgives poor results.4.1.1 Relinearization of the ETEAs shown later, linearization of the ETE shows promise in dramatically re-ducing the computational cost without sacrificing asymptotically higher or-der accuracy in the error estimate. However, for a mesh with a particularlength scale, it is possible that the linearization error is large in comparisonto the other terms dictating the accuracy of the ETE. Here, the plausibility ofrelinearizing the ETE is explored. This can be used whenever the lineariza-744.2. Properties of ETE Discretizationtion error is large in the error estimate, and the ETE can be solved againusing the linearization at the new corrected solution as needed. Anotherway to identify this is when the second order solution is far away from theexact or higher order solution. Since the linearization error behaves likethe square of the second order discretization error, this quantity may not besmall in such cases. This is typical when using turbulence model which maynot have all features resolved. In our results, large variations in the turbu-lence variable can be observed between second and fourth order solutionson typical meshes for turbulent airfoil cases. Relinearization of the ETE isperformed once using the corrected solution from the original linearizedETE error estimate.Another possible application of ETE relinearization may be an instancewhere a sequence of progressively higher order accurate error estimates isdesired. If the meshes are in the asymptotic region, it is possible to con-struct a sequence of corrected solutions accurate to order p → 2p → 22p...using only an order p primal solution and solutions of linear systems usingthe higher order discretization matrices and forcing term on the right. Inother words, an order p corrected solution can in principle be obtained bya low order nonlinear solution and dlog2 pe − 1 solutions of linear algebraicsystems, which can be less computationally expensive than a full nonlin-ear solve to order p. For example, an 8th order accurate solution could beobtained by solving the 2nd order primal problem, followed by a linear al-gebraic system for the 4th order discretization, and then a linear algebraicsystem for the 8th order discretization. The 2-4-6 sequence is shown inFigure 4.5 for a 1D viscous Burgers’ equation example. This sequence waschosen instead of 2-4-8 due to solver restrictions, but it can be seen thatthe iteratively corrected solution is 6th order accurate instead of 4th orderaccurate if (2, 6, 6)L were used.4.2 Properties of ETE DiscretizationIn this section, we wish to generalize to other orders (p, q, r) the previousobservation that the (2, 2, 4) scheme gives only second order accurate errorestimates while the (2, 4, 4) scheme gives a fourth order error estimate.4.2.1 Exact Discrete Residual Source TermLet Lp and Lˆq be the discrete derivative operators for the primal equation(to order p) and the ETE (to order q) respectively. For the discussion, we754.2. Properties of ETE DiscretizationCorrecteddiscretizationerror10−2 10−110−1410−1210−1010−810−610−410−2100h (2,0,0)(2,4,4) linearized(2,6,6) linearized(2,4,4),(2,6,6) relinearizedsecond orderfourth ordersixth orderFigure 4.5: Accuracy of relinearization for 1D Burgers’ equation.764.2. Properties of ETE Discretizationconsider the case when the operators are linear or linearized, where L = Lˆ,although numerical results show that these properties hold for nonlinearproblems also. Here, unlike the discrete residual, the operators do not in-clude the source term. We have thenLpUp = S. (4.14)The ETE operator is defined asLˆqpqr := Lq(Up + pqr)− LqUp, (4.15)and the estimate is obtained by solvingLˆqpqr = −(LrUp − S). (4.16)From Equation 4.15, the exact error satisfiesLˆqp := Lq(Up + p)− LqUp= Lq(Up + p)− S − (LqUp − S)= LqU − S − (LqUp − S)= τq − (LqUp − S) (4.17)which states that estimating the ETE source term using a reconstructionto order q of the control volume averages of the exact solution minus thediscrete residual using reconstruction to order q of the converged primal so-lution accurate to order p will yield the exact error estimate. This allowsfor some insight on the accuracy of the error estimate. We will first demon-strate that r ≥ p is necessary in obtaining an error estimate that is evenasymptotically correct, and that r = p does not give anything useful.Case 1 : r < pEquation 4.16 can be written asLˆqpqr = −(LrUp − S)= −Lr(Up − Ur)− (LrUr − S)= −Lr(Up − Ur) (4.18)where Ur satisfies the discrete equations to order r. Then,774.2. Properties of ETE Discretizationpqr = − (Lq)−1 Lr(Up − Ur)||pqr|| = O(1)O(hmin(p,r))= O(hmin(p,r)). (4.19)For r < p this is O(hr), which is asymptotically larger than ||p|| = O(hp),for Lq nonsingular. This says that the error estimate is asymptotically largerthan the exact error estimate, which is undesirable.Case 2 : r = pIn this case Equation 4.16 becomesLˆqpqr = −(LrUp − S) = 0 (4.20)which means the error estimate is identically zero. As before the exact erroris ||p|| = O(hp). Therefore, the only schemes that are candidates in obtain-ing accurate error estimates are those with r > p, for which ||p|| = O(hp),and ||pqr|| = O(hp). However, this does not guarantee that the error es-timate will be asymptotically more accurate or even more accurate by aconstant factor. We seek, therefore, schemes for which p − pqr is a higherorder quantity, which is possible with the remaining possibility of r > p,with additional conditions.4.2.2 Expected Order of Accuracy for the Error EstimateNote that when q = r, Equation 4.16 becomesLˆqpqq = −(LqUp − S). (4.21)Using Equation 4.15 this can be written asLq(Up + pqq)− LqUp + LqUp − S = 0 (4.22)whenceLq(Up + pqq)− S = 0. (4.23)Comparing this to Equation 4.14 at order q, this meansUp + pqq = Uq. (4.24)784.2. Properties of ETE DiscretizationHence we havep − pqq = q||p − pqq|| = O(hq) (4.25)Therefore our claim is that in this context, choosing p < q = r guaranteesa higher order accurate error estimate; we will verify this experimentally aswell.The accuracy of the error estimate for unstructured nonsmooth meshescan be summarized from these results as||p − pqr|| = O(hp+(q−p)δqr) (4.26)for the nonlinear ETE, and||p − pqr|| = O(hmin(2p,p+(q−p)δqr)) (4.27)for the linearized ETE. The extra 2p term is needed for the linearized version,since terms on the order of the square of discretization error, were neglected.The expressions for uniform meshes, as discussed in [10], are O(hmin(p+q,r))and O(hmin(2p,p+q,r)), respectively. A major difference is that for uniformmeshes, p, q, and r can be arbitrarily increased, keeping r > p, withoutdetrimental effects on the accuracy of the error estimate, whereas this isnot true for the nonsmooth case, due to the fact that the accuracy of theresidual is not increased as r increases. For example, by Equation 4.26, a(2, 4, 4) scheme produces an error estimate that is fourth order accurate,but increasing the discretization order of the ETE only to (2, 6, 4) leads toan error estimate that is only second order accurate. Recall that this wasconfirmed in 1D in Figure 3.23.4.2.3 Equivalence of Defect Correction and Higher OrderScheme when q = rEquation 4.24 shows that solving the nonlinear ETE with q = r and applyingthe error estimate as a defect correction yields the same result as solving thehigher order primal problem. This result says that the difference in errorestimate, p − pqq, is independent of p, although the individual terms dodepend on p. This independence is observed numerically using the fullynonlinear error flux, and within linearization error for the linearized ETE.The terms of the left side of Equation 4.25 are individually of size O(hp),and the term of the right is O(hq), a higher order term if p < q.794.3. Compressible Inviscid Flow4.2.4 Equivalence of Nonlinear ETE and Order Ramping whenq = rA common approach, as in [23], for obtaining an initial guess for a higherorder scheme is to use the converged solution from a lower order problem.The approach taken in this work in solving the fully nonlinear ETE to higherorder is exactly equivalent to this. Consider solving the primal problemusing a second order scheme and using the converged solution to initializea fourth order scheme. The implicit pseudo time-stepping for the primalproblem is given by Equation 2.16, with initial conditions U0 = U2, whilefor the ETE with T (Up) = R(Up), Equation 4.6 becomes(1k+∂∂(R(U2 + n)−R(U2)))∆n = −(R(U2 + n)−R(U2) +R(U2)),(4.28)where ∆n = n+1 − n. Hence,(1k+∂R(U2 + n)∂)∆n = −R(U2 + n). (4.29)Here the discrete residuals are evaluated to order q, the ETE discretizationaccuracy. The initial condition for the ETE is 0 = 0. If there is only onesolution that maps R to 0, then it can be shown by induction on n for Equa-tion 4.29 and Equation 4.6 that Un = U2 + n, meaning that at every stepthe two methods are equivalent, up to round-off error.4.3 Compressible Inviscid FlowWe proceed with showing the results for more complex test cases. First,the compressible Euler equations as a model of inviscid flow are consideredbefore proceeding with more sophisticated models. In two dimensions thevector of conserved variables consists of density ρ, velocity v =(u v)T ,and total energy E, which is related to pressure P by the equation of statefor a calorically perfect gasP = (γ − 1)(E − 12ρ|v|2). (4.30)804.3. Compressible Inviscid FlowName Value Name ValuerI 2 γ 1.4rO 3 UrI rIMIργ−12IMI 2 r√x2 + y2ρI 1Table 4.1: Summary of values used for the supersonic vortex test case.In these equations, γ is the ratio of specific heats cpcv , a constant. The Eulerequations can then be written as∂t ρρvE+∇ · ρvρvv + P1v(E + P ) = 0. (4.31)4.3.1 Supersonic VortexA purely supersonic flow around a quarter annulus is used as a test case, asstudied by [2]. The exact solution isρuvP =ρI[1 + γ−12 M2I(1− r2Ir2)] 1γ−1yUrI/r2−xUrI/r2ργ/γ (4.32)defined on Ω ={(r, θ) ∈ [rI , rO]×[0, pi2]}. The test case parameters aresummarized in Table 4.1. The ETE will not be explicitly written out buttakes the form of Equation 2.6. The solution and error estimates are shownin Figure 4.6, and the asymptotic accuracy is shown in Figure 4.7. It canbe seen that as before, the (2, 2, 4) scheme is only second order accurate forthe error estimate, while (2, 4, 4) and (2, 4, 4)L gives fourth order accuracy.The meshes used in this section are quite isotropic, but not unreasonableconsidering the inviscid flow physics. The test cases in the following sectionswill investigate the ETE performance on anisotropic families of meshes.814.3. Compressible Inviscid Flow4.05 6.30 8.55(a) Energy-0.0500 0.0000 0.0500(b) (2, 0, 0) error-0.0500 0.0000 0.0500(c) (2, 4, 4) error estimate-0.0500 0.0000 0.0500(d) (2, 4, 4)L error estimate824.3. Compressible Inviscid Flow-0.0073 -0.0021 0.0032(e) (2, 4, 4) error in error-0.0073 -0.0021 0.0032(f) (2, 4, 4)L error in errorFigure 4.6: Exact error, error estimate, and difference for the energy variableusing the (2, 4, 4) discretization for the supersonic vortex test case.10−2 10−110−610−410−21004.21412121length scale||ǫp−ǫpqr|| (4,0,0)(2,4,4)(2,4,4), Linear(2,0,0)(2,2,4)Figure 4.7: Convergence of the error estimate for the supersonic vortexproblem.834.4. Compressible Viscous Flow4.3.2 Gaussian BumpAnother classical test case is the Gaussian bump. The computational domainis defined asΩ ={(x, y) ∈ [−1.5, 1.5]× [0, 0.8] : y > 116e−25x2}(4.33)The freestream Mach number is M = 0.5. This is a good test case for theEuler equations, and a measurement of the entropy s is used as a proxy foraccuracy. This is convenient since the entropy error can be computed de-spite the nontrivial geometry of the problem, without necessarily requiringthe knowledge of the exact solution. The entropy of the exact solution iszero since the flow is inviscid with no external thermodynamic interactions.Hence, the entropy that is computed can be used as a measure of error,computed ass =PP∞(ρ∞ρ)γ− 1 (4.34)The nondimensionalization constants are ρ∞ = 1 and P∞ = 1/γ. The plotsof the entropy using only primal solutions, and the solution corrected by alinearized ETE error estimate are shown in Figure 4.8, and the asymptoticbehaviors are shown in Figure 4.9. The Richardson extrapolated orders ofaccuracy are shown in Table 4.2. For the finest meshes, the (2, 0, 0) entropyis second order, as expected, and the (4, 0, 0) and (2, 4, 4)L schemes giveorders of accuracy that are quite close to each other and quite close to fourthorder. As expected from the analysis of Section 4.2, the (2, 4, 4) values foreach mesh are identical to those obtained from the (4, 0, 0) scheme, andthe results match the expectation that (2, 4, 4)L would give results close to(2, 4, 4).4.4 Compressible Viscous FlowWe consider next the compressible laminar Navier-Stokes equations as amodel of viscous fluid flow. Detailed discussions can be found in the litera-ture [13]. In two dimensions the vector of conserved variables is the sameas the Euler equations. The Navier-Stokes equations contains the viscousand heat transfer terms, which can in all be written as∂t ρρvE+∇ · ρvρvv + P1− τv(E + P )− τ · v + q = 0 (4.35)844.4. Compressible Viscous Flow-0.0001 -0.0000 0.0001(a) (2, 0, 0) entropy-0.0001 -0.0000 0.0001(b) (4, 0, 0) entropy-0.0001 -0.0000 0.0001(c) (2, 4, 4)L entropyFigure 4.8: Entropy error from primal solutions and ETE corrected solutionfor the Gaussian bump test case.854.4. Compressible Viscous Flow10−210−710−610−510−410−3||s|| 2h = N−1/2 (2,0,0)(2,4,4) Linear(4,0,0)Figure 4.9: Accuracy of entropy norm for the Gaussian bump test case.N ,||s||2 (2, 0, 0) (2, 4, 4) (2, 4, 4)L (4, 0, 0)1599 1.14177e-4 6.46338e-5 6.09974e-5 6.46338e-56470 3.27943e-5 1.29657e-5 1.29734e-5 1.29657e-525750 8.16123e-6 5.62541e-7 5.65209e-7 5.62541e-7Order 2.0 4.5 4.5 4.5Table 4.2: Accuracy of entropy (error) for the different schemes for thelaminar Euler equations test case on the Gaussian bump geometry. Theorders of accuracy are calculated using the finest meshes considered.864.4. Compressible Viscous Flowwhere the heat flux is proportional to the gradient of temperature as givenby Fourier’s lawq =(qxqy)= −K∇T, (4.36)and the stress tensor isτ =(τxx τxyτyx τyy)= −23µ(∇ · v)1+ µ(∇v +∇vT ). (4.37)The same equation of state of Equation 4.30 applies. The solution algorithmis mostly the same as the previous cases. One difference is that for viscouscases, better conditioning was observed if the reconstruction is unweighted,or setting β = 0 in Equation 2.31, as per [60], whereas the inviscid cases useβ = 1. Furthermore, it was found that linear systems show better conver-gence behavior for these viscous cases using the Quotient Minimum Degree(QMD) ordering [42] instead of RCM ordering. The specific command linearguments for these test cases are summarized in Appendix C.4.4.1 Manufactured SolutionA good starting test case can be obtained by the method of manufacturedsolutions. The exact solution was selected to beρuvP =1 + 0.1 sinpix sinpiy2 + 0.1 sinpixsinpiy2 + 0.1 sinpix sinpiy (4.38)and the exact solution for the conserved variables is determined from com-binations of these. The results are presented for the energy component ofthe conserved variables, whose error was the largest in magnitude for thesetest cases. The convergence of the error estimate for a sequence of meshesis shown in Figure 4.10.4.4.2 NACA 0012 AirfoilIt is also possible to perform some of these comparisons for flows withouta known analytical exact solution. We consider the application of the ETEapproach to laminar flow around a NACA 0012 airfoil [59], with flow con-ditions of Reynolds number Re = 5000, Mach number M = 0.5, and angle874.4. Compressible Viscous Flow10−2 10−110−810−610−410−21004.214.112.31length scale||ǫp−ǫpqr|| (4,0,0)(2,4,4)(2,4,4), Linear(2,0,0)Figure 4.10: Convergence of the error estimate in nondimensionalized en-ergy for the 2D Navier-Stokes equations using manufactured solutions.of attack α = 0◦. This test case demonstrates the applicability of the ETE torealistic flow problems.The viscous drag coefficient (CD,Viscous) values for a sequence of meshesare shown in Table 4.3. Again the (2, 4, 4) values for each mesh are identicalto those obtained from the (4, 0, 0) scheme. This sequence of meshes doesnot have a constant refinement of the length scale. To perform a Richardsonextrapolation based estimate of the order of accuracy, the following equa-tions are solved, as described in [24]:p =1log(h2/h3)∣∣∣∣log ∣∣∣∣U1 − U2U2 − U3∣∣∣∣+ q(p)∣∣∣∣ (4.39)q(p) = log((h2/h3)p − ζ(h1/h2)p − ζ)(4.40)ζ = sgn(U1 − U2U2 − U3)(4.41)If the length scale is reduced by the same factor, Equation 4.40 simply eval-uates to zero. The results for the primal discretizations are close to what884.4. Compressible Viscous FlowN ,CD,Viscous (2, 0, 0) (2, 4, 4) (2, 4, 4)L (4, 0, 0)7440 0.030939 0.032211 0.032220 0.03221113120 0.031673 0.032571 0.032574 0.03257129760 0.032402 0.032743 0.032741 0.03274352480 0.032670 0.032779 0.032777 0.032779Order 1.8 3.3 3.3 3.3Table 4.3: Accuracy of CD,Viscous for the different schemes for the laminarNavier-Stokes test case on the NACA 0012 geometry.is expected, but is slightly under the respective nominal orders. Since theRichardson extrapolation techniques assumes that the coarsest mesh in thesequence is in the asymptotic region, it is possible that adding one morelevel of mesh resolution will make the orders of accuracy close to what isexpected. Nevertheless, from the point of view of the ETE methodology, theobjective has been achieved in that the results are extremely close to thefourth order primal discretization at a fraction of the computational cost.That is, there is no particular reason to expect that the solution correctedby the error estimate would outperform the primal fourth order scheme interms of accuracy.4.4.2.1 Efficiency ConsiderationsNext, we examine the efficiency of using the ETE in terms of computationaltime for a given level of error for the 2D Navier-Stokes equations. We com-pare the overall time required for a level of error resulting from solving onlythe primal problem, versus the time required for the same level of errorresulting from the defect corrected primal solution using the ETE. For thelinearized ETE, Equation 4.13 is solved as a single linear system. The per-formance tests were conducted using a single core on a typical workstationwith an Intel i7-4790 processor at 3.60GHz. It is emphasized that the cur-rent development can be completely parallelized as in the primal problem,but the efficiency comparisons are performed serially. Averaged timing re-sults for the manufactured solution case, with freestream initial conditions,an initial guess that is the typical for flow simulations, and an idealizedgood initial guess of the continuous exact solution projected on the discretespace, are shown in Figure 4.11. From these plots, it can be determinedhow much time is required for a certain level of accuracy, and the converserequirement. This visualization encapsulates the linearization error that is894.4. Compressible Viscous Flow10−2 10−1 100 101 10210−710−610−510−410−310−210−1Time (s)Error in E using exact IC (4,0,0) (2,4,4) (2,4,4), Linear (2,0,0)10−2 10−1 100 101 10210−710−610−510−410−310−210−1Error in E using freestream ICTime (s)Figure 4.11: Comparison of solution error plotted against overall computa-tional time required for the 2D Navier-Stokes equations using the manufac-tured solution test case, using freestream (left) and continuous exact (right)initial conditions. The (4, 0, 0) scheme did not converge for the finest meshusing freestream initial conditions.10−3 10−2 10−10.03050.0310.03150.0320.03250.0330.0335 C D,ViscousLength scale, N−0.5100 101 102 1030.03050.0310.03150.0320.03250.0330.0335Time (s) (4,0,0) (2,4,4) (2,4,4), Linear (2,0,0)Figure 4.12: Comparison of an output functional (CD,Viscous) plotted againstmesh length scale and overall computational time required for 2D Navier-Stokes equations for the NACA 0012 test case.904.4. Compressible Viscous Flowmade for the linearized ETE and shows the resulting performance.A few conclusions can be drawn. First, it can be observed that (4, 0, 0),(2, 4, 4), and linearized (2, 4, 4) schemes all give similar error magnitudesand are asymptotically better than (2, 0, 0), which supports the efficiencymotivation for using higher order methods. Furthermore, it can be seen thatthe linearized (2, 4, 4) scheme outperforms the alternatives, while the non-linear (2, 4, 4) and the (4, 0, 0) primal problem are comparable, and relativeperformance changes depending on choice of initial conditions, for example.The incremental time that is needed to solve the ETE can be determined bythe horizontal distance between the (2, 0, 0) points and the correspondingpoints on the ETE data. It can be seen that the extra time that is neededfor the linearized (2, 4, 4) scheme is never more than another (2, 0, 0) pri-mal solve. For cases where the primal solve takes longer, the linearized ETEtakes an even smaller fraction of the primal time. For example, for the finestmesh using freestream conditions, the linearized (2, 4, 4) scheme takes onlyapproximately an additional 10% of the primal time, but for an answer thatis asymptotically higher order. The reduction of error for this particularmesh size is a factor of more than 700, and increases as the mesh is refined,owing to the error estimate being higher order accurate and computation-ally affordable. Finally, the linearized (2, 4, 4) scheme gives asymptoticallythe same results as the (4, 0, 0) higher order primal scheme, but in an overalltime that is a factor of two or more faster. The gains are more apparent usingfreestream initial conditions than the idealized continuous exact solution asthe initial guess.For the NACA 0012 test case, whose results are plotted in Figure 4.12,there is an appreciable advantage in using the linearized (2, 4, 4) schemeonce again. For a given level of error the linearized (2, 4, 4) scheme givesthe same answer in a shorter time by about a factor of three than the (4, 0, 0)primal scheme. Conversely, for a given computational time budget, the lin-earized (2, 4, 4) scheme gives an answer that is much closer to the true value.4.4.2.2 Robustness ConsiderationsWe have demonstrated gains in efficiency in using the linearized ETE. Forthe (4, 0, 0) scheme in Figure 4.11, there is no data point for the finest meshusing freestream conditions, since it failed to converge in that case, whichmotivates our discussion about potential advantages in robustness in solv-ing the ETE instead of solving a higher order primal problem. As discussed,solving the linearized ETE requires only solving the lower order primal prob-lem and a single linear algebraic system. This should be more robust than914.4. Compressible Viscous Flow(a) Elliptic mesh in which the wake region isnot properly resolved.(b) Hyperbolic mesh in which the boundarylayer region is not properly resolved.Figure 4.13: Particularly unsuitable meshes which converge for the lin-earized (2, 4, 4) scheme but not (4, 0, 0) or (2, 4, 4).solving the nonlinear ETE, which as we have shown is equivalent to solv-ing the higher order primal problem with a lower order primal initial guess,which in turn should be more robust than solving the higher order primalproblem with freestream initial conditions.Shown in Figure 4.13 are two particularly poor meshes for viscous flowaround a NACA 0012 airfoil, where one is an elliptic mesh with the wakeregion poorly resolved, and the other is a hyperbolic mesh with the wallboundary layer anisotropy poorly captured. In both of these cases, the lin-earized ETE converges to an answer (albeit one that may not be very accu-rate), whereas the nonlinear ETE and higher order primal do not converge.It is also possible, as in the manufactured solution test case previously, thatfor some cases only the linearized ETE and nonlinear ETE converge butnot the higher order problem. This is just equivalent to cases where forthe higher order problem using a better initial guess converges but usingfreestream does not. We observe thus a hierarchy of increasing robustnessin the three schemes.924.5. Simulation of Turbulent Flows4.5 Simulation of Turbulent FlowsWith advances in both hardware and numerical algorithms throughout thelast several decades, much progress has been made to simulate more realis-tic aerodynamic models, notably in the simulation of turbulent flows, a char-acteristic of the chaotic behavior inherent to the nonlinearity of the Navier-Stokes equations. The use of the RANS equations to model the mean flow,in conjunction with a turbulence model, has been widely studied, includ-ing some models that have found much success in predicting observation.We take our ETE approach further by applying the methodology to turbu-lent flow. It is emphasized once again that in this context, even though themodeling errors can potentially dominate, the study of discretization erroris also important, since the fidelity of a particular turbulence model cannotbe reliably evaluated without an adequate understanding of discretizationerror.An abundance of turbulence models exists in the literature, with vary-ing degrees of complexity and success in different applications. A shortoverview can be found in [102]. The purpose in the current work is not toreview or develop new approaches in existing turbulence models. Rather,the applicability of the ETE methodology to complex models in CFD mod-eling is demonstrated through one particular turbulence model: the nega-tive variant of the Spalart-Allmaras (SA-neg) [5] one-equation model [103],sometimes referred to as the RANS-SA model. The application to this modelreinforces the idea that the ETE approach for accurately estimating dis-cretization error can be applied for complex model equations, and by ex-tension, even more complex ones that might be developed in the future.One major difference between this model and the previous simpler modelsis that the source term depends on the solution and gradients. This needs tobe taken into account when the ETE and its linearization is formulated. In-stances of the implementation of this turbulence model and can be found inthe literature [35, 55, 60, 72]. A brief outline of the SA-neg model equationsis presented, starting with the model equations∂tρρvEρν˜+∇ ·ρρvv + P1− τv(E + P )− τ · v + cp(µPr +µTPrT)∇Tρvν˜ + 1σ (µ+ µ′fnρν˜)∇ν˜ = S (4.42)934.5. Simulation of Turbulent FlowswithS =000ρ(P− D) + µ′σ cb2ρ∇ν˜ · ∇ν˜ − 1σ (ν + µ′ν˜)∇ρ · ∇ν˜ , (4.43)where µ′ = 1000 is a scaling factor to make the turbulence variable com-parable in magnitude to the other solution quantities [27], which improvesthe conditioning and hence convergence rate of the iterative linear systemsolver. In Equation 4.43, P and D are the eddy viscosity production anddestruction terms, defined asP ={cb1(1− ft2)S˜ν˜, ν˜ ≥ 0cb1(1− ct3)Sν˜, ν˜ < 0 (4.44)D =µ′(cw1fw − cb1κ2ft2)( ν˜d)2, ν˜ ≥ 0−µ′cw1(ν˜d)2, ν˜ < 0. (4.45)The vorticity S and its modified versions are defined asS˜ =S + S, S ≥ −cv2SS +S(c2v2S + cv3S)(cv3 − 2cv2)S− S, S < −cv2S (4.46)S =12∑i∑j(∂vi∂xj− ∂vj∂xi)21/2 (4.47)S =µ′ν˜fv2κ2d2. (4.48)Other functions in the model are defined asfv1 =χ3χ3 + c3v1(4.49)fv2 = 1− χ1 + χfv1(4.50)χ =µ′ν˜ν(4.51)944.6. Zero Pressure Gradient Flat PlateName Value Name Value Name Valueσ 0.66 cb1 0.1355 ct3 1.2κ 0.41 cb2 0.622 ct4 0.5cv1 7.1 cw1 cb1κ2 +1+cb2σ cn1 16cv2 0.7 cw2 0.3 Pr 0.72cv3 0.9 cw3 2.0 PrT 0.9Table 4.4: Empirical values for the Spalart-Allmaras turbulence model.fn =1, ν˜ ≥ 0cn1 + χ3cn1 − χ3 , ν˜ < 0 (4.52)ft2 = ct3e−ct4χ2(4.53)r = min(µ′ν˜S˜κ2d2, 10)(4.54)g = r + cw2(r6 − r) (4.55)fw = g(1 + c6w3g6 + c6w3). (4.56)The turbulent eddy viscosity isµT ={µ′ρν˜fv1, ν˜ ≥ 00, ν˜ < 0. (4.57)The parameters used in these equations are summarized in Table 4.4.4.6 Zero Pressure Gradient Flat PlateFirst, we consider turbulent flow on a 2D turbulent flat plate geometry,which is a classical verification benchmark test case. Our results will becompared to those of CFL3D [64] and FUN3D [17] legacy codes from theNASA Turbulence Modeling Resource [94]. These reference codes havedata for this test case for the finest mesh with 208896 cells. The nega-tive Spalart-Allmaras turbulence model of Equation 4.42 is used in a do-main Ω = {(x, y) ∈ [−0.33, 2] × [0, 1]}, where the flow parameters include954.6. Zero Pressure Gradient Flat PlateN ,CD (2, 0, 0) (2, 4, 4)L (4, 0, 0)2040 0.003617 0.002373 0.0023238160 0.002980 0.002846 0.00284932640 0.002872 0.002862 0.002862Order 2.6 4.9 5.4Table 4.5: Accuracy of CD for the different schemes for the flat plate testcase.the Reynolds number Re = 5 × 106, the Mach number M = 0.2, and theturbulence variable scaled as ν˜/ν∞ = 3.0. The geometry and boundary con-ditions are shown in Figure 4.14. To evaluate the accuracy of the differentschemes, the drag coefficient is monitored. In this case, only the viscousdrag is nonzero, so this is calculated from the wall shear stress τw asCD = CD,Viscous =τw12ρ|v∞|2(4.58)where v∞ is the freestream velocity. Equation 4.58 is integrated along theflat plate to get the total drag coefficient. In Figure 4.15, it can be seen thatthe turbulence variable is dramatically improved with the (2, 0, 0) solutionis corrected by the (2, 4, 4)L ETE error estimate, bringing it rather close tothe (4, 0, 0) primal solution. It should also be noted that both the (2, 4, 4)Land (4, 0, 0) schemes capture the peak values quite well in comparison to thereference solutions which have about 26 times as many degrees of freedom.These results demonstrate the superior accuracy of higher order methodsover the (2, 0, 0) solution. A singularity in the solution is produced at theleading edge, and some difficulties showing the nominal order of accuracyof the reference codes have been reported [94]. Nevertheless, we are able toshow order of accuracy results in Table 4.5, plotted in Figure 4.16, namelythat the linearized (2, 4, 4) is an asymptotic improvement to the second or-der primal scheme in such a way that the corrected solution is very close tothe fourth order primal solution. Figure 4.17 shows the skin friction coeffi-cient distribution and µT/µ∞ at x = 0.97 for the primal and ETE correctedschemes, for a relatively coarse mesh of 8160 quadrilaterals, to adequatelyshow the difference between the schemes. It can be seen that the (2, 4, 4)Lscheme is extremely close to the (4, 0, 0) scheme.964.7. Flow over a NACA 0012 AirfoilXY-1 -0.5 0 0.5 1 1.5 2 2.5-0.500.511.522.5Figure 4.14: Geometry and boundary conditions for the flat plate test case.4.7 Flow over a NACA 0012 AirfoilThe NACA 0012 airfoil geometry is used next as a test case for turbulentflow. The flow parameters are Re = 6 × 106,M = 0.15, and ν˜/ν∞ = 3.0.The angle of attack is α = 10◦. Computations were performed with thefinest mesh having 100352 control volumes. No-slip adiabatic boundaryconditions are used on the surface of the airfoil, and characteristic boundaryconditions are used at the farfield at about 500 chords away, where inflow oroutflow is determined based on the local velocity variables. The turbulencevariable, along with the ETE error estimate and correction using the ETEwith (2, 4, 4)L are plotted in Figure 4.18. It can be seen from Figure 4.19that the (2, 0, 0) solution is much further from the (4, 0, 0) solution thanthe (2, 4, 4)L corrected solution is, warranting our decision to relinearizethe ETE once more. The corrected solution using this error estimate from(2, 4, 4)L gives a result that is remarkably close to the fourth order solution,with only two additional linear system solves.The lift coefficient (CL) functional for the different schemes are shownin Table 4.6. For this set of data, it can be seen that the results are not quiteasymptotic, especially for the coarsest mesh using (4, 0, 0). Nevertheless, itcan be seen that for the finest mesh considered, the (2, 4, 4)L is discerniblycloser to the (4, 0, 0) result than (2, 0, 0). Next the profiles of the CP and Cf974.7. Flow over a NACA 0012 Airfoil(a) CFL3D reference (b) FUN3D reference(c) (2, 0, 0) (d) (2, 4, 4)L(e) (4, 0, 0)Figure 4.15: Plot of the turbulence variable µ′ν˜ for N = 8160 and compar-ison with reference solutions with N = 209825. It is qualitatively evidentthat the (2, 0, 0) is the least accurate, and the ETE corrected solution is veryclose to both (4, 0, 0) and the reference solutions. 984.7. Flow over a NACA 0012 Airfoil0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.042.22.42.62.833.23.43.6x 10−3C Dh = N−1/2 (2,0,0)(2,4,4) Linear(4,0,0)CFL3DFUN3DFigure 4.16: Accuracy of CD for the flat plate test case.994.7. Flow over a NACA 0012 Airfoil0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00x0.00200.00250.00300.00350.00400.00450.00500.00550.0060Cf(2,0,0), N = 8,160(2,0,0), N = 130,560(4,0,0), N = 8,160(2,4,4)L, N = 8,160Reference(a) Cf , FUN3D reference0 25 50 75 100 125 150 175 200μT/μ∞0.00000.00250.00500.00750.01000.01250.01500.01750.0200y(2,0,0), N = 8,160(4,0,0), N = 8,160(2,4,4)L, N = 8,160FUN3D referenceCFL3D reference(b) µT /µ∞ at x = 0.97Figure 4.17: Functionals for the flat plate test case.1004.7. Flow over a NACA 0012 AirfoilXY-0.16 0.1 0.36 0.62 0.88 1.14(a) (2, 0, 0)XY-0.17 -0.05 0.07 0.19 0.31 0.43 0.55(b) Error estimate from (2, 4, 4)LXY-0.16 0.1 0.36 0.62 0.88 1.14(c) Corrected solution using (2, 4, 4)LXY-0.16 0.1 0.36 0.62 0.88 1.14(d) (4,0,0)Figure 4.18: Turbulence variable ν˜, error estimates, and corrected solutions.The error estimate is large in some parts of the domain, warranting the needfor ETE relinearization.1014.7. Flow over a NACA 0012 AirfoilXY-0.2 -0.124 -0.048 0.028 0.104 0.18(a) Difference between (2, 0, 0) and (4, 0, 0)XY-0.2 -0.124 -0.048 0.028 0.104 0.18(b) Difference between (2, 4, 4)L correctedand (4, 0, 0)Figure 4.19: Estimate of error in the turbulence variable ν˜ for the second or-der primal solution and the linearized ETE corrected solution. This demon-strates that the first linearization error is large and can potentially dominate.functionals are examined, which are defined asCP =P − P∞12ρ|v∞|2(4.59)Cf =τw12ρ|v∞|2. (4.60)It can be seen from Figure 4.20 that for the finest mesh considered, the(2, 4, 4)L and (4, 0, 0) are closer to the reference solutions than (2, 0, 0),which can be seen in the insets in the figures.1024.7. Flow over a NACA 0012 Airfoil0.0 0.2 0.4 0.6 0.8 1.0x−6−5−4−3−2−101Cp(2,0,0), N = 100,352(4,0,0), N = 100,352(2,4,4)L, N = 100,352Reference(a) Coefficient of pressure distribution, CFL3D reference0.0 0.2 0.4 0.6 0.8 1.0x−0.010−0.0050.0000.0050.0100.0150.0200.0250.030Cf(2,0,0), N = 100,352(4,0,0), N = 100,352(2,4,4)L, N = 100,352Reference(b) Skin friction coefficient distribution, FUN3D referenceFigure 4.20: Functional accuracy for the NACA 0012 test case. The refer-ences use N = 230529.1034.7. Flow over a NACA 0012 AirfoilCFL3D 1.0909FUN3D 1.0983N ,CL (2, 0, 0) (2, 4, 4)L (4, 0, 0)6272 1.0781 1.0822 0.582825088 1.0779 1.0635 1.0952100352 1.0882 1.0913 1.0905Table 4.6: Accuracy of CL for the NACA 0012 test case, along with referencevalues.104Chapter 5ETE Methodology Applied toUnsteady ProblemsFor many applications of interest the model equations are inherently un-steady. In this chapter we shall examine how the error estimation method-ology by means of the ETE can apply to two classifications of unsteadyproblems: those that are inherently unsteady and those that reach a time-periodic steady state. In the results, we will investigate the error estimateaccuracy in the former case for the solution variables to demonstrate theETE method, and in the latter case we compare typical aerodynamic func-tionals using corrected solutions as a proxy to deduce the accuracy in theerror estimates.5.1 Time-Dependent Residual EvaluationRecall that we have formulated the primal problem in Equation 2.34 andETE in Equation 2.35 after discretizing in space. The objective is to solvethese equations by choosing an accurate source term T that approximatesthe true residual as the source term for the ETE. For steady problems, it wasevident that careful evaluation of the spatial part of the ETE source termwas critical to the accuracy of the error estimate; we expect that discretiza-tion choices for the time-dependent term in unsteady problems would beimportant in an analogous way. We give here a survey of several estimatesof the ETE source term T before discussing and showing results for thefinite-difference in time ETE source term that gives an accurate estimate.5.1.1 Source Term for Exact Error EstimateStarting from Equation 2.40, it is instructive to first see what the targetresidual source term for the ETE is. To see what discrete ETE source termwould be required for the corrected solution V to be the exact solutionaverage U , we simply make this substitution in Equation 2.40 to obtain1055.1. Time-Dependent Residual EvaluationTn+1 =Un+1 − Unk− fq(Un+1)− fp(Un+1) + fq(Un+1) (5.1)for backward Euler andTn+1 + Tn2=Un+1 − Unk+12(−fq(Un+1)− fp(Un+1) + fq(Un+1)−fq(Un)− fp(Un) + fq(Un) (5.2)for Crank-Nicolson from Equation 2.41. These are the ETE source terms thatgive the exact solution after the ETE correction, but these are not practicallycomputable since the exact solution U is, of course, unknown. For practicalETE source terms, the objective can then be thought of as finding an accurateapproximation of this exact unsteady ETE source term. Several possibilitiesare examined as approximations to the ETE source term and the implicationsfor the error estimate are investigated.5.1.2 T = 0 Gives the Trivial Zero EstimateThe error estimate is zero in the simplest case of using zero for the ETEsource term. To see this, consider the backward Euler discretization, andobserve that Equation 2.40 becomesV n+1 − V nk= fq(Vn+1) + fp(Un+1)− fq(Un+1). (5.3)When the residual source T = 0, having an uncorrected solution V = Uwill satisfy Equation 5.3, because this just reduces to Equation 2.36 for theprimal solution U , hence is zero. A similar argument can be made for theCrank-Nicolson scheme.5.1.3 Higher Order Error Estimate in Space OnlyFor either backward Euler or Crank-Nicolson, using Tn = fq(Un) − fp(Un)with q > p will give higher order accuracy in space only. For example, usingp = 2, q = 4, pt = qt = 2, this gives an identical error in the error estimate asthe primal error using p = 4, pt = 2, where the time accuracy is unchanged.To see this, Equation 2.41 for the Crank-Nicolson discretization becomesV n+1 − V nk=12(fq(Vn+1) + fq(Vn)) (5.4)1065.1. Time-Dependent Residual Evaluationwhich means the corrected solution satisfies the space discretization to orderq but time accuracy is unchanged, meaning ||U −V || = O(hq) +O(k2). Thisis a direct consequence of what was concluded in Equation 4.24 for steadyproblems when q = r, namely that the corrected solution is identical to aprimal solution to order q. A diffusion test case in 2D is used, as∂tu−∆u = s, on Ω× ΩT = {(x, y, t) ∈ [0, 1]2 × [0, 1]} (5.5)with Dirichlet boundary conditions u|∂Ω = 0. The effect of using this ETEsource term in shown in Figure 5.1 for q = 4, and the expected order ofaccuracy of the resulting error estimate has been verified in both space andtime. In this case the difference in error is dominated by the time discretiza-tion, hence the smooth nature of the uniform time discretization prevailsin the profile. From a standard truncation error analysis, the leading orderterm can be shown to depend on the derivatives of the solution as a smoothsignal, where the ETE source term gives a discretization error estimate withan accuracy of O(h4) +O(k2), and a clear sinusoidal pattern can be seen inthe plot for the error in error estimate. For cases with errors that are notfully dominated by the time discretization, as we see later when higher or-der accuracy in both space and time is achieved, we would get error patternsthat resemble the smooth spatial derivatives of the solution only if the meshhad smoothly varying geometrical features. In our context of unstructuredmeshes with no underlying smooth mesh function, the mesh irregularitycreates no discernible signal in the error patterns as derivatives of solutions.If spatial higher order accuracy is sufficient, this estimate of the ETEsource term can be used. However, as seen in the steady case, without lin-earizing the ETE for nonlinear problems, the cost of solving the primal prob-lem and ETE with p = 2, q = r = 4 is not very competitive. For the unsteadycase, direct linearization in the same way is not efficient, since the Jacobianwould not be constant in time and a different matrix and preconditionerswould need to be evaluated at each time step, which is prohibitive. Hencewe look for alternative ETE source terms that give an error estimate that ishigher order accurate in time, which would be more promising candidatesas efficient schemes.5.1.4 MEA Source TermAnother approach for the unsteady ETE source term is motivated by modi-fied equation analysis (MEA) [114]. We consider this approach because itis a common approach for smooth meshes [10]. First we observe that after1075.1. Time-Dependent Residual Evaluation(a) Exact error (b) Error estimate (c) Error in error estimateFigure 5.1: The diffusion problem using p = 2, q = 4, pt = qt = 2 andfq(U)− fp(U) as the ETE source.discretizing in space, a system of ODEs, Equation 2.34, is obtained. Withsufficient smoothness in the right hand side when a sequence of smallerh and k is taken, the standard MEA approach is used to exchange temporalderivatives in the error in favor of spatial derivatives. The approach is to usethe leading order term in a truncation error analysis as the ETE source term.To demonstrate this method, we start with the backward Euler discretiza-tion, Equation 2.36, and expand in a Taylor series about Un, collecting termsin powers of k. Using the shorthand u = Un, f = f(Un), ut = ∂tUn, andsimilarly for higher derivatives, we haveu + kut + k2utt/2− uk+O(k2) = f(u + kut +O(k2)), (5.6)from which we findut +k2utt = f + kfuut +O(k2). (5.7)From this we know that ut = f +O(k), and hence utt = ft +O(k) = fuut +O(k) = fuf + O(k) from the chain rule, where fu is the Jacobian matrix attime level n. Substituting and rearranging, this givesut − f = k2fuf +O(k2). (5.8)The ETE source term we use is1085.1. Time-Dependent Residual EvaluationTn+1 = −k2(fUf)|n+1, (5.9)and using this ETE source term will give second order accuracy in time. Ina similar manner, using Crank-Nicolson it can be shown thatut − f = k212(fuuf2 + f2uf) +O(k4) (5.10)The ETE source term used is thusTn+1 = − k212(fUUf2 + f2Uf)|n+1Tn = − k212(fUUf2 + f2Uf)|n. (5.11)For nonlinear problems, the Hessian fUU , a third rank tensor, is nonzeroand an inner product with f is taken twice, denoted by f2. One approach isto directly compute the Hessian and compute the inner products. Anotherapproach that avoids the explicit computation of a third rank tensor is tofirst observe thatfUUf2 + f2Uf = (fUf)Uf. (5.12)Practically, this can be computed as a Fre´chet derivative in the direction off ,(fUf)(U + δU)− (fUf)(U)δf (5.13)which can be evaluated as(fUf)(U + f(U)δ)− (fUf)(U)δ(5.14)as a finite-difference, for small δ ≈ 10−8.To demonstrate this method, consider a simple ODE system for u(t) =[ v(t) w(t) ]Tdudt=ddt[vw]=[ −2/(vw4)1/(vw)](5.15)with initial conditions u(0) = [ 1 1 ]T . The exact solution u(t) = [ e−2t et ]Tsatisfies the ODE system. On the domain t ∈ [0, 1], the forcing term on the1095.1. Time-Dependent Residual Evaluationright is smooth, unlike the PDE discretization on an unstructured mesh. Theprimal and ETE problems are discretized with pt = 2, qt = 2.It can be seen in Figure 5.2 that the primal solution is second orderaccurate, and using the MEA ETE source term the error estimate is fourthorder accurate.0 0.2 0.4 0.6 0.8 10123v(t),w(t) v(t) exactw(t) exactv(t) computedw(t) computed0 0.2 0.4 0.6 0.8 1−6−4−202x 10−3tPrimal error v(t) error exactw(t) error exactv(t) error computedw(t) error computed(a) ODE system solution and error.10−2 10−110−1010−810−610−410−2k Primal errorError in error estimatesecond orderfourth order(b) Norms of the primal error and error in error.Figure 5.2: Accuracy of error estimate using the MEA source for a simpleODE system.We then turn to applying this for a PDE discretization for unstructuredmeshes, which gives a system of ODEs after discretizing in space. For alinear problem, we have fUU = 0, and fUη = f(η) for any admissible η.Then Equation 5.11 becomesTn+1 = −k212f3(Un+1)Tn = −k212f3(Un) (5.16)with f3 indicating composition of f three times. For a (2, 4, 4; 2, 2, 4) schemefor diffusion however, we do not get an asymptotically correct error estimateusing MEA, as seen in Figure 5.3. The qualitative difference between thisand the model ODE system is the source term for the ETE on the right,rather than the scheme itself being unstable. For the residual right handside of a second order scheme for diffusion, the asymptotically large noisedominates any underlying signal as a sequence of finer unstructured meshesis used, as was seen in the steady case. Hence it is not surprising that1105.1. Time-Dependent Residual Evaluationrepeated composition of f (three times for Crank-Nicolson) for diffusionusing a second order primal scheme gives poor error estimates. The normof this quantity at the initial step is shown in Figure 5.4.(a) Exact error (b) Error estimate (c) Error in error estimateFigure 5.3: The diffusion problem using p = 2, q = 4, pt = qt = 2 and MEAETE source.5.1.5 Computing the Time Derivative in −(∂tu˜+R(u˜)) Directlyby Finite-DifferenceThe method that gives accurate error estimates, which we demonstrate inthe results section, is to apply finite-differencing in time, which we refer toas the FD method. As in the steady case, the spatial residualR is computedusing the reconstruction of the primal solution to order r, which is readilycomputable since solving the primal problem gives the solution at each timestep, Un, n = 0, 1, ...NT . The time derivative ∂tu˜ can be computed usinga finite-difference in time of the primal solution at nearby time steps. Forinterior nodes in time a centered finite-difference operator Drt to order rt isused, and the ETE source isTn = −DrtUn + fr(Un). (5.17)In this study, rt = 4 is kept fixed, which is sufficiently accurate for all testcases, and abbreviate the notation as (p, q, r; pt, qt). In the first and lastfew time steps, there is an insufficient amount of data in one time directionto use a fully centered finite-difference approximation in time. One-sideddifferences of sufficient accuracy are used in those cases.1115.2. Expected Error Estimate Accuracy Using FD ETE Source10−2100101102103hinitial norm of MEA source term uniform meshunstructured meshFigure 5.4: Norms of the MEA source term ||k2f3(U)/12|| at the initial time,keeping k/h approximately fixed.To gain more insight into the error behavior as a function of time, weconsider estimating the norm ||e||2 :=(1|Ω|∫Ω e2dx)1/2as t increases. Dis-cretely we compute this as||||2 :=(∑i|Ωi|2i|Ω|)1/2. (5.18)Figure 5.5 demonstrates this for the exact error, and error estimates us-ing the MEA and FD source terms for the ETE. We observe that the MEAerror estimate is not unstable in the sense of unbounded growth in time,but rather the MEA source term is just not asymptotically correct for un-structured meshes due to the repeated application of the residual operator,whereas using the FD source term for the ETE gives an error estimate thatmatches the exact error quite well, as explored further in the results.5.2 Expected Error Estimate Accuracy Using FD ETESourceTo demonstrate the expected accuracy of the error estimate using the FDETE source term, we start from Equation 2.38, expanding the Crank-Nicolson1125.2. Expected Error Estimate Accuracy Using FD ETE Source0 2 4 6 8 1010−810−610−410−2100102t||ǫ||2 Exact errorError estimate (FD ETE source)Error estimate (MEA ETE source)Figure 5.5: Plot of the evolution in time of the L2 norm of the exact error,and error estimates using the MEA and FD source terms for the ETE.discretization in time to obtain the leading order error term1k(Un+1 − Un − k2(f(Un) + f(Un+1)))= −k212∂tttUn +O(k3). (5.19)In the FD approach, we have DrtUn = ∂tUn+O(krt) in Equation 5.17. ThenEquation 2.41 becomesV n+1 − V nk=12(fq(Vn)− fq(Un) + fq(V n+1)− fq(Un+1)) + Un+1 − Unk+12(−DrtUn + fr(Un)−DrtUn+1 + fr(Un+1)). (5.20)For q = r, this simplifies toV n+1 − V nk=12(fq(Vn) + fq(Vn+1)) +Un+1 − Unk+12(−DrtUn −DrtUn+1). (5.21)We expand about u := Un to getV n+1 − V nk=12(fq(Vn) + fq(Vn+1)) + ut +k2utt +k26uttt +12(−ut −(ut + kutt +k22uttt)). (5.22)1135.3. Co-Advancing the ETEThis can be written asV n+1 − V nk=12(fq(Vn) + fq(Vn+1))− k212uttt, (5.23)which means the FD ETE source term corrects the primal scheme by theterm −k212 uttt. In the same way, the accuracy of this corrected solution canbe determined by performing a truncation error analysis. Equation 5.23 canbe rearranged, and using Equation 5.19 for V , which gives1k(V n+1 − V n − k2(fq(Vn) + fq(Vn+1)))+k212uttt =−k212∂tttVn +k212∂tttUn = −k212∂tttn = O(hpk2) +O(k4) (5.24)since = O(hp) +O(k2). The accuracy of the forcing term is O(hq), whichsuggests that V is accurate to O(hpk2) +O(k4) +O(hq). A similar analysiscan be performed for different discretizations in time, from which we expectthe accuracy of the corrected solution and hence the error estimate to be||U − V || = O(hpkqt) +O(kpt+qt) +O(hq). (5.25)The assumption here is that the solutions have sufficient smoothness in timeso that higher time derivatives remain bounded. Equation 5.25 gives muchinsight into the expected accuracy of the error estimate and corrected solu-tion. The cross term containing both h and k is not at all obvious. In theresults, we confirm this result numerically, and in particular we can see thata scheme with (2, 4, 4; 2, 2) leads to a fourth order accurate error estimate,or equivalently a fourth order accurate corrected solution in both space andtime.5.3 Co-Advancing the ETEThe primal problem and ETE can of course be solved sequentially, in thesense of advancing to the final time for the primal problem, then advanc-ing to the final time for the ETE. However, whereas for steady problems theETE source term is calculated from a single evaluation of the converged pri-mal solution via p-truncation error, for unsteady problems the entire historywould need to be stored in order to compute the time-dependent ETE sourceterm if this sequential approach were taken. For large unsteady problemswith many time steps, the memory required to evaluate this ETE source termbecomes prohibitive.1145.4. Accuracy Comparisons for Error in SolutionIn fact, this large storage requirement is one of the largest difficultieswhen using unsteady adjoint methods. Even in some recent advances incheckpointing [105], there are still some difficulties for problems with manytime steps.An approach for the unsteady ETE that can reduce this memory require-ment is to advance in time both the primal problem and ETE concurrently,so that only a few local solutions and residuals in time need to be stored.One way of achieving this is to advance the primal problem forward at leastrt + 1 steps, with the ETE lagging. This is required for the FD approachof the previous section when using a centered difference in time. Once asufficient number of primal steps is taken, the primal problem and ETE canbe advanced forward together one step at a time, computing the ETE sourceterm Equation 5.17, until the primal problem reaches the final time. The al-gorithm terminates with the ETE advancing until the final time. A schematicof this is shown in Algorithm 1. A data structured that is used to achieve thisis a double-ended queue (deque) [63], which is implemented in many mod-ern programming languages, allowing push and pop operations from eitherend. This is a natural data structure for this application, because it sup-ports computationally inexpensive enqueue and dequeue operations, whichis needed as solutions at time steps are computed in the problem problem,used for residual source term calculation, and discarded. This is an inter-esting approach, because in principle it is possible to parallelize this in timedue to the delay in time step information. In the current implementation, itis only parallel in the mesh processing and linear algebra operations throughPETSc.5.4 Accuracy Comparisons for Error in SolutionThe results for using the FD method to compute the ETE source term arepresented. The chosen test cases are relatively well behaved in time. Casesthat exhibit chaotic behavior are not addressed in this work, but an outlookin the context of the ETE approach is discussed in Chapter 6.In this section, the asymptotic accuracy of errors and error estimatesevaluated by taking the difference between the exact and computed so-lutions at the final time is used before examining functionals in the nextsection. These results serve two purposes: first to verify the expected ac-curacy given by Equation 5.25 by varying the choices of p, q = r, pt, andqt, and to demonstrate the practicality of solving the ETE to get a higherorder accurate error estimate using models of compressible flow for a par-1155.4. Accuracy Comparisons for Error in SolutionAlgorithm 1 Co-advance pseudocode for solving the unsteady primal andETE1: for n=1,...nOffset do2: U ← NEXTPRIMALSTEP3: ENQUEUE(U)4: end for5: for n=nOffset+1,...nSteps do6: U ← NEXTPRIMALSTEP7: ENQUEUE(U)8: DEQUEUE9: T ← COMPUTERESETESOURCE(U)10: ← NEXTETESTEP(T )11: end for12: for n=nSteps+1,...nSteps+nOffset do13: DEQUEUE14: T ← COMPUTERESETESOURCE(U)15: ← NEXTETESTEP(T )16: end forticular combination, (2, 4, 4; 2, 2), obtaining fourth order accuracy in bothspace and time, in support of the theory in the previous sections.The advection, diffusion, and Navier-Stokes equations are used, and re-finement studies are performed by dividing both k and h by two, therebykeeping k/h constant for each test case. If the ansatz of Equation 5.25holds, the leading order terms are c1h4 + c2h2k2 + c3k4, for some constantsc1, c2, and c3. Fixing k/h and subsequently dividing each by two duringrefinement will eventually show the desired asymptotic behavior regardlessof the relative magnitudes of the constants c1, c2, and c3. The question iswhether the magnitudes are close enough so that the asymptotic range canbe reached. One way of getting some insight into the relative magnitudesof the constants is to fix h and taking a very small k, and vice versa. Byfixing h and taking a very small k, c1 can be approximated, and by fixing kand taking a very small h, c3 can be approximated. It was found that theseconstants are close in order of magnitude.5.4.1 DiffusionWe first consider the unsteady diffusion problem in 2D, defined in Equation5.5. The exact solution is u(x, y, t) = sin(pix) sin(piy) sin(t) and the source1165.4. Accuracy Comparisons for Error in Solutionterm is s(x, y, t) = sin(pix) sin(piy)(cos(t) + 2pi2 sin(t)). Plots of the exacterror, estimate, and difference for this case are shown in Figure 5.6. Theerror estimate and exact error are very close to each other, and for thisspace and time resolution the error in the error estimate is approximatelythree orders of magnitude smaller.The asymptotic behavior of the lower order primal error, the error inerror, and higher order primal error for diffusion are shown in Figure 5.7. Itcan be observed that the ETE approach gives results that are quite close towhat is obtained from a fully higher order primal discretization.(a) Exact error (b) Error estimate (c) Error in error estimateFigure 5.6: The diffusion problem using (2, 4, 4; 2, 2) FD time derivative ETEsource, with h ≈ 0.0059, k = 0.025.5.4.2 AdvectionNext, the linear advection problem is considered, as∂tu+ c · ∇u = 0, on Ω× ΩT = {(x, y, t) ∈ [0, 1]2 × [0, 0.2]} (5.26)with c =(1 0)T , initial conditions u(x, y, 0) = e−λ((x−x0)2+(y−y0)2), wallconditions u(x, 0, t) = u(x, 1, t) = 0, inflow at u(0, y, t), outflow at u(1, y, t),and parameters λ = 500, x0 = 0.4, and y0 = 0.5. A constant velocity c withonly one nonzero component was selected to simplify the formulation ofboundary conditions in this test case. Since there is no particular directionalignment of the mesh, we do not expect fortuitous cancellations comparedto the full 2D model. We will verify later that the same conclusions hold forthe more difficult Navier-Stokes model. The exact solution for the advection1175.4. Accuracy Comparisons for Error in SolutionCorrecteddiscretizationerror10−2 10−110−1010−910−810−710−610−510−410−310−2k Diffusion exact error (lower order primal)Diffusion error in errorDiffusion exact error (higher order primal)second orderfourth orderFigure 5.7: A comparison of lower order primal error, error in error usingETE, and higher order primal error for diffusion using (2, 4, 4; 2, 2), withrefinement in both k and h, k/h ≈ 4.2 fixed.case is u(x, y, t) = e−λ((x−x0−t)2+(y−y0)2). Plots of the exact error, estimate,and error in error for this case using (2, 4, 4; 2, 2) are shown in Figure 5.8.It can be seen that the error estimate is remarkably close to the exact error.Another observation is that the error in the error estimate is approximatelyan order of magnitude smaller than the error itself for this choice of h andk, and the difference compared to the exact error increases asymptoticallysince the error in error estimate is fourth order as opposed to second orderfor the error.The lower order primal error, the error in error, and higher order primalerror are compared as a sequence of h and k is taken. These results foradvection are shown in Figure 5.9. It can be observed that the ETE approachgives results that are quite close to what is obtained from a fully higher orderprimal discretization, which would require a higher order discretization inboth space and time.For this test case, we also varied the choices of p, q = r, pt, and qt, andthe resulting norms of the exact primal error and the error in error estimatesare shown in Figure 5.10. It can be seen that the accuracy of the errorestimate matches Equation 5.25 quite well, and some combinations usedhere were chosen so that the h and k cross term is dominant. As predicted,1185.4. Accuracy Comparisons for Error in Solution(a) Exact error (b) Error estimate (c) Error in error estimateFigure 5.8: The advection problem using (2, 4, 4; 2, 2) FD time derivativeETE source, with h ≈ 0.0059, k = 0.0025.Correcteddiscretizationerror10−3 10−210−610−510−410−310−210−1k Advection exact error (lower order primal)Advection error in errorAdvection exact error (higher order primal)second orderfourth orderFigure 5.9: A comparison of lower order primal error, error in error usingETE, and higher order primal error for advection using (2, 4, 4; 2, 2), withrefinement in both k and h, k/h ≈ 0.42 fixed.the (1, 3, 3; 2, 1) scheme gives close to a second order accurate error estimatewhile (1, 3, 3; 1, 2) gives third order precisely because of the difference in thiscross term, even though the computational work is approximately the same,a result that is difficult to deduce without careful examination.1195.4. Accuracy Comparisons for Error in SolutionCorrecteddiscretizationerror10−3 10−210−510−410−310−210−1k (1,−,−;2,−) primal error(1,−,−;1,−) primal error(2,−,−;1,−) primal error(2,−,−;2,−) primal error(1,3,3;2,1) error in error (2)(1,3,3;1,2) error in error (3)(2,3,3;1,2) error in error (3)(2,3,3;2,1) error in error (3)(2,4,4;2,2) error in error (4)first ordersecond orderthird orderfourth orderFigure 5.10: Plots of the norms of exact primal error and error in error es-timate for the advection problem using different discretization orders. Thevalue in the parentheses indicates the expected accuracy of the error esti-mate from Equation 5.25.5.4.3 Navier-Stokes EquationsThe Navier-Stokes equations are described in Section 4.4. The test case weselect has a translating vortex exact solution, described in [67], computedasu(x, y, t)a∞=u∞a∞− G2pia∞y¯e−Br¯/2v(x, y, t)a∞=v∞a∞+G2pia∞x¯e−Br¯/2T (x, y, t)T∞= 1− G2(γ − 1)8pi2Ba2∞e−Br¯2ρ(x, y, t)ρ∞=(T (x, y, t)T∞)1/(γ−1)P (x, y, t)P∞=(T (x, y, t)T∞)γ/(γ−1)(5.27)on a domain Ω× ΩT = {(x, y, t) ∈ [0, 1]2 × [0, 0.1]}, with1205.5. Accuracy Comparisons for Error in FunctionalsName Value Name Valuex0 0.4 u∞ 1.01y0 0.4 v∞ 1.01B 500 P∞ 1/γG 100 T∞ 1/γρ∞ 1 a∞√γP∞/ρ∞ = 1Table 5.1: Summary of constants used for the translating vortex test case.x¯ = x− x0 − u∞ty¯ = y − y0 − v∞tr¯ =(x¯2 + y¯2)1/2. (5.28)This is an exact solution to the inviscid Euler equations, but we can usethis for the Navier-Stokes model by applying the method of manufacturedsolutions. The source term that results from this exact solution is used inthe Navier-Stokes test case.The nondimensionalized solution uses the constants shown in Table 5.1.Plots of the exact error, estimate, and error in error for this case are shownin Figure 5.11. The accuracy of the primal solution and error estimate for allof these 2D test cases is summarized in Figure 5.12. As in the steady case,for each data point, the vertical distance between the lower order primalerror and error in the error estimate can be interpreted as a measure of howgood the error estimate is, and is also the amount by which the solutionwill improve if the error estimate is used as a correction, which followsfrom Equation 2.22. It can be seen that the primal solution is second orderaccurate, while the error estimate is fourth order accurate. These casesdemonstrate that a higher order accurate error estimate is possible withoutdiscretizing both space and time to that full order.5.5 Accuracy Comparisons for Error in FunctionalsHaving explored the higher order accuracy properties of the unsteady ETEapproach for solutions at the final time, we demonstrate the extension totime-periodic problems. Here we indirectly measure the order of accuracyof the ETE error estimate, through functionals often encountered in appli-cation problems, calculated from the ETE corrected solution. We note that1215.5. Accuracy Comparisons for Error in Functionals(a) Exact error (b) Error estimate (c) Error in error estimateFigure 5.11: Energy variable for the Navier-Stokes test case using(2, 4, 4; 2, 2) FD time derivative ETE source, with h ≈ 0.0059, k = 0.0025.Correcteddiscretizationerror10−3 10−210−710−610−510−410−310−210−1k Navier−Stokes exact error (lower order primal)Navier−Stokes error in errorNavier−Stokes exact error (lower order primal)second orderfourth orderFigure 5.12: A comparison of lower order primal error, error in error usingETE, and higher order primal error for Navier-Stokes using (2, 4, 4; 2, 2), withrefinement in both k and h, k/h ≈ 0.42 fixed.1225.5. Accuracy Comparisons for Error in Functionalswhile adjoint methods for functional correction are constrained to targetfunctionals instead of the solution variables themselves, the current ETE ap-proach is by no means restricted to the examination of functionals, as thesolution variables were examined in the previous section. To proceed, weexamine test cases where a periodic steady state in time exists. To performanalysis on functionals J(t) of the solution at periodic steady state, at timesit will be convenient to consider the signal in the frequency domain. This isdone by computing the Fourier Transform, in the continuous sense, by thepairsJ (ω) = Ĵ(t) =∫ ∞−∞J(t)e−iωtdtJ(t) = J (ω) ∧= 12pi∫ ∞−∞J (ω)eiωtdω (5.29)and discretely asJ [r] = Ĵ [n] =NT−1∑n=0J [n]e−2piirn/NTJ [n] = J [r] ∧= 1NTNT−1∑r=0J [r]e2piirn/NT . (5.30)We can take advantage of some well known properties of this duality inthe analysis of the results. For instance, Parseval’s theorem [82] says thatthe following integral norm can be taken, in either the time or frequencydomain, in the continuous sense as∫ ∞−∞|J(t)|2dt = 12pi∫ ∞−∞|J (ω)|2dω (5.31)and discretely asNT−1∑n=0|J [n]|2 = 1NTNT−1∑r=0|J [r]|2. (5.32)Accuracy of the unsteady functionals will sometimes need to be consideredin the frequency domain. In many instances quantities of interest are relatedto signal frequencies in a natural way. In the vortex shedding case, we willsee that the Strouhal number is directly related to the shedding frequency.Higher order accuracy of unsteady functionals has not been extensivelystudied in the literature. To analyze the results, we summarize several plau-sible derived scalar quantities of the signal that we will use to determine the1235.5. Accuracy Comparisons for Error in Functionalsorder of accuracy of the primal and ETE corrected functionals, from whichthe accuracy of the primal and ETE corrected solutions can be inferred:• Pointwise values.This has been demonstrated for the solution variables in the previoussection, where the solution at the end time is used to monitor theorder of accuracy. However, one disadvantage as we will see for moredifficult test cases is that it is possible for a phase error to arise fromthe discretization. Monitoring pointwise values may be ineffective forcases with varying phases, causing the solution at a fixed time to notnecessarily converge to a value over a particular sequence of meshes.• Maximum amplitude in time.To alleviate the phase problem, we seek other characteristic measuresof the signal that are independent of phase. One such quantity toconsider is taking the maximum amplitude in the time domain (orsimilarly, peak-to-peak values).• Energy-like integrated quantities.Integrating the square of the modulus of the signal gives physically ameasure of the energy, or power when divided by the signal length.The accuracy of this measure will be examined in the results. Othersimilar measures such as the root-mean-square (RMS) in the time do-main, or the power spectral density (PSD) in the frequency domaincan also be used, and some of these are equivalent, at least in thecontinuous domain, through Parseval’s theorem.• Quantities derived from a one-frequency least-squares fit.As mentioned, the single frequency is related to some output quan-tities in a natural way for some test problems. The amplitude andfrequency can be obtained by fitting the data to a sinusoidal ansatzwith a single frequency, asJ(t) = a0 + a1 cos(ωt) + b1 sin(ωt) = a0 +A cos(ωt− ϕ). (5.33)The amplitude A and phase ϕ can be determined from these coeffi-cients as A =√a21 + b21, ϕ = tan−1(b1a1).1245.5. Accuracy Comparisons for Error in FunctionalsThese procedures to obtain a proxy to monitor the asymptotic accuracy areplausible, and in the following test cases we will consider several of thesemeasures. The same model PDEs with slightly different inputs and parame-ters are used to illustrate the functional accuracy.5.5.1 DiffusionThe diffusion problem is posed here as∂tu−∆u = s, on Ω× ΩT = {(x, y, t) ∈ [0, 1]2 × [0, 10]} (5.34)with Dirichlet boundary conditions u|∂Ω = 0. The time-dependent sourceterm iss = sin(pix) sin(piy)(ν cos(νt) + 2pi2 sin(νt)), ν = 2pi (5.35)and the initial condition is u(x, y, 0) = 0. The exact solution is thenu(x, y, t) = sin(pix) sin(piy) sin(νt). (5.36)The time step and length scale are reduced by approximately a factor of two,with ratio fixed at k/h ≈ 4.2 over a sequence of three unnested refinements.The solution, error estimate, and error in error for the medium refinementat the end time are shown in Figure 5.13. We will examine the convergenceof two quantitiesJ1(t) =∫Ωudx =4pi2sin(νt), (5.37)a linear functional, and the 2-norm functionalJ2(t) =(1|Ω|∫Ωu2dx)1/2=12| sin(νt)| (5.38)which is nonlinear. These functionals are shown in Figure 5.14, with themaximum detectable frequency related to the number of samples by theNyquist-Shannon sampling theorem [98]. The results demonstrate the sameconclusions for both smooth functionals.In the frequency domain the solution is very sharp about the forcing fre-quency. This smoothing of high frequencies can be confirmed by a classicalvon Neumann analysis [28, 30] on the diffusion problem. For a sequence ofthree spatial and temporal refinements, Figure 5.15 shows the convergence1255.5. Accuracy Comparisons for Error in Functionals0 101 -0.00059 -0.00041 -0.00023 -5E-05(a) solution u(x, y, Tf )0 101 1E-05 0.00018 0.00035 0.00052(b) error estimate e(x, y, Tf )0 101 -2.05E-05 -1.35E-05 -6.5E-06(c) error in error estimateFigure 5.13: Solution and errors at the final time for the diffusion test case.of the functionals using only the second order primal solution, and usingthe ETE corrected solution.To determine the order of accuracy in the absence of the exact value, weuse Richardson extrapolation techniques [87, 88]. The order of accuracyfrom solutions on three meshes are estimated in the same way as before, perEquation 4.39. The results for measures of the functional are summarizedin Table 5.2. The integrals of the functional data in time are performedusing composite Simpson’s rule for quadrature, which is sufficient to detectfourth order accuracy. Indeed the expected results are observed, where thesecond order primal discretization gives second order accurate functionals,and using the ETE corrected solutions gives fourth order accuracy.5.5.2 AdvectionNext, we consider the linear advection problem∂tu+ c · ∇u = s, on Ω× ΩT = {(x, y, t) ∈ ×[0, 1]2 × [0, 10]}(5.39)with s = e−b((x−x0)2+(y−y0)2) sin(νt), b = 100, x0 = y0 = 0.5, ν = 10, and c =(1 0)T . The problem has initial conditions u(x, y, 0) = 0 , wall boundaryconditions u(x, 0, t) = u(x, 1, t) = 0, inflow at u(0, y, t), and outflow atu(1, y, t). We examine the performance of the ETE method on these twoproblems using a sequence of three meshes, again fixing k/h ≈ 4.2, with theresults for the medium refinement shown in Figure 5.16, where the primalsolution is plotted along with the error estimate.1265.5. Accuracy Comparisons for Error in Functionals0 1 2 3 4 5 6 7 8 9 10−0.6−0.4−0.200.20.40.6t J1 primalJ1 correctedJ2 primalJ2 corrected(a) time domain0 2 4 6 8 10 12 14 16 18 2000.050.10.150.20.250.30.350.40.45frequency (Hz)single−sided amplitude spectrum J1 primalJ1 correctedJ2 primalJ2 corrected(b) frequency domainFigure 5.14: Functionals J1, J2 for the diffusion test case.1275.5. Accuracy Comparisons for Error in Functionals10−1.9 10−1.8 10−1.7 10−1.6 10−1.5 10−1.40.4050.40510.40520.4053k 0.0820.08205 0.08210.08215max J(t) using primalmax J(t) using corrected(1/|ΩT|)∫ΩT|J(t)|2dt using primal(1/|ΩT|)∫ΩT|J(t)|2dt using correctedFigure 5.15: Accuracy of measures of primal and corrected diffusion func-tional.(h, k) (2, 0, 0; 2, 0, 0) (2, 4, 4; 2, 2, 4) corrected(0.01180, 0.05) 0.4050567 0.4052989(0.005906, 0.025) 0.4052256 0.4052856(0.002956, 0.0125) 0.4052704 0.4052848Order 1.91 4.00(a) max J(t)(h, k) (2, 0, 0; 2, 0, 0) (2, 4, 4; 2, 2, 4) corrected(0.01180, 0.05) 0.0820359 0.0821335(0.005906, 0.025) 0.0821039 0.0821282(0.002956, 0.0125) 0.0821221 0.0821279Order 1.90 3.98(b) 1|ΩT |∫ΩT|J(t)|2dtTable 5.2: Orders of accuracy for measures of functionals in the diffusiontest case, using primal solution only and ETE corrected solution.1285.5. Accuracy Comparisons for Error in Functionals0 101 -0.13 -0.085 -0.04 0.005 0.05 0.095(a) solution u(x, y, Tf )0 101 -0.0019 -0.0007 0.0005 0.0017(b) error estimate e(x, y, Tf ).Figure 5.16: Solution and errors at the final time for the advection test case.We investigate the effect of different initial conditions for the ETE, com-paring the usual e(x, y, 0) = 0 with the modified e(x, y, 0) = 0.001ex, whichis a prototype for the following test cases where the ETE initial conditionmay not be known but is taken as zero arbitrarily. The magnitude is smallcompared to the solution, but any adverse effects of this magnitude will beobservable since the error in error is asymptotically smaller as the spatialand time scale is reduced. The functional examined for the order of accu-racy calculations isJ(t) = (u, g) =∫Ωugdx, g = 2pi2 sin(pix) sin(piy). (5.40)Figure 5.17 shows this functional in the time and frequency domains com-puted from the primal and ETE corrected solutions. The influence of theperturbed initial condition diminishes as t increases. The convergence ofthe measures of the functional are shown in Figure 5.18. In Table 5.3, it canbe seen that the order of accuracy of the measures of the periodic functionalis once again as expected, even when the initial conditions for the ETE areperturbed.1295.5. Accuracy Comparisons for Error in Functionals0 1 2 3 4 5 6 7 8 9 10−0.1−0.0500.050.10.15t J primalJ correctedJ corrected with perturbed ETE IC(a) time domain0 2 4 6 8 10 12 14 16 18 2000.010.020.030.040.050.060.070.080.09frequency (Hz)single−sided amplitude spectrum at periodic steady state J primalJ correctedJ corrected with perturbed ETE IC(b) frequency domain.Figure 5.17: Functional for the advection test case using perturbed initialconditions for the ETE.1305.5. Accuracy Comparisons for Error in Functionals10−1.9 10−1.8 10−1.7 10−1.6 10−1.5 10−1.40.0750.080.085k 0.0030.00320.0034max J(t) using primalmax J(t) using corrected(1/|ΩT|)∫ΩT|J(t)|2dt using primal(1/|ΩT|)∫ΩT|J(t)|2dt using correctedFigure 5.18: Accuracy of measures of primal and corrected advection func-tional using perturbed initial conditions for the ETE.(h, k) (2, 0, 0; 2, 0, 0) (2, 4, 4; 2, 2, 4) corrected(0.01181, 0.05) 0.0784193 0.0815988(0.005886, 0.025) 0.0806404 0.0813898(0.002957, 0.0125) 0.0811908 0.0813736Order 2.01 3.69(a) max J(t)(h, k) (2, 0, 0; 2, 0, 0) (2, 4, 4; 2, 2, 4) corrected(0.01181, 0.05) 0.0030662 0.0033219(0.005886, 0.025) 0.0032411 0.0033015(0.002957, 0.0125) 0.0032854 0.0033002Order 1.98 3.95(b) 1|ΩT |∫ΩT|J(t)|2dtTable 5.3: Orders of accuracy for measures of the functional in the advectiontest case using perturbed ETE initial conditions.1315.5. Accuracy Comparisons for Error in FunctionalsName ValueDomain {(x, y, t) ∈ [0, 30]× [0, 40]× [0, 1000]}Center (10,20)Diameter 1M 0.2ReD = ρUD/µ 200Table 5.4: Flow parameters in vortex shedding test case.5.5.3 Navier-Stokes EquationsWe now turn to examining unsteady functionals for the Navier-Stokes equa-tions as a model for viscous compressible flow, which were presented inEquation 4.35. For the ETE, the initial conditions are taken to be zero arbi-trarily, in hopes that higher order accuracy can still be achieved at periodicsteady state, as demonstrated for the simpler model equations.This ETE approach for the Navier-Stokes equations has been applied pre-viously for a manufactured solution test case of a translating vortex. Theoscillatory nature of this problem renders itself to an examination of the pe-riodic steady state, when u(t + $) = u(t) for some period $. A classicalexample commonly used for benchmark testing is the phenomenon of vonKa´rma´n vortex shedding for flow around a circular cylinder of diameter D,a more difficult test case. The selected flow parameters are summarizedin Table 5.4. The left and right boundaries are inflow and outflow, respec-tively. The top and bottom boundaries are imposed as symmetric walls, anda no-slip condition is imposed on the cylinder surface.Over a sequence of three meshes, the time step to length scale ratio isfixed at k/h ≈ 2.0, and the solution and error estimate for the x-momentumvariable for the medium refinement can be seen in Figure 5.19. It was ob-served that if we solve only the primal problem by discretizing in space tofourth order and using the fourth order backward differentiation formula(BDF4), our scheme does not converge, but our ETE approach, which re-quires only second order discretizations in time, did not have difficulties.We use the lift coefficient due to pressure as the functional, defined asCL =1ρ |v∞|2 c/2∣∣∣∣∮ΓP nˆdσ∣∣∣∣ sin θ, (5.41)with c being the reference length, and θ the angle v∞ makes with the resul-1325.5. Accuracy Comparisons for Error in Functionals-0.066 -0.009 0.048 0.105 0.162 0.219(a) x-momentum solution-0.167 -0.106 -0.045 0.016 0.077 0.138(b) x-momentum error estimateFigure 5.19: Solution and errors at the final time for the vortex sheddingtest case.1335.5. Accuracy Comparisons for Error in Functionals(h, k), ω (2, 0, 0; 2, 0, 0) (2, 4, 4; 2, 2, 4) corrected(0.5096, 1) 0.2308491 0.2356796(0.2550, 0.5) 0.2370824 0.2391180(0.1274, 0.25) 0.2386945 0.2388040Order 1.95 3.45Table 5.5: Orders of accuracy for the angular frequency of the functional inthe vortex shedding test case.tant force ∮ΓP nˆdσ. (5.42)The integrals are taken over the cylinder boundary Γ.This functional is plotted in Figure 5.20, where the presence of someovershoots can be seen. This is not unusual in some applications of theCrank-Nicolson method [79], and not necessarily a deficiency of the ETEmethod. This is likely due to boundary condition implementations thatare not totally non-reflecting even in the primal problem, and the Crank-Nicolson scheme only very weakly damps these oscillations, which persisteven after many periods. Using non-reflecting boundary conditions, as in[15, 83, 92, 93], much smoother results can be seen in Figure 5.21.Nevertheless, we can use the original data and infer the fundamentalfrequency from the computed lift coefficient through the fitting method dis-cussed in Section 5.5. These frequency values are summarized in Table 5.5.To proceed further, we can nondimensionalize this quantity to arrive at theStrouhal number, which we use as a proxy to measure higher order accuracy.We make a comparison with the empirical relationshipSt =(ω/2pi)D|v∞| = 0.198(1.0− 19.7ReD)(5.43)obtained in the literature from experimental data [68]. Our results areshown in Figure 5.22. From a validation point of view, it is quite remarkablethat the obtained results are so close to the experimental values, consideringthe fact that experimental error exists in the empirical formula, as well asother forms of error present in our discretization and model. For reference,other numerical computations in the literature are also included in Figure5.22 [22, 34, 48, 70]. The observation here is that using the ETE approach,we can approach our grid-converged Strouhal number at a higher order ratecompared to just using the second order primal solution.1345.5. Accuracy Comparisons for Error in Functionals0 100 200 300 400 500 600 700 800 900 1000−0.8−0.6−0.4−0.200.20.40.60.8tLift coefficient (CL) primalcorrected(a) time domain0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.50.60.7frequency (Hz)single−sided amplitude spectrum at periodic steady state J primalJ corrected(b) frequency domainFigure 5.20: CL for the vortex shedding case.1355.5. Accuracy Comparisons for Error in FunctionalsFigure 5.21: CL for the vortex shedding case using non-reflecting boundaryconditions.10−0.6 10−0.4 10−0.2 1000.1750.180.1850.190.1950.2kStrouhal number using primalusing correctedexperimental modelHarichandan, RoyMeneghini et al.Ding et al.Braza et al.Figure 5.22: Accuracy of Strouhal number computed from the frequency ofthe lift coefficient functional using primal and ETE corrected solutions. Val-ues from the literature for the same ReD = 200 are included for comparison.1365.6. Comparisons to the Higher Order Primal Problem5.6 Comparisons to the Higher Order PrimalProblemFor the steady results of the previous chapter, the computational times takenfor a fixed level of error for the corrected solution using linearized (2, 4, 4)ETE versus the (4, 0, 0) primal solution were compared, and it was foundthat the linearized (2, 4, 4) ETE took less aggregate time than the higher or-der primal problem. We attempt to perform a similar comparison in the con-text of unsteady problems. The computational time comparison between the(2, 4, 4; 2, 2, 4) ETE approach versus solving the higher order primal problem(4, 0, 0; 4, 0, 0) for the 2D diffusion problem is shown in Figure 5.23. The re-sults of the current chapter indicate that they exhibit the same higher orderaccuracy. It can be seen that the ETE approach is marginally faster com-pared to the higher order primal problem. The ETE approach gives more0.02 0.04 0.06 0.08 0.110−1100101102103time stepWall time (s) ETE with lower order primalHigher order primalLower order primalFigure 5.23: Comparison of the efficiency of the ETE approach (combinedtime) versus solving only the primal problem. A diffusion test case is usedwith fixed ratio of time step to mesh length scale, and implicit multistage(Runge-Kutta) schemes were used for time integration.information than just a higher order primal solution, however. As shown,the ETE approach computes a higher order accurate error estimate for alower order accurate solution, and correcting this solution using this errorestimate gives a higher order accurate solution. This is useful in some ap-1375.7. Smoothly Varying Time Stepsplications, but we may alternatively wish for an error bound on the lowerorder solution, which is not available in the case of solving only the higherorder primal problem. A comparison between the cost of the ETE approachagainst the cost of the higher order primal and the lower order primal maybe more fair, in which case it can be shown a priori that the ETE approachis less expensive.Analogous to the steady case, the unsteady ETE approach has robustnessadvantages, in addition to efficiency considerations, over the higher orderprimal problem. The solution of a fourth order primal problem requires dis-cretization to fourth order in both space and time. From the point of view ofstability, it can be shown that A-stable multistage methods can have arbitrar-ily high order, whereas this is not the case for multistep methods [18, 32];there exist so-called order barriers beyond which linear multistep methodsthat satisfy some notion of stability cannot be constructed. Therefore, it isdesirable to keep time discretization orders low in order to avoid difficultieswith stability. For example, in the BDF family of schemes, the unstable re-gion grows into the left half of the complex plane as order increases, as seenin Figure 5.24. Therefore, by using lower order time discretization schemeswe are able to stay further away from the Dahlquist stability barriers. In theexample of the (2, 4, 4; 2, 2) scheme, since the ETE approach only requiresdiscretizing in time to second order, it is more robust than a higher orderprimal solve, which would require discretizing to fourth order in time.5.7 Smoothly Varying Time StepsThe conclusions about mesh smoothness carries over analogously to thetime discretization. For example, one can expect that using a nonsmoothmesh in space and nonsmooth time steps, both the space and time discretiza-tions of the ETE would need to be higher order, or (2, 4, 4; 2, 4, 4), in orderfor a higher order accurate error estimate to be obtained. We have shownthat when uniform time steps are used, some savings can be made in thetime discretization of the ETE, even though the spatial mesh is nonsmooth.A natural question would be whether the uniform time step requirementcan be relaxed to the time steps varying smoothly and still retain the higherorder accuracy in the error estimate, analogous to the case of the general-ization of uniform to smooth meshes in space. Indeed this was observed tobe the case. For example, using uniform meshes, a (2, 2, 4; 2, 2, 4) schemeproduces error estimates that are fourth order accurate in both space andtime. The accuracy of error estimates using uniform meshes in space with1385.7. Smoothly Varying Time Steps−10 0 10 20 30−20−1001020Re(z)Im(z)Stable1 2 3 4 5 6Figure 5.24: Stability region of the BDF family of time discretizationschemes with different orders of accuracy. The exterior of the curves arestable.1395.7. Smoothly Varying Time StepsCorrecteddiscretizationerror10−2 10−110−1510−1010−5100h (2,0,0;2,0,0) Uniform k(2,2,4;2,2,4) Uniform k(2,0,0;2,0,0) Smooth k(2,2,4;2,2,4) Smooth k(2,0,0;2,0,0) Perturbed k(2,2,4;2,2,4) Perturbed ksecond orderfourth orderFigure 5.25: Error accuracy with uniform in space, time step variedsmoothly and randomly.smooth and nonsmooth time step variations for the 1D diffusion problem isshown in Figure 5.25. It can be seen that as expected, in the case where thetime steps are smoothly expanding, the higher order accuracy of the errorestimate is preserved, whereas this is not the case if perturbed time steps areused. This generalization to nonuniform time steps would be quite usefulin practice, since it can serve to provide flexibility and guidance in selectingnonuniform adaptive time steps, for example.140Chapter 6Concluding Remarks6.1 SummaryThe research gap that this work addresses is establishing a framework inobtaining asymptotically accurate and robust discretization error estimatesusing the error transport equation, for general unstructured meshes withoutnecessarily smooth geometric features.The deficiencies in general discretization choices for the primal prob-lem, ETE, and ETE source term were demonstrated in the context of steadyproblems, where the condition of requiring higher order discretization ofthe ETE and evaluation of the ETE source term to the same order, for ahigher order accurate error estimate, was established. The drawback of theadded cost of higher order discretization can be offset by linearizing theETE. For the Navier-Stokes test cases examined, the linearized ETE methodwas able to provide several fold reduction in computation time in the com-bined primal problem with ETE versus the higher order primal only, for afixed error tolerance. This is in addition to robustness advantages versusthe higher order primal, since no root finding is needed, either by pseudo-transient continuation or otherwise; only a linear algebraic system needs tobe solved. This has significant ramifications in practice, both in terms ofobtaining a higher order accurate corrected solution faster than the obviousway of solving the higher order primal problem, but also readily providesa reliable error estimate for the lower order primal solution, which wouldnot be otherwise available. For steady-state problems, the methodology hasbeen demonstrated for models of compressible flow, including the RANS-SAmodel, as well as some simplifications.For unsteady problems, the ETE methodology was applied, with thetime-dependent ETE source term computed using the time derivative of theprimal solution using finite-differences in time, while keeping a lower ordertime discretization for both the primal problem and ETE. This gives higherorder accurate error estimates in both time and space, which is remark-able because this does not require discretizations in both space and time tobe higher order. We also considered some measures of output functionals,1416.2. Future Workwhich are signals in time, to monitor higher order accuracy in the correctedsolution or functionals using the error estimate. As shown, the examinationof functionals is not a restrictive condition; our approach by constructionfirst aims for higher order accuracy in the error estimate of solution vari-ables, and output functionals are computed as postprocessing steps. Time-periodic functionals were analyzed and measures of their accuracy using theETE corrected solution were found to be higher order accurate.The relevance of this work is embodied in the fact that the ETE method-ology can be applied to any primal model equations in an analogous way,and many of the conclusions drawn would apply.6.2 Future WorkHaving demonstrated the research gap in accurate and robust discretiza-tion error estimation using unstructured meshes and how the ETE approachcan address many of these requirements, there are several logical next stepsthat can be investigated further. Much of the proposed further work focusesaround extensions to more complex model equations, in attempt to confirmour claim that the ETE methodology can be applied independent of the cho-sen model equations. In these extensions of the methodology, the generalapproach should be first to establish whether there is a meaningful interpre-tation of higher order accuracy. Once that framework is in place, the ETEapproaches investigated in this thesis can be applied in a similar manner.6.2.1 Extensions to 3DA natural extension is to apply the ETE methodology to 3D problems on un-structured meshes. From the point of view of formulating and solving theETE in tandem with the primal equations, extending the methodology from2D to 3D problems is likely to be conceptually straightforward. To demon-strate, a prototypical example is shown using the 3D Poisson equation,−∆u = s, on Ω = {(x, y, z) ∈ [0, 1]3} (6.1)with Dirichlet boundary conditions u|∂Ω = 0, and source terms = 3pi2 sin(pix) sin(piy) sin(piz), (6.2)so that the exact solution isu = sin(pix) sin(piy) sin(piz). (6.3)1426.2. Future Work00.5100.5100.51 xy z00.20.40.60.81(a) (2,0,0) solution00.5100.5100.51 xy z−8−7−6−5−4−3−2−10x 10−3(b) exact discretization error1436.2. Future Work00.5100.5100.51 xy z−8−7−6−5−4−3−2−10x 10−3(c) (2,2,4) error estimate00.5100.5100.51 xy z0123456x 10−5(d) error in errorFigure 6.1: Results for a 3D Poisson problem using a uniform mesh.1446.2. Future Work10−2 10−110−710−610−510−410−310−210−1hCorrected discretization error (2,0,0)(2,2,4)second orderfourth orderFigure 6.2: Accuracy of the primal solution and error estimate for a 3DPoisson problem, on uniform meshes.The solution, exact discretization error, error estimate, and error in erroris shown in Figure 6.1. Applying the ETE approach with (2,2,4) using auniform mesh, the same conclusion can be obtained for the accuracy of theerror estimate as in 2D, as seen in Figure 6.2. It is expected, then, that the(2,2,4) scheme will not be fourth order accurate if nonsmooth meshes wereused, and a (2,4,4) scheme will be needed.6.2.2 PDEs of Higher OrderIt is important for the method to be able to handle PDEs with higher deriva-tives, which arise quite commonly in practice. For example, the biharmonicequation∆∆u = 0 (6.4)is a fourth order PDE that arises in the modeling of linear elasticity of solidmechanics or Stokes flow. In the finite-volume discretization of PDEs with1456.2. Future Workorder d derivatives, the fluxes have derivatives of order d − 1. Therefore,it is reasonable to require that the discretization be polynomials of degreeat least d − 1, hence being order d accurate. In terms of the ETE, the re-quirements are the same for nonsmooth meshes as we have seen for thesecond derivative, only that the minimal discretization orders for all of p, q,and r would need to be at least d. Further work can confirm this claim, andalso investigate whether accuracy is retained if the equations are recast as asystem of first order equations.6.2.3 Turbulence Modeling ExtensionsThe ETE methodology applied to the Spalart-Allmaras (SA-neg) turbulencemodel has been shown in the results as a proof-of-concept for more difficultmodel equations in CFD. The error estimation approach can be extendedto other approaches to turbulence, such as Large Eddy Simulation (LES)[100], or even direct numerical simulations (DNS) to resolve all scales.Through this process, it is expected that intrinsically chaotic systems willbe encountered, where small perturbations to the data do not necessarilyremain small. Therefore, one direction of further study is to investigate theextent to which the ETE methodology can be applied in these cases. To il-lustrate, the Lorenz equations are a classical system of ODEs used to studysuch chaotic behavior, which are formulated as: xyzt= σ(y − x)x(ρ− z)− yxy − βz (6.5)which are stable [54] if the parameters satisfy the conditionρ < σσ + β + 3σ − β − 1 . (6.6)The values that are used to exhibit chaos are σ = 10, β = 8/3, ρ = 28, andthe values σ = 10, β = 8/3, ρ = 5 are used for stable behavior. The set ofODEs are solved using the forward Euler method, which has nominal firstorder accuracy. The time trajectory of the solution in both cases is shownin Figure 6.3. The discretization error at the final time is estimated usingreference solutions computed with an extremely small time step, and theseare shown in Figure 6.4. It can be seen that while the stable system hasa discretization error that is close to the nominal first order accuracy, nopointwise order of accuracy can be meaningfully discerned for the chaotic1466.2. Future Worksystem. This highlights the difficulties that will be encountered when theETE approach is used in attempt to estimate the discretization error. Theeffectiveness of other approaches such as averaging will remain to be seen.6.2.4 Weak SolutionsThe extension of the unstructured mesh ETE methodology to problems withweak solutions requires careful investigation. Some results for uniformmeshes exist in the literature [10]; it is documented that for solutions withshocks, the fully nonlinear problem needs to be solved for accurate error es-timates for smooth meshes. This is due to the fact that linearization assumesthat the error and its derivatives are small in comparison to the solution, butnear a shock this no longer holds. An example would be using the fully non-linear (2,2,4) scheme to obtain a fourth order accurate error estimate. Whilethis is still economical to solve versus a fully nonlinear fourth order primalproblem, for nonsmooth meshes a fully nonlinear (2,4,4) scheme would beneeded, drastically diminishing any efficiency advantages. Furthermore, thevirtue of higher order accuracy for solutions with shocks may be in questionaltogether, due to the fact that the error terms in the Taylor series expan-sion depend on higher and higher derivatives of the solution, which areunbounded near a discontinuity. The performance of the unstructured meshETE methodology for classical techniques for discontinuous data should alsobe examined. This current work has shown the effectiveness on handlingnonlinear model equations, but performance for nonlinear schemes is notso clear, which typically is the case regardless of the linearity of the modelequations. More sophisticated schemes such as the use of a Godunov-typeflux [45] for entropy satisfying solutions, and the use of limiters [109] arejust a few examples that fall in this category.1476.2. Future Work−20 −10 010 20 30−500500102030405060xyz(a) σ = 10, β = 8/3, ρ = 281 1.5 22.5 3 3.5 44.502460123456xyz(b) σ = 10, β = 8/3, ρ = 5Figure 6.3: Lorenz solution using initial conditions that differ by a 0.1%perturbation. The solid marker indicates the solutions at the final time. Forthe nonchaotic system, the solutions at the final time are close, while forthe chaotic system the solutions at the final time as well as trajectories arenoticeably different.1486.2. Future Work0.0005 0.001 0.001410−610−510−410−310−210−1100101102hDiscretization error σ = 10, ρ = 8/3, β = 5σ = 10, ρ = 8/3, β = 28First orderFigure 6.4: Discretization error for stable and chaotic Lorenz system.149Bibliography[1] Adjerid, S., and Baccouch, M., 2010. “Asymptotically exact a posteri-ori error estimates for a one-dimensional linear hyperbolic problem”.Applied Numerical Mathematics, 60(9), pp. 903–914.[2] Aftosmis, M., Gaitonde, D., and Tavares, T. S., 1995. “Behavior of lin-ear reconstruction techniques on unstructured meshes”. AIAA Jour-nal, 33, pp. 2038–2049.[3] Ainsworth, M., and Oden, J., 1997. “A posteriori error estimation infinite element analysis”. Computer Methods in Applied Mechanics andEngineering, 142(1–2), pp. 1–88.[4] Alcrudo, F., and Garcia-Navarro, P., 1993. “A high-resolutionGodunov-type scheme in finite volumes for the 2D shallow-waterequations”. International Journal for Numerical Methods in Fluids,16(6), pp. 489–505.[5] Allmaras, S. R., Johnson, F. T., and Spalart, P., 2012. “Modificationsand clarifications for the implementation of the Spalart-Allmaras tur-bulence model”. In Proceedings of the Seventh International Confer-ence on Computational Fluid Dynamics. ICCFD7-1902.[6] Alnæs, M. S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A.,Richardson, C., Ring, J., Rognes, M. E., and Wells, G. N., 2015. “TheFEniCS project version 1.5”. Archive of Numerical Software, 3(100).[7] Arnold, D. N., Brezzi, F., Cockburn, B., and Marini, L. D., 2002.“Unified analysis of discontinuous Galerkin methods for elliptic prob-lems”. SIAM Journal on Numerical Analysis, 39(5), pp. 1749–1779.[8] Babuska, I., and Strouboulis, T., 2001. The Reliability of the FEMethod. Oxford Press.[9] Baker, A. H., Falgout, R. D., Kolev, T. V., and Yang, U. M., 2012. “Scal-ing hypre’s multigrid solvers to 100,000 cores”. In High-PerformanceScientific Computing. Springer, pp. 261–279.150Bibliography[10] Banks, J., Hittinger, J., Connors, J., and Woodward, C., 2012. “Nu-merical error estimation for nonlinear hyperbolic PDEs via nonlinearerror transport”. Computer Methods in Applied Mechanics and Engi-neering, 213-216, pp. 1–15.[11] Barth, T. J., and Frederickson, P. O., 1990. “Higher order solution ofthe Euler equations on unstructured grids using quadratic reconstruc-tion”. In Proceedings of the Twenty-Eighth AIAA Aerospace SciencesMeeting. AIAA Paper 90-0013.[12] Barth, T. J., and Larson, M. G., 2002. A posteriori error estimatesfor higher order Godunov finite volume methods on unstructuredmeshes. Tech. rep., Complex Applications III, R. Herbin and D. Kro-ner (Eds), HERMES Science Publishing Ltd.[13] Batchelor, G., 2000. An Introduction to Fluid Dynamics. CambridgeMathematical Library. Cambridge University Press.[14] Batill, S. M., 1994. Experimental uncertainty and drag measurementsin the national transonic facility. Technical Report 4600, NASA.[15] Bayliss, A., and Turkel, E., 1982. “Far field boundary conditionsfor compressible flows”. Journal of Computational Physics, 48(2),pp. 182–199.[16] Bermejo-Moreno, I., Bodart, J., Larsson, J., Barney, B. M., Nichols,J. W., and Jones, S., 2013. “Solving the compressible Navier-Stokes equations on up to 1.97 million cores and 4.1 trillion gridpoints”. In Proceedings of the International Conference on HighPerformance Computing, Networking, Storage and Analysis, ACM,pp. 62:1–62:10.[17] Biedron, R. T., Carlson, J.-R., Derlaga, J. M., Gnoffo, P. A., Ham-mond, D. P., Jones, W. T., Kleb, B., Lee-Rausch, E. M., Nielsen, E. J.,Park, M. A., Rumsey, C. L., Thomas, J. L., and Wood, W. A., 2017.FUN3D manual: 13.1. NASA TM-2016-219580.[18] Bijl, H., Carpenter, M. H., and Vatsa, V. N., 2001. “Time integrationschemes for the unsteady Navier-Stokes equations”. In Proceedings ofthe Fifteenth AIAA Computational Fluid Dynamics Conference. AIAAPaper 2001-2612.[19] Black, F., and Scholes, M., 1973. “The pricing of options and corpo-rate liabilities”. Journal of Political Economy, 81(3), pp. 637–654.151Bibliography[20] Blinchikoff, H., and Krause, H., 2001. Filtering in the time and fre-quency domains. The Institution of Engineering and Technology.[21] Brandt, A., 1977. “Multi-level adaptive solutions to boundary-valueproblems”. Mathematics of Computation, 31(138), pp. 333–390.[22] Braza, M., Chassaing, P., and Ha Minh, H., 1986. “Numerical studyand physical analysis of the pressure and velocity fields in the nearwake of a circular cylinder”. Journal of Fluids Mechanics, 163, pp. 79–130.[23] Burgess, N. K., and Mavriplis, D. J., 2012. “High-order discontinuousGalerkin methods for turbulent high-lift flows”. In Seventh Interna-tional Conference on Computational Fluid Dynamics.[24] Celik, I., Ghia, U., Roache, P. J., et al., 2008. “Procedure for es-timation and reporting of uncertainty due to discretization in CFDapplications”. Journal of Fluids Engineering-Transactions of the ASME,130(7).[25] Celik, I., and Hu, G., 2004. “Single grid error estimation using errortransport equation”. Journal of Fluids Engineering, 126(5), pp. 778–790.[26] Celik, I., and Karatekin, O., 1997. “Numerical experiments on appli-cation of Richardson extrapolation with nonuniform grids”. Journalof Fluids Engineering, 119(3), pp. 584–590.[27] Ceze, M., and Fidkowski, K. J., 2013. “Pseudo-transient continuation,solution update methods, and CFL strategies for DG discretizationsof the RANS-SA equations”. In Proceedings of the Twenty-First AIAAComputational Fluid Dynamics Conference. AIAA Paper 2013-2686.[28] Charney, J. G., Fjörtoft, R., and von Neumann, J., 1950. “Numer-ical integration of the barotropic vorticity equation”. Tellus, 2(4),pp. 237–254.[29] Christlieb, A. J., Macdonald, C. B., and Ong, B. W., 2010. “Parallelhigh-order integrators”. SIAM Journal on Scientific Computing, 32(2),pp. 818–835.[30] Crank, J., and Nicolson, P., 1947. “A practical method for numer-ical evaluation of solutions of partial differential equations of the152Bibliographyheat-conduction type”. Mathematical Proceedings of the CambridgePhilosophical Society, 43(1), pp. 50–67.[31] Cuthill, E. H., and McKee, J. M., 1969. “Reducing the bandwidthof sparse symmetric matrices”. In Proceedings of the 24th NationalConference of the Association for Computing Machinery, pp. 157–172.[32] Dahlquist, G. G., 1963. “A special stability problem for linear multi-step methods”. BIT Numerical Mathematics, 3(1), pp. 27–43.[33] Davisson, C., and Germer, L., 1928. “Reflection of electrons by acrystal of nickel”. Proceedings of the National Academy of Sciences,14(4), pp. 317–322.[34] Ding, H., Shu, C., Yeo, K., and Xu, D., 2007. “Numerical simulation offlows around two circular cylinders by mesh-free least square-basedfinite difference methods”. International Journal for Numerical Meth-ods in Fluids, 53(2), pp. 305–332.[35] Diskin, B., Thomas, J. L., Rumsey, C. L., and Schwöppe, A., 2016.“Grid-convergence of Reynolds-averaged Navier–Stokes solutions forbenchmark flows in two dimensions”. AIAA Journal.[36] Dongarra, J. J., Hempel, R., Hey, A. J., and Walker, D. W., 1993. Aproposal for a user-level, message passing interface in a distributedmemory environment. Tech. rep., Oak Ridge National Laboratory.[37] Dutt, A., Greengard, L., and Rokhlin, V., 2000. “Spectral deferred cor-rection methods for ordinary differential equations”. BIT NumericalMathematics, 40(2), pp. 241–266.[38] Ferziger, J. H., 1989. Estimation and reduction of numerical error.ASME Winter Annual Meeting, Dec.[39] Ferziger, J. H., and Peric, M., 1996. “Further discussion of numericalerrors in CFD”. International Journal for Numerical Methods in Fluids,23(12), pp. 1263–1274.[40] Fortune, P., 1996. “Anomalies in option pricing: The Black-Scholesmodel revisited”. New England Economic Review, pp. 17–41.[41] Fraysse, F., de Vicente, J., and Valero, E., 2012. “The estimation oftruncation error by τ -estimation revisited”. Journal of ComputationalPhysics, 231(9), pp. 3457–3482.153Bibliography[42] George, A., and Liu, J. W.-H., 1981. Computer Solution of Large SparsePositive Definite Systems. Prentice-Hall, Englewood Cliffs, N.J.[43] George, J. A., 1971. Computer implementation of the finite elementmethod. Tech. Rep. STAN–CS–71–208, Computer Science Depart-ment, Stanford University.[44] Giles, M. B., and Pierce, N. A., 2000. “An introduction to the adjointapproach to design”. Flow, Turbulence, and Combustion, 65(3–4),pp. 393–415.[45] Godunov, S. K., 1959. “A finite difference method for the numeri-cal computation of discontinuous solutions of the equations of fluiddynamics”. Matematicheskie Sborniki, 47(8–9), pp. 357–393.[46] Goldberg, D., 1991. “What every computer scientist should knowabout floating-point arithmetic”. ACM Computing Surveys (CSUR),23(1), pp. 5–48.[47] Griewank, A., and Walther, A., 2000. “Algorithm 799: revolve: animplementation of checkpointing for the reverse or adjoint mode ofcomputational differentiation”. ACM Transactions on MathematicalSoftware (TOMS), 26(1), pp. 19–45.[48] Harichandan, A., and Roy, A., 2010. “Numerical investigation of lowReynolds number flow past two and three circular cylinders usingunstructured grid CFR scheme”. International Journal of Heat andFluid Flow, 31(2), pp. 154–171.[49] Hay, A., and Visonneau, M., 2004. “CFD uncertainty analysis forturbulent flows using error equation method”. In Workshop on CFDUncertainty Analysis Lisbon.[50] Hay, A., and Visonneau, M., 2006. “Error estimation using the errortransport equation for finite-volume methods and arbitrary meshes”.International Journal Computational Fluid Dynamics, 20(7), pp. 463–479.[51] Herrmann, L. R., 1976. “Laplacian-isoparametric grid generationscheme”. Journal of the Engineering Mechanics Division, 102(5),pp. 749–907.154Bibliography[52] Hestenes, M. R., and Stiefel, E., 1952. “Method of conjugate gradi-ents for solving linear systems”. Journal of Research of the NationalBureau of Standards, 49(6), December, pp. 409–436.[53] Hirsch, C., 1988. Numerical Computation of Internal and ExternalFlows. Wiley.[54] Hirsch, M. W., Smale, S., and Devaney, R., 2003. Differential Equa-tion, Dynamical Systems, & An Introduction to Chaos, second ed. Aca-demic Press.[55] Hoshyari, S., 2017. “A higher-order unstructured finite volume solverfor three-dimensional compressible flows”. Master’s thesis, The Uni-versity of British Columbia.[56] Huynh, H., 2013. “High-order space-time methods for conservationlaws”. In Proceedings of the Twenty-First AIAA Computational FluidDynamics Conference. AIAA Paper 87-1184.[57] Ilinca, C., Zhang, X., Trépanier, J.-Y., and Camarero, R., 2000. “Acomparison of three error estimation techniques for finite-volume so-lutions of compressible flows”. Computer Methods in Applied Mechan-ics and Engineering, 189(4), pp. 1277–1294.[58] Illingworth, C., 1950. “Some solutions of the equations of flow of aviscous compressible fluid”. Proceedings of the Cambridge Philosophi-cal Society, 46(3), pp. 469–478.[59] Jacobs, E. N., Ward, K. E., and Pinkerton, R. M., 1933. The character-istics of 78 related airfoil sections from tests in the variable-densitywind tunnel. Tech. Rep. No. 460, National Advisory Committee forAeronautics.[60] Jalali, A., 2017. “An adaptive higher-order unstructured finite volumesolver for turbulent compressible flows”. PhD thesis, The Universityof British Columbia.[61] Jameson, A., 1988. “Aerodynamic design via control theory”. Journalof Scientific Computing, 3(3), pp. 233–260.[62] Kelley, C. T., and Keyes, D. E., 1998. “Convergence analysis ofpseudo-transient continuation”. SIAM Journal on Numerical Analy-sis, 35(2), pp. 508–523.155Bibliography[63] Knuth, D. E., 1997. The Art of Computer Programming, Volume 1 :Fundamental Algorithms, third ed. Addison Wesley Longman Pub-lishing Co., Inc., Redwood City, CA, USA.[64] Krist, S. L., Biedron, R. T., and Rumsey, C. L., 1998. CFL3D user’smanual (version 5.0). NASA TM-1998-208444.[65] LeVeque, R., 1992. Numerical Methods for Conservation Laws. Lec-tures in Mathematics ETH Zürich, Department of Mathematics Re-search Institute of Mathematics. Springer.[66] LeVeque, R., 2007. Finite Difference Methods for Ordinary and Par-tial Differential Equations: Steady-State and Time-Dependent Prob-lems. Other Titles in Applied Mathematics. Society for Industrial andApplied Mathematics.[67] Masatsuka, K., 2013. I do like CFD, Vol.1, second ed. Katate Masat-suka.[68] Massey, B., and Ward-Smith, J., 2006. Mechanics of Fluids. Taylorand Francis Group.[69] Mehta, U. B., 1996. “Guide to credible computer simulations of fluidflows”. Journal of Propulsion and Power, 12(5), pp. 940–948.[70] Meneghini, J., Saltara, F., Siqueira, C., and Ferrari Jr, J., 2001. “Nu-merical simulation of flow interference between two circular cylin-ders in tandem and side-by-side arrangement”. Journal of Fluids andStructures, 15(2), pp. 327–350.[71] Moore, G. E., 1965. “Cramming more components onto integratedcircuits”. Electronics, 38(8).[72] Nguyen, N. C., Persson, P.-O., and Peraire, J., 2007. “RANS solutionsusing high order discontinuous Galerkin methods”. In Proceedings ofthe Forty-Fifth AIAA Aerospace Sciences Meeting. AIAA Paper 2007-914.[73] Nishikawa, H., 2010. “Beyond interface gradient: A general principlefor constructing diffusion schemes”. In Proceedings of the FortiethAIAA Fluid Dynamics Conference and Exhibit. AIAA paper 2010-5093.156Bibliography[74] Nogueira, X., Cueto-Felgueroso, L., Colominas, I., Gomez, H., Navar-rina, F., and Casteleiro, M., 2009. “On the accuracy of finite volumeand discontinuous Galerkin discretizations for compressible flow onunstructured grids”. International Journal for Numerical Methods inEngineering, 78(13), pp. 1553–1584.[75] Oberkampf, W. L., and Roy, C. J., 2010. Verification and validation inscientific computing. Cambridge University Press.[76] Oberkampf, W. L., Trucano, T. G., and Hirsch, C., 2010. Verifica-tion, validation, and predictive capability in computational engineer-ing and physics. Tech. rep., Sandia National Laboratories. SandiaReport SAND 2003-3769.[77] Ollivier-Gooch, C. F., 1998–2017. GRUMMP — Generation and Re-finement of Unstructured, Mixed-element Meshes in Parallel. http://tetra.mech.ubc.ca/GRUMMP.[78] Ollivier-Gooch, C. F., and Van Altena, M., 2002. “A high-order ac-curate unstructured mesh finite-volume scheme for the advection-diffusion equation”. Journal of Computational Physics, 181(2),pp. 729–752.[79] Østerby, O., 2003. “Five ways of reducing the Crank-Nicolson oscilla-tions”. BIT Numerical Mathematics, 43(4), pp. 811–822.[80] Owen, S., 1998. “A survey of unstructured mesh generation technol-ogy”. In Proceedings of the 7th International Meshing Roundtable,pp. 239–267.[81] Paige, C. C., and Saunders, M. A., 1975. “Solution of sparse indefinitesystems of linear equations”. SIAM journal on numerical analysis,12(4), pp. 617–629.[82] Parseval des Chênes, M.-A., 1806. “Mémoire sur les séries etsur l’intégration complète d’une équation aux différences partielleslinéaires du second ordre, à coefficients constants”. Mémoires présen-tés par divers savants, Académie des Sciences, Paris, 1, pp. 638–648.[83] Poinsot, T., and Lele, S., 1992. “Boundary conditions for direct sim-ulations of compressible viscous flows”. Journal of ComputationalPhysics, 101(1), pp. 104–129.157Bibliography[84] Qin, Y., and Shih, T. I.-P., 2003. “A method for estimating grid-induced errors in finite-difference and finite-volume methods”. InProceedings of the Forty-First AIAA Aerospace Sciences Meeting.AIAA Paper 2003-0845.[85] Richardson, L. F., 1911. “The approximate arithmetical solution by fi-nite differences of physical problems involving differential equations,with an application to the stresses in a masonry dam”. PhilosophicalTransactions of the Royal Society of London A: Mathematical, Physicaland Engineering Sciences, 210(459-470), pp. 307–357.[86] Richardson, L. F., 1922. “Weather prediction by numerical pro-cess”. Quarterly Journal of the Royal Meteorological Society, 48(203),pp. 282–284.[87] Richardson, L. F., and Gaunt, J. A., 1927. “The deferred ap-proach to the limit. Part I. Single lattice. Part II. Interpenetratinglattices”. Philosophical Transactions of the Royal Society London A,226, pp. 299–361.[88] Roache, P. J., 1994. “A method for uniform reporting of grid refine-ment studies”. Journal of Fluids Engineering, 116(3), pp. 405–413.[89] Roe, P. L., 1981. “Approximate Riemann solvers, parameter vectors,and difference schemes”. Journal of Computational Physics, 43(2),pp. 357–372.[90] Romberg, W., 1955. “Vereinfachte numerische integration”. Det Kon-gelige Norske Videnskabers Selskab Forhandlinger, Trondheim, 28(7),pp. 30–36.[91] Roy, C. J., 2010. “Review of discretization error estimators in scien-tific computing”. In Proceedings of the Forty-Eighth AIAA AerospaceSciences Meeting. AIAA Paper 2010-126.[92] Rudy, D. H., and Strikwerda, J. C., 1980. “A nonreflecting outflowboundary condition for subsonic Navier-Stokes calculations”. Journalof Computational Physics, 36(1), pp. 55–70.[93] Rudy, D. H., and Strikwerda, J. C., 1981. “Boundary conditions forsubsonic compressible Navier-Stokes calculations”. Computers & Flu-ids, 9(3), pp. 327–338.158Bibliography[94] Rumsey, C., 2014. Turbulence modeling resource, NASA Langley Re-search Center. http://turbmodels.larc.nasa.gov.[95] Saad, Y., and Schultz, M. H., 1986. “GMRES: A generalized minimalresidual algorithm for solving nonsymmetric linear systems”. SIAMJournal of Scientific and Statistical Computing, 7(3), July, pp. 856–869.[96] Salari, K., and Knupp, P., 2000. Code verification by the method ofmanufactured solutions. Tech. rep., Sandia National Laboratories.Sandia Report SAND 2000-1444.[97] Schrödinger, E., 1926. “An undulatory theory of the mechanics ofatoms and molecules”. Physical Review, 28(6), p. 1049.[98] Shannon, C. E., 1949. “Communication in the presence of noise”.Proceedings of the IRE, 37(1), pp. 10–21.[99] Slotnick, J., Khodadoust, A., Alonso, J., Darmofal, D., Gropp, W.,Lurie, E., and Mavriplis, D., 2014. CFD vision 2030 study: a path torevolutionary computational aerosciences. Technical Report 218178,NASA.[100] Smagorinsky, J., 1963. “General circulation experiments with theprimitive equations: I. the basic experiment”. Monthly weather re-view, 91(3), pp. 99–164.[101] Smith, B., Balay, S., Brown, J., Dalcin, L., Karpeev, D., Knepley, M.,and Zhang, H., 2011. PETSc home page. http://www.mcs.anl.gov/petsc.[102] Spalart, P. R., 2000. “Strategies for turbulence modelling and simula-tions”. International Journal of Heat and Fluid Flow, 21(3), pp. 252–263.[103] Spalart, P. R., Allmaras, S. R., et al., 1994. “A one equation turbulencemodel for aerodynamic flows”. La Recherche Aérospatiale, pp. 5–21.[104] Stetter, H., 1978. “The defect correction principle and discretizationmethods”. Numerische Mathematik, 29(4), pp. 425–443.[105] Stumm, P., and Walther, A., 2008. Towards the economical computa-tion of adjoints in PDEs using optimal online checkpointing. Techni-cal Report DFG-SPP 1253-15-04, Deutsche ForschungsgemeinschaftSchwerpunktprogramm.159Bibliography[106] Succi, S., 2001. The lattice Boltzmann equation: for fluid dynamicsand beyond. Oxford university press.[107] Tyson, W. C., S´wirydowicz, K. K., Derlaga, J. M., Roy, C. J., andde Sturler, E., 2016. “Improved functional-based error estimationand adaptation without adjoints”. In Proceedings of the Fifty-FourthAIAA Aerospace Sciences Meeting. AIAA Paper 2016-3809.[108] van der Vorst, H. A., 1992. “Bi-CGStab — a fast and smoothly con-verging variant of Bi-CG for the solution of nonsymmetric linear sys-tems”. SIAM Journal of Scientific and Statistical Computing, 13(2),Mar., pp. 631–644.[109] van Leer, B., 1979. “Towards the ultimate conservative differencescheme. V. A second-order sequel to Godunov’s method”. Journal ofComputational Physics, 32(1), pp. 101–136.[110] Vassberg, J. C., DeHaan, M. A., and Sclafani, T. J., 2003. “Grid gener-ation requirements for accurate drag predictions based on overflowcalculations”. In Proceedings of the Sixteenth AIAA ComputationalFluid Dynamics Conference. AIAA Paper 2003-4124.[111] Venditti, D. A., and Darmofal, D. L., 2000. “Adjoint error estimationand grid adaptation for functional outputs: Application to quasi-one-dimensional flow”. Journal of Computational Physics, 164(1), Oct.,pp. 204–227.[112] Veress, Á., and Rohács, J., 2012. Application of Finite Volume Methodin Fluid Dynamics and Inverse Design Based Optimization, Finite Vol-ume Method - Powerful Means of Engineering Design. InTech.[113] Wang, Z. J., Fidkowski, K., Abgrall, R., Bassi, F., Caraeni, D., Cary,A., Deconinck, H., Hartmann, R., Hillewaert, K., Huynh, H., Kroll,N., May, G., Persson, P.-O., van Leer, B., and Visbal, M., 2013. “High-order CFD methods: Current status and perspective”. InternationalJournal for Numerical Methods in Engineering, 72(8), pp. 811–845.[114] Warming, R., and Hyett, B., 1974. “The modified equation approachto the stability and accuracy analysis of finite-difference methods”.Journal of Computational Physics, 14(2), pp. 159–179.[115] Yan, G., and Ollivier-Gooch, C., 2015. “Accuracy of discretizationerror estimation by the error transport equation on unstructured160meshes”. In Proceedings of the Fifty-Third AIAA Aerospace SciencesMeeting. AIAA Paper 2015-1264.[116] Yan, G., and Ollivier-Gooch, C., 2015. “Accuracy of discretiza-tion error estimation by the error transport equation on unstruc-tured meshes - nonlinear systems of equations”. In Proceedings ofthe Twenty-Second AIAA Computational Fluid Dynamics Conference.AIAA Paper 2015-2747.[117] Yan, G., and Ollivier-Gooch, C., 2016. “Discretization error estima-tion by the error transport equation on unstructured meshes - appli-cations to viscous flows”. In Proceedings of the Fifty-Fourth AIAAAerospace Sciences Meeting. AIAA Paper 2016-0836.[118] Yan, G., and Ollivier-Gooch, C., 2016. “Towards higher order dis-cretization error estimation by error transport using unstructuredfinite-volume methods for unsteady problems”. In Proceedings of theNinth International Conference on Computational Fluid Dynamics.ICCFD9-174.[119] Yan, G., and Ollivier-Gooch, C., 2017. “Applications of the unsteadyerror transport equation on unstructured meshes”. In Proceedings ofthe Twenty-Third AIAA Computational Fluid Dynamics Conference.AIAA Paper 2017-4112.[120] Yan, G., and Ollivier-Gooch, C., 2017. “On efficiently obtaininghigher order accurate discretization error estimates for unstructuredfinite-volume methods using the error transport equation”. Journalof Verification, Validation, and Uncertainty Quantification, 2(4).[121] Yan, G., and Ollivier-Gooch, C., 2017. “Towards higher order dis-cretization error estimation by error transport using unstructuredfinite-volume methods for unsteady problems”. Computers and Flu-ids, 154, pp. 245–255.[122] Zhang, X., Trépanier, J.-Y., and Camarero, R., 2000. “A posteriorierror estimation for finite-volume solutions of hyperbolic conserva-tion laws”. Computer Methods in Applied Mechanics and Engineering,185(1), pp. 1–19.161Appendix ARichardson ExtrapolationThe details of the Richardson extrapolation procedure is presented for ascalar exact solution u. Using three successive refinements, the solutionsU1(h), U1(h2 ), U1(h4 ) are obtained, which satisfyu = U1(h) + c0hp + c1hp+1 + c2hp+2 + c3hp+3 +O(hp+4)u = U1(h2)+ c0hp2p+ c1hp+12p+1+ c2hp+22p+2+ c3hp+32p+3+O(hp+4)u = U1(h4)+ c0hp4p+ c1hp+14p+1+ c2hp+24p+2+ c3hp+34p+3+O(hp+4) (A.1)We define an extrapolated solution using the two coarsest solutions, asU2(h) :=2p2p − 1U1(h2)− 12p − 1U1(h), (A.2)which can be verified to satisfyu = U2(h) + c˜0hp+1 +O(hp+2), (A.3)an asymptotically more accurate solution.A.1 Known pIf the order of accuracy p is known, then the error on the coarsest mesh canbe approximated asc0hp ≈ 2p2p − 1(U1(h2)− U1(h)). (A.4)Furthermore, with three solutions we can perform elimination using thesolutions on the second and third meshes to getU2(h2):=2p2p − 1U1(h4)− 12p − 1U1(h2)(A.5)162A.2. Unknown pwhich satisfiesu = U2(h2)+ c˜0hp+12p+1+O(hp+2). (A.6)Using the two extrapolated solutions, we can perform extrapolation oncemore to getU3(h) :=2p+12p+1 − 1U2(h2)− 12p+1 − 1U2(h) (A.7)so thatu = U3(h) + cˆ1hp+2 +O(hp+3). (A.8)A.2 Unknown pIf the order of accuracy is unknown, at least three mesh resolutions areneeded. Subtracting consecutive pairs of solutions, we haveU1(h2)− U1(h) ≈ c0hp 2p − 12pU1(h4)− U1(h2)≈ c0hp2p2p − 12p. (A.9)Hence p can be computed asp ≈ log[(U1(h2)− U1(h)) / (U1 (h4 )− U1 (h2 ))]log 2(A.10)and the error on the coarsest mesh can be computed.163Appendix BCFD Model EquationsThe simplified model equations of advection, diffusion, Burgers, and Eu-ler tested in this thesis are useful for testing numerical schemes, but theycan also be considered as simplifications of the compressible Navier-Stokesequations under various assumptions. We start from the compressible Navier-Stokes equations∂t ρρvE+∇ · ρvρvv + P1− τv(E + P )− τ · v + q = 0. (B.1)If the flow is incompressible, and pressure contributions are negligible, weget∇ · v = 0 (B.2)∂v∂t+∇ · vv = ν∇2v. (B.3)Combining these gives the viscous Burgers’ equation∂v∂t+ v · ∇v = ν∇2v. (B.4)Under the incompressible flow assumption, the energy equation can be de-coupled from the mass and momentum equations, as∂E∂t+ v · ∇E = ∇ · (τ · v)−∇ · q. (B.5)This can be written in terms of the temperature using the thermodynamicrelations E = ρcpT, and q = −K∇T, giving∂T∂t+ v · ∇T = 1ρcp∇ · (τ · v) + kρcp∇2T. (B.6)A further assumption of zero velocity v = 0 gives the diffusion equation∂T∂t= κ∇2T, (B.7)164Appendix B. CFD Model Equationswhich is the diffusion model used in this thesis. Separately, if we assume ofzero viscosity and thermal conductivity in Equation B.6, this would give theadvection equation∂T∂t+ v · ∇T = 0. (B.8)Starting again from the compressible Navier-Stokes equations, Equation B.1,with an inviscid assumption, we get the Euler equations∂t ρρvE+∇ · ρvρvv + P1v(E + P ) = 0 (B.9)As before, with a further incompressible assumption and negligible pressurecontributions, the inviscid Burgers’ equation is obtained:∂v∂t+ v · ∇v = 0. (B.10)165Appendix CANSLib Command LineOptionsAn outline of command line arguments used for the in-house ANSLib solveris provided here. The command line interface allows for customization ofarguments for different test cases as well as options specific to PETSc [101],the external numerical linear algebra library. Support exists for reading thepotentially long list of arguments through a file, allow for batch processing.Additional details for options for solving the primal problem, especially forthe Spalart-Allmaras turbulence model test cases is discussed in [55]. Forthese cases, additional arguments are needed for the primal problem, suchas additional geometry and mesh information, whose details are describedin [60].For the 2D cases discussed in the current thesis, the arguments for solv-ing the ETE for these cases are similar to the other model problems. Somecommands common throughout different test cases are first highlighted be-fore listing the specifics for each test case. For parallel processing, the Mes-sage Passing Interface (MPI) protocol [36] is invoked by appending$ mpiexec -np < number of processors >before the solver executable. To use the error transport equation (ETE) inconjunction with the primal problem, the command line option-pqr <order of: primal , ETE , reconstruction for p-TE>is used. For ETE linearization,-NewtonErroris invoked. If the linearized ETE should be solved using a single GMRESiteration, the option-OneLinearis included in addition. Another feature is co-advancing the ETE with theprimal problem for unsteady flow. The option166C.1. Steady Problems-coadvanceis used for this. Examples of the list of command line arguments are in-cluded here for the main test cases considered in the results.C.1 Steady ProblemsC.1.1 Inviscid Supersonic Vortex./ Driver-physics Roe2D-f <mesh file >-mesh_type c-reconConserved-mach 2-SSV-newtonits 100-pqr 244C.1.2 Inviscid Gaussian Bump./ Driver-physics Roe2D-f <mesh file >-mesh_type c-reconConserved-mach 0.5C.1.3 NACA 0012 Laminar Navier-Stokes./ Driver-physics RoeVisc2D-f <mesh file >-mesh_type c-pqr 244-mach 0.5-angle 0-reynolds 5000-reconConserved-no_distance_weight167C.1. Steady Problems-C 0.1-show_DC-NewtonError-OneLinearC.1.4 Flat plate with SA-neg Turbulence Model./ Solver-mesh_type c-f <mesh file >-physics RoeTurbSA2D-mach 0.2-reynolds 5000000-no_distance_weight-pseudolts_fixed-C 0.1-a 2-pc_factor_levels 3-pc_factor_mat_ordering_type qmd-pqr 244-NewtonError-OneLinear-show_DCC.1.5 NACA 0012 airfoil with SA-neg Turbulence Model./ Solver-physics RoeTurbSA2D-f <mesh file >-a 2-mesh_type c-mach 0.15 #flow parameters: Mach-angle 10 #angle of attack-max_iter 100-jacobian_type recanalytic-B <bdry file >-fec <fec file >-no_distance_weight-pseudolts_fixed-pc_factor_mat_ordering_type qmd-pc_factor_levels 4-reynolds 6000000168C.2. Unsteady Problems-C 0.1-ita_target_residual 1e-6-show_DC-pqr 244-NewtonError-OneLinear-sol ~/ Downloads/n0012_449 -129-a2 -backup-pc_type lu-ksp_final_residual-output_type tecplotC.2 Unsteady ProblemsC.2.1 Diffusion./ Driver-f <mesh file >-physics Poisson-g t-final_time 0.1-i s-C 0.4-mesh_type c-ts_type cn-ts_final_time 0.1-ts_monitor-ts_exact_final_time matchstep-ts_write_solution-pqr 244-te-time_step 0.01C.2.2 Advection./ Driver-f <mesh file >-physics Advection-g t-final_time 10-i s-C 0.4169C.2. Unsteady Problems-mesh_type c-ts_type cn-ts_final_time 10-ts_monitor-ts_exact_final_time matchstep-ts_write_solution-pqr 244-time_step 0.05-post_process_every_step-periodic-coadvanceC.2.3 Translating Vortex./ Driver-f <mesh file >-physics RoeVisc2D-g t -final_time 0.2-i s-C 0.4-mesh_type c-ts_type cn-ts_final_time 0.2-ts_monitor-ts_exact_final_time matchstep-ts_write_solution-pqr 244-unsNSvortex-reconConserved-mach 2-reynolds 5000-ETE-time_step 0.02-coadvanceC.2.4 Vortex sheddingmpiexec -np 8 ./ Driver-f <mesh file >-physics RoeVisc2D-g t-final_time 1000170C.2. Unsteady Problems-i s-C 0.4-mesh_type c-ts_type cn-ts_final_time 1000-ts_monitor-ts_exact_final_time matchstep-ts_write_solution-pqr 244-vortexshed-reconConserved-mach 0.2-reynolds 200-ETE-time_step 0.125-const_visc-show_DC-post_process_every_step-force_no_write-coadvance-output_type tecplot-Bdry_Cond 5171
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical estimation of discretization error on unstructured...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical estimation of discretization error on unstructured meshes Yan, Kai Kin Gary 2018
pdf
Page Metadata
Item Metadata
Title | Numerical estimation of discretization error on unstructured meshes |
Creator |
Yan, Kai Kin Gary |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | A numerical estimation of discretization error is performed for solutions to steady and unsteady models for compressible flow. An accurate and reliable estimate of discretization error is useful in obtaining a more accurate defect corrected solution, as well as a tight uncertainty bound as error bars. The error estimation procedure is performed by solving an auxiliary problem, known as the error transport equation (ETE), solved on the same mesh as the original model equations. Unlike unsteady adjoint methods for error in functionals, the ETE requires only one other set of equations to be solved, agnostic to the choice and number of output functionals, including common aerodynamic quantities such as lift or drag. Furthermore, co-advancing the ETE in time only requires the storage of local solutions in time and not the entire history, reducing memory requirements. This method of error estimation is performed in the context of higher order finite-volume methods on unstructured meshes. Approaches based on solving the ETE can be found in the literature for uniform or smooth meshes, but this has not been well studied for unstructured meshes. Such meshes necessarily have nonsmooth geometric features, which create many difficulties in accurate error estimation. These difficulties in accurate discretizations of the ETE are investigated, including the discretization of the ETE source term, which is critical to error estimate accuracy. It was found that the proposed schemes by the ETE approach can be more efficient and robust compared to solving the higher order problem. The choices of discretization schemes need to be made carefully, and these results demonstrate how it is possible, along with justification, to obtain asymptotically accurate, efficient, and robust error estimates that can be used with vast possibilities of model equations in practice. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-03-12 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0364237 |
URI | http://hdl.handle.net/2429/64792 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_may_yan_kaikingary.pdf [ 20.51MB ]
- Metadata
- JSON: 24-1.0364237.json
- JSON-LD: 24-1.0364237-ld.json
- RDF/XML (Pretty): 24-1.0364237-rdf.xml
- RDF/JSON: 24-1.0364237-rdf.json
- Turtle: 24-1.0364237-turtle.txt
- N-Triples: 24-1.0364237-rdf-ntriples.txt
- Original Record: 24-1.0364237-source.json
- Full Text
- 24-1.0364237-fulltext.txt
- Citation
- 24-1.0364237.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0364237/manifest