Applied Science, Faculty of
Mechanical Engineering, Department of
DSpace
UBCV
Medina Lino, Ivonne Leonor
Master of Applied Science - MASc
University of British Columbia
We present a GPU implementation of non-Newtonian Hele-Shaw flow that models the displacement of Herschel-Bulkley fluids along narrow eccentric annuli. This flow is characteristic of many long-thin flows that require extensive calculation due to an inherent nonlinearity in the constitutive law. A common method of dealing with such flows is via an augmented Lagrangian algorithm, which is often painfully slow. Here we show that such algorithms, although involving slow iterations, can often be accelerated via parallel implementation on graphic processor units (GPUs). Indeed, such algorithms explicitly solve the nonlinear aspects only locally on each mesh cell (or node), which makes them ideal candidates for GPU. Combined with other advances, the optimized GPU implementation takes ≈ 2.5% of the time of the original algorithm.
GPU computing of yield stress fluid flows in narrow gapsbyIvonne Leonor Medina LinoB.Sc., National Autonomous University of Mexico, 2013M.Sc., The University of Nottingham, 2016A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mechanical Engineering)The University of British Columbia(Vancouver)July 2023© Ivonne Leonor Medina Lino, 2023The following individuals certify that they have read, and recommend to theFaculty of Graduate and Postdoctoral Studies for acceptance,the thesis entitled:GPU computing of yield stress fluid flows in narrow gapssubmitted by Ivonne Leonor Medina Lino in partial fulfillmentof the requirements for the degree of Master of Applied Sciencein Mechanical EngineeringExamining Committee:Ian Frigaard, Professor, Mechanical Engineering Department, UBCSupervisorDana Grecov, Professor, Mechanical Engineering Department, UBCSupervisory Committee MemberCarl Ollivier-Gooch, Professor, Mechanical Engineering Department, UBCSupervisory Committee MemberiiAbstractWe present a GPU implementation of non-Newtonian Hele-Shaw flow that modelsthe displacement of Herschel-Bulkley fluids along narrow eccentric annuli. Thisflow is characteristic of many long-thin flows that require extensive calculationdue to an inherent nonlinearity in the constitutive law. A common method of deal-ing with such flows is via an augmented Lagrangian algorithm, which is oftenpainfully slow. Here we show that such algorithms, although involving slow iter-ations, can often be accelerated via parallel implementation on graphic processorunits (GPUs). Indeed, such algorithms explicitly solve the nonlinear aspects onlylocally on each mesh cell (or node), which makes them ideal candidates for GPU.Combined with other advances, the optimized GPU implementation takes ≈ 2.5%of the time of the original algorithm.iiiLay SummaryThe study of non-Newtonian fluids is a very dynamic field of science with a lot ofscientific and technological applications. Some examples of non-Newtonian fluidare cement, mud, mayonnaise and blood. We studied an annular cementing modelthat has applications in the oil industry. Our contribution to the model consists inthe optimization and parallelization of the algorithm. We used the graphic proces-sor unit (GPU) technology better known as a graphic card to parallelize the code.We found that the original model was suitable for parallelization, and we exploitthis fact to optimize the model. We analyzed the original model and found a newalgorithm that reduces the execution time. The new version of the code and algo-rithm takes just 2.5% of the original time. This is a major advantage given that theold version could take weeks to be executed and now the times have been reducedto some hours for the most computationally challenging cases.ivPrefaceThe research of this dissertation was conducted by the author, Ivonne Leonor Med-ina Lino under the supervision of Professor Ian Frigaard. The original code wasdeveloped by Professor Frigaard and group. Dr. Mariana Carrasco modified andimproved the model and the code. She kindly provided her version to the authorwho parallelized, optimized and developed a new code. The results of this researchhave been presented in the CAIMS/SCMAI Conference 2022 and 2021. Some ofthis results have been submitted for publication as: Ivonne Medina Lino, Mari-ana Carraco-Teja and Ian Frigaards. GPU computing of yield stress fluid flowsin narrow gaps. Ian Frigaard was the supervisor on this project and was involvedthroughout the project in concept formation and manuscript composition. MarianaCarrasco-Teja also contributed providing a version of the original code and withthe manuscript edits.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Parallel Computing . . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Primary cementing displacement flows . . . . . . . . . . . . . . . 92 Model framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.1 Computational overview . . . . . . . . . . . . . . . . . . . . . . 152.1.1 Finding the streamfunction . . . . . . . . . . . . . . . . . 163 Computational improvements . . . . . . . . . . . . . . . . . . . . . 193.1 GPU parallelization . . . . . . . . . . . . . . . . . . . . . . . . . 19vi3.2 Algorithmic improvements . . . . . . . . . . . . . . . . . . . . . 213.3 Theoretical Peak Performance Calculation . . . . . . . . . . . . . 213.4 Scalability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.5 Implementation and novelty . . . . . . . . . . . . . . . . . . . . 253.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.7 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294 Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45A Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52A.1 Optimized algorithm pseudocode . . . . . . . . . . . . . . . . . . 52A.2 Convergence Analysis . . . . . . . . . . . . . . . . . . . . . . . 54viiList of FiguresFigure 1.1 Description of the primary cementing process by stage. . . . . 3Figure 1.2 (a) Scheme of transversal view of a well. (b) Scheme of across-section of the well that describes the eccentric annulargeometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 1.3 Draw of the eccentric annulus model. . . . . . . . . . . . . . 5Figure 1.4 The GPU Devotes more transistors to data processing. Takenfrom [48] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8Figure 3.1 Test of computational performance between GPU vs. CPU us-ing GeForce RTX 2080 Ti. . . . . . . . . . . . . . . . . . . . 20Figure 3.2 Case A concentration colormap, computed in parallel. Se-quence of times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59s, from left to right. The white lines show the streamlines φ =0, 0.1, 0.2, ..., 1. . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 3.3 Case B concentration colormap, computed in parallel. Se-quence of times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59s, from left to right. The white lines show the streamlines φ =0, 0.1, 0.2, ..., 1. . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 3.4 Case C concentration colormap, computed in parallel. Se-quence of times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59s, from left to right. The white lines show the streamlines φ =0, 0.1, 0.2, ..., 1. . . . . . . . . . . . . . . . . . . . . . . . . . 31viiiFigure 3.5 Case D concentration colormap, computed in parallel. Se-quence of times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59s, from left to right. The white lines show the streamlines φ =0, 0.1, 0.2, ..., 1. . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 3.6 Case E concentration colormap, computed in parallel. Se-quence of times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59s, from left to right. The white lines show the streamlines φ =0, 0.1, 0.2, ..., 1. . . . . . . . . . . . . . . . . . . . . . . . . . 33Figure 3.7 Case G, computed in parallel. Sequence of time: 2042.14 s.,5105.26 s., 7147.33 s., 1010.59 s. from top to bottom. Thewhite lines show the streamlines φ = 0, 0.1, 0.2, ..., 1. . . . . . 34Figure 4.1 Case A: a) computed solution; b) analytical advected interface;c) comparison. . . . . . . . . . . . . . . . . . . . . . . . . . 38Figure 4.2 Case E in optimized parallel at t = 4333; (b) Case E in nonop-timized serial at t = 4333; (c) Case F in optimized parallel at t= 4982; (d) Case F in non-optimized serial at t = 4982. . . . . 39Figure 4.3 Case A: a) L2 error between optimized parallel and non-optimizedserial codes; b) annular timing benchmark; red line: non-optimizedserial code; blue line: optimized parallel code. . . . . . . . . . 39Figure 4.4 Case B: a) L2 error between optimized parallel and non-optimizedserial codes; b) annular timing benchmark; red line: non-optimizedserial code; blue line: optimized parallel code. . . . . . . . . . 40Figure 4.5 Case E: a) L2 error between optimized parallel and non-optimizedserial codes; b) annular timing benchmark; red line: non-optimizedserial code; blue line: optimized parallel code. . . . . . . . . . 40ixList of Symbolsγ˙ . . . . . . . . . . . . . . . . Shear rateτY . . . . . . . . . . . . . . Yield stressκ . . . . . . . . . . . . . . . Consistency factorn . . . . . . . . . . . . . . . . Power-law indexRe . . . . . . . . . . . . . . Reynolds numberA . . . . . . . . . . . . . . . Cross-sectional areau . . . . . . . . . . . . . . . Velocity vectoru . . . . . . . . . . . . . . . . Radial velocityv . . . . . . . . . . . . . . . . Azimuthal velocityw . . . . . . . . . . . . . . . Axial velocityβ . . . . . . . . . . . . . . . Angle of inclinationµ . . . . . . . . . . . . . . . Fluid viscosityg . . . . . . . . . . . . . . . . Gravity of Earthza . . . . . . . . . . . . . . . Annulus surface heightzbh . . . . . . . . . . . . . . Bottom hole heightr . . . . . . . . . . . . . . . . Pipe radiusxri . . . . . . . . . . . . . . . Inner radius of the annulusro . . . . . . . . . . . . . . . Outer radius of the annuluse . . . . . . . . . . . . . . . . Eccentricity of the annulustˆ = π rˆ∗awˆ∗ t . . . . . . . . . Timeβ . . . . . . . . . . . . . . . Angle of inclinationwˆ . . . . . . . . . . . . . . . Mean velocitylˆ . . . . . . . . . . . . . . . . Length of the casinglˆ . . . . . . . . . . . . . . . . Length of the casingQˆ = π(rˆ20− rˆ2i ) . . . Imposed flow raten j . . . . . . . . . . . . . . . Fluid j power-law indexm j = 1/n j . . . . . . . Fluid j inverse power-law indexτY, j = τˆY, j/τˆ∗ . . . . Fluid yield stressτˆ∗} . . . . . . . . . . . . . Viscous stress scaleˆ˙γ∗ = 2wˆ∗/(rˆ0− rˆi) Strain rateκ j . . . . . . . . . . . . . . . Fluid j consistencySt∗ . . . . . . . . . . . . . . Stokes numberxiAcknowledgmentsThis research has been carried out at the University of British Columbia, supportedfinancially by Schlumberger and NSERC through CRD project No. 514472-17.This support is gratefully acknowledged. I.M. gratefully acknowledges supportfrom the University of British Columbia via a university FYF award.I would like to kindly and deeply thank my supervisor Ian Frigaard for hisenormous support from the beginning of the project and all his insightful commentsand recommendations. I am specially thankful with Dr. Mariana Carrasco for allher kindness and relentless disposition to teach me and answer all my questions andher total wiliness of collaboration, without her technical help this project wouldn’tbe possible.Finally, I want to express my love and appreciation to all my family who sup-ported me during all this years, specially my two beautiful and beloved kids Ed-uardo and Sofia for whom I work everyday.xiiDedicationTo my beloved mother Marı´a del Carmen Lino Ramı´rez. (March 3, 1948 – Febru-ary 10, 2021)xiiiChapter 1IntroductionThe aim of this thesis is the development of a parallel implementation using GraphicProcessing Unit (GPU) technology of a yield stress fluids in narrow annular gapsmodel. Yield stress fluids are non-Newtonian fluids that occur commonly in nat-ural and industrial settings [2]. These fluids do not deform unless the deviatoricstress exceeds a particular threshold, the yield stress. Often, yielding behaviouris accompanied by other rheological behaviours. However, here we consider onlysimple yield stress fluids [23], as typified by rheological models such as the Bing-ham, Casson and Herschel-Bulkley fluids.1.1 Literature reviewThe process in the oil industry, which consists of depositing cement between theformations and the casing in the annulus of a wellbore, is called primary cement-ing, described in [27], outlined in [34] and explained in more detail in [37]. Theprincipal purpose of primary cementing is to provide zonal isolation in water, gas,and oil wells to exclude gas or water in one part from oil to another part of the well.Hydraulic isolation must be placed between the formation and the casing. At thesame time, the cement sheath has to be free of fluid channels.Primary cement is very important to perform a safe operation in the well. Itavoids drilling hazards. Failing to place proper primary cementing can lead to welloperational problems, dangerous environmental accidents, and early well abandon-1ment.It turns out an important question about this process is: is it turbulent or laminarflow the best regime to cement a well? In the literature review we find that [27]and [7] suggest turbulent regime should be use preferably. In [14, 31, 33] a guidein laminar mud displacement is developed. A model base simulation is used in[3, 52]. Some analytical result on the topic are in [11, 51, 53]. Simulations that arecurrently used in the industry can be found in [5, 26, 49, 61]. Some of them usesophisticate mechanical models, [3, 38]. A description of local mechanical effectsare included in [1, 24, 73, 73]. Description of defects that can appear during thecement process are described in [39, 40, 63, 69]. Leakage factors are studied in[17, 70, 71].In the next paragraphs we will briefly described the theory behind a Non-Newtonian flow thorough an eccentric annulus. Different flow regimes and nondi-mensionalization are described as well as algorithms and numerical solution. Af-terwards, concepts of parallel computing are introduced, and in the last part, weexplain the procedures that we implemented to optimize and parallelise the code.Primary cement should be planned carefully and executed flawlessly. Thereis just one opportunity to pour the slurry between the borehole and the casing. Ashort time after placement, the cement hardens and seals the wellbore. Avoidingleakages, mixing between cement and mud, and correct centralization is critical toperform a successful completion and well drilling [27].We show 1.1, a schematic of the primary cementing process. This may bedescribed as follows:After an exploration process, drilling engineers plan a trajectory that maxi-mizes wellbore exposure to the pay zones and goes ahead with the drilling. Infigure 1.1 a) it is shown the drilling of a new stage of the well. The next step is toremove the drill pipe as in b). In c) The steel casing it is inserted. After the casingis set, a spacer fluid like water is pumped and after that lead and tail slurry can beseen in e). In the last step, the mud and water are displaced and replaced by thecement as in f). Finally, cement cover all the annulus and the operation is ended aswe can see in g).2a) Drill new stage of welle) Pumplead & tailslurryf) Displace mud inannulusg) End of operationd) Pumpspacerfluidb) Remove drillpipec) InsertsteelcasingFigure 1.1: Description of the primary cementing process by stage.We use the model as in [34] and [4] to define our variables. It is used in a curvi-linear coordinate system. The following figure 1.2(a) where rˆ measures the lengthfrom the center of the casing, θ is the azimuthal coordinate, and ξˆ is the depthfrom the bottom-hole upwards. The inclination of the well is denoted β (ξˆ ). Thecross-section of the well as in figure 1.2(b), with inner radius rˆi(ξˆ ), representingthe outer radius of the casing and outer radius rˆo(ξˆ ) representing the inner radiusof the drilled hole. The configuration is assumed to be eccentric and the radii areknown. The distance between the two centers of the cylinders, eˆ(ξˆ ), is given. Itis assumed that the cylinders do not touch, i.e. eˆ(ξˆ ) < 2dˆ(ξˆ ) and the variation inthe geometry of the model with ξˆ is slow. Another assumption is that the casingremains stationary, and the lower part of the well is found on the narrow side of theconfiguration. [10] and [64] explore these assumptions in more detail.3β(ξ) ξθr(ξ)ir(ξ)0e(ξ)(a) (b)Figure 1.2: (a) Scheme of transversal view of a well. (b) Scheme of a cross-section of the well that describes the eccentric annular geometry.At each depth ξˆ , the mean radius rˆa(ξˆ ) and the mean half-gap width dˆ(ξˆ ) aredefined by:rˆa(ξˆ )≡ 12 [rˆo(ξˆ )+ rˆi(ξˆ )], dˆ(ξˆ )≡12[rˆo(ξˆ )− rˆi(ξˆ )]. (1.1)The model also can be explained with figure 1.3. The draw in the left siderepresents the annular displacement with an eccentric disposition. We can see themud and the cement flowing upwards. In the right side, we can see the draw ofthe unwrapped annulus into a Hele-Shaw displacement cell. The fluid moves inthe same direction, in the same space but now in the Cartesian system. Calculationof turbulent flows for the Newtonian version of this model and non-linearity canbe found in [12, 22, 42, 44, 59]. An analysis of the mass transport and dispersioneffects are included in [21, 62, 65, 74].A yield stress rheology introduces nonlinearity in the deviatoric stress tensorand typically a singularity in the effective viscosity. One often wants to resolve theyielding behaviour exactly. Methods for doing so were first pioneered by Glowin-4Figure 1.3: Draw of the eccentric annulus model.ski in the 1970s [67], as outlined in the reviews [25, 58]. Although there havebeen many advances in meshing/discretization, implementations and applications(e.g. [16, 57, 68]), the basic method for computing the exact yield stress flow re-tains its original framework, as follows. First, the flow problem on each time step isrewritten as a velocity (v) minimization. The yield stress means that the minimiza-tion is non-differentiable for flows that contain unyielded regions. To work aroundthis, the strain rate tensor is replaced by a new variable (q). The constraint aris-ing from this relaxation is then enforced using a Lagrange multiplier (λ ), so thatthe original convex minimization is now a Lagrangian saddle point problem. Fornumerical stability the Lagrangian is augmented with a penalty term that vanishesat the solution. The augmented Lagrangian problem is solved using an iterativeUzawa algorithm, which has 3 steps, repeated iteratively until convergence. Step 1is an elliptic linear problem for v, solved over the whole flow domain. The secondstep is a nonlinear problem for q, solved locally on each element. The third step isa simple projection/update of λ .The inherent interest for GPU parallelization comes in the structure of the basicalgorithm described above. In solving a problem on N cells, the dominant costmight be thought to arise from the linear algebra problem of step 1. However,5this same linear problem is solved repeatedly within the Uzawa iteration. For amoderate sized problem, one can use e.g. a QR decomposition, after which thestep 1 problem requires O(N) operations at each iteration. For very large problemsthe linear step 1 might require an iterative solution but good initial guesses andpreconditioners are likely available to accelerate. Thus, although step 1 can be asignificant cost timewise, it is a standard problem and not the focus here.In step 2 the problem is more interesting. In a very simple application such as aBingham fluid, this nonlinear problem is a quadratic equation solved algebraically.In such cases step 2 is computationally trivial. However, for more complex rheo-logical laws and/or for modeling techniques common in long-thin geometries, aswe treat here, the local nonlinear problem in step 2 must itself be solved numer-ically. We then enter a regime where step 2 is (by far) the dominant cost com-putationally. However, the beauty of the underlying mathematical framework isthat step 2 can be solved locally on each cell, making this the perfect problemfor a highly parallel method such as GPU programming. In this thesis, we indeeddemonstrate the efficacy of such an implementation.The class of flows that we consider are flows with small aspect ratios. Theseoccur widely, in both thin film flows and Hele-Shaw cell applications. Ours isa version of the latter and it is motivated by the industrial process of primary ce-menting, used to complete oil & gas wells prior to production [41]. In the process asequence of fluids (mainly yield stress fluids) are pumped upwards along a narroweccentric annulus formed by the borehole and the outside of a steel casing insertedinto the well. In practice the geometry varies due to changes in eccentricity, that area result of both geomechanical strains on the borehole and due to geological anddrilling induced defects (e.g. washouts). Thus, we have both a complex geometryand a complex rheology to contend with.To illustrate different applications that lead to similar mathematical structure,see for example the food processing flows modelled in [20], or some of the ice-sheet flows reviewed in [60], where analogous “yielding” comes from frictionalbehaviour at the base of the ice sheet. Indeed, disregarding yield stress applica-tions, it is a common modelling technique to use asymptotic methods to reducethe momentum balance to a simplified shear flow, so that more complex physi-cal features can be added and resolved. This can include multi-phase scenarios,6heat transfer, species transport etc., as well as certain boundary layer flows. Themain point is that the narrow dimension has a physical model that should be solved(numerically) prior to being able to compute the flow in the other 2 dimensions.1.2 Parallel ComputingFor the last decades, one important method to improve the performance of comput-ers has been to increase the speed of the processor clock operation. In recent years,however, due to different limitations in the fabrication of integrated circuits, themanufacturers have looked for alternatives to this traditional source of increasedcomputational power. After 2015, manufacturers began offering processors withtwo processors and, in the following years they add up to eight so far. A differ-ent technology called Graphic Processing Unit enters the market in the late 1980s.In the beginning, these devices had the purpose of displaying three-dimensionalgraphics. They were designed to produce a colour for every pixel on the screenusing programmable arithmetic units known as pixel shaders. The arithmetic wascontrolled by the programmer entirely, which means the input could be any data.At the beginning, there were some difficulties due to the programming lan-guages to interact with the GPU. In November 2005, NVIDIA unveiled the firstGPU to be built with NVIDIA’s CUDA architecture [56],[48],[46],[47]. To reachthe maximum number of developers possible, the CUDA language developed theindustry standard of C language and used a small number of keywords. It was alsoincluded a compiler for this new CUDA C. All these were designed to facilitate thegeneral-purpose computing on GPUs.From that time a huge variety of application have been developed in CUDAC, like in medical imaging, computational fluid dynamics, environmental scienceamong others [13],[19], [56],[45].The tremendous horsepower and very hight memory bandwidth delivered by aGPU is illustrated in Figure 3.1. For this benchmark, it is used a model GeForceRTX 2080 Ti. The graph estimates the peak performance on the GPU and theCPU, in Gigafloating-point operations per second (GFLOPS). We can see that ifwe have more operations to calculate, the difference between the GPU and the CPUperformance is bigger.7The architecture of the GPU is devoted to data processing, and more transistorsare used to do calculations. This is shown in the left image of figure Figure 1.4 andis the main difference with the CPU where data catching and flow control is themain task and reason behind a better capability of the GPU [48].The increase in speed comes at a cost. The access to the memory in the archi-tecture of a GPU like in the right figure in Figure 1.4 works like a bottleneck. Datamust be sent from the CPU to the GPU before the calculations and then retrievedfrom it afterwards. The memory access is slower in the GPU than in the CPU andbeside the data has to pass through a PCI Express bus from the GPU to the CPU,causing the bottleneck of the information.Due to this drawbacks, we have to apply two criteria to decide if an applicationis suitable to be parallelizable:1. It is computationally intensive. The time spends on transferring data to anfrom a GPU memory is significantly less than the time in computation.And 2. The computation can be broken down into hundreds or thousands ofindependents unit of work, which means that they can be massively parallelizable.If the program does not satisfy these criteria, it can run slower.Figure 1.4: The GPU Devotes more transistors to data processing. Takenfrom [48]A GPU is especially well-suited to solve problems with data-parallel computa-8tion, i.e. the same operation or function is executed at the same time or in parallel.This allows the memory access latency to be hidden with high arithmetic intensityinstead of the big data cache. Applications that process large data sets can use adata-parallel programming model to speed up computations.A multi-threaded program is partitioned into a block of threads that execute in-dependently from each other. Therefore a GPU with more Cuda cores will executecode in less time than a GPU with fewer multiprocessors.This programming model is called scalable and allows the GPU architecture tobe scalable to a more multiprocessors and memory partitions.1.3 Primary cementing displacement flowsAn example of the above is given by our work on cementing displacement flows.The model considered here was derived by [3], and we outline it below. It assumesthe flow is laminar, uses scaling arguments to reduce the Navier-Stokes equationsto a local shear flow, and then averages across the narrow annular gap to derivean elliptic equation for the stream function (essentially a Hele-Shaw approach).The mathematical foundation of this model was developed by [51], who put themodel into its variational framework and established key properties of the functionsinvolved. Ideal cementing flows involve a steady state displacement, in which adisplacement front advances along the well at a uniform pumping speed. Thatsuch solutions exist and how to predict them and incorporate them into processdesigns for vertical wells was studied by [51, 53]. Using the same variationalformulation, [52] developed an augmented Lagrangian algorithm for the cementingdisplacement flow model. This numerical solution is used to study the effects ofincreasing eccentricity, yield stress and how the interface can become unsteady.The early 2000’s also saw the rise of horizontal wells as a dominant architec-ture for oil and gas wells in many parts of the world. The change in well orientationcauses the annulus to be more eccentric and buoyancy effects, due to fluid densitydifferences, promote stratification. Carrasco-Teja et al. developed analyses, basedon the Bittleston et al. model [3], specifically designed for horizontal well ce-menting [11]. They were able to show that steady displacements may still arise inhorizontal wells, even with significant density differences, and that the axial length9of the interface would grow as the density difference increased [11].Although these results are reassuring, in many situations a long slumping inter-face is simply not desirable, e.g. because such stratified interfaces tend to destabi-lize. An industry practice developed to aid horizontal well displacements has beento rotate the casing during displacement, with the idea that the induced Couettecomponent will drag the fluids around in the azimuthal direction. This was firstmodelled by [9, 10]. Although it is relatively easy to model for Newtonian fluids[9], for non-Newtonian fluids the nonlinearity results in a local closure problemthat must be resolved as a 2D flow, i.e. the Poiseuille and Couette components ofthe flow do not decouple. When yield stress fluid are present this local closureproblem again is computed using an augmented Lagrangian method, becomingpainfully slow to compute.More recently, this modeling approach has been extended to more complexcementing scenarios. For example, turbulent, transitional and mixed regimes aremodelled [34–36, 54]. Such models must cover all flow regimes, which meansthat the yield stress dominated behaviour at slow flow rates is augmented by othernonlinear behaviours at high flow rates, eventually becoming turbulent. The aug-mented Lagrangian algorithm is again used computationally. Both the rotating cas-ing scenarios and dealing with mixed flow regimes require progressively expensivemodels for the local closure problem, i.e. due to transitions between regimes orother parametric variations. In all cases where the local closure problem is expen-sive, a GPU implementation is progressively attractive.This thesis is divided as follows. In Chapter 2 we briefly describe the theo-retical framework of the annular displacement flow model. In the next section wedescribe the software and hardware that we used to parallelize our model and thecriteria that we used to choose them. The fourth section is related to the modifica-tion of the original algorithm that we developed to optimize the code, as well as thenew parallel code. We explain in detail this mathematical modification. We alsoshow the pseudocode with its three variants. We included a subsection with con-vergence analysis. The fifth section is devoted to show our results in serial and inparallel and the corresponding comparison. We included error analysis and physi-cally relevant examples to validate the algorithm’s change and the parallelization.10Chapter 2Model frameworkThe approach uses scaling arguments to reduce the Navier-Stokes equations locallyto a shear flow. The main consideration is thatannular gap ≪ annular circumference ≪ axial length-scale, (2.1)together with some requirements that the geometry vary slowly along the well. Thereduced equations are then averaged across the annular gap, in the same way that aHele-Shaw model is derived.The annular gap half-width H(φ ,ξ ) is given by:H(φ ,ξ ) = H¯(ξ )[1+ e(ξ )cosπξ ] (2.2)∂∂φ[Hv¯]+∂∂ξ[Hraw¯] = 0. (2.3)As the gap-averaged flow is two-dimensional, the velocity components are bestrepresented via a stream function, automatically satisfying the incompressibilitycondition above.∂Ψ∂φ= Hraw¯,∂Ψ∂ξ=−Hv¯. (2.4)Determining a field equation for Ψ is addressed in detail in [3]. It starts withthe observation that the streamlines are parallel to the modified pressure gradientof the flow, i.e. the pressure gradient minus the mean static pressure gradient. This11allows for a local re-orientation of the coordinates along the streamline. In theselocal coordinates the flow is basically a Poiseuille flow through a plane channel.Cementing fluids are characterized as Herschel-Bulkley fluids, with 3 parameters:yield stress τY , consistency κ and power-law index n. For a plane Poiseuille flowthe relationship between the modified pressure gradient G, and the area flow rate|∇aΨ|, is given by:|∇aΨ|=0, χ ⩽ 0,H2+1/nκ1/n(2+1/nχ1+1/n(χ+ τY/H)2[χ+(2+1/n)τY(1+1/n)H], χ > 0.(2.5)The variable χ measures the excess of the pressure gradient G over the criticalvalue τY/H at which we transition from flow to no-flow. The operator∇a is definedas:∇a =(1ra∂∂φ,∂∂ξ).To understand (2.5), note that GH gives the wall shear stress in the annular gap.If GH ≤ τY then the fluid does not deform and consequently there is no flow,i.e. |∇aΨ|= 0. When χ =G−τY/H > 0, there is a flow and, due to the nonlinearityof the fluid, the flow rate is no longer linearly related to the pressure gradient. Thisnonlinearity is described by the second part of (2.5), which is an implicit equationfor χ(|∇aΨ|). The properties of χ(|∇aΨ|) determine the function space in which aweak solution may be found, as discussed in [51]. Note that when the rheologicalparameters are: n = 1, τY = 0, we recover a Newtonian fluid behaviour.The flow law (2.5) can be used to either define an equation for the pressure orfor the stream function. As the pressure gradient is not determined in regions where|∇aΨ| = 0, the latter is the preferred formulation. Indeed determining regionswhere |∇aΨ|= 0 has practical relevance for the cementing process, often indicatingdefective cementing. The stream function satisfies∇a ·S =− f , (2.6)12where S =[raχ(|∇aψ|)|∇aψ| +raτYH|∇aψ|]∇aψ, ⇔ |S|> raτYH ,|∇aψ|= 0, ⇔ |S| ≤ raτYH .(2.7)The term f represents the influence of buoyancy forces:f = ∇a ·(raρ(c¯)cosβSt∗,raρ(c¯)sinβ sinπφSt∗)= ∇a · f˜, (2.8)where St∗ is the Stokes number, representing the balance of viscous to buoyancystresses, as defined in Nomenclature; see also [3]. This term arises from cross-differentiating to eliminate the pressure, which results in taking derivatives of thestatic pressure gradient components.Equation (2.7) is solved for Ψ at each time, in order to derive the velocityfield. The flows considered are however displacement flows. To model the differentfluids a concentration vector c¯ is introduced, representing the gap-averaged volumefraction of each fluid. Each individual fluid component c¯k satisfies∂∂ t[Hrac¯k]+∂∂φ[Hv¯c¯k]+∂∂ξ[Hraw¯c¯k] = 0, k = 1,2, .. (2.9)Here we assume only 2 fluids, with fluid 2 displacing fluid 1. Conservation ofvolume means that c¯1+ c¯2 = 1. Note that in some formulations, the above equationis supplemented by a diffusion term on the right-hand side which can represent bothmolecular diffusivity and dispersion effects.The fluid concentrations advance in time. Time dependency may also enterin the boundary conditions, e.g. the amount of each fluid that is pumped into theannulus. Coupling of c¯ andΨ is two-way. The concentrations are used to define the(mixture of) rheological properties τY , κ and n, from those of the individual fluids.The density ρ(c¯) also depends on the fluids, and hence the buoyancy comes fromgradients in ρ(c¯). In the other direction, the fluid velocities determine transport ofc¯.We impose the boundary conditions of symmetry at the wide and narrow sides13for the concentration equation:∂ c¯k∂φ= 0, φ = 0,1. (2.10)At the bottom of the annulus, inflow concentrations are specified and at the top anoutflow condition is used. For the elliptic second-order equation (2.7), on the wideside of the annulus we impose:Ψ(0,ξ , t) = 0, (2.11)and on the narrow side:Ψ(1,ξ , t) = Q(t), (2.12)where Q(t) is an appropriately scaled flow rate. In the axial direction we imposeDirichlet conditions given by:Ψ(φ ,0, t) =Ψ0(φ , t), Ψ(φ ,Z, t) =ΨZ(φ , t), (2.13)with Ψ0 and ΨZ specified. These boundary conditions can vary for different op-erational settings, or for slight variants of the model. For example, when the fullannulus is modelled, φ is extended to [0,2] and all variables are made periodic inφ .142.1 Computational overviewWe now describe what is the typical computational challenge of simulating a pri-mary cementing displacement flow. The model (2.4)-(2.13) is solved over a spa-tial domain of size 1× Z. Since oil and gas wells tend to be very long, Z ∼5× 102− 104, and even Z ∼ 102 for a laboratory experiment; equation (2.1). Asthe aim is to displace the drilling mud and replace it with the cement slurry, in anysimulation the fluids must travel the length of the well Z, and the transport equation(2.9) is advanced over a total time T ∼ Z, given that the velocities have been scaledto be O(1). At each time step both the stream function equation (2.7) and the trans-port equation (2.9) must be solved. Equation (2.9) imposes a CFL constraint onthe time step used, but is otherwise relatively cheap. In more complex models ofthe mass transport diffusive effects may also be included (molecular and turbulentdiffusivities, Taylor dispersion, etc.; see [34]). Nevertheless, in the dimensionlesssetting the (physically scaled) diffusivities remains small enough that an explicitmethod can be used without compromising the time step restrictions, i.e. the CFLcondition can be considered the main time step restriction.The above represent the unmovable constraints of the computational problem,in terms of number of time steps that need to be executed, for any given spatialdiscretization. The spatial discretization is determined by accuracy and resolutionrequirements. Here we will assume a simple rectangular mesh. Although one canelongate the mesh spacing in ξ to help reduce the impact of Z ≫ 1, this comes witha loss of resolution. Regarding azimuthal resolution, the point of the annular Hele-Shaw approach is to resolve 2D effects around the displacement front (secondaryflows), so that Nphi = 10−100 azimuthal meshpoints would be sensible. Assuminga mesh dφ = 1/Nphi in the azimuthal direction, we will frequently want dξ ∼ dφ , toavoid distorting the mesh too much. This might also be a reason for not selectingNphi too large. Other factors in selecting dξ include the fact that for industrialsettings, the eccentricity, mean radius and gap width may all vary considerablyalong a well and one wishes to capture these variations. On the other hand, asthe Hele-Shaw approach involves averaging across the annular gap, resolving flowfeatures on the gap length-scale is the limit of model validity. Thus having balancedall consideration, the number of meshpoints is N ∼ O(ZN2phi), with potentially a15large constant of proportionality (reflecting the mesh aspect ratio).The second-order equation (2.7) is elliptic and nonlinear. Broadly speaking,at each time step we must solve this nonlinear problem of size N, and the usualmethods involve repeatedly resolving a sequence of linear problems, each of sizeN. Although of size N, these elliptic problems are sparse and the matrices aregenerally banded with width 3Nφ . For now lets assume that (2.7) is solved in atime ∝ AALGNaALG on each time step, where AALG & aALG(≥ 1) will depend on thealgorithm used and its implementation. Lets now consider the total time requiredto simulate a displacement, Tcomp. We suppose an end time T ∼ Z and that a CFLcondition imposed, which lead to:Tcomp ∝ [AALGNaALG +O(N)]ZNphi.The O(N) term comes from the cost of solving (2.9), assumed cheap. The maincost at each time step is in resolving (2.7), which we discuss below and focus onimproving in this thesis. Optimistically, we might be able to ensure that aALG = 1,as we will discuss, which leaves AALG to be improved upon. Nevertheless, we seethat at best we end up with:Tcomp ∝ AALGZ2N3phi.Although we work on reducing AALG in this thesis, the other terms above remainimmovable and Z ≫ 1 is a serious constraint to rapid computations in realistic wellgeometries.2.1.1 Finding the streamfunctionAlthough elliptic, the streamfunction equation (2.7) is not defined in regions wherethe velocity is zero, i.e. where the gradient of Ψ vanishes. Instead of the classicalformulation, we adopt the weak formulation of (2.7) as a variational inequality, asis discussed at length in [51]. This leads naturally to a minimization problem forthe streamfunction (at each time step), formally written as:minΦ∈V¯0(Ω){J(Φ)}, (2.14)16where Ω = (0,1)× (0,Z). The functional only depends on the gradient of thevelocity: q=∇aΦ, and can be split as follows. We have homogenized the boundaryconditions for Ψ, along φ = 0,1, and at the ends ξ = 0,Z, by settingΨ=Ψ∗+Φ, Ψ∗ ∈ V¯ (Ω), Φ ∈ V¯0(Ω), (2.15)with suitably constructed Ψ∗. In the continuous setting we takeV0(Ω) = {Φ ∈C∞(Ω) : Φ(0,ξ ) = 0, Φ(1,ξ ) = 0},V (Ω) = {Φ ∈C∞(Ω) : Φ(0,ξ ) = 0, Φ(1,ξ ) = 1},and the test spaces V¯0(Ω) and V¯ (Ω) are the closures of V0(Ω) and V (Ω), respec-tively, with respect to the W 1,1+n(Ω) norm. Hence, X = L1+n(Ω), and n wouldbe the smallest value of n in Ω, (which may change with time if the fluid mix).Nevertheless, practically speaking we work with finite dimensional spaces for thenumerical solution (when n = 1 is equivalent).This method replaces the minimization with a saddle point problem, that resultsfrom introducing an additional variable for the velocity gradients. The Lagrangiansaddle point problem is also typically augmented with a stabilizing term. Thus, theaugmented Lagrangian problem concerns the functional Lr:Lr(Φ,q,µ) = J(q)+∫Ωµ · (∇aΦ−q)dΩ+ r2∫Ω|∇aΦ−q|2dΩ (2.16)The saddle point of Lr, ∀r > 0, is found using an iterative Uzawa algorithm, inwhich successive iterates solve for one variable.{φ n,pn,λ n} are found numerically.For example, holding {pn−1,λ n−1} fixed, we solve for Φn. The problem forΦn minimizes Lr(v,pn−1,λ n−1) with respect to v, which we see is a linear ellip-tic problem (a Poisson equation). The discrete matrix does not change on eachiteration, which means that this problem can generally be solved cheaply.The second problem, to define pn, holds {Φn,λ n−1} fixed, and minimizesLr(Φn,q,λ n−1) with respect to q. After some algebra, we find that ∇aΨ∗ + qshould be parallel to the vector (λ n − 1+ r∇a(Ψ∗+Φn)− f˜). If x = |λ n − 1+17r∇a(Ψ∗+Φn)− f˜| then (∇aΨ∗+q) has magnitude θx, where θ ≥ 0 minimizesM(θ) =12∫ (θx)20χ(s1/2)s1/2ds+τYH|θ |x+ r2(θx)2−θx2. (2.17)This local minimization has two cases. First, if x≤ τY/H then M(θ) is minimizedby θ = 0, with solution qn = −∇Ψ∗, that represents an unyielded flow. The othercase represents yielded flow, when x > τY/H. In this second case we simply dif-ferentiate M(θ) and equate to zero to obtain the minimum:χ(θx)+τYH+ rθx− x = 0 (2.18)Recall that χ is defined by (2.5), Which is an implicit definition of χ as a functionof the areal flow rate |∇aΨ|. Thus, the problem for qn at each iteration is time con-suming, i.e. we must solve (2.18) iteratively and each evaluation of χ would requirethe implicit function χ in (2.5) to be inverted. This is a significant contribution tothe total computational time.18Chapter 3Computational improvementsWe refer to the algorithm explained above as the serial algorithm. We explore theperformance of this code, which is predictably slow, and improve it. In order toaccelerate the code, we have implemented two strategies: GPU parallelization andalgorithmic improvements.3.1 GPU parallelizationGenerally, we have to apply two criteria to decide if an application is suitable to beparallelized via GPU:1. It is computationally intensive. The time spends on transferring data to aGPU memory is significantly less than the time in computation.2. The computation can be broken down into hundreds/thousands of indepen-dent units of work, which means that they can be massively parallelized.If the program does not satisfy both criteria, it may run slower on a GPU. A GPUis especially well-suited to solve problems with data-parallel computation, i.e. thesame operation or function is executed at the same time or in parallel. This allowsthe memory access latency to be hidden with high arithmetic intensity instead ofthe big data cache. Applications that process large data sets can use a data-parallelprogramming model to speed up computations. A multi-threaded program is parti-tioned into a block of threads that execute independently of each other. Therefore,1910 4 10 6 10 8Matrix size0100200300400500Calculation Rate (GFLOPS)GPUCPUFigure 3.1: Test of computational performance between GPU vs. CPU usingGeForce RTX 2080 Ti.a GPU with more CUDA cores will execute code in less time than a GPU withfewer multiprocessors.Turning to our annular model, the mathematical operations are varied. We willprovide a profile of the main computational costs below, but overall the main costcomes from the augmented Lagrangian algorithm outlined in 2.1.1. The key pointis that the most expensive step is solved locally (i.e. independently) on each finitevolume cell, allowing us to identify each with one thread on the GPU and mak-ing the parallelization process very straightforward. Other parts of the algorithmare either not suitable for parallelization, or are not costly enough to warrant thisattention.203.2 Algorithmic improvementsIn addition to the code parallelization, we modified the way the nonlinear problem(2.18) is solved. Briefly, note that the physical definition of χ is as the excess ofthe pressure gradient over its limiting value, given byχ = G− τYH. (3.1)Thus, when χ > 0 the fluid is flowing locally, meaning that |∇aΨ∗+qn|> 0. Alsonote that the function χ(θx) increases monotonically, where θx = |∇aΨ∗+ qn|.With a simple change of variable, we can rewrite (2.18) not as an equation for θ ,but rather as an equation for χ , i.e.f (χ) = χ+ rHm+2κm(m+2)χm+1(χ+ τY/H)2[χ+(m+2)τY(m+1)H]−(x− τYH)= 0.(3.2)The function f (χ) is evidently monotonically increasing, and considering that x >τYH we also see that f (0)< 0. It is not hard to establish that (3.2) has a single zero,which needs to be found numerically. Note however that this eliminates one levelof nested iteration in defining χ . Having found the solution to (3.2), we can use(2.5) to find θx explicitly, and as |∇ψ∗+qn|= θx, then qn is found.There are many options that can be used to solve (3.2), e.g. bisection method,regula falsi, Ridder’s method. Although nonlinear, the properties of (2.5) are ana-lyzed in depth in [52], which makes it relatively easy to find bounds for the solutionand iterate reliably. We shall refer to using (3.2) in place of (2.18) as the optimizedalgorithm. For all results presented in the thesis the bisection method has beenused.3.3 Theoretical Peak Performance CalculationThe Theoretical Peak Performance is given by the multiplication of:(#Processors)x(Cores per processor)x(clock speed) x(FMAunit). In this regard, Figure 3.1 is a test of computational performance. The cal-culation in Figure 3.1 is the calculation of the rate given by (Number of floating21Table 3.1: Percentage of time by function#Processors Cores perproces-sorsclock speed (cy-cles per second)GHz(GPU) FMAunit = 2 (Fusemultiply-add)or FLOPS(CPU)Total(GFLOPS)GPU 68 64 7 2 60928CPU 4 8 3.9 2 249.6point calculations)/(Time of host/GPU)*(1e9), i.e. for every size of the matrix, thecode calculates the number of floating point calculations and the time that it takes todo the matrix-matrix multiplication. Notice that this time is a dependent variable,the independent variable is the size of the matrix. The CPU performance reachesa plateau before the size of the matrix 107 in about 75 GFLOPS, lower than the249.6 GFLOPS that is the theoretical peak. We observe that the GPU theoreticalpeak is very high in comparison with the maximum rate of less than 450 GFLOPS,this is for the last calculation or the largest size of the matrix. We anticipate that toreach a maximum in the calculation rate we should increase the size of the matrixto a much larger one and consistent with the increment in time. We did some ex-periments with larger matrices, but the theoretical peak is approximately 135 timesthe maximum rate of Figure 3.1. This means that depending on the time of thematrix-matrix multiplication, this calculation could be so time-consuming that itis impractical to attempt it. Having said this the purpose of Figure 3.1 was not toshow the peak of performance but rather to compare the performance between theCPU and the GPU given an intensive and common operation in scientific compu-tation. It is worth mentioning that the time to reach the peak is not infinite butimpractical to attempt. However, we would argue that it is a false discussion toreduce the analysis to a simple comparison between a CPU and a GPU and wewould like to elaborate on this matter. In the code that we uploaded one can noticethat there are four different versions of the code. The first one is the original one,which uses the resources of the CPU, the second one is the modified algorithmthat also uses just the resources of the CPU but the third and fourth ones are ver-sions of the code that were changed to run in parallel, being the last one, the one22Table 3.2: Percentage of time by functionAlgorithm CPU resources GPU resourcesOriginal ♢Original with extra algorithm ♢Parallelized algorithm ♢ ♢Parallelized algorithm with extra algorithm ♢ ♢that also includes an extra modification in the algorithm. As we describe abovethe algorithm itself was modified in some parts to use the resources of the GPU,although it still uses the resources of the CPU for the rest of the code. Hence, wecan assure that the third and four versions of the algorithm use the same and bothtogether the CPU and GPU. Therefore, we do not want to prove that the GPU alonehas a better performance but rather that those versions that uses the GPU and theCPU together are more efficient than those that only use the CPU. Lets emphasizehere that we are talking about the global purpose of the article in contrast with theparticular analysis that we did for Figure 3.1 where we want to specifically discussthe difference between one and the other. Furthermore, it might be that these twodifferent discussions lead to confusion or a false discussion, given the similarity ofthe topics, concepts, and definitions between them. For the sake of clarity, we addTable 3.2 that explains visually what we just argued.3.4 ScalabilityGiven the broad sense of the concept of scalability, we will address this topic usingdifferent definitions. The first definition of scalability is related to the size of theproblem and the time it takes to be executed, there are many possibilities to makethis calculation and we will present here some of them. When we want to analyzethe size of the problem, the number of processors that are used to calculate thealgorithm become relevant. The parallelization in a GPU uses different kinds ofmemories, like share memory, global memory, etc., and arithmetic processor unitsare set in different levels, like threads, blocks, warps, SM (Streaming Multiproces-sors), etc. Hence if we want to anticipate the scalability of our application, we canuse Amdahl’s law that relates the number of the processor, with the proportion of23Table 3.3: Amdahl’s law calculationTotal time of the code in serial 41394 s.Total time of the code to be parallelized 40115 s.% Time of the program running in serial 3.09 %The number of multiprocessors 4352S 32.13XTable 3.4: Percentage of time by functionMatrix size 10 100 1000 10000 100000 1000000Time (s) 0.36 0.363 0.367 0.37 0.373 0.375the code in parallel, and give us an approximation of the speed upS =1f + (1− f )p, (3.3)where f is part of the program running in serial, p is the number of processors.When we use a profile over our application, we obtain the results in Table 3.3.We equivalently can use the values of Table 3.5; however, the definitions fora GPU using this law has few studies and other laws related like Gustafson lawaccept different analysis and interpretations. In general, it is still an active topic ofresearch as we can see in [15].Another study and definition related to scalability it is a time measurement vs.the size of the problem. An example is the calculation of the solution of the linearalgebra problem Ax=b with different sizes of matrices. Figure 3.1 is a similar testbut here we are just taking the time directly, and a few cases for clarity.In practical terms, the time increase linearly tending to be constant. The matrixsize increases constantly by the order of magnitude. We anticipate this kind ofscaling for our algorithm. We observe this tendency in Figure 4.3 b) in the blueline that represents the parallel version and linear for the same case in Figure 4.4and Figure 4.5. This behavior in the last to curves can be explained given theincrement in the number of operations in these cases.Regarding the speed-up performance, as we described above, we changed dif-ferent parts of the code and implemented different functions to improve the speed-2410 1 10 2 10 3 10 4 10 5 10 6mesh size0.360.3620.3640.3660.3680.370.3720.3740.376time (s)Scalabilityup. We also reduced and optimized the code taking advantage of the original al-gorithm. All these strategies along with a detailed analysis and profiling of theoriginal code produced a notable improvement of the run time. The speed-up alsodepends on the model of the GPU, there are many examples in the literature e.g.in this link a company made a study on this matter: [43]. In this study they useda GeForce GTX 1070 with 1920 cuda cores in total, it claims that the speed-upfor this model was 200X. Note that the model that we use for this work has 4352Cuda cores. In [18] the speed-up was 58X. In [32] the speed-up was 45X. Thesefew examples show that it is relatively common to find simulations with 40X ofspeed up or above when they use a single GPU as we did. I do not include anyresearch with GPU cluster as we do not use it but for these cases,but in those cases,the speed-up tens to be is higher.3.5 Implementation and noveltyA deeper explanation of the process to use the GPU and parallelization for the codewill be given. One idea that is implicit in the thesis is that given the nature of theGPU coding, the change of the code in serial to the code in parallel using a GPU im-plies a change in the algorithm itself. Not all parallelization implies the same task.25There are some cases where the parallelization of a code just requires changing thenumber of cores that the code uses and it is very automatic. The programmer doesnot need to make major changes to the code. In this case, the parallelization usingthe GPU meant a massive change in the algorithm. We did not use libraries withprefabricated code. We implement new functions that use primitive kernels to useparallelization. In that sense, instead of using some of the solvers available in theCuda toolkit we developed our parallel algorithms for these important functions inthe original code, and we use fundamental concepts of parallel and concurrent pro-gramming. We use concepts of concurrency and parallelism at their fundamentallevels to modify the algorithm. In terms of Flynn’s taxonomy, we used multipleinstructions, and multiple data technics to allow our new algorithm to take advan-tage of the GPU, in that sense we use a different paradigm of programming forthis code. All the code suffered modifications but we will describe some particularlines of this new functions. These new algorithms and implementation in parallelamong other kernels (functions that run in the GPU) are the next:attributes(global) subroutine AdvanceLambdad,this function uses the thread and blocks (here just a basic definition not all themodification of the algorithm):x = (blockIdx%x-1) * TILE_DIM + threadIdx%x,y = (blockIdx%y-1) * TILE_DIM + threadIdx%y.Available from the CUDA cores to be able to calculate the variable related tothe calculation of lambda and hence fundamentally modifying the augmented la-grangian algorithm, this function was written as a global function that can be calledusing a function in the host (CPU) and also from the GPU. This modification usedthe heterogeneous parallel computing technique for the implementation.attributes(global) subroutine RHS(lambdan_d,fij_d,pn_d),this subroutine calculates the right-hand side of the linear algebra problem as weassemble the Finite Volume method for solving the advection equation. Using thesame parallelization techniques as above all the assembly happens in parallel.attributes(device) function hgp_d(x),26hgp is a small function that calculates the annular gap half-width and uses theeccentricity of the annulus, this function is used intensively. We changed its im-plementation using the concept of pure function on the GPU (using the deviceattribute), which means that the result of the function does not need to go back-and-forth between the CPU and the GPU during the use of the function. The calculationhappens in the GPU and stays there saving the time to transfer the information.attributes(global) subroutine extractuvwParallel(deltat_t),this function modifies the Flux Correct Transport algorithm using the threads andblocks to calculate in parallel all the velocities defined, across all the domain inthe algorithm. It also uses the concurrency function syncthread() that assures allthe calculation are correct. Similarly, as when one uses share memory in the im-plementation, the memory is allocated per thread block and it is shared betweenthe threads. This speeds up the calculation. In the literature, we can find that thiscapability gives a memory latency of 100X lower than uncached memory latency.However, this can create incorrect calculations or deadlocks if the threads are notsynchronized. To solve this problem, we used of syncthreads() function to avoidsrace condition. A race condition means that the thread runs logically in parallelbut does not execute physically at the same time. Lets say for example that wehave thread A and thread B and lets assume A and B are not in the same warp, ifB has not finished writing its result when A tries to read it, we have a race condi-tion. Given the complexity of this global subroutine, the syncthreads() assure us toobtain correct results. Here is the line where we use this technic in the algorithm:max(dabs(v_d(x+1,y+1)), maxv),call syncthreads().3.6 ResultsBefore presenting results of representative cases, we give an overview of the rel-ative computational costs. We have measured the runtime of the entire code andprofiled runtimes of the different parts of the code. Apart from the initializationsteps in the code, the main subroutines that are called on each time step are asfollows:27Table 3.5: Percentage of time by functionFunction time(sec.) on 124x256 mesh %Advancep 376.50 97.2Advancelambda 9.41 2.42AdvanceOneTimeStep 1.34 0.34Advanceu 0.13158 0.035• advanceonetimestep: this calculates the concentration of each fluidinside the annulus, advancing one time step using an explicit (flux correctedtransport) method.• advanceu: this calculates a new value of Φ from the saddle point problem.• advancep: this calculates a new value of q from the saddle point problem.• advancelambda: this calculates a new value of λ .As we have noted in the prior discussion, advanceonetimestep, advancepand advancelambda represent operations that naturally scale linearly with thesize of the mesh. In general, advanceu involves solving a linear Poisson problemon the mesh. However, the Poisson problem remains unchanged for successivetime steps. Thus, here we have performed a QR decomposition at the start ofthe computation and used this repeatedly on each time step to solve the Poissonproblem. This ensures that advanceu will scale linearly with the size of themesh. This method is suitable for the moderate sized meshes that we use, but someform of iterative method would be needed for very large problems.Table 3.4 shows the percentage of computational time taken by each functionfor a typical computation using yield stress fluids and with a mesh size 124x256.We see that advancep represents 97% of the runtime. Although each functionshould scale linearly with time, the nested nonlinear iteration performed on eachmesh cell is the dominant cost. This explains our motivation for targeting our GPUimplementation for this operation.280 0.5 1100.0200.0300.0400.0500.0600.0700.0800.0900.01000.00 0.5 1 0 0.5 1 0 0.5 100.10.20.30.40.50.60.70.80.9ConcentrationFigure 3.2: Case A concentration colormap, computed in parallel. Sequenceof times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59 s, from left to right.The white lines show the streamlines φ = 0, 0.1, 0.2, ..., 1.3.7 Case studiesWe have developed and analyzed six different cases to illustrate typical rheologicalparameters and operational variables. Our cases are greatly simplified comparedto the ranges found in operations, by having only 2 fluids present and in that theannulus is uniform. The fluid properties are specified in table 3.6. All cases con-sider inner and outer radii: rˆi = 0.15m and rˆo = 0.175m, and a length lˆ = 1000m.Each annulus has a constant standoff (eccentricity). The annuli are all vertical,which would relate to a surface or intermediate casing operation. For simplicitythe flow rate is fixed at Qˆ = 0.002525 m3/s. Even though simplified, each case isdescribed by 15 parameters in the dimensional setting, which can be reduced to 11dimensionless parameters by scaling the equations and variables. The dimensional290 0.5 1100.0200.0300.0400.0500.0600.0700.0800.0900.01000.00 0.5 1 0 0.5 1 0 0.5 100.10.20.30.40.50.60.70.80.9ConcentrationFigure 3.3: Case B concentration colormap, computed in parallel. Sequenceof times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59 s, from left to right.The white lines show the streamlines φ = 0, 0.1, 0.2, ..., 1.parameters are described and discussed in [34].In table table 3.7 we listed the dimensionless parameters used in our experi-ments. Note that ρ2 = 1 in all the cases due to scaling, and the scaled length of thewell is Z = 1958.86. We use Nξ = 1024 meshpoints in the axial direction through-out. Cases A-C use Nφ = 64 and the remainder use Nφ = 32. All computationsstart with the annulus completely full with displaced fluid (colored yellow below),and with the interface at the bottom (ξ = 0) of the annulus.Cases A & B are Newtonian. The densities are the same for both fluids. Inprinciple the augmented-Lagrangian procedure is not needed for these linear cases.These are included for bench-marking purposes. In Figure 3.2 we show resultsfrom case A at 4 successive times through the displacement flow. The streamlines(white) are plotted over the colormap of the fluid concentrations. For this base300 0.5 1100.0200.0300.0400.0500.0600.0700.0800.0900.01000.00 0.5 1 0 0.5 1 0 0.5 100.10.20.30.40.50.60.70.80.9ConcentrationFigure 3.4: Case C concentration colormap, computed in parallel. Sequenceof times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59 s, from left to right.The white lines show the streamlines φ = 0, 0.1, 0.2, ..., 1.Table 3.6: Dimensional fluid properties: rˆi = 0.15m, rˆo = 0.175m, lˆ =1000m, β = 0, Qˆ = 0.002525 m3/s. Fluid 1 (displaced fluid), Fluid 2(displacing fluid).Case e κˆ1=1.0(Pa.sn)κˆ2=1.0(Pa.sn)n1 n2 τˆY,1(Pa)τˆY,2(Pa)ρˆ1(kg/m3)ρˆ2(kg/m3)A 0.2 1.0 1.0 1.0 1.0 0 0 1200 1200B 0.5 1.0 1.0 1.0 1.0 0 0 1200 1200C 0.5 1.0 1.0 0.6 0.5 0 0 1200 1200D 0.5 0.02 0.01 0.7 1.0 7.05 2.35 1440 1800E 0.5 0.02 0.03 0.7 1.0 4.79 2.35 1440 1800F 0.4 1.0 1.0 1.0 0.6 0.3 1.0 1200 1300G 0.4 1.0 1.0 0.7 1.0 4.79 2.35 1440 1440310 0.5 1100.0200.0300.0400.0500.0600.0700.0800.0900.01000.00 0.5 1 0 0.5 1 0 0.5 100.10.20.30.40.50.60.70.80.9ConcentrationFigure 3.5: Case D concentration colormap, computed in parallel. Sequenceof times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59 s, from left to right.The white lines show the streamlines φ = 0, 0.1, 0.2, ..., 1.Table 3.7: Dimensionless parameters: ρ2 = 1, Z = 1958.86. Fluid 1 (dis-placed fluid), Fluid 2 (displacing fluid).Case St∗ e κ1 κ2 n1 n2 τY,1 τY,2 ρ1A 0.233 0.2 1.0 1.0 1.0 1.0 0.0 0.0 1.0B 0.233 0.5 1.0 1.0 1.0 1.0 0.0 0.0 1.0C 0.135 0.5 1.0 0.8123 0.6 0.5 0.0 0.0 1.0D 0.143 0.5 0.0120 0.0336 0.7 1.0 0.9880 0.3290 0.8E 0.143 0.5 0.0176 0.0492 0.7 1.0 0.6705 0.3290 0.8F 0.0806 0.4 0.0773 0.0336 1.0 0.6 0.2899 0.9664 0.9G 0.2421 0.3 0.4142 0.7729 0.7 1.0 0.4628 0.2421 1.0320 0.5 1100.0200.0300.0400.0500.0600.0700.0800.0900.01000.00 0.5 1 0 0.5 1 0 0.5 100.10.20.30.40.50.60.70.80.9ConcentrationFigure 3.6: Case E concentration colormap, computed in parallel. Sequenceof times: 2042.14 s, 5105.26 s, 7147.33 s, 1010.59 s, from left to right.The white lines show the streamlines φ = 0, 0.1, 0.2, ..., 1.case, the only effect is due to the eccentricity (e = 0.2), which causes the fluid toflow faster in the wider part of the annulus. From an initially horizontal interfacewe observe that the interface elongates along the wide side of the annulus andexits. Case B is a similar flow but with larger eccentricity (e= 0.5); see Figure 3.3.The behaviour is qualitatively the same as in Figure 3.2, but as we increase theeccentricity the fluid in the wider part of the annulus increases and that in thenarrower part decelerates, so that the interface elongates in each frame.Note that for these two cases there is no density difference. The 2 fluids rhe-ologies are identical, and we are viewing a pure dispersion problem. The degreeof spreading of the concentration profile at the interface is entirely indicative ofnumerical diffusion/dispersion in solving the concentration equation (2.9). Unlikein more complex situations, the streamlines are entirely parallel and there is no sec-33Figure 3.7: Case G, computed in parallel. Sequence of time: 2042.14 s.,5105.26 s., 7147.33 s., 1010.59 s. from top to bottom. The white linesshow the streamlines φ = 0, 0.1, 0.2, ..., 1.34ondary flow close to the interface. It can be observed that there is relatively littlespreading of the interface with the current algorithm. Note that since the fluids up-stream and downstream are identical, the streamlines should remain parallel acrossthe interface.Case C represents the first displacement flow; see Figure 3.4. The fluids arenow shear-thinning (n1 = 0.5 and n2 = 0.6) and the eccentricity is still e = 0.5.Recall that the flow rate imposed is identical for all cases, to allow comparison.Due to the shear-thinning effect both fluids flow faster on the wide side, relativeto the Newtonian case B, i.e. faster flow shears the fluid more, reducing its effec-tive viscosity. We observe this in that the streamlines are now denser on the wideside than in case B. The consequence of the faster velocities is that the interfacesadvances faster on the wide side, making the displacement less effective than thatof case B. Note too that fluid 2 has smaller n and hence is effectively less viscousthan fluid 1, which may contribute to the finger-like phenomenon observed. Thestreamlines now distort as they cross the displacement front, only becoming par-allel upstream and downstream of the front. We see that the wide side interfacedevelops a shock-like profile.A different non-Newtonian phenomenon is introduced in case D, where bothfluids are shear thinning, n1 < n2 < 1, and both have a yield stress, that is largerin the displaced fluid, as is often the case for drilling mud and cement slurry; seeFigure 3.5. The eccentricity is still e = 0.5, but there is a significant density dif-ference, ρ2 > ρ1, with the density of the displacing fluid being approximately thatof cement. Similar to the previous example, the shear-thinning of the fluids leadsto a faster flow in the wide side of the annulus, but this effect is countered by thebuoyancy force coming from the density difference. The result is a more efficientdisplacement. The distortion around the interface is more visible, and we see sig-nificant dispersion ahead of the front close to φ = 0, the wide side of the annulus.It’s been shown that the numerical method is able to preserve the sharpness of thefront, so the dispersion is due to secondary flows that emerge. On the wide side,near φ = 0, there is flow ahead of the front and on the narrow side, around φ = 0,the secondary flow is behind the front. This secondary flow, which acts to flattenthe interface, is due to buoyancy.Note that ahead of the front, in the displaced fluid, in the narrowest part of the35annulus, there are no streamlines, |∇aΨ| ∼ 0; see (2.5). Interestingly, far behind thefront, when fully in fluid 2, there is also static fluid on the narrow side. However,buoyancy induces weak secondary flows on the narrow side, sufficient to movefluid 1 away from φ = 1 in fluid 1 and towards φ = 1 in fluid 2. This allows thefront to move along the narrow side, albeit slowly. It is interesting to compare withcase C, where there is also poor narrow side displacement. The slow movement ofthe front on the narrow side is comparable between cases C and D. However, caseC has no density difference and no yield stress, thus the fluids are always movingon the narrow side, and not just close to the front. Evidently, both cases C and D arepoor displacements. In the preceding figures the final time shown corresponds toapproximately 1 annular volume pumped. During a cementing operation, an excessof ∼20% annular volume may be pumped, but continuous pumping of cement issimply not possible, due to environmental waste disposal and cost concerns. Thus,displacement simulations such as cases C and D show wells that will have verypoor cement coverage over significant lengths, likely leading to leakage and givingconcern for corrosion, i.e. the cement does not protect the narrow side of the casing.Case E is similar to case D, but with a reduced yield stress in fluid 1 and anincreased consistency in fluid 2. Both effects lead to a better displacement, albeitstill poor, as we see that there remains a static channel along the narrow side ofthe annulus; see Figure 3.6. The residual channel of fluid 1 is narrower than incase D. In order to displace a yield stress fluid fully, one needs a combination offeatures: reduced yield stress, increased viscosity on fluid 2, significant buoyancyand reduced eccentricity. Finding the parametric ranges in which displacementscan be effective is a key engineering challenge, more so when the well geometriesare not as simple as those here. Case F illustrates a situation where a the yieldstress fluid 1 is fully displaced.36Chapter 4ComparisonsCases A and B involve only passive advection, i.e. fluids 1 and 2 are identical. Thissituation admits an analytical solution for the streamfunction and hence velocity;see [53]. The computed concentration may then be compared with the advectedinterface. For identical Newtonian fluids the velocity depends only on φ and isgiven by:w¯(φ) =H2(φ)∫ 10 H3(φ) dφ(4.1)with H(φ) = 1+ ecosπφ . This integral is readily evaluated and the (analytical)position of the interface is simply advected from an initially flat initial condition.The comparison is made in Figure 4.1 for case A.Figure 4.2 shows comparisons of computed results using both the non-optimizedserial and optimized parallel codes, for cases E and F. The results are very close.Displacement flows are advective, hence small differences in solution will flowdownstream and can be amplified in doing so. Such differences in numerical solu-tion remain small in cases where there is some underlying stability to the solution,as here. However, this might not always hold as many cementing displacementflows are unstable.To give a quantitative comparison of the difference in solutions, we calculatedthe L2 norm between the results from the non-optimized serial code and optimizedparallel code, for meshes: 8x8, 8x16, 16x16, 16x32, 32x32, 32x64 64x64 and64x128. We chose these specific size of meshes because the grids in the parallel37Figure 4.1: Case A: a) computed solution; b) analytical advected interface;c) comparison.version of the code should be a power of 2, to match the number of threads in theGPU that run the functions in parallel.As we can verify in Figure 4.3(a) and Figure 4.4(a), the difference is less than10−13 for all the meshes, for the Newtonian fluid computations. For the non-Newtonian cases this increases to ≲ 10−11, as shown in Figure 4.5(a) for caseE. Although the algorithm is the same for Newtonian and non-Newtonian cases,in general the non-Newtonian flows require more iteration, which allows for largererror in convergence. As commented with regard to Figure 4.2, small numericaldifferences are transported with the displacement flow and hence the discrepanciesbetween methods will grow with time in most simulations, and more so in responseto the flow complexity.For the same 3 cases, A, B and E, Figure 4.3(b), Figure 4.4(b) and Figure 4.5(b)show the computing time taken for various meshes. The red line denotes thenon-optimized serial code and the blue line is the optimized parallel code. Bothcodes runtime increase with the meshsize. However, the decrease in computingtime using the optimised parallel code appears to be more significant for the non-38(a)0 0.5 1.0100.0200.0300.0400.0500.0600.0700.0800.0900.01000.0(b)0 0.5 1.0(c)0 0.5 1.0(d)0 0.5 1.000.10.20.30.40.50.60.70.80.9ConcentrationFigure 4.2: Case E in optimized parallel at t = 4333; (b) Case E in nonopti-mized serial at t = 4333; (c) Case F in optimized parallel at t = 4982; (d)Case F in non-optimized serial at t = 4982.Figure 4.3: Case A: a) L2 error between optimized parallel and non-optimized serial codes; b) annular timing benchmark; red line: non-optimized serial code; blue line: optimized parallel code.39Figure 4.4: Case B: a) L2 error between optimized parallel and non-optimizedserial codes; b) annular timing benchmark; red line: non-optimized se-rial code; blue line: optimized parallel code.Figure 4.5: Case E: a) L2 error between optimized parallel and non-optimizedserial codes; b) annular timing benchmark; red line: non-optimized se-rial code; blue line: optimized parallel code.Newtonian simulations compared to the Newtonian simulations, as we might ex-pect.In table 4.1 we show the time for all the cases ran in serial and in parallel with a32X1024 mesh. In table 4.2 we show the time of the optimized algorithm in serialand the optimized algorithm in parallel. We notice the time reduction for all casesfor the parallel version. We can also see that the time of the optimized algorithm islower than the original algorithm. The best results are obtained with the code with40the optimized algorithm and in parallel.Table 4.1: Benchmark for cases: in parallel and in serialBenchmark mesh32×1024Serial ParallelCase A 150 min. 55 s 86m 29s.Case B 189 min. 61 s. 70m 8 s.Case C 1197 min. 46 s. 78m 49s.Case D 2160 min. 20 s. 162m 25s.Case E 5755 min. 30 s. 177 min 1 s.Case F 6324 min. 52 s. 101 min. 6sCase G 1194 min. 55 s. 131 min 40 s.Table 4.2: Benchmark for cases: original algorithm and optimised algorithmBenchmarktime mesh32×1024Optimized algorithmserialOptimized algorithmparallelCase A 140 min. 3 s. 46 min. 43 s.Case B 139 min 5 s. 46 min. 45 s.Case C 140 min 29 s. 61 min. 54 s.Case D 589 min. 36 s. 67 min. 17 s.Case E 593 min. 5 s. 71 min. 58 s.Case F 360 min. 5 s. 69 min. 58 s.Case G 274 min. 25 s. 80 min. 50 sTo quantify the improvement we compared the average times of these cases,over the 3 different improvement strategies. We calculated the percentage of thetime take with respect to the original serial algorithm. The time in parallel is≈ 13%of the time in serial. The optimized algorithm in serial takes≈ 4.45% of the time ofthe original code. The optimized algorithm in parallel takes ≈ 2.5% of the originalcode. As we discussed before, the specific improvements will also depend on themesh size and rheological parameters.Noticeable here is that the very long runtimes (cases E & F, Non-Newtonianfluids) improved the most. As we know, Non-Newtonian simulations require morecomputational resources. Hence, the optimization becomes more relevant and valu-able in the light of the real-time simulations. Running 10-15 simulations of 1-hour41duration to enhance the design of a well is feasible, but not when each takes 100hours. The most important point is that the improvements bring the simulationsinto the range of a single desktop machine with a good GPU, as is available topracticing engineers.42Chapter 5ConclusionsIn this study, we implemented a faster version of the annular displacement flowssimulations first presented in [50]. We parallelized the original code using a GPU,and modified the algorithm where it had the most expensive computational cost.This resulted in an average reduction of the computing time to around 2.5% ofthe original cost. The cost reduction appears to vary significantly with the fluidrheologies: perhaps unsurprising as these change the nature of the non-linearities.The study serves as a example of a class of fluid flow problems that may beamenable to significant GPU acceleration. It is fairly commonplace to use scal-ing arguments to reduce 3D Navier-Stokes problems to simpler situations wherethe flow occurs in a thin film, slowly varying channel or similar geometry. Suchmodels are amenable to analysis and often are fast enough to compute, such thatthey may serve as a digital twin for process design or control. In these modelsthere is typically a reduced direction in which the equations are integrated analyt-ically (across a boundary layer, or thin film, etc). In dealing with complex fluids,these reduced models are not always analytically tractable and their computationcan quickly become the main cost. Here we have show that both GPU and smartconsideration of the reduced models, can both lead to significant savings in com-putational cost.Among the limitation of this work we can find the hardware that we used torun all the simulations. We used a specific model of GPU during all the timethat we run the code and we kept using the same one to be consistent throughout43all our analysis; However, we can upgrade the GPU and the CPU and we willgain velocity. The code in parallel does not need to be modify to run in a newGPU and we will improve its performance. Another limitation is that we did notparalellize all the functions of the code, just those ones that we identify as theones that slowed down the time of simulation the most. We can parallelized thosesecondary functions and obtain a better speed-up to the simulation. This last taskcould be part of a future work. Also we can do more benchmarks between all thedifferent versions of the code. As we described in the results, we obtained fourdifferent versions of the code with different algorithms or level of parallelization.At the end of our work we released all our code with a public license as a big pieceof software. Hence, detail comparisons of errors and benchmarks e.g. per case andper different meshes could be very useful to improve and upgrade this software.44Bibliography[1] M. Allouche, I. Frigaard, and G. Sona. Static wall layers in the displacementof two visco-plastic fluids in a plane channel. Journal of Fluid Mechanics,424:243–277, 2000. → pages 2[2] N. J. Balmforth, I. A. Frigaard, and G. Ovarlez. Yielding to stress: recentdevelopments in viscoplastic fluid mechanics. Annual Review of FluidMechanics, 46:121–146, 2014. → pages 1[3] S. Bittleston, J. Ferguson, and I. Frigaard. Mud removal and cementplacement during primary cementing of an oil well–laminar non-newtoniandisplacements in an eccentric annular hele-shaw cell. Journal ofEngineering Mathematics, 43(2-4):229–253, 2002. → pages 2, 9, 11, 13[4] F. J. Bittleston SH and F. IA. Mud removal and cement placement duringprimary cementing of an oil well- Laminar non-Newtonian displacements inan eccentric annular Hele-Shaw cell., volume 43. 2002. → pages 3[5] M. Bogaerts, C. Azwar, M. Bellabarba, M. Dooply, A. Salehpour, et al.Wellbore cementing: an integral part of well integrity. In OffshoreTechnology Conference. Offshore Technology Conference, 2015. → pages 2[6] B. Bradie. A friendly introduction to numerical analysis. Pearson EducationIndia, 2006. → pages 55[7] J. W. Brice Jr, B. Holmes, et al. Engineered casing cementing programsusing turbulent flow techniques. Journal of Petroleum Technology, 16(05):503–508, 1964. → pages 2[8] R. Burden and J. Faires. Numerical differentiation & integration numericaldifferentiation i. Numerical analysis, pages 174–184, 2011. → pages 55[9] M. Carrasco-Teja and I. Frigaard. Displacement flows in horizontal, narrow,eccentric annuli with a moving inner cylinder. Physics of Fluids, 21(7):073102, 2009. → pages 1045[10] M. Carrasco-Teja and I. Frigaard. Non-newtonian fluid displacements inhorizontal narrow eccentric annuli: effects of slow motion of the innercylinder. Journal of Fluid Mechanics, 653:137–173, 2010. → pages 3, 10[11] M. Carrasco-Teja, I. Frigaard, B. Seymour, and S. Storey. Viscoplastic fluiddisplacements in horizontal narrow eccentric annuli: stratification andtravelling wave solutions. Journal of Fluid Mechanics, 605:293, 2008. →pages 2, 9, 10[12] S. Y. Chung, G. H. Rhee, and H. J. Sung. Direct numerical simulation ofturbulent concentric annular pipe flow: Part 1: Flow field. InternationalJournal of Heat and Fluid Flow, 23(4):426–440, 2002. → pages 4[13] S. Cook. CUDA programming: a developer’s guide to parallel computingwith GPUs. Newnes, 2012. → pages 7[14] M. Couturler, D. Guillot, H. Hendriks, F. Callet, et al. Design rules andassociated spacer properties for optimal mud removal in eccentric annuli. InCIM/SPE International Technical Meeting. Society of Petroleum Engineers,1990. → pages 2[15] F. De´vai. The refutation of amdahl’s law and its variants. In Transactions onComputational Science XXXIII, pages 79–96. Springer, 2018. → pages 24[16] Y. Dimakopoulos, G. Makrigiorgos, G. C. Georgiou, and J. Tsamopoulos.The PAL (penalized augmented lagrangian) method for computingviscoplastic flows: A new fast converging scheme. Journal ofNon-Newtonian Fluid Mechanics, 256:23–41, 2018. → pages 5[17] M. B. Dusseault, R. E. Jackson, and D. Macdonald. Towards a road map formitigating the rates and occurrences of long-term wellbore leakage.University of Waterloo Waterloo, ON, Canada, 2014. → pages 2[18] P. Eastman and V. S. Pande. Efficient nonbonded interactions for moleculardynamics on a graphics processing unit. Journal of computationalchemistry, 31(6):1268–1272, 2010. → pages 25[19] R. Farber. CUDA application design and development. Elsevier, 2011. →pages 7[20] A. Fitt and C. Please. Asymptotic analysis of the flow of shear-thinningfoodstuffs in annular scraped heat exchangers. Journal of EngineeringMathematics, 39(1):345–366, 2001. → pages 646[21] A. C. Fowler, A. C. Fowler, and A. Fowler. Mathematical models in theapplied sciences, volume 17. Cambridge University Press, 1997. → pages 4[22] J. B. Freund, S. K. Lele, and P. Moin. Compressibility effects in a turbulentannular mixing layer. part 1. turbulence and growth rate. Journal of FluidMechanics, 421:229–267, 2000. → pages 4[23] I. Frigaard. Simple yield stress fluids. Current Opinion in Colloid &Interface Science, 43:80–93, 2019. → pages 1[24] I. Frigaard, S. Leimgruber, and O. Scherzer. Variational methods andmaximal residual wall layers. Journal of Fluid Mechanics, 483:37–65, 2003.→ pages 2[25] R. Glowinski and A. Wachs. On the numerical simulation of viscoplasticfluid flow. In Handbook of Numerical Analysis, volume 16, pages 483–717.2011. → pages 5[26] A. Gregatti, G. Carvalho, T. Piedade, G. Campos, et al. Cementing in frontof saline zone and carbonate reservoir with presence of co 2. In OTC Brasil.Offshore Technology Conference, 2015. → pages 2[27] N. Guillot D. Well Cementing. Schlumberger Educational Services, SugarLand, 2006. → pages 1, 2[28] A. IK1* and G. S2. Local convergence for a regula falsi-type method underweak convergence. Journal of Applied & Computational Mathematics, 4:3,2015. → pages 55[29] Ishwaranand. Fortran - bisection method, 2020. URLhttps://www.ishwaranand.com/2020/08/bisection-method.html. → pages 52[30] Ishwaranand. Fortran - regula falsi method, 2020. URLhttps://www.ishwaranand.com/2020/08/regula-falsi-method.html. → pages52[31] A. Jamot. De´placement de la boue par le latier de ciment dans l’espaceannulaire tubage-paroi d’un puits. Revue Assoc. Franc. Techn. Petr, 224(27-37):352, 1974. → pages 2[32] B. Jiang, G. A. Woodell, and D. J. Jobson. Novel multi-scale retinex withcolor restoration on graphics processing unit. Journal of Real-Time ImageProcessing, 10:239–253, 2015. → pages 2547[33] C. F. Lockyear, D. F. Ryan, M. M. Gunningham, et al. Cement channeling:how to predict and prevent. SPE Drilling Engineering, 5(03):201–208,1990. → pages 2[34] A. Maleki and I. Frigaard. Primary cementing of oil and gas wells inturbulent and mixed regimes. Journal of Engineering Mathematics, 107(1):201–230, 2017. → pages 1, 3, 10, 15, 30[35] A. Maleki and I. Frigaard. Rapid classification of primary cementing flows.Chemical Engineering Science, page 115506, 2020. → pages[36] A. Maleki and I. A. Frigaard. Turbulent displacement flows in primarycementing of oil and gas wells. Physics of Fluids, 30(12):123101, 2018. →pages 10[37] A. Maleki Zamenjani. Annular displacement flows in turbulent and mixedflow regimes. PhD thesis, University of British Columbia, 2019. → pages 1[38] R. McLean, C. Manry, W. Whitaker, et al. Displacement mechanics inprimary cementing. Journal of Petroleum Technology, 19(02):251–260,1967. → pages 2[39] F. I. Moyers-Gonzalez MA. Kinematic instabilities in two-layer eccentricannular flows, part 1: Newtonian fluids. J Eng Math, (62):103–131, 2008.→ pages 2[40] F. I. Moyers-Gonzalez MA. Kinematic instabilities in two-layer eccentricannular flows, part 2: Newtonian fluids. J Eng Math, (65):25–52, 2009. →pages 2[41] E. Nelson. Well cementing. Schlumberger Educational Services, 1990. →pages 6[42] N. Nikitin, H. Wang, and S. Chernyshenko. Turbulent flow and heat transferin eccentric annulus. Journal of Fluid Mechanics, 638:95–116, 2009. →pages 4[43] D. A. A. NIKOLA VASILEV. Speed up of numeric calculations using agraphics processing unit (gpu), 2020. URL chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/http://www.eurorisksystems.com/documents/speed up of numeric calculations using GPU.pdf. → pages 25[44] J. Nouri and J. Whitelaw. Flow of newtonian and non-newtonian fluids in aneccentric annulus with rotation of the inner cylinder. International Journalof Heat and Fluid Flow, 18(2):236–246, 1997. → pages 448[45] C. NVidia. Cuda c programming guide, version 9. 1. NVIDIA Corp, 2018.→ pages 7[46] C. Nvidia. Nvidia cuda installation guide for linux. Nvidia Corporation,(10.1):43, 2019. → pages 7[47] C. Nvidia. Cusparse library. Nvidia Corporation, (10.1):358, 2019. →pages 7[48] C. NVIDIA®. NVIDIA® CUDA C, programming guide. NVIDIA®Corporation, 120(8):265, 2018. → pages viii, 7, 8[49] V. Osayande, I. Isgenderov, M. Bogaerts, I. Kurawle, P. Florez, et al. Beatingthe odds with channeling in erd well cementing: solution to persistentchanneling problems. In SPE Russian Oil and Gas Exploration &Production Technical Conference and Exhibition. Society of PetroleumEngineers, 2014. → pages 2[50] S. Pelipenko and I. Frigaard. Mud removal and cement placement duringprimary cementing of an oil well–part 2; steady-state displacements. journalof Engineering Mathematics, 48(1):1–26, 2004. → pages 43[51] S. Pelipenko and I. Frigaard. On steady state displacements in primarycementing of an oil well. J. Eng. Math, 46:1–26, 2004. → pages 2, 9, 12, 16[52] S. Pelipenko and I. Frigaard. Two-dimensional computational simulation ofeccentric annular cementing displacements. IMA Journal of AppliedMathematics, 69(6):557–583, 2004. → pages 2, 9, 21[53] S. Pelipenko and I. Frigaard. Visco-plastic fluid displacements innear-vertical narrow eccentric annuli: prediction of travelling-wave solutionsand interfacial instability. Journal of Fluid Mechanics, 520:343–377, 2004.→ pages 2, 9, 37[54] A. Renteria, A. Maleki, I. Frigaard, B. Lund, A. Taghipour, and J.-D.Ytrehus. Displacement efficiency for primary cementing of washout sectionsin highly deviated wells. In SPE Asia Pacific Oil and Gas Conference andExhibition, page 137. Society of Petroleum Engineers, 2018. → pages 10[55] C. Ridders. A new algorithm for computing a single root of a realcontinuous function. IEEE Transactions on Circuits and Systems, 26(11):979–980, 1979. doi:10.1109/TCS.1979.1084580. → pages 5249[56] J. Sanders and E. Kandrot. CUDA by example: an introduction togeneral-purpose GPU programming. Addison-Wesley Professional, 2010.→ pages 7[57] P. Saramito. A damped Newton algorithm for computing viscoplastic fluidflows. Journal of Non-Newtonian fluid mechanics, 238:6–15, 2016. →pages 5[58] P. Saramito and A. Wachs. Progress in numerical simulation of yield stressfluid flows. Rheologica Acta, 56(3):211–230, 2017. → pages 5[59] R. Sawko. Mathematical and computational methods of non-newtonian,multiphase flows. 2012. → pages 4[60] C. Schoof and I. Hewitt. Ice-sheet dynamics. Annual Review of FluidMechanics, 45:217–239, 2013. → pages 6[61] H. Su, Y. Zheng, J. Qi, F. Gao, Y. Guo, S. Taoutaou, R. Wang, Y. An,H. Guo, et al. Cementing optimization through an enhanced ultrasonicimaging tool. In SPE/IATMI Asia Pacific Oil & Gas Conference andExhibition. Society of Petroleum Engineers, 2015. → pages 2[62] S. Taghavi and I. Frigaard. Estimation of mixing volumes in buoyantmiscible displacement flows along near-horizontal pipes. The CanadianJournal of Chemical Engineering, 91(3):399–412, 2013. → pages 4[63] S. W.-B. M. F. I. Taghavi SM, AlbaK. Miscible displacement flows innear-horizontal ducts at low atwood number. J Fluid Mech 696, (696):175–214, 2012. → pages 2[64] P. Tardy and S. Bittleston. A model for annular displacements of wellborecompletion fluids involving casing movement. Journal of Petroleum Scienceand Engineering, 126:105–123, 2015. → pages 3[65] G. I. Taylor. The dispersion of matter in turbulent flow through a pipe.Proceedings of the Royal Society of London. Series A. Mathematical andPhysical Sciences, 223(1155):446–468, 1954. → pages 4[66] L. N. Trefethen and D. Bau III. Numerical linear algebra, volume 50. Siam,1997. → pages 54[67] R. Tre´molie`res, J.-L. Lions, and R. Glowinski. Elsevier, 1976. → pages 550[68] T. Treskatis, A. Roustaei, I. Frigaard, and A. Wachs. Practical guidelines forfast, efficient and robust simulations of yield-stress flows withoutregularisation: A study of accelerated proximal gradient and augmentedLagrangian methods. Journal of Non-Newtonian Fluid Mechanics, 262:149–164, 2018. → pages 5[69] T. Watson et al. Surface casing vent flow repair-a process. In CanadianInternational Petroleum Conference. Petroleum Society of Canada, 2004. →pages 2[70] T. L. Watson, S. Bachu, et al. Identification of wells with high co2-leakagepotential in mature oil fields developed for co2-enhanced oil recovery. InSPE Symposium on Improved Oil Recovery. Society of Petroleum Engineers,2008. → pages 2[71] T. L. Watson, S. Bachu, et al. Evaluation of the potential for gas and co2leakage along wellbores. SPE Drilling & Completion, 24(01):115–126,2009. → pages 2[72] H. P. William. Numerical recipes. The art of scientific computing, 745,2007. → pages 55[73] M. Zare, A. Roustaei, and I. Frigaard. Buoyancy effects on micro-annulusformation: Density stable displacement of newtonian–bingham fluids.Journal of Non-Newtonian Fluid Mechanics, 247:22–40, 2017. → pages 2[74] J. Zhang and I. Frigaard. Dispersion effects in the miscible displacement oftwo fluids in a duct of large aspect ratio. Journal of Fluid Mechanics, 549:225–251, 2006. → pages 451AppendixA.1 Optimized algorithm pseudocodeIn this section, the modification of the algorithm is described in more detail. Thetheoretical deduction can be seen in 3.2. We obtained the final equation in terms ofχ given by section 3.2. Because this is a polynomial equation of χ , we need to usea root finder to obtain the final solution. We can observe that even before using theroot finder, we have significantly optimized the algorithm, avoiding all the implicitcalculations that were previously needed. The particular implementation of thealgorithm is shown step by step in Algorithm 1We used three root finders: Ridder’s, regula falsi and bisection methods ([55],[29], [30]). The structure of these three algorithms is different, but we programthem in a way that all the steps and implementation are the same except for step 8.In, Algorithm 1 we describe the algorithm from the definition of the initial con-ditions until the result is found it. In line 2, we define the function f (ξ ) recallingsection 3.2. In line 4, we define the endpoints, tolerance, and maximum iterationsexactly as the original algorithm. In line 5, it starts the new algorithm, and we cansee the first step of the root finder, in this case, the bisection method. In line 5, weassure that the endpoints meet the intermediate value theorem conditions. Fromsteps 6 to 16, we enter the bisection method loop using an Nmax as a criterion tostop the loop and finish the algorithm. In lines 13 and 15, we obtain the last valueof c, which is the numerical approximation of the root. This value is passed as52the solution of, eq. (2.18) used in the original algorithm. The other two variationsof the codes have the same steps and implementations, except step 8. The onlydifference between these two algorithms is the calculation of the point c that rep-resents the root approximation. In Algorithm 2 and Algorithm 3, we show themodification of this step for regula falsi and Ridder’s correspondingly.53A.2 Convergence AnalysisLastly, we calculate the theoretical speed of the three implementations. Eventhough the differences between them could be slight, if we take into account thatthis part of the code is used intensively, when this small-time is added hundreds orthousands of times, it can make a considerable difference. Therefore, we do thefollowing numerical analysis to decide which algorithm could lead us to obtain thebest improvement. We will use the concept of the floating-point operation as a firstapproximation and then the definition of rate of convergence to make a more the-oretical approach. Using the definition in [66] for floating-point operations (flops)for algorithm 1, we can see that step 8 has two flops. Algorithm 2 has four flopsfor that step, and algorithm 3 has nine flops. This step in the algorithm is the mostrelevant because it is the one that makes the difference between the root finders.Table A.1: Floating point operationsAlgorithm Floating point operations (flops) in step 8Bisection 2Regula falsi 4Ridder’s 8When the code is running, we can choose between two stopping criteria, thefirst one is a fixed number of iterations and the second one is tolerance to conver-gence. If we decide to use the first one, the number of flops determines the timedifference between the algorithms. In this case, Ridder’s would be the slowest asit has more flops to calculate. Furthermore, Bisection is the fastest as it has thesmallest number of flops. On the other hand, if we decide to use the tolerance toconvergence, we need to know the value of the asymptotic error constant O(βk).We explain the process briefly to calculate this. To predict which algorithm will54run faster an important factor is the theoretical rate of convergence [8]. For thebisection method, given a domain [a,b] ∈ R we have that for each iteration con-tained in the interval {ak,bk},then bk−ak = 2−k(b−a) and we can bound the errorby ek = bk − ak. And as r = 1 we have that the convergence is linear with theasymptotic error constant equals to 1/2.Following analogous concepts and procedure for the Ridder’s method, we canfind for example in [72] that O(βk) =√2≈ 1.4142.In the case of regula falsi, the rate of convergence is linear, similarly to the bi-section method ([28]). In this case, the asymptotic error constant depends on whichpoint remains fixed. In [6] we can be found an important result about O(βk). Theyshow that |O(βk)|< 1. This result gives us the last information for our comparisonthat we summarize in table A.2Table A.2: Benchmark for modification of the algorithm.Algorithm Asymptotic error constant O(βk)Bisection 1/2Regula falsi ≤ 1Ridder’s 1.4142This analysis tells us that if we choose the tolerance as stopping criteria, a goodoption would be the Ridder’s algorithm, but if we use a fixed number of iterations,the bisection would be better. One last factor to take into account is the stabilityof the algorithm. Bisection converges for all cases; however, Ridder’s and regulafalsi do not. In practice, this fact is more precise, and therefore we chose bisectionto do the benchmark that we performed for the new algorithm.55
Thesis/Dissertation
eng
Mechanical Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Frigaard, Ian
GPU computing of yield stress fluid flows in narrow gaps
Text