Thermal Analysis of a Lengthy Section of a Gas Pipeline on Permafrost

Numerical solution of the problem was implemented in Frost 3D Universal software, with a computational mesh of 58.5 million nodes, to predict ground thaw under the influence of the pipeline. A 2-year simulation of the ground thaw around the pipeline was performed. The computation in Frost 3D Universal software on a quad-core CPU took approximately 70 hours and used 17 GB of RAM. We had already simulated this model before with less detailed discretization: the mesh consisted of 22.3 million nodes and this wasn’t enough to enable us to factor in small elements such as thin heat insulators. Using an NVIDIA Titan graphics accelerator, the prediction of ground thaw over a 20-year period took about 7 hours.

The model of a lengthy section of pipeline

The model of a lengthy section of a gas pipeline

The specific nature of the current problem – simulation of the ground thaw formation along a 500-meter section of the oil pipeline, 1.2 meters in diameter – rendered a large computational mesh necessary. This section is characterized by a complex geological-lithological structure of grounds, which also includes layers of ice deposits. A total of 21 ground typologies with various thermophysical properties were revealed in this section alone. Altitude differences in the daylight surface of 35 meters were also taken into account.

Continue reading

Parallel Computing GPU vs. CPU. Technical Details, Comments

The article “Nvidia Tesla K20 vs Intel Core i7: Speed Comparison”, published on the company's website, and, at the request of many readers, additional technical information regarding the computations performed on the NVIDIA Tesla K20 and NVIDIA Geforce GTX 560 Ti (see the original post on Nvidia web site) has been placed here on the corporate blog.

The first benchmark tests on the heat problem with phase transition (Stefan problem) calculation were executed for a sphere with a radius of 10 meters for 364 days. The summary table below gives the computational time taken by the corresponding graphics accelerator.

Number of cells along each axis X, Y, Z 100 150 200 250 300 350
Number of nodes (millions) 1,03 3,44 8,12 15,81 27,27 43,24
Time in seconds for NVIDIA Tesla K20 7 52 212 637 1625 3394
Time in seconds for NVIDIA GeForce GTX 560 Ti 10 75 311 923

Continue reading

Algorithms and Mathematical Notes of Creating Cylinder Shaped Voxel Meshes

In the range of applied problems for various industries, for example, architecture, fine arts, animation (movies), meshes of complex configuration may be of great demand. Let us show, how it is possible to describe a configurationally complex computational mesh, consisted of voxels (bubble shaped), which are located on a sheet and twisted into cylinder.

Computational mesh,which is constituent of voxels, situated on a sheet and twisted into cylinder

Let us describe a mathematical model used at the design of algorithms of cylinder shaped voxel meshes.


There is a sheet with voxels (bubbles of cylinder shape) located on its surface. Then, it is assumed that voxels are situated on the sheet according to the template, demonstrated in the picture 1.

Voxels location. Template

Picture 1– Voxel location template

Continue reading

Application of Static Adaptive Partitioning of the Computational Domain


Many software packages for numerical computations allow users to use a static adaptive (hereby referred to as adaptive) step in the construction of an orthogonal hexahedral structured computational mesh. This means that informed users can employ their experience to get a more accurate computation without significantly increasing the computation time by specifying the areas of the computational domain in which, in their opinion, it is necessary to apply more detailed partitioning (use a smaller spatial step) as compared to the rest of the computational domain.

When properly used, adaptive partitioning of the computational domain is a powerful tool in numerical computations to increase accuracy. However, when the above option is overused, the computational time can increase dramatically without necessarily altering the accuracy of the computation to any significant degree. In this article, we describe the theoretical advantages and disadvantages of using adaptive partitioning of the computational domain, and also give two examples for numerical computations of thermal fields in soils. In the first example, application of the adaptive step is appropriate; this is not the case, however, in the second.

Continue reading

How to Eliminate Blank Space Error While Building a New NVIDIA CUDA Runtime Project in Microsoft Visual Studio

Assume the developer has Microsoft Visual Studio 2008, 2010 or 2012 and NVIDIA CUDA Toolkit 5.5 installed on the PC. If the Windows username contains a blank space, then building a new NVIDIA CUDA 5.5 Runtime project might fail. Typically, this error is either “Could not setup the environment for Microsoft Visual Studio using …” referring to the files of vcvars32.bat or vsvars32.bat, or “stumbling” over a blank space in the name of the user’s home directory (see the screenshots below).

