Search

Close to solid walls, at high Reynolds numbers, fluids may develop steep gradients which require a fine mesh for an accurate simulation of the turbulent boundary layer. An often used cure is to use a wall model instead of a fine mesh, with the drawback that modeling is introduced, leading to possibly unstable numerical schemes. In this paper, we leave the modeling aside, take it for granted, and propose a new set of provably energy stable boundary procedures for the incompressible Navier-Stokes equations. We show that these new boundary procedures lead to numerical results with high accuracy even for coarse meshes where data is partially obtained from a wall model. Funding: Vetenskapsradet, SwedenSwedish Research Council [2018-05084 VR]; Swedish e-Science Research Center (SeRC)

Numerical solvers of initial boundary value problems will exhibit instabilities and loss of accuracy unless carefully designed. The key property that leads to convergence is stability, which this t ...

We consider accurate and stable interpolation procedures for numerical simulations utilizingtime dependent adaptive meshes. The interpolation of numerical solution valuesbetween meshes is considered as a transmission problem with respect to the underlying semidiscretizedequations, and a theoretical framework using inner product preserving operatorsis developed, which allows for both explicit and implicit implementations. The theory issupplemented with numerical experiments demonstrating practical benefits of the new stableframework. For this purpose, new interpolation operators have been designed to be used withmulti-block finite difference schemes involving non-collocated, moving interfaces. Funding:National Research Foundation of South AfricaNational Research Foundation - South Africa [89916]; Vetenskapsradet, SwedenSwedish Research Council [2018-05084_VR]

The influence of different boundary conditions on the spectral properties of the incompressible Navier-Stokes equations is investigated. By using the Fourier-Laplace transform technique, we determine the spectra, extract the decay rate in time, and investigate the dispersion relation. In contrast to an infinite domain, where only diffusion affects the convergence, we show that also the propagation speed influence the rate of convergence to steady state for a bounded domain. Once the continuous equations are analyzed, we discretize using high-order finite-difference operators on summation-by-parts form and demonstrate that the continuous analysis carries over to the discrete setting. The theoretical results are verified by numerical experiments, where we highlight the necessity of high accuracy for a correct description of time-dependent phenomena. Funding agency: The Swedish e-Science Research Centre (SeRC)

We use the behavior of the $$L_{2}$$ L 2 norm of the solutions of linear hyperbolic equations with discontinuous coefficient matrices as a surrogate to infer stability of discontinuous Galerkin spectral element methods (DGSEM). Although the $$L_{2}$$ L 2 norm is not bounded in terms of the initial data for homogeneous and dissipative boundary conditions for such systems, the $$L_{2}$$ L 2 norm is easier to work with than a norm that discounts growth due to the discontinuities. We show that the DGSEM with an upwind numerical flux that satisfies the Rankine–Hugoniot (or conservation) condition has the same energy bound as the partial differential equation does in the $$L_{2}$$ L 2 norm, plus an added dissipation that depends on how much the approximate solution fails to satisfy the Rankine–Hugoniot jump.

We develop a novel numerical scheme for the simulation of dissipative quantum dynamics following from two-body Lindblad master equations. All defining continuum properties of the Lindblad dynamics, hermiticity, positivity and in particular trace conservation of the evolved density matrix are preserved. The central ingredient is a new spatial difference operator, which not only fulfils the summation by parts (SBP) property but also implements a continuum reparametrization property. Using the time evolution of a heavy-quark anti-quark bound state in a hot thermal medium as an explicit example, we show how the reparametrization neutral summation-by-parts (RN-SBP) operator preserves the continuum properties of the theory., 34 pages, 7 figures, open-access code available via https://doi.org/10.5281/zenodo.3744460

We investigate the spatial operator in the incompressible Navier–Stokes, Oseen and Stokes equations and show how to avoid singularities associated with null spaces by choosing specific boundary conditions. The theoretical results are derived for a general form of energy stable boundary conditions, and applied to a few commonly used ones. The analysis is done on a system that simultaneously covers the nonlinear incompressible Navier–Stokes, the Oseen and the Stokes equations. When the spectrum of the spatial operator is investigated, we restrict the analysis to the Oseen and Stokes equations. The continuous analysis carries over to the discrete setting by using the summation-by-parts framework.

In this paper we study boundary conditions for linear hyperbolic systems of equations and the corresponding dual problem. In particular, we show that the primal and dual boundary conditions are related by a simple scaling relation. It is also shown that the weak dual problem can be derived directly from the weak primal problem. Based on the continuous analysis, we discretize and perform computations with a high-order finite difference scheme on summation- by-parts form with weak boundary conditions. It is shown that the results obtained in the continuous analysis lead directly to stability results for the primal and dual discrete problems. Numerical experiments corroborate the theoretical results.

High frequency errors are always present in numerical simulations since no difference stencil is accurate in the vicinity of the $$\pi $$π-mode. To remove the defective high wave number information from the solution, artificial dissipation operators or filter operators may be applied. Since stability is our main concern, we are interested in schemes on summation-by-parts (SBP) form with weak imposition of boundary conditions. Artificial dissipation operators preserving the accuracy and energy stability of SBP schemes are available. However, for filtering procedures it was recently shown that stability problems may occur, even for originally energy stable (in the absence of filtering) SBP based schemes. More precisely, it was shown that even the sharpest possible energy bound becomes very weak as the number of filtrations grow. This suggest that successful filtering include a delicate balance between the need to remove high frequency oscillations (filter often) and the need to avoid possible growth (filter seldom). We will discuss this problem and propose a remedy.

GPU-accelerated computing is becoming a popular technology due to the emergence of techniques such as OpenACC, which makes it easy to port codes in their original form to GPU systems using compiler directives, and thereby speeding up computation times relatively simply. In this study we have developed an OpenACC implementation of the high order finite difference CFD solver ESSENSE for simulating compressible flows. The solver is based on summation-by-part form difference operators, and the boundary and interface conditions are weakly implemented using simultaneous approximation terms. This case study focuses on porting code to GPUs for the most time-consuming parts namely sparse matrix vector multiplications and the evaluations of fluxes. The resulting OpenACC implementation is used to simulate the Taylor-Green vortex which produces a maximum speed-up of 61.3 on a single V100 GPU by compared to serial CPU version.

We study an incompletely parabolic system in three space dimensions with stochastic boundary and initial data. We show how the variance of the solution can be manipulated by the boundary conditions, while keeping the mean value of the solution unaffected. Estimates of the variance of the solution is presented both analytically and numerically. We exemplify the technique by applying it to an incompletely parabolic model problem, as well as the one-dimensional compressible Navier–Stokes equations. Funding agencies: European Commission [ACP3-GA-2013-605036]

To impose boundary conditions, data at the boundaries must be known, and consequently measurements of the imposed quantities must be available. In this paper, we consider the two most commonly used inflow boundary conditions with available data for internal flow calculations: the specification of the total temperature and total pressure. We use the energy method to prove that the specification of the total temperature and the total pressure together with the tangential velocity at an inflow boundary lead to well-posedness for the linearized compressible Euler equations. Next, these equations are discretized in space using high-order finite-difference operators on summation-by-parts form, and the boundary conditions are weakly imposed. The resulting numerical scheme is proven to be stable and the implementation of the corresponding nonlinear scheme is verified with the method of manufactured solutions. We also derive the spectrum for the continuous and discrete problems and show how to predict the convergence rate to steady state. Finally, nonlinear steady-state computations are performed, and they confirm the predicted convergence rates.

In the hyperbolic community, discontinuous Galerkin approaches are mainly applied when finite element methods are considered. As the name suggested, the DG framework allows a discontinuity at the element interfaces, which seems for many researchers a favorable property in case of hyperbolic balance laws. On the contrary, continuous Galerkin method obtained from a straightforward discretisation of the weak form of the PDEs appear to be unsuitable for hyperbolic problems. To remedy this issue, stabilization terms are usually added and various formulations can be found in the literature. There exists still the perception that continuous Galerkin methods are not suited to hyperbolic problems, and the reason of this is the continuity of the approximation. However, this perception is not true and the stabilization terms can be removed, in general, provided the boundary conditions are suitable. In this paper, we deal with this problem, and present a different approach. We use the boundary conditions to stabilize the scheme following a procedure that are frequently used in the finite difference community. Here, the main idea is to impose the boundary conditions weakly and specific boundary operators are constructed such that they guarantee stability. This approach has already been used in the DG framework, but here we apply it with a continuous Galerkin scheme. No internal dissipation is needed even if unstructured grids are used. Further, we point out that we do not need exact integration, it suffices if the quadrature rule and the norm in the differential operator are the same, such that the summation-by-parts (SBP) property is fulfilled meaning that a discrete Gauss Th. is valid. This contradicts the perception in the hyperbolic community that stability issues for pure Galerkin scheme exist. In numerical simulations, we verify our theoretical analysis., 28 pages, 10 figures

We introduce the notion of a transmission problem to describe a general class of problems where different dynamics are coupled in time. Well-posedness and stability are analysed for continuous and discrete problems using both strong and weak formulations, and a general transmission condition is obtained. The theory is applied to the coupling of fluid-acoustic models, multi-grid implementations, adaptive mesh refinements, multi-block formulations and numerical filtering. Funding agencies: Swedish Meteorological and Hydrological Institute (SMHI)

Multigrid schemes for high order finite difference methods on summation-by-parts form are studied by comparing the effect of different interpolation operators. By using the standard linear prolongation and restriction operators, the Galerkin condition leads to inaccurate coarse grid discretizations. In this paper, an alternative class of interpolation operators that bypass this issue and preserve the summation-by-parts property on each grid level is considered. Clear improvements of the convergence rate for relevant model problems are achieved. Funding agencies: VINNOVA, the Swedish Governmental Agency for Innovation Systems [2013-01209]

Temporal error bounds for the wave equation expressed on second order form are investigated. We show that, with appropriate choices of boundary conditions, the time and space derivatives of the error are bounded even for long times. No long time bound on the error itself is obtained, although numerical experiments indicate that a bound exists.

A one-dimensional steady-state advection-diffusion problem using summation-by-parts operators is investigated. For approximating the second derivative, a wide stencil is used, which simplifies implementation and stability proofs. However, it also introduces spurious, oscillating, modes for all mesh sizes. We prove that the size of the spurious modes is equal to the size of the truncation error for a stable approximation and hence disappears with the convergence rate. The theoretical results are verified with numerical experiments. Funding agencies:This project was funded by the Swedishe-science Research Center (SeRC). Thefunding source had no involvement in thestudy design, collection and analysis ofdata, or in writing and submitting thisarticle

The aim of the study was to describe the reoperation rates and postoperative complications associated with different methods of osteosynthesis in all extra-articular, closed fractures of the proximal and middle phalanges operated on in the Department of Hand Surgery at Södersjukhuset beween 2010-2014, and to describe the associated patient demographics.This study included all the relevant operations, which comprised operations on 181 fractures in 159 patients (82 male, 77 female), median and mean age = 43 (range = 14-88 years). The clinical records and radiographs were examined retrospectively. A logistic regression analysis was performed on the reoperation data to determine which method of osteosynthesis was the most important descriptor for reoperation, and whether the fracture type was a significant confounder.Forty-seven patients (26%) were reoperated on, mainly due to finger stiffness. The reoperation rates were 25% after K-wire, 15% after screws, and 42% after plate fixation. The unadjusted reoperation rate after plate fixation was significantly higher than for the other methods, but not after adjustment for fracture complexity. The proximal phalanx was affected in 88% of the fractures, and 77% were located in the fourth or fifth finger. Falls, animal-related, and sports injuries were the most frequent causes of injuries.Open reduction with plate fixation was associated with a higher reoperation rate, but this method was also used for the more complex fractures. Plate fixation for phalangeal fractures often entails a need for later tenolysis and plate removal. More aggressive mobilisation regimes might be indicated to prevent adhesion problems.

Many physical phenomena can be described mathematically by means of partial differential equations. These mathematical formulations are said to be well-posed if a unique solution, bounded by the gi ...

We consider a stochastic analysis of non-linear viscous fluid flow problems with smooth and sharp gradients in stochastic space. As a representative example we consider the viscous Burgers’ equation and compare two typical intrusive and non-intrusive uncertainty quantification methods. The specific intrusive approach uses a combination of polynomial chaos and stochastic Galerkin projection. The specific non-intrusive method uses numerical integration by combining quadrature rules and the probability density functions of the prescribed uncertainties. The two methods are compared in terms of error in the estimated variance, computational efficiency and accuracy. This comparison, although not general, provide insight into uncertainty quantification of problems with a combination of sharp and smooth variations in stochastic space. It suggests that combining intrusive and non-intrusive methods could be advantageous.

This paper outlines the development of a parallel three-dimensional hybrid finite volume finite difference capability. The specific application area under consideration is modeling the trailing vortices shed from the wings of aircraft under transonic flight conditions. For this purpose, the Elemental finite volume code is employed in the vicinity of the aircraft, whereas the ESSENSE finite difference software is employed to accurately resolve the trailing vortices. The former method is spatially formally second-order, and the latter is set to sixth-order accuracy. The coupling of the two methods is achieved in a stable manner through the use of summation-by-parts operators and weak imposition of boundary conditions using simultaneous approximation terms. The developed hybrid solver is successfully validated against an analytical test case. This is followed by demonstrating the ability to model the flowfield, including trailing vortex structures, around the NASA Common Research Model under transonic flow conditions. The interface treatment is shown to describe the intersecting vortices in a smooth manner. In addition, insights gained in resolving the vortices include violation of underlying assumptions of analytical vortex modeling methods. Funding agencies: National Aerospace Centre of the University of Witwatersrand, Johannesburg; South African Research (SARChI) Chair in Industrial CFD - Department of Science and Technology; South African Research (SARChI) Chair in Industrial CFD - National Research Foundation

We consider a stochastic analysis of the non-linear viscous Burgers’ equation and focus on the comparison between intrusive and non-intrusive uncertainty quantification methods. The intrusive approach uses a combination of polynomial chaos and stochastic Galerkin projection. The non-intrusive method uses numerical integration by combining quadrature rules and the probability density functions of the prescribed uncertainties. The two methods are applied to a provably stable formulation of the viscous Burgers’ equation, and compared. As measures of comparison: variance size, computational efficiency and accuracy are used.

We derive analytic solutions to the scalar and vector advection equation with variable coefficients in one spatial dimension using Laplace transform methods. These solutions are used to investigate how accuracy and stability are influenced by the presence of discontinuous wave speeds when applying high-order-accurate, skew-symmetric finite difference methods designed for smooth wave speeds. The methods satisfy a summation-by-parts rule with weak enforcement of boundary conditions and formal order of accuracy equal to 2, 3, 4 and 5. We study accuracy, stability and convergence rates for linear wave speeds that are (a) constant, (b) non-constant but smooth, (c) continuous with a discontinuous derivative, and (d) constant with a jump discontinuity. Cases (a) and (b) correspond to smooth wave speeds and yield stable schemes and theoretical convergence rates. Non-smooth wave speeds [cases (c) and (d)], however, reveal reductions in theoretical convergence rates and in the latter case, the presence of an instability.

The coupling of the compressible and incompressible Navier-Stokes equations is considered. Our ambition is to take a first step towards a provably well posed and stable coupling procedure. We study a simplified setting with a stationary planar interface and small disturbances from a steady background flow with zero velocity normal to the interface. The simplified setting motivates the use of the linearized equations, and we derive interface conditions such that the continuous problem satisfy an energy estimate. The interface conditions can be imposed both strongly and weakly. It is shown that the weak and strong interface imposition produce similar continuous energy estimates. We discretize the problem in time and space by employing finite difference operators that satisfy a summation-by-parts rule. The interface and initial conditions are imposed weakly using a penalty formulation. It is shown that the results obtained for the weak interface conditions in the continuous case, lead directly to stability of the fully discrete problem.

We present a modified formulation of the dual time-stepping technique which makes use of two derivatives in pseudo-time. This new technique retains and improves the convergence properties to the stationary solution. When compared with the conventional dual time-stepping, the method with two derivatives reduces the stiffness of the problem and requires fewer iterations for full convergence to steady-state. In the current formulation, these positive effects require that an approximation of the square root of the spatial operator is available and inexpensive. Funding agencies: Linkoping University; Swedish Governmental Agency for Innovation SystemsVinnova [2013-01209]; VINNOVAVinnova

Affordable, high order simulations of turbulent flows on unstructured grids for very high Reynolds number flows require wall models for efficiency. However, different wall models have different accuracy and stability properties. Here, we develop a kinetic energy stability estimate to investigate stability of wall model boundary conditions. Using this norm, two wall models are studied, a popular equilibrium stress wall model, which is found to be unstable and the dynamic slip wall model which is found to be stable. These results are extended to the discrete case using the Summation by parts (SBP) property of the discontinuous Galerkin method. Numerical tests show that while the equilibrium stress wall model is accurate but unstable, the dynamic slip wall model is inaccurate but stable., Accepted in Computers and Fluids

Unresolved gradients produce numerical oscillations and inaccurate results. The most straightforward solution to such a problem is to increase the resolution of the computational grid. However, this is often prohibitively expensive and may lead to ecessive execution times. By training a neural network to predict the shape of the solution, we show that it is possible to reduce numerical oscillations and increase both accuracy and efficiency. Data from the neural network prediction is imposed using multiple penalty terms inside the domain.

Artificial neural networks together with associated computational libraries provide a powerful framework for constructing both classification and regression algorithms. In this paper we use neural networks to design linear and non-linear discrete differential operators. We show that neural network based operators can be used to construct stable discretizations of initial boundary-value problems by ensuring that the operators satisfy a discrete analogue of integration-by-parts known as summation-by-parts. Our neural network approach with linear activation functions is compared and contrasted with a more traditional linear algebra approach. An application to overlapping grids is explored. The strategy developed in this work opens the door for constructing stable differential operators on general meshes.

Most computational fluid dynamics simulations are, at present, performed in a body-fixed frame, for aeronautical purposes. With the advent of sharp manoeuvre, which may lead to transient effects originating in the acceleration of the centre of mass, there is a need to have a consistent formulation of the Navier–Stokes equations in an arbitrarily moving frame. These expressions should be in a form that allows terms to be transformed between non-inertial and inertial frames and includes gravity, viscous terms, and linear and angular acceleration. Since no effects of body acceleration appear in the inertial frame Navier–Stokes equations themselves, but only in their boundary conditions, it is useful to investigate acceleration source terms in the non-inertial frame. In this paper, a derivation of the energy equation is provided in addition to the continuity and momentum equations previously published. Relevant dimensionless constants are derived which can be used to obtain an indication of the relative significance of acceleration effects. The necessity for using computational fluid dynamics to capture nonlinear effects remains, and various implementation schemes for accelerating bodies are discussed. This theoretical treatment is intended to provide a foundation for interpretation of aerodynamic effects observed in manoeuvre, particularly for accelerating missiles. Funding agencies: DRDB [KT466921, KT528944, KT470887]

The effect of an inaccurate geometry description on the solution accuracy of a hyperbolic problem is discussed. The inaccurate geometry can for example come from an imperfect CAD system, a faulty mesh generator, bad measurements or simply a misconception. We show that inaccurate geometry descriptions might lead to the wrong wave speeds, a misplacement of the boundary conditions, to the wrong boundary operator and a mismatch of boundary data. The errors caused by an inaccurate geometry description may affect the solution more than the accuracy of the specific discretization techniques used. In extreme cases, the order of accuracy goes to zero. Numerical experiments corroborate the theoretical results.

We analyze the extension of summation-by-parts operators and weak boundary conditions for solving initial boundary value problems involving second derivatives in time. A wide formulation is obtained by first rewriting the problem on first order form. This formulation leads to optimally sharp fully discrete energy estimates that are unconditionally stable and high order accurate. Furthermore, it provides a natural way to impose mixed boundary conditions of Robin type, including time and space derivatives. We apply the new formulation to the wave equation and derive optimal fully discrete energy estimates for general Robin boundary conditions, including nonreflecting ones. The scheme utilizes wide stencil operators in time, whereas the spatial operators can have both wide and compact stencils. Numerical calculations verify the stability and accuracy of the method. We also include a detailed discussion on the added complications when using compact operators in time and give an example showing that an energy estimate cannot be obtained using a standard second order accurate compact stencil. Funding agencies: Swedish Research Council [621-2012-1689]

We make an initial investigation into the temporal efficiency of a fully discrete summation-by-parts approach for unsteady flows. As a model problem for the Navier–Stokes equations we consider a two-dimensional advection–diffusion problem with a boundary layer. The problem is discretized in space using finite difference approximations on summation-by-parts form together with weak boundary conditions, leading to optimal stability estimates. For the time integration part we consider various forms of high order summation-by-parts operators and compare with an existing popular fourth order diagonally implicit Runge–Kutta method. To solve the resulting fully discrete equation system, we employ a multi-grid scheme with dual time stepping.

We construct accurate central difference stencils for problems involving high frequency waves or multi-frequency solutions over long time intervals with a relatively coarse spatial mesh, and with an easily obtained bound on the dispersion error. This is done by demonstrating that the problem of constructing central difference stencils that have minimal dispersion error in the infinity norm can be recast into a problem of approximating a continuous function from a finite dimensional subspace with a basis forming a Chebyshev set. In this new formulation, characterising and numerically obtaining optimised schemes can be done using established theory.

Stochastic problems governed by nonlinear conservation laws are challenging due to solution discontinuities in stochastic and physical space. In this paper, we present a level set method to track discontinuities in stochastic space by solving a Hamilton-Jacobi equation. By introducing a speed function that vanishes at discontinuities, the iso-zeros of the level set problem coincide with the discontinuities of the conservation law. The level set problem is solved on a sequence of successively finer grids in stochastic space. The method is adaptive in the sense that costly evaluations of the conservation law of interest are only performed in the vicinity of the discontinuities during the refinement stage. In regions of stochastic space where the solution is smooth, a surrogate method replaces expensive evaluations of the conservation law. The proposed method is tested in conjunction with different sets of localized orthogonal basis functions on simplex elements, as well as frames based on piecewise polynomials conforming to the level set function. The performance of the proposed method is compared to existing adaptive multi-element generalized polynomial chaos methods. Funding agencies: Research Council of Norway through CLIMIT program [244035/E20 CONQUER]; U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research [DE-SC0006402]; NSF [CMMI-1454601]

We study hyperbolic and incompletely parabolic systems with stochastic boundary and initial data. Estimates of the variance of the solution are presented both analytically and numerically. It is shown that one can reduce the variance for a given input, with a specific choice of boundary condition. The technique is applied to the Maxwell, Euler, and Navier–Stokes equations. Numerical calculations corroborate the theoretical conclusions.

The suitability of a discretization method is highly dependent on the shape of the domain. Finite difference schemes are typically efficient, but struggle with complex geometry, while finite element methods are expensive but well suited for complex geometries. In this paper we propose a provably stable hybrid method for a 2D advection–diffusion problem, using a class of inner product compatible projection operators to couple the non-conforming grids that arise due to varying the discretization method throughout the domain.

In this paper we generalize results regarding the order of accuracy of finite difference operators on summation-by-parts (SBP) form, previously known to hold on uniform grids, to grids with arbitrary point distributions near domain boundaries. We give a definite proof that the order of accuracy in the interior of a diagonal norm based SBP operator must be at least twice that of the boundary stencil, irrespective of the grid point distribution near the boundary. Additionally, we prove that if the order of accuracy in the interior is precisely twice that of the boundary, then the diagonal norm defines a quadrature rule of the same order as the interior stencil. Again, this result is independent of the grid point distribution near the domain boundaries.

We construct fully discrete stable and accurate numerical schemes for solving partial differential equations posed on non-simply connected spatial domains. The schemes are constructed using summation-by-parts operators in combination with a weak imposition of initial and boundary conditions using the simultaneous approximation term technique. In the theoretical part, we consider the two-dimensional constant coefficient advection equation posed on a rectangular spatial domain with a hole. We construct the new scheme and study well-posedness and stability. Once the theoretical development is done, the technique is extended to more complex non-simply connected geometries. Numerical experiments corroborate the theoretical results and show the applicability of the new approach and its advantages over the standard multiblock technique. Finally, an application using the linearized Euler equations for sound propagation is presented.

Corrigendum to “On the relation between conservation and dual consistency for summation-by-parts schemes”[J. Comput. Phys. 344 (2017) 437–439]

This is a short response to some statements by Wang and Kreiss in “Convergence of Summation-by-Parts Finite Difference Methods for the Wave equation” (Wang and Kreiss in J Sci Comput 71(1):219–245, 2017) that questions results in our paper “On the order of accuracy for difference approximations of initial-boundary value problems” (Svard and Nordstrom in J Comput Phys 218(1):333–352, 2006). We show that our results still stand.

The two dimensional advection–diffusion equation in a stochastically varyinggeometry is considered. The varying domain is transformed into a fixed one andthe numerical solution is computed using a high-order finite difference formulationon summation-by-parts form with weakly imposed boundary conditions. Statistics ofthe solution are computed non-intrusively using quadrature rules given by the probabilitydensity function of the random variable. As a quality control, we prove that thecontinuous problem is strongly well-posed, that the semi-discrete problem is stronglystable and verify the accuracy of the scheme. The technique is applied to a heat transferproblem in incompressible flow. Statistical properties such as confidence intervals andvariance of the solution in terms of two functionals are computed and discussed. Weshow that there is a decreasing sensitivity to geometric uncertainty as we graduallylower the frequency and amplitude of the randomness. The results are less sensitiveto variations in the correlation length of the geometry. Funding agencies: European Union [ACP3-GA-2013-605036]