Build errors of a new NVIDIA CUDA 5.5 Runtime project in Microsoft Visual Studio 2012

Build error of a new NVIDIA CUDA 5.5 Runtime project in Microsoft Visual Studio 2008

Build errors of a new NVIDIA CUDA 5.5 Runtime project in Microsoft Visual Studio 2012. Error D8022

More details on these errors can be found in the build log if the verbose logging is on: use --verbose command line option of NVCC (CUDA C/C++ compilation driver).

To fix the problem you may choose any of the following ways:

  1. Create a new user account without a blank space in the name (be careful with the previous user data!),
  2. Set the user environment variable TEMP to a path that doesn’t contain a space (the folder at this path must exist), e.g. C:\Temp,
  3. Insert into the text file %CUDA_PATH%/bin/nvcc.profile, where %CUDA_PATH% may look like C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v5.5, the line “TEMP=path_without_space”, where path_without_space is the path that doesn’t contain a blank space (the folder at the path must exist), e. g. C:\Temp.

CUDA Project Template for Microsoft Visual Studio 2008

The Microsoft Visual Studio 2008 IDE has a lot of standard templates for a New Visual C++ Project. But when using such parallel programming tool as NVidia CUDA Toolkit 5.0 (without any add-ins) one has to set every new C++ CUDA project manually.

Fortunately, the routine of setting becomes automated after installing the free CUDA VS Wizard 2.9. In few clicks one may create an empty CUDA project that compiles on Win32 configuration. To compile on x64 (see item 2) and to achieve purely cosmetic effects (items 3 and 4) complete the following steps.

1. Find the IDE installation folder, e.g. C:/Program Files (x86)/Microsoft Visual Studio 9.0. Further it is denoted as $(VSInstallDir).

2. Open the file $(VSInstallDir)/VC/VCWizards/CUDA/CUDAWinApp/Scripts/1033/default.js in a text editor. Find all occurrences of the path «$(CUDA_PATH)\lib» and replace them with «$(CUDA_PATH)/lib/$(PlatformName);».

3. Add the line «CUDAWinApp.vsz| |CUDAWinApp|1|A sample CUDA application.| |6777| |CUDAWinApp» to the text file $(VSInstallDir)/VC/vcprojects/vc.vsdir.

4. Correct the text of the hint «TODO: Wizard Description.», changing it to «A sample CUDA application.», in the file $(VSInstallDir)/VC/vcprojects/CUDA/CUDAWinApp.vsdir

The transfer of boundary conditions on an orthogonal hexahedral mesh

In the simulation of physical processes and phenomena overall there is a problem of numerical solutions of differential problems in partial derivatives. One of the methods for numerical solutions of mathematical physics equations is the approach, based on the finite difference approximation. However, the major drawback of this method is the need to use orthogonal hexahedral mesh. While solving problems in practice, sometimes a complex geometric configuration of the computational domain (for example, see Figure 1) is used, and therefore, the drawback, mentioned above, is enough critical.

Picture 1 – The example of a complex surface as a part of the computational domain

Due to the aforesaid, there is quite a natural problem of the geometrical configuration approximation of computational domain by cell faces of the given orthogonal hexahedral mesh.

Let us turn to a more detailed problem statement. Let the computational mesh in 3D space be orthogonal and hexahedral. Namely, let be a set of mesh divisions along axis arranged in ascending order, similar to it let us introduce the ordered sets of mesh divisions along and respectively. According to designations, the following set of points is the set nodes of computational mesh . Geometric configuration that requires approximation by mesh faces, is defined by triangulated surface. Let be a set of triangles. This set is a part of surface triangulation, and its quantity is .With the numerical solution of mathematical physics problem, besides geometric approximation of surface by faces of mesh cells there is a problem of adequate transfer of triangulated surface area on cells faces, involved in its approximation. Thus, by means of set of triangles it is necessary to determine the set of faces of mesh cells, that approximate the triangulated surface, and match the transferred area with each of such face.

Let us give the algorithm, which solves the set above problem

1) Let us consider the mesh, which is dual to the initial one, i.e. such a mesh the nodes of which are cells centers of the initial computational mesh. Let us store a 3D array of real numbers, elements of which correspond to the cells of dual mesh. First of all we initialize elements of array by zeros.

Continue reading

Quasilinear Heat Equation in Three Dimensions and Stefan Problem in Permafrost Soils in the Frame of Alternating Directions Finite Difference Scheme

You can read the original article at the World Congress Engineering web site.


Abstract — The quasilinear heat equation with thermal conductivity and heat capacity depending on the temperature field in three spatial dimensions is studied in application to the phase transition problem in permafrost soils. The conditions under which the alternating directions Douglas – Rachford finite difference scheme retains numerical stability are explicitly formulated. The comparison with the known analytical similarity solution to the Stefan problem in one spatial dimension is performed.

Index Terms—quasilinear heat equation, Stefan problem, finite differences, alternating directions scheme, numerical stability.


Since the early formulations of the alternating directions implicit methods (ADI methods) [1], [2], they have been tremendously developed and found a vast number of applications [3], [4]. Nevertheless, serious difficulties are encountered with the use of these methods in application to problems with complex geometries [5] and/or nonlinear equations of mathematical physics [6].

Although the schemes of the ADI methods are proved to be efficient and economic with respect to time consumption and, in most cases, unconditionally stable, they exhibit some disadvantages:

1) Their finite-difference formulations permit to consider only rectangular spatial domains (due to commutativity requirements imposed on the factorized and split operators) [7];

2) The application of the ADI schemes to the problems with Neumann and Robin boundary conditions that are varying in time encounters serious problems due to the necessity of evaluation or approximation of these boundary conditions at the intermediate steps of the scheme [8];

3) When applied to the solution of nonlinear heat equations, the operators constituting an ADI scheme do not commute, thus leading to the loss of unconditional stability of the scheme [6].

The first of the above disadvantages can be overcome either by the use of finite elements methods in conjunction with operator splitting techniques, or by domain decomposition techniques. The second and the third disadvantages present an important problem for the successful application of the ADI scheme. To the best of our knowledge, no complete stability analysis for an ADI scheme applied to the nonlinear heat equation in a three-dimensional spatial domain is available in the literature, thus motivating this work.

Another motivation for the present work is the application of ADI scheme to the modeling of heat transfer in large scale environmental systems (e.g., large areas of permafrost soils) which, in the case of purely explicit finite-differences schemes, imposes stiff constraints on the time-step in order to guarantee the numerical stability. At the same time, implementation of implicit schemes can often lead to much greater computing expenses than that of explicit schemes, especially for the problems with rapidly changing coefficients in complex geometries and substantially nonhomogeneous meshes. Thus, in modeling of heat transfer in large scale systems the necessity of making an optimal choice between explicit and implicit schemes arises. In case of finite-elements method, applied to modeling of processes in permafrost soils, the analysis of numerical stability appears to be so complex, that the stability criterion is often established empirically [9].

In the present paper we discuss the application of the ADI Douglas – Rachford scheme to the solution of Stefan problem in porous permafrost soils. The paper is organized as follows: next section contains the problem formulation and some assumptions that will be used in the proof of numerical stability of the ADI scheme while section 3 exposes the proof itself. Section 4 presents some numerical results and is followed by Conclusions.

Continue reading

3D geological modeling based on boreholes data


3D geological modeling is a very actual issue nowadays in building development, environmental assessment of soil (variably-saturated porous medium) pollution, assessment of mineral deposits, etc. There are different approaches to solve this problem by means of modern software designed for simulation in geology [1-3]. The most frequently used method is that of reconstructing geological model. This method is based on information about the levels of geological horizons occurrence received from the results of drilling [4-6]. The implementation of this method itself may have some peculiarities.

In this article an alternative approach for 3D geological model creation is being proposed. It is based on the following:
1) Surface triangulation of site topology
2) Automatic cross-section generation
3) Segment height interpolation for each layer of geological model.
This approach allows both simplify and accelerate 3D geological model creation while maintaining acceptable 3D site building accuracy.

The proposed method consists of six basic steps described below. The following information on boreholes is considered as given data: 1) borehole coordinates; 2) seamark; 3) capacity of geological horizons.

Continue reading