**TABLE OF CONTENTS**

1\. *V.L.Averbukh, A.Yu.Baidalin, P.A.Vasev, A.I.Zenkov, D.R.Ismagilov,
D.V.Manakov, D.S.Perevalov, M.R.Shagubakov*

PARALLEL COMPUTATION VISUALIZATION
PROBLEMS..........................................\...7

2\. *A.L. Adrianov*

ASYMPTOTIC MODEL OF SHOCK-SHEAR LAYER
INTERACTION...................................7

3\. *A.V. Alekseev, O.E. Vlasova, A.V. Gichuk, V.V. Evdokimov, Е.A
Irinichev,S.V. Mzhachikh, AA. Nuzhdin, N.P. Pleteneva, R.M. Shagaliev,
H. Scott*

THE RESULTS OF NUMERICAL STUDIES OF VARIOUS APPROXIMATION METHODS FOR
THE 2D TRANSPORT EQUATION OVER ANGULAR
VARIABLES....................................8

4\. *T.N.Alekseeva, S.G.Volkov, S.E.Kuratov, O.V.Olkhov*

\"MEDUZA\" COMPUTATIONS OF RICHMYER-MESHKOV INSTABILITY WITH INCLUSION
OF VORTEX
COMPONENT....................................................................................\...9

5\. *Raymond E. Alcouffe*

A PARALLEL MULTIGRID METHOD FOR INVERSION OF THE DIFFUSION OPERATOR IN
NEUTRONICS
APPLICATIONS..................................................................................9

6\. *А.А.Anishchenko, A.M.Grigorjev, N.M.Mikiytchuk, S.N.Nikonov,
E.E.Santalova, O.N.Tazetdinov*

NUCLEAR MATERIAL CONTROL AND ACCOUNTING SYSTEM \"ACCORD-2000\":
PRINCIPAL
CONCEPTS..........................................................................................10

7\. *A.Yu. Artemiev, Yu.G.Bartenev, Yu.A.Bondarenko, V.A.Erzunov,
M.A.Kolesov, A.V.Lomtev, A.S.Maksimov, A.I.Panov, A.I.Prokofiev,
M.D.Romanova, A.P.Trubitsin, N.V.Frolova, T.V.Tsareva,*

*E.B.Shchannikova, Y.Dandass, D.E.Nielsen, A.Skjellum, V.Valsalam, W.Wu*

THE LIBRARY OF SOLVERS FOR SPARSE LINEAR SYSTEMS. CAPABILITIES, USE,
DEVELOPMENT
PROSPECTS..................................................................................11

8\. *D.A. Afremov, A.V. Kuzin, V.E. Men'shikov, G.S. Mingaleeva, Yu.V.
Mironov, M.Yu. Chernov*

COMPARATIVE COMPUTATIONS FOR SOME MODEL PROBLEMS ON GENERAL-CIRCUIT

THERMOHYDRAULIC
CODES................................................................................11

9\. *Afremov D.A., Mironov Y.V., Radkevich V.E., Zhuravleva J.V.*

A TECHNIQUE FOR statistical analysis of thermohydraulic COMPUTATION
uncertainties.................................................................................................12

10\. *R.A.Barabanov, V.I.Budnikov, O.I.Butnev, V.I.Delov, O.K.Loginova,
V.A.Pronin, V.V.Sadchikov*

AUTOMATIC CORRECTION OF 2D MESH FRAGMENTS IN D CODES WITH MIXED
CELLS................................................................................................................13

11\. *R.A.Barabanov, O.I.Butnev, S.G.Volkov, B.M.Zhogov, O.G.Loginova,
V.A.Pronin, V.V.Sadchikov*

MEDUZA COMPUTATIONS OF INSTABILITY DEVELOPMENT AT GAS-GAS INTERFACE WITH
CONTACT LINE SEPARATION IN MIXED
CELLS...............................................14

12\. *V.A.Batalov, V.F.Volgin, I.D.Sofronov*

SOFTWARE DEVELOPMENT AND OPERATION TECHNOLOGY AT VNIIEF COMPUTATIONAL
DIVISION..................................................................................15

13\. *S.M.Bakhrakh, I.P.Bakumova, P.N.Nizovtsev, E.V.Shuvalova*

RELAXATION MODEL FOR SHEAR STRENGTH IN EULERIAN
FORMULATION...............16

14\. *S.M.Bakhrakh, N.A.Volodina, I.V.Kuz'mitsky, M.N.Leontyev,
K.V.Tsiberev*

COMPUTATIONS FOR LOW-VELOCITY MODE OF HE DETONATION DEVELOPMENT...\...16

15\. *S.M.Bakhrakh, N.A.Volodina, P.N.Nizovtsev, Spiridonov V.F.*

NUMERICAL SIMULATION OF INITIAL PERTURBATION GROWTH WITH OBLIQUE-ANGLE
IMPACT OF METAL
SLABS.........................................................................17

16\. *S.M.Bakhrakh, S.V.Velichko, V.F.Spiridonov*

CRASH-FREE TECHNOLOGY FOR CONTINUUM FLOW COMPUTATIONS USING \"LEGAK\"
CODE ON MULTYPROCESSOR COMPUTER
SYSTEMS................................................18

17\. *V.A. Bakhtin, N.A. Konovalov, V.A. Kryukov*

OpenMP Fortran EXTENSION FOR DISTRIBUTED-MEMORY
SYSTEMS............................18

18\. *V.V.Bashurov, N.V.Samodolova*

SPH CODE FOR SIMULATING 3D ELASTIC-PLASTIC FLOWS USING THE METHOD OF
SMOOTHED
PARTICLES........................................................................................18

19\. *V.V.Bashurov, F.A.Pletenev, A.V.Ryzhikh*

CODES FOR VISUALIZING MULTIDIMENSIONAL DATA ON ARBITRARY
GRIDS............19

20\. *V.V.Bashurov, V.P.Bashurin, O.A.Vinokurov*

A CODE FOR COMPUTATION OF POLISHING IN 3D
GEOMETRY..................................20

21\. *S.A.Belkov, C.V.Bondarenko, O.A.Vinokurov, G.G.Kochemasov,
L.S.Mkhitaryan, L.F.Potapkina, N.A.Ryabikina*

\"MIMOZA-ND\" SIMULATIONS FOR EXPERIMENTS WITH A SPHERICAL BOX ON LASER
FACILITY
\"ISKRA-5\".............................................................................................20

22\. *S.P.Belyaev, L.I.Degtyarenko, A.L.Stadnik, Yu.V.Yanilkin*

A CODE FOR 3D HYDRODYNAMICS WITH FINE-GRAIN POINTWISE
PARALLELIZATION..............................................................................................21

# 23. *G. Ben-Dor, E.I. Vasilev, L.F. Henderson, and T.Elperin*

# THE JET STREAM EFFECT IN PSEUDOSTATIONARY MACH REFLECTION.....................21

24\. *V.L.Bogdankevich*

COMPUTATION PARALLELIZATION USING MPI 1.3 MESSAGE PASSING INTERFACE
WITHIN INTELLIGENT SOFTWARE SHELL ShIPRW FOR MODELING NUCLEAR REACTORS
.........................................................................................................................22

25\. *S.N.Boldyrev, S.A.Sukov, M.V.Yakobovsky*

USE OF HIERARCHICAL GRAPH PARTITION METHODS FOR NUMERICAL SOLUTION OF
APPLICATION PROBLEMS USING IRREGULAR GRIDS ON MULTIPROCESSOR
SYSTEMS............................................................................................................23

26\. *A.I.Bochkov, V.A.Shumilin*

PARALLELIZATION OVER DIRECTIONS FOR THE CALCULATION OF THE 2D TRANSPORT
EQUATION WITH SATURN-3 CODES USING THE OpenMp
INTERFACE........................\...24

# 27. *T.N.Bronina, Irina A.Gasilova, O.V.Ushakova*

# TECHNOLOGIES FOR GENERATION OF 3D STRUCTURED GRIDS.................................25

28\. *Yu.N. Brylev, N.F. Podderyugina, I.F. Podlivaev*

EFFICIENCY OF USING POWERFUL COMPUTERS IN CONSTRUCTION OF COMPUTATIONAL
METHODS AND TREATMENT OF COMPUTED DATA FOR ELECTROMAGNETIC WAVE
IONOSPHERE REFLECTION WITH ACCOUNT FOR NONLINEAR HEATING
.........................................................................................26

29\. *A.N.Bykov, A.S.Zhdanov*

COMBINED MODEL OpenMP & MPI FOR PARALLEL COMPUTATIONS IN GAS DYNAMICS
PROGRAM...........................................................................................................26

30\. *M.A.Vagin, A.V.Kosterin, K.A.Potashev*

PROBLEMS OF LIQUID POLLUTIONS IN
SOILS..........................................................27

31\. *E.I.Vasilyev, D.V.Tsitser*

FLOW FORMATION AND DEVELOPMENT IN A SHOCK TUBE AT DIAPHRAGM OPENING

USING A
LEAF...................................................................................................\...28

32\. *E.I.Vasilyev, A.V.Zenovich, A.V.Karpov*

AUTO-OSCILLATORY MODES OF SUPERSONIC FLOW ABOUT A DOUBLE
WEDGE......\...28

33\. *B.A.Voinov*

A CODE SIMULATING FRACTURE GROWTH IN ELASTIC-PLASTIC
MEDIUM.................29

34\. *E.V.Voskresensky, T.F.Mamedova*

METHODS FOR STABILIZING PROGRAM
MOTIONS...................................................29

35\. *K.G.Gainullin, V.A.Zhmailo, Yu.F.Kirianov*

NUMERICAL ALGORITHM FOR THE CALCULATION OF 3D EQUATIONS DESCRIBING THE
PROPAGATION OF MAGNETO HYDRODYNAMIC WAVES IN
IONOSPHERE...................30

36\. *S.M.Garina, N.V.Zmitrenko, M.E.Ladonkina, N.G.Proncheva,
V.F.Tishkin*

NUMERICAL SIMULATION OF GRAVITATIONAL MIXING USING MULTIPROCESSOR
COMPUTER
SYSTEMS...........................................................................................31

37\. *A.E. Godin*

GENERALIZED ASSIGNMENT
PROBLEM...............................................................\...31

######## 38. *A.I. Golovanov, S.A. Kuznetsov*

CALCULATIONS OF LARGE ELASTIC-PLASTIC DEFORMATIONS IN RATES WITH
ADDITIVE AND MULTICATIVE DEFORMATION
EXPANSIONS....................................32

39\. *A.I.Golubev, T.G.Sysoyeva*

A SCHEME OF SPLITTING FOR NUMERICAL SOLUTION OF THE VLASOV KINETIC
EQUATION.........................................................................................................32

40\. *A.N.Grebennikov, N.A.Krutko, G.G.Farafontov*

THE INFORMATION AND REFERENCE SYSTEM FOR NUCLEAR DATA
NDX..................33

41\. *Guk N.A*., *Obodan N.I.*

OPTIMAL COOLING OF HEATED THIN-WALL SYSTEMS IN THERMAL SHIELD
BREAKDOWN......................................................................................................34

42\. *N.N.Degtyarenko G.M.Yeliseyev*

COMPUTATION OF PHOTON ABSORPTION CROSS SECTIONS IN MIXTURES USING DATA
FOR
COMPONENTS..............................................................................................35

43\. *Yu.N.Deryugin, A.A.Bazin, O.M.Velichko, N.F.Gavrilov, V.V.Gorev,
I.V.Gorev, G.N.Gubkova, D.K.Zelensky, A.I.Panov, E.I.Ryabov, V.I.Seliv,
L.I.Sizova, B.P.Tikhomirov, A.V.Kosterin, A.N.Chekalin, V.M.Entov*

NUMERICAL SOLUTION OF 3D GROUNDWATER IMPURITY TRANSPORT PROBLEMS ...36

44\. *L.V.Dmitrieva, I.I.Karpenko, N.V.Korepova, V.G.Morozov*

NUMERICAL SIMULATION OF THE FLEGMATIZED OCTOGEN SAMPLE RESPONSE TO THE
SHOCK CONSIDERING THE DETONATION KINETICS IN THE D
CODES..................37

45\. *E.N.Donskoy, A.N.Zalyalov*

INCLUSION OF BREMSSTRAHLUNG IN GAMMA-RAY TRANSPORT
PROBLEMS.............38

46\. *S.P. Egorshin, T.V. Rezvova*

A PARALLEL ROUTINE FOR THE GLOBAL SCALING OF GRID VALUES IN \"LEGAK\"
CODES................................................................................................................39

47\. *G.S Elizarov, V.V. Karatanov, V.V. Korneev, V.K Levin*

DESIGN OF MULTIPROCESSOR
SYSTEMS................................................................40

48\. *N.V. Elkina, V.D. Levchenko*

THE METHOD OF LOCAL SPACE-TIME DECOMPOSITION IN PROBLEMS OF SHOT LASER
PULSE-PLASMA INTERACTION
SIMULATIONS.........................................................41

49\. *A.S.Zapasskaya*

IMPLEMENTATION OF THE MIXED PARALLELIZATION BASED ON MPI AND OPEN MP
STANDARDS FOR 2D STATIONARY TRANSPORT
EQUATION......................................42

50\. *V.V.Zmushko*

\"ISOMIM\", A CODE FOR GRAPHICALLY PROCESSING 2D COMPUTATION RESULTS IN
\"MIMOZA\"
COMPLEX...........................................................................................43

51\. *K.M.Zuyev*

INTEGRABILITY OF (SU(n),S(U(k~1~)×U(k~2~)×U(k~3~))) PAIR IN SMALL
DIMENSIONS..............43

52\. *R.T.Islamov*

BRANCH COMPUTER CODES FOR PROBABILISTIC ANALYSIS OF SAFETY AND
UNCERTAINTY OF DETERMINISTIC
MODELS..........................................................44

53\. *T.Istomin, O.Panin*

THE PARALLEL PROGRAMMING SYSTEM BASED ON TYPE ALRORITHMIC
STRUCTUERS......................................................................................................45

# 54. *I.M.Kanayev*

# ON LOWERING THE COMPLEXITY OF DIGITAL IMAGES...........................................46

55\. *D.V.Kovalenko, D.P.Sidorov, S.A.Fedosin*

EMPLOYMENT OF DISTRIBUTED COMPUTERS FOR LARGE NUMBER
FACTORIZATION.................................................................................................47

56\. *A.A. Kovalishin, N.I Laletin, N.V.Sultanov*

APPLICATION OF MSPS AND MSH METHODS IS AN EFFICIENT AND USABLE WAY T0
PARALLELIZE NUCLEAR
COMPUTATIONS..............................................................49

57\. *N.A.Konovalov, V.A.Kryukov, N.V.Podderyugina, Yu.L.Sazanov*

PARALLELIZATION OF COMPUTATIONS IN DVM
SYSTEM........................................50

58\. *P.K.Konosavsky, K.A.Soloveychik*

PARALLELIZATION OF ALGORITHMS EOR NUMERICAL GROUNDWATER DYNAMICS
MODELS............................................................................................................51

59\. *A.I.Korotkiy, I.A.Tsepelev*

3D SIMULATION OF DIRECT AND INVERSE RALEIGH-BERNARD AND RALEIGH-TAYLOR
PROBLEMS.........................................................................................................52

60\. *V.M. Kosenkov*

COMPUTER SIMULATION OF NEAR-BREAST AREA DYNAMICS OF WATER WELL UNDER
ELECTRICAL DISCHARGE IN
WATER.....................................................................54

61\. *V.M.Kyashkin, S.V.Rod'kin, O.B.Tomilin, S.A.Fedosin*

ON USING ERROR BACK PROPAGATION NEURON NETWORKS IN SIMULATIONS OF
MULTI-MATERIAL
SYSTEMS.................................................................................55

62\. *E.F.Lelikova, L.I.Rubina, O.N.Ul'yanov, M.A.Chashchin*

PARALLEL COMPUTATION TECHNOLOGIES FOR THE PROBLEM OF RADIATION
TRANSPORT........................................................................................................58

63\. *N.V.Linich*

A REMOTE DATA VISUALIZATION
SYSTEM.............................................................59

64\. *Yu.B. Malykhanov*

SOLVING A MANY-ELECTRON PROBLEM USING THE METHODS OF MINIMIZING
FUNCTIONS OF SEVERAL
VARIABLES..................................................................\...60

65\. *V.S.Popov, S.A.Stepanenko, A.A.Kholostov*

HARDWARE ARCHITECTURE IN THE INTERPROCESSOR COMMUNICATION SYSTEM OF MP-X
MULTIPROCESSOR
SYSTEM..........................................................................60

66\. *S.V.Rebrov, O.G.Sin'kova, A.L.Stadnik, V.P.Statsenko,
Y.V.Yanilkin*

STUDYING TURBULENCE CHARACTERISTICS BY DIRECT 3D NUMERICAL SIMULATION OF
GRAVITATIONAL
MIXING..............................................................................61

67\. *N.M.Rud'ko, S.I.Skrypnik*

TECHNIQUES FOR NUMERICAL SOLUTION OF A 2D NONLINEAR HEAT CONDUCTION
EQUATION IN COMPLEX GEOMETRY
ROBLEMS......................................................62

68\. *A.A.Sadovoi, N.V.Sokolova*

MODELING OF ISOTROPIC DESTRUCTION OF ELASTIC-PLASTIC MEDIA WITH MIMOZA
CODE.................................................................................................................62

69\. *A. Sazonov, A.N. Salnikov*

A PROTOTYPE APPLICATION DEVELOPMENT SYSTEM FOR PROCESSING EXPERIMENTAL
DATA USING HETEROGENEOUS MULTIPROCESSOR SYSTEMS...........63

70\. *S.S. Sokolov*

TECHNIQUES FOR SOLVING TIME-DEPENDENT ELASTIC-PLASTIC PROBLEMS USING
IRREGULAR POLYHEDRAL LAGRANGIAN
GRIDS....................................................64

71\. *E.A.Solovyeva, A.M.Stenin*

STUDYING THE EFFECT OF SOME ARTIFICIAL VISCOSITY FORMS ON THE \"ENTROPY
TRACE\" NATURE IN NUMERICAL SOLUTION OF GAS DYNAMICS
EQUATIONS.............64

72\. *S.A.Stepanenko*

TOPOLOGICAL REDUNDANCY IN MULTIPROCESSOR
ENVIRONMENTS.......................66

73\. *S.A.Stepanenko*

THE EFFICIENCY OF MULTIPROCESSOR ENVIRONMENTS WITH GLOBAL
COMMUNICATIONS..............................................................................................67

74\. *E.G. Sukhov*

PARALLEL COMPUTATIONS WITH TOEPLITZ
MATRICES...........................................67

75\. *L.N.Sukhova*

ON COMPUTATIONS OF TWO-COMPONENT FLOWS IN COMPLEX-GEOMETRY
CHANNELS.........................................................................................................68

76\. *I.T.Usmanov*

A MULTI-GRID METHOD IMPROVEMENT FOR HYDRODYNAMIC SIMULATIONS OF
OILFIELDS WITH BOTTOMHOLE PRESSURES
SPECIFIED...........................................68

77\. *Filippov N.A.*

THE CHANNEL-TEMPORARY DIAGRAM OF FAST CALCULATIONS OF INDIVIDUAL
SUMS..................................................................................................................69

78\. *A.F. Khairullin, O.B. Khairullina*

CONSTRUCTION OF OPTIMAL MESHES IN MULTIPLY CONNECTED COMPLEX-TOPOLOGY
DOMAINS ON MULTIPROCESSOR
COMPUTERS.......................................70

79\. *P.A. Shamanaev*

ON LOCAL REPRODUCIBILITY OF NONLINEAR DIFFERENTIAL EQUATION
SYSTEM.....71

80\. *V.E.Shemarulin*

STRUCTURE OF 3D STATIONARY ISOBARIC DOUBLE WAVES IN PERFECT
INCOMPRESSIBLE
FLUID......................................................................................72

81\. *Yu.V.Yanilkin, A.V.Gorodnichev, S.P.Belyaev, L.I.Degtyarenko,
A.L.Stadnik, N.A.Khovrin, A.V.Volgin*

NUMERICAL SIMULATIONS OF 2D GAS-DYNAMIC AND ELASTIC-PLASTIC FLOWS ON
ADAPTIVELY EMBEDDING FRAGMENTED COMPUTATIONAL
GRID............................73

**PARALLEL COMPUTATION VISUALIZATION PROBLEMS**

***V.L.Averbukh, A.Yu.Baidalin, P.A.Vasev, A.I.Zenkov, D.R.Ismagilov,
D.V.Manakov, D.S.Perevalov, M.R.Shagubakov***

# (IMM of RAS UrB, Ekaterinburg)

Both problems of representing parallel computation results proper and
issues concerning parallel software visualization should be referred to
a set of parallel computation visualization problems. The latter
traditionally includes issues of developing visual languages for
parallel programming, visual debuggers and tools for parallel code
performances measuring, debugging and adjustment.

RAS UrB Institute for Mechanics and Mathematics (IMM) performs studies
and developments concerning the whole complex of visualization problems
to represent results of computations using parallel codes.

Visualization of parallel computation results is related to problems of
data amounts to be visualized, however, more significant problems are
caused by computational model complexities and, hence, by the complexity
of the problem of abstract data imaging and adequately interpreting such
data by users.

RAS UrB IMM developed a number of specialized systems that support
visualization of computations using specific computational models and
carefully consider each problem and parallel implementation specifics.
In our opinion, interesting results have been obtained when representing
4D sets that occur during simulation of chemical reactions of splitting
on a parallel computer. Another line is to develop tools providing
on-line visualization of parallel computations. Such tools, in
particular, are required to support initial stages of simulation, when
the user needs to force the computation process. Such tools are also
required to create tools for parallel software visualization, for
example, visual debuggers for debugging the correctness and efficiency
of computations.

Developments on parallel software visualization are, first of all,
related to DVM parallel programming system developed at RAS M.V.Keldysh
IAM. DVM possesses a smart system of debugging data acquisition that
allows debugging the correctness and performance of massively parallel
codes in text format.

Nowadays, a pilot version of the system for visual representation of
data on DVM-codes performances is implemented. The system gives the user
a rich set of various capabilities to visualize complex temporal
characteristics of code execution process. Design and development of a
full-value environment for visual programming, including tools to
support DVM-code designing, and a visual correctness debugger are
planned. It is supposed to use 3 dimensions, animation and elements of
virtual reality among the methods of representing DVM system's program
entities.

**ASYMPTOTIC MODEL OF SHOCK-SHEAR LAYER INTERACTION**

A.L. Adrianov

(Academician M.F.Reshetnev SAA, Krasnoyarsk)

A new efficient model of shock-shear layer interaction is suggested. A
time-independent analytical solution to the problem is obtained. It is
essential that the shock is therewith replaced with a surface of strong
gas-dynamic discontinuity, which ordinary oblique shock conditions and
differential conditions including medium viscosity and heat conduction
effects are imposed on /1-3/. In addition, the boundary effect following
the shock /4/ which is insignificant to the problem is excluded. The
conditions allow the initial boundary problem for Navier-Stokes
equations to be reduced to the initial value problem for a stiff
ordinary differential equation system with a complex right-hand side. An
efficient numerical method for solving the system without its
normalization is suggested.

The influence of Reynolds number, initial incident shock intensity,
boundary (initial) condition type, etc. on the \"refraction\" mechanism
of the shock in the shear layer is numerically studied. The model is
used to computationally prove the existence of a perturbed flow
\"self-similarity zone\" with respect to Reynolds number that appears at
significant shock intensities, that is when pressure forces begin to
prevail over viscous forces and, hence, the concept of the gas-dynamic
approximation \"works\". This physical fact is also confirmed with more
complex calculations using the complete model of time-dependent
Navier-Stokes equations.

References

1.  Bai Shi-i. *Introduction to the theory of compressible fluid flow*.
    Moscow, *Izdatelstvo Inostrannoy Literatury Publishers*, 1962. 411p.

2.  Adrianov A.L., Starykh A.L., Uskov V.N. *Interference of stationary
    gas-dynamic discontinuities*. Novosibirsk, *Nauka Publishers*, 1995.
    180 p.

3.  Adrianov A.L. *Differential relations on a shock in viscous gas at
    large Reynolds numbers. Proceedings of Seminar \"Computer
    Simulations in Mechanics\"*. RAS SB VTsK. Krasnoyarsk, 1996. 18 p.
    Deposited at VINITI 04/01/96, No. 1052, V-96.

4.  Adrianov A.L. *On model shock curvature in nonuniform flow.
    Vychislitelnyje Tekhnologii*. 2000, V.5, No.6. Pp. 3-14.

**THE RESULTS OF NUMERICAL STUDIES OF VARIOUS APPROXIMATION METHODS FOR
THE 2D TRANSPORT EQUATION IN ANGULAR VARIABLES**

***A.V. Alekseev, O.E. Vlasova, A.V. Gichuk, V.V. Evdokimov, Е.A
Irinichev,S.V. Mzhachikh, AA. Nuzhdin, N.P. Pleteneva, R.M. Shagaliev,
H. Scott***

(RFNC-VNIIEF; LLNL, USA)

The paper was prepared based on the data resulting from the contract
between the Sarov Open Computing Center and LLNL to study various
approximation methods for 2D transport equation written in cylindrical
coordinates in angular variables.

The traditional methods like Sn schemes for the phase space
discretization in angular variables require a very fine angular grid.
Otherwise the well known beam effect occurs.For the system with
optically thin emitting domains the requirements for the number of
intervals in the angular variables are so severe that the CPU and main
memory time cost of the computations can grow dramatically. One of the
solutions for this problem may be the development of more efficient
computational methods for the multidimensional transport equation in
angular variables.

Within the contract, a collection of benchmarks was developed, on the
one hand, relatively simple in terms of physical properties and geometry
specification and on the other hand sufficiently representative to
reflect real methodology problems. In addition, the resulting benchmarks
are relatively complicated for the traditional computational methods for
2D transport equation in angular variables.

Four schemes are examined for the approximation of the transport
equation in angular variables. These include the traditional method for
the transport equation using the schemes like DS~n~ method, improved
scheme meeting the requirements such as accurate calculation of momentum
integrals or ensuring the diffusion limit, two versions of WWD scheme
that scales the outgoing flow through the phase cell faces.

These schemes were implemented in the test TRAST routine used for a
series of numerical studies. The paper discusses the results of studies
for the approximation methods of interest in terms of the grid solution
accuracy, the conclusions are made of how they are efficient on the
given class of applications, the potential for further method
improvement is examined.

**MEDUZA COMPUTATIONS OF RICHMYER-MESHKOV INSTABILITY WITH INCLUSION OF
VORTEX COMPONENT**

***T.N.Alekseeva, S.G.Volkov, S.E.Kuratov, O.V.Olkhov***

(RFNC-VNIIEF, Sarov)

The paper briefly describes the methodology for inclusion of the
possibility to account the vortex component in the Lagrangian-Eulerian
code MEDUZA. The well-known algorithm \"vortex-in-cell\" served as a
basis for the modification to MEDUZA. The obtained codes were used for
the test computations of Richtmyer-Meshkov instability development for
two problem types: interface pulsed loading with one-mode perturbation
in the plane case and instability growth with convergence of cylindrical
shell. The agreement with experimental data and results of computations
with other codes is good.

**A PARALLEL MULTIGRID METHOD FOR INVERSION OF THE DIFFUSION OPERATOR IN
NEUTRONICS APPLICATIONS**

***Raymond E. Alcouffe***

(Transport Methods Group, LANL; USA)

We seek to invert the diffusion operator encountered in the solution of
the neutral particle transport equation using the diffusion synthetic
acceleration method (DSA). In the single processor regime, the most
effective methods for inverting the multidimensional diffusion operator
employs the multigrid method either by itself or as a preconditioner for
a Krylov method. With the advent of the desire to solve very large
problems consisting of spatial domain sizes where the number of mesh
cells is on the order of 1 billion, it is necessary to use a parallel
computer architecture. Thus in order to use the multigrid method in
these circumstances, we must make the algorithm parallel. This must be
done in as efficient a manner as possible in order to make the method
viable in the DSA context of solving the transport equation. In this
paper we present our method of making the multigrid V-cycle parallel by
examining each of the constituents, relaxation, condensation, and
interpolation where the method of parallelization is spatial domain
decomposition. With spatial domain decomposition we use a message
passing paradigm through MPI to move the requisite data among the
processors. The diffusion operator is discretized on the spatial domains
as a seven-point, finite differenced form on a regular orthogonal mesh.
We present a study of the efficiency of solution of the diffusion
problem using the multigrid method on up to 2048 processors. The problem
is chosen from neutronics to show the scaling behavior of the method on
selected parallel architectures and we compare this with other methods
of solution. We then show how the solution efficiency impacts the
solution of the transport equation in the DSA context.

**NUCLEAR MATERIAL CONTROL AND ACCOUNTING SYSTEM \"ACCORD-2000\":
PRINCIPAL CONCEPTS**

***А.А.Anishchenko, A.M.Grigorjev, N.M.Mikiytchuk,***

***S.N.Nikonov, E.E.Santalova, O.N.Tazetdinov***

# (RFNC-VNIIEF, Sarov)

Principal concepts of ACCORD-2000 system that provides computerized
control and accounting of nuclear materials are presented. The system is
intended for its use at various-type enterprises of nuclear complex,
including those operating bulk-form materials. The model of distributed
multicomponent items by Microsoft-DCOM is taken as the technology basis.
The system has a three-layer architecture: a database, a business-logic
level, and a user's interface. ACCORD-2000 is a flexible, easily
modifiable and adjustable system. Use of simple text files of item
descriptions allows the system to be adapted to any enterprise specifics
using an ordinary text editor.

Particular attention is paid to the problem of data protection. Data
access is performed only via the stored procedures. Capabilities of
adjusting user's rights to access any items of the system (containers,
accounting items, operations, documents, etc.) are provided for. Control
of a person's data access rights is implemented on a server and cannot
be violated by client actions. To receive data, client programs call
special ActiveX objects that are executed on the server and responsible
for data protection in the system.

The system provides execution of several dozens of various accounting
operations, it is open for new operations to be included.

The system's preparations for printing a large set of paper accounting
and reporting documents in Ms Word 2000 format have been implemented.
The mechanism of including any new documents has been developed.

ACCORD-2000 uses barcode readers and weighting equipment to enter data
on containers. Codes for a number of devices have been implemented,
including those for TRAKKER ANTARES-2481 stationary reader and TRAKKER
ANTARES-2420, and JANUS-2010 portable readers.

The system has been installed and is successfully operated now at some
RF Minatom's enterprises. New versions for other enterprises are being
developed.

**THE LIBRARY OF SOLVERS FOR SPARSE LINEAR SYSTEMS. CAPABILITIES, USE,
DEVELOPMENT PROSPECTS**

***A.Yu.Artemiev, Yu.G.Bartenev, Yu.A.Bondarenko, V.A.Erzunov,***

***M.A.Kolesov, A.V.Lomtev, A.S.Maksimov, A.I.Panov, A.I.Prokofiev,***

***M.D.Romanova, A.P.Trubitsin, N.V.Frolova, T.V.Tsareva,***

***E.B.Shchannikova, Y.Dandass, D.E.Nielsen, A.Skjellum, V.Valsalam,
W.Wu***

(RFNC-VNIIEF, Sarov; USA)

The paper deals with the purpose of the library, its comparison with
other similar libraries; we present the description of its functional
capabilities.

1\. The collection of iterative methods is implemented to solve large
algebraic systems resulting from different approximations of
differential equations on regular and irregular grids, on sequential and
parallel distributed and shared memory computers.

2\. The collection of direct solvers is implemented for small symmetric
and asymmetric linear systems resulting from different approximations of
differential equations on regular and irregular grids, on sequential and
parallel shared memory computers.

3\. The tools are implemented for renumbering of unknowns, permutation
of equations, splitting a linear system into subsystems for more
efficient solving and parallelization.

4\. The above features are implemented as portable libraries with a
unique interface oriented to distributed and non-distributed data
representation that can be used within Windows and Linux operating
systems.

5\. An easy in operation and flexible interface was designed to access
these libraries from C++, C, and Fortran codes. Efficiently implemented
features of the libraries have been and are being introduced (as they
become ready) to existing programs.

6\. The tools are designed for numerical studies of linear system
examples.

Generally, the development of the parallel linear algebra library not
only facilitates the shakedown of high performance computer systems but
also allows to use the hardware and software of the system to be used in
the most efficient way.

The paper is based on those prepared for publication in VANT Review:

1\. Implemented methods for solving sparse linear systems in the linear
algebra library. Bartenev et al.

2\. Graphics shell and research tools of the library of sparse linear
system solvers. Tsareva et al.

**COMPARATIVE COMPUTATIONS FOR SOME**

**MODEL PROBLEMS ON GENERAL-CIRCUIT**

**THERMOHYDRAULIC CODES**

***D.A. Afremov, A.V. Kuzin, V.E. Men'shikov, G.S. Mingaleeva,***

***Yu.V. Mironov, M.Yu. Chernov***

(NIKIET, Moscow)

The paper presents the results of comparative calculations for some
simple model thermohydraulic problems using domestic and foreign codes.

The computations were with American codes MELCOR 1.8.3 and
RELAP5/mod.3.2, the German code ATHLET 1.1 cycle D and domestic general-
circuit codes RATEG, KORSAR and BAGIRA.

The following model problems have been considered:

-   Gravity separation of phases in originally homogeneous two-phase
    medium;

-   Testing mass balance conservation equations;

-   Testing departure nucleate boiling models of the general-circuit
    thermohydraulic codes;

-   Studying internal properties of computational models for boiling
    water and steam-water mixture used in the modern general-circuit
    thermohydraulic codes.

The computations were performed for a small sample from available
experiments, as at this step the principal goal was assessment of
quality of the hydrodynamic models allowing physical interpretation of
the computed data rather than complete verification of domestic and
foreign codes.

**A TECHNIQUE FOR statistical analysis of thermohydraulic COMPUTATION
uncertainties**

***Afremov D.A., Mironov Y.V., Radkevich V.E., Zhuravleva J.V.***

(N. A. Dollezhal Research and Development Institute

of Power Engineering NIKIET)

The technique for statistical analysis of thermohydraulic computation
uncertainties has been implemented by NIKIET. The technique makes it
possible to find the confidence interval for simulations and determine
areas of further research for efficient enhancement of accuracy of
computations.

In the general case, sources of computational errors are technological
deviations in facility structural elements, assumptions made in facility
nodalization, and errors of closing relations used in the code. Casual
deviations in accordance with the established distribution functions are
superimposed on all possible uncertainties. The necessary number of
variance calculations is determined using methods of non-parametric
statistics and depends on given confidence probabilities and analysis
reliability. The results of the variance calculations are used to
determine confidence intervals for the key parameters and the closing
relations whose dispersion makes the major contribution to the total
computational error.

The technique was used to analyze errors in RELAP5/MOD3.2 and MELCOR
1.8.3 simulations. The standard version of the RELAP5/MOD3.2 does not
allow introduction of correction factors, so separate subprograms of the
code have been appropriately corrected. Code MELCOR 1.8.3 has not been
modified as the varied parameters were accessible through input data.

**AUTOMATIC CORRECTION OF 2D MESH FRAGMENTS**

**IN D CODES WITH MIXED CELLS**

***R.A.Barabanov, V.I.Budnikov, O.I.Butnev, V.I.Delov,***

***O.K.Loginova, V.A.Pronin, V.V.Sadchikov***

(RFNC-VNIIEF, Sarov)

A method and code were previously developed /3/ for global correction of
Lagrangian grid fragments to improve the technology of computations with
D codes on modern high performance massively parallel computers in
particular to avoid alert conditions. The scaling of values to a new
grid is performed based on mass, momentum and internal energy
conservation laws. The computations with the routine for automatic
global correction of the grid fragments showed how promising the
approach is for increasing the level of abend-free computations with D
codes. This is especially seen when the grids with a great number of
cells are used and vortex flows occur in the domains to be computed.
Currently this method and routine are generalized and used for boundary
cells of computational domains. However only one material with the
greatest volume was put into the new cell containing several materials
when scaling the values. The frequent use of the scaling routine with
this approach could lead to serious material unbalance and the damage of
real material configurations.

Once the method of concentrations was implemented in 2D codes it is
possible to place several materials into the cell according to the
volume concentrations and proceed with the computations using mixed
cells. The next natural step was to develop the method and routine to
scale the values of the previous grid to the new one keeping in mind the
mixed cells in the grid fragments to be corrected. The scaling is based
on the introduction of embedded adaptive grids, that is each mixed cell
divides into smaller nonintersecting cells with only one material
specified in each of them. The finer is the adaptive grid in mixed
cells, the smaller is the material unbalance during the scaling; however
this increases the run time of the scaling routine. To avoid the mass
unbalance in the scaling process, the method is proposed for the
correction of the material masses. The test computations of some
problems showed that keeping in mind the mixed cells in the scaling
domains better retains the interface shape of the materials and
considerably reduces the mass unbalance.

We take as a test problem the instability evolution when the shock wave
passes through the sinusoidally perturbed interface of a heavy and light
gas. The main objective of the computations was to verify whether the D
codes are accurate and abend-free on this class of problems for the
nonlinear growth phase of perturbations. Since no analytical solutions
exist for this phase the computational results were compared with
similar data obtained with MEDUZA-P codes /4/. The comparison
demonstrated a good agreement of the results.

1\. Sofronov I.D., Dmitriev N.A., Dmitrieva L.V., Malinovskaya E.V. *The
computational method for 2D gas dynamics non-stationary problems in
Lagrange variables*. Keldysh Institute, *preprint* N 59, 1976.

2\. Sofronov I.D., Delov V.I., Dmitrieva L.V. *D codes for the
calculation of multidimensional solid mechanics problems in Lagrange
variables on regular grids. VANT. Series: Mat. Mod.Fiz. Proc*.1999. No.
4, p.42-50.

3\. Budnikov V.I., Delov V.I., Loginova O.K. *The method and routine for
the global automatic correction of the 2D grid inside the computational
domains in D codes. VANT. Series: Mat. Mod.Fiz. Proc*.2001. No. 3,
p.49-59.

4\. Barabanov R.A., Butnev O.I., Volkov S.G., Voronin B.L., Zhogov B.M.,
Pronin V.A., Sofronov I.D. *Parallelization of a 2D gas dynamics problem
on unstructured grids on the distributed memory MP-3 system. VANT.
Series: Mat. Mod.Fiz. Proc*.2000. No. 2, p.3-9.

**MEDUZA COMPUTATIONS OF INSTABILITY DEVELOPMENT AT GAS-GAS INTERFACE
WITH CONTACT LINE SEPARATION IN MIXED CELLS**

***R.A.Barabanov, O.I.Butnev, S.G.Volkov, B.M.Zhogov,***

***O.G.Loginova, V.A.Pronin, V.V.Sadchikov***

(RFNC-VNIIEF, Sarov)

Numerical simulations of Richtmyer-Meshkov instability development have
been attracting much attention of application mathematicians in the past
30 years. This is primarily because of the practical significance of
understanding of this phenomenon, especially in studying the initial
phase of the turbulent mixing development. Besides, simulations of the
nonlinear phase of the instability development are a good test for any
technique from the standpoint of its possibilities to describe severe
jet and vortex flows. The recent development of computer equipment and
parallel computation technology allowed us to transfer to a more
comprehensive study of the Richtmyer-Meshkov instability.

The irregular explicit free-Lagrangian code MEDUZA developed by
RFNC-VNIIEF Computational Division showed itself to advantage in
solution of gas dynamics problems with geometrically complex initial and
boundary conditions and with severe shear flows. The major features of
the code are specification of all grid quantities at cell centers, a
variable difference template for differential equation numerical
integration, and abilities to change the grid topology during the
solution. For one-domain model of problem solution, mixed cells appear
at media interfaces, the solution in which is sought using the
multi-material approach without explicit interface separation. During
the problem solution the neighborhood is re-specified on the basis of
the well-known Dirichlet principle (principle of metric proximity), some
points are discarded or new ones are added, and the grid functions are
re-calculated from the old mesh to the new one by their superposition.

The technique simulates the interfaces using mixed cells, i.e. the ones
composed of several (typically two) materials. The paper discusses the
algorithm for construction of separated contact lines in a mixed
polygonal nonconvex cell containing an arbitrary number of mixtures
using the volume concentration field of each material.

Results of test convergence computations for a plane shock wave incident
on a rectangular region filled with heavy gas are presented. The
obtained results are compared with experimental results and results of a
computation without the contact line separation in mixed cells. In
addition, the comparison is made to the results of the computations with
regular Lagrangian code D with automatic invocation of instruction
ISPVEL for global correction of a grid fragment using the method of
concentrations. The results of the computations are in a good agreement
both with each other and with the experimental results.

The code is parallelized by the following principles: geometry
decomposition with partial data overlap, using the library of MPI
communications, communication minimization, separation of point sets
adjacent to processor interfaces. (It is the computation of these points
that primarily allows making the computation and communication
concurrent through asynchronous communications.) The paper presents
measured parallelization efficiency for the problem of Richtmyer-Meshkov
instability development with the number of processors up to 32.

**SOFTWARE DEVELOPMENT AND OPERATION TECHNOLOGY AT VNIIEF COMPUTATIONAL
DIVISION**

***V.A.Batalov, V.F.Volgin, I.D.Sofronov***

(RFNC-VNIIEF, Sarov)

The computational division of VNIIEF calculates basic classified
applications. The full series of computations requires a wide spectrum
of codes for multidimensional non-stationary problems of computational
physics. These include gas dynamics, heat conduction, strength
consideration, transport of charged particles, reaction kinetics,
detonation etc.

The special nature of problems to be solved requires the development of
designated computational methods and computer codes. Since the early
1960s of the previous century all codes were written by the
mathematicians themselves, that is the methods were developed, the codes
were written and those were used to calculate all the items designed at
VNIIEF.

High requirements for computer codes were the point of honor of our
institution. The administration of the computational division introduced
the method and code development technology for massive production
computations. The paper examines each issue of the technology.

The required normative documents are changed , if necessary. For
example, an ensemble of documents was approved in 2001 similar to that
developed earlier accounting for more severe requirements for the
computational accuracy, halt-free runs, use of PCs, multiprocessor
systems.

The technology development allowed the teams to use and adopt the
experience of other teams This factor seemed to be the most effective
incentive for improving the software quality, the maximum advantage was
taken from the principles and modules, that is the valuable experience
gained by various development teams was widely spread.

References

1.  *Integrated software documentation system. National standards*. Nos.
    19-0001-77 through 19-0909-77.Moscow.1978

2.  Myers G. *Art of software testing*.1983.

3.  Lipayev V.V. *Software quality*. Moscow. 1983.

4.  Fox G. *Software and its development*. Moscow. 1985.

5.  Lipayev V.V. *Software development management*. Moscow. 1993.

6.  Bezhanova M.M., Kolker N.I. *Applications packages. Preparation and
    acceptance testing. Siberian Division of Russian Academy of
    Sciences*. Novosibirsk. 1987.

**RELAXATION MODEL FOR SHEAR STRENGTH**

**IN EULERIAN FORMULATION**

***S.M.Bakhrakh, I.P.Bakumova, P.N.Nizovtsev, E.V.Shuvalova***

(RFNC-VNIIEF, Sarov)

The paper presents a model describing the relaxation of the shear
strength both in Lagrangian and Eulerian formulation. For the
generalization to the flows in Eulerian variables, the rules are
established for scaling all quantities describing the strength. This
allowed calculation of a broad class of problems with the relaxation
model on Lagrangian-Eulerian and Eulerian grids.

The relaxation model yields good results in the pressure range 10-65 GPa
and deformation rates 10^2^ - 10^5^ s^\_1^ and describes experimentally
measured instability growth fairly well.

The model allows a correct description of the shock propagation. In
particular, a moderate intensity shock wave clearly demonstrates the
damping of the elastic precursor and progressive conversion of the
elastic precursor to the plastic wave is modeled.

**COMPUTATIONS FOR LOW-VELOCITY MODE**

**OF HE DETONATION DEVELOPMENT**

***S.M.Bakhrakh, N.A.Volodina, I.V.Kuz'mitsky, M.N.Leontyev,
K.V.Tsiberev***

(RFNC-VNIIEF, Sarov)

The paper presents results of 2D computations for octogen high explosive
initiation through impact by a steel spherical fragment and compares the
computations to VNIIEF experiments.

The computations were performed within kinetic model \"Ignition And
Growth\" implemented in program complex LEGAK. The kinetics had been
developed by Tarver and employees of US Lawrence Livermore National
Laboratory using Zeldovich-Neumann-During model.

HE and EP properties are described using equation of state JWL (after
the initial letters of authors Jones -- Wilkins -- Lee). At VNIIEF,
I.V.Kuz'mitsky was one of the first to use the EOS.

High explosive 12 cm in diameter, 7 cm in thickness was placed into a
steel casing. HE was loaded by a steel ball 1.43 cm in diameter
accelerated up to fixed belocity 1.6 km/s. The intensity of the action
on HE was varied by varying thickness of the cover, an aluminum plate,
that covered the HE.

The computations involved plates 0.13 cm, 0.14 cm, 0.15 cm, 0.2 cm, 0.3
cm, and 2.0 cm thick on Eulerian grid with 0.01-cm spacing. Thus, the
problem dimension was 1.1\~1.2 million computational cells. The
computations were performed on multiprocessor distributed-memory
computer.

The computations resulted in the following HE initiation threshold:
thickness of the Al plate that covered the HE was 0.14 cm.

The next computational step was estimation of average indicator plate
velocity on 0.2...0.7 gauge length. The right point of the plate surface
on the system axis of symmetry was used. For the 0.2-cm-thick HE
covering plate, the average plate velocity was \~0.36 km/s, which is in
quite a good agreement with the experimental result: 0.40 km/s; for the
0.3-cm-thick covering plate, this was \~0.35 km/s, which is also in a
good enough agreement with the experimental result: 0.39 km/s. For
2-cm-thick shield, the computation of the indicator plate yielded
\~0.17-0.19 km/s velocity, which is about two times as small as the
experimental indicator plate velocity of 0.38 km/s.

We can offer the following explanation for the difference in the
indicator plate velocity between the problem computation and the 2.0-cm
Al plate experiment. The detonation development pattern can be
significantly affected by HE elastic-plastic properties. (Recall that
the computations were conducted in the hydrodynamic approximation.) For
example, 1D computations which included HE elastic-plastic properties
demonstrated a significant increase in the indicator plate loading
pressure. Presently additional studies of the suggestion are underway.
On the whole it can be inferred that our computations for the
experiments within the kinetic model \"Ignition And Growth\" have
reproduced the pattern of constant indicator plate velocity on the
0.2...0.7 gauge length obtained in its experimental recording in the
pressure range, where the hydrodynamic approximation works with a good
accuracy.

**NUMERICAL SIMULATION**

**OF INITIAL PERTURBATION GROWTH**

**WITH OBLIQUE-ANGLE IMPACT OF METAL SLABS**

***S.M.Bakhrakh, N.A.Volodina, P.N.Nizovtsev, Spiridonov V.F.***

(RFNC-VNIIEF, Sarov)

Initial perturbation growth with oblique-angle impact of metal slabs is
simulated numerically. It is known that with supersonic impacts in
hydrodynamic approximation growth of initial perturbations is observed.

If strong slabs are impacted, the problem becomes more complicated.

When metal layers are impacted at an oblique angle, intensive shear
deformations will evolve over the contact area, near-interface layers of
metals will be heated to high temperatures and form shaped jets. These
effects cause distortions of the interface profile. Regular waves and
non-symmetric distorted waves occur.

Oblique-angle impacts under subsonic conditions, i.e.
U~points\ of\ contact~ \< C~sound~, have been studied in detail
recently. Under such loading conditions, a shaped jet continuously forms
at the impact point, if pressure in the vicinity of this point exceeds
the metal strength.

If C~sound~\<U~point\ of\ contact~\<U~critical~ , detached oblique shock
waves occur and perturbations grow, as in case of subsonic impact.

With U~point\ of\ contact~\>U~critical~ attached oblique shock waves
settle at the contact point and jet formation is not observed. Earlier,
evolution of perturbations in such a case was thought to be impossible
because of the absence of their main generator, the shaped jet.

One of the VNIIEF divisions performed experiments on impacting aluminum
slabs. It was expected to obtain an essentially prompt interruption of
perturbation growth when passing to jet-free conditions of slab impacts.
That is, as soon as U~point\ of\ contact~=U~critical~, no perturbations
should appear. It was observed during experiments, however, that
amplitude of perturbations decreased monotonically as the area of the
contact increased further. Such data have been obtained for the first
time.

States of material interfaces under such supersonic loading conditions
are poorly studied. So, we performed numerical simulation of experiments
with varying impact conditions for strong slabs. Computational and
experimental data obtained are in a good agreement.

**CRASH-FREE TECHNOLOGY FOR CONTINUUM FLOW COMPUTATIONS USING \"LEGAK\"
CODE ON MULTYPROCESSOR COMPUTER SYSTEMS**

# ***S.M.Bakhrakh, S.V.Velichko, V.F.Spiridonov***

## (RFNC-VNIIEF, Sarov)

The paper describes major principles of LEGAK techniques and LEGAK code
complex intended for computations of time-dependent continuum flows. The
ideas of the code system implementation on distributed-memory
multiprocessor systems are discussed. Data on the parallelization
performance and crash resistance for 2D problem computations are given.

**OpenMP Fortran EXTENSION**

**FOR DISTRIBUTED-MEMORY SYSTEMS**

######### **V.A. Bakhtin, N.A. Konovalov, V.A. Kryukov**

(MSU, Moscow)

####### OpenMP Standard is an extension of low-level shared-memory parallelism models to the Fortran 77, Fortran 90, C and C++ languages. OpenMP is built-in into the programming languages in the form of some directives describing the program parallelism. An unquestionable merit of OpenMP is the high-level parallelism model that ensures portability and sufficient efficiency of parallel programs for multiprocessor shared-memory systems.

####### A demerit of OpenMP is the lack of tools for data localization control. State-of-the-art massively parallel systems are a network of computational nodes. The multiprocessors can therewith be only nodes of the network, which significantly restricts the OpenMP applicability, as the program in this Standard can only be executed on one network node.

####### The report deals with the proposals on OpenMP extension to Fortran by adding directives for data allocation and access to remote data which would allow using OpenMP on distributed-memory systems.

**SPH CODE FOR SIMULATING 3D ELASTIC-PLASTIC FLOWS**

**USING THE METHOD OF SMOOTHED PARTICLES**

# ***V.V.Bashurov, N.V.Samodolova***

## (RFNC-VNIIEF, Sarov)

SPH code basing on the algorithm of smoothed particles is used to solve
time-dependent 3D elastic-plastic problems. The code uses algorithms for
calculating the continuity equation, equations of state and artificial
viscosity factors taken from MIMOZA code complex. The paper describes
verification of new types of artificial viscosity proposed in papers by
various authors via experiments with steel rods impacting a barrier in a
wide range of velocities (from 200 m/s to 3 km/s).

Computation results are compared to experimental data and results
obtained using other codes. Based on the results of numerical and
experimental data comparison, the scheme of approximating the continuity
equation has been updated.

Numerical computations show that the new artificial viscosity type gives
good enough results (as compared to experiments) only for problems with
subsonic velocities (about 0.2-0.5 km/s). At higher velocities (more
than 1 km/s), the new viscosity turns to be ineffective and causes
non-physical penetration of one material's particles into another
material.

Besides, a new algorithm of initial data computation has been used
during numerical simulations that converts regular data specified in
MIMOZA code complex into an irregular set of points used by SPH code.
The new algorithm significantly improves the symmetry and smoothness of
data mapped.

**CODES FOR VISUALIZING MULTIDIMENSIONAL DATA**

**ON ARBITRARY GRIDS**

***V.V.Bashurov, F.A.Pletenev, A.V.Ryzhikh***

# (RFNC-VNIIEF, Sarov)

Graphical processing is a basis for studying results of 2D and 3D
problem computations. Currently, it is difficult to imagine all
numerical simulation stages, including initial geometry computations and
initial data preprocessing, analysis of results during computation
process and making decisions concerning the tactics of computations, as
well as final processing of computation results and output data
representation, without graphical processing of data.

The paper presents ZOOM code developed on the base of algorithms by
VNIIEF Mathematics Division for visual representation of 2D and 3D data
determined using various grid types (both regular and irregular) and
data specified in separate nodes not connected by the grid.

Information about capabilities of the code developed under WINDOWS
environments is given, data viewing capabilities and formats of input
values are described.

The problem of processing results of numerical computations performed on
a multi-parallel supercomputer Argum operating under Linux is discussed.
Graphical processing using ARGUM control computer is offered as one of
possible approaches. Users obtain graphical processing results on their
personal computers using Exceed code package which allows a PC operating
under *Microsoft Windows* to have access to a Linux system and execute
X-applications.

The paper is concluded by the description of prospects for further
development of visualization codes using various platforms and under
various operating systems.

**A CODE FOR COMPUTATION OF POLISHING IN 3D GEOMETRY**

***V.V.Bashurov, V.P.Bashurin, O.A.Vinokurov***

(RFNC-VNIIEF, Sarov)

A code has been developed for computation of the polishing of complex 3D
profile surfaces with abrasive in a fluid flow produced by a rotating
glazer pressed with a given effort onto the surface to be polished.

The code includes a block for computing Navier-Stokes equations for the
fluid dynamic description, a block for computing the fluid-entrained
abrasive particle motion, a block for statistical computation of the
surface loading with the abrasive particles, and a block for mass
removal in multiple action of the abrasive particles on the surface
being polished.

The Navier-Stokes equations are computed by an implicit scheme on a
regular rectangular grid.

The abrasive particle motion computation includes all forces acting on
the particles, including Stokes force and elastic response forces of the
glazer and the surface being polished when they contact with the
abrasive particles.

The polished surface loading is computed with the method of statistical
averaging of the action performed by the model abrasive particles
randomly introduced to the space between the glazer and the surface.

The fatigue failure of loaded material and destroyed material removal
from the surface are calculated by the model of material defect buildup
under action of tensile stresses in excess of the threshold values.

**MIMOZA-ND SIMULATIONS FOR EXPERIMENTS WITH**

**A SPHERICAL BOX ON LASER FACILITY \"ISKRA-5\"**

***S.A.Belkov, C.V.Bondarenko, O.A.Vinokurov, G.G.Kochemasov,***

***L.S.Mkhitaryan, L.F.Potapkina, N.A.Ryabikina***

(RFNC-VNIIEF, Sarov)

The analysis and numerical simulations of the experimental studies on
hot and dense plasma physics on laser facility \"ISKRA-5\" required
development of 2D methods that would include features of powerful laser
radiation-material interaction, nonequilibrium X radiation generation
and propagation inside spherical boxes.

The ISKRA-5 experiments demonstrated that experimental and computed data
could be compared in such important areas, as:

1\. Effect of asymmetry in initial geometry of targets on their
compression and neutron generation.

2\. Effect of asymmetry of X radiation flux on the target compression
and neutron generation.

3\. X radiation propagation inside a spherical box.

The paper discusses setting up the experiments and their analysis,
development of and modifications to 2D programs of radiation gas
dynamics within complex MIMOZA-ND with using the spectral diffusion
approximation for X radiation transport in laser plasma. Issues of
parallel computations for this problem class are discussed.

Also, the paper compares data from ISKRA-5 experiments that studied
compression of asymmetric targets with thermonuclear DT fuel in
symmetric X-ray field as well as of symmetric targets in asymmetric
conditions proceeding in spherical converter box with numerical
simulations by 2D complex MIMOZA-ND. Results of simulations for X
radiation propagation in the spherical box of target LABIRINT are also
presented.

**A CODE FOR 3D HYDRODYNAMICS**

**WITH FINE-GRAIN POINTWISE PARALLELIZATION**

# ***S.P.Belyaev, L.I.Degtyarenko, A.L.Stadnik, Yu.V.Yanilkin***

## (RFNC-VNIIEF, Sarov)

TREK++ code complex that has been put into operation is intended for
solving complex 3D problems using parallel distributed-memory systems.

Its computational modules implement fine-grain pointwise irregular
parallelization that can provide highly efficient parallel computations.
So the complex is of a composite arrangement, where a parallel system's
architecture is taken into account not only at the level of service
programs (utilities). This is also true for computational modules having
their special arrangement as well that allows:

-   efficient combination of a problem computation by a processor and
    interprocessor communications;

-   selection of data transfer organization and size that are most
    efficient with respect to computer;

-   dynamic balancing.

The following methodical and production problems have been solved using
codes for multicomponent gas dynamics computations in Eulerian
variables:

-   a problem of a steel projectile penetrating a barrier made of steel
    that impacts the barrier at velocity u=4km/s and at angles 30°, 45°,
    60°, and 90°;

-   a problem of two gases with density ratio 3 mixed at a plane
    interface under constant gravity force (constant gravitational
    acceleration) that is solved by direct numerical simulation;

-   a spherical-symmetry point explosion on a rectangular grid.

All the computations performed show good enough accuracy of TREK++
Eulerian hydrodynamic code.

Measurements performed using 1 to 16 processors of the special computer
complex (SCC) show high enough performance (89.5%) of fine-grain
irregular pointwise parallelization implemented in TREK++ .

# **THE JET STREAM EFFECT IN PSEUDOSTATIONARY**

# **MACH REFLECTION**

***G. Ben-Dor, E.I. Vasilev, L.F. Henderson, and T.Elperin***

(Volgograd State University)

The jet stream effect is present in many phenomena of scientific
interest, such as the explosion of a shaped charge; the impact at small
angles of two explosively driven plates; the diffraction of shock waves
over ramps in gases.

The jetting effect in the Mach reflection of a shock wave is studied
under Euler basic equations. These equations were solved numerically
using a W‑modification of the Godunov method which has a second order
accuracy both in space and in time. Curvilinear moving grids and the
shock tracking procedure were used in the calculations.

It is shown that there is no correspondence between the type of Mach
reflection (single, transitional or double) and the occurrence of
jetting. Furthermore there are two kinds of jetting: *strong jetting*
that occurs when there is a branch point on the ramp surface where the
streamlines divide into an upstream moving *jet* and a downstream moving
*slug*; and *weak jetting* that has no branch point and may occur at
small and large values of the ramp angle. The width of the jet for
transitional- and double-Mach reflections is determined by the angle of
the Mach stem at the shock triple point (also called the Mach node or
3-shock node). Strong jetting is unstable and its primary instability is
in the jet itself. The main reason for the jet instability is that the
self-similar Mach number was nearly sonic in the inlet region of the
jet. The contact discontinuity is also unstable, but its instability is
secondary with respect to the jet instability. Two new types of shock
reflection patterns are identified in the dual-solution-domain. They are
a system comprising a 3-shock node followed by a 4-shock node, and
another which seems to be intermediate between the previous system and a
triple-Mach reflection. An approximate criterion for the
jetting ↔ no-jetting transition is presented. It is derived by an
analysis of the system of Euler equations for self-similar flow, and has
a simple geometrical interpretation. An adaptive moving grid allows us
to detect the structure of the interaction of the second Mach shock with
the contact discontinuity in a double-Mach reflection. This interaction
results in an alternating cascade of rarefaction and compression waves
with decreasing intensity along the contact discontinuity.

**COMPUTATION PARALLELIZATION USING MPI 1.3 MESSAGE**

**PASSING INTERFACE WITHIN INTELLIGENT**

**SOFTWARE SHELL ShIPRW FOR MODELING NUCLEAR REACTORS**

# 

# ***V.L.Bogdankevich***

## (RSC \"Kurchatov Institute\", Moscow)

Two basically different conditions are possible that occur with an
attempt to run several computational modules concurrently. The first of
them occurs when computational modules are different and independent on
each other. The second one consists in running the same module several
times. Both of them are readily supported by ShIPWR shell and standard
methods for solving the problem can be developed. The paper describes
methods for operation under such conditions implemented using WMPI 1.3
program libray.

**USE OF HIERARCHICAL GRAPH PARTITION METHODS FOR NUMERICAL SOLUTION OF
APPLICATION PROBLEMS USING**

**IRREGULAR GRIDS ON MULTIPROCESSOR SYSTEMS**

# ***S.N.Boldyrev, S.A.Sukov, M.V.Yakobovsky***

## (Institute of Computer Simulation, Moscow)

The method of geometric parallelism is a classical method for
parallelizing a lot of spatial-grid numerical algorithms, when a
computational domain is partitioned into a set of sub-domains depending
on the number of processors used in computations. If distributed-memory
systems are used, it is necessary to provide distribution of grid nodes
(computational cells) such that allows maximum possible even
distribution of the total amount of calculations required over
processors, on the one hand, and minimization of interprocessor data
communications, on the other hand, thereby reducing total run times for
each of the processors.

For irregular grids, the problem of rational decomposition of a
computational domain covering grid over processors is not a trivial one
and requires special methods. A graph is constructed for the combination
\"numerical algorithm-computational grid\"; the nodes of the graph
reflect amounts of calculations in various computational cells and its
edges reflect information communication between them, i.e. the necessity
of data exchange during the computation process. The problem formulated
is, hence, reduced to partitioning the graph into the specified number
of approximately equal portions (domains) with the requirement of
simultaneous minimization of the number of edges cut along adjacent
boundaries of those domains.

When solving the problem for large-size graphs, such as real
computational grids, it becomes reasonable to apply hierarchical
methods. With their use, one obtains a sequence of approximations to the
original graph in the order of descending sizes and only the last one
(i.e. the smallest graph) is subject to partitioning. The resultant
partition is then projected onto the original graph through the whole
sequence constructed and, thus, the desired solution to the problem is
found. In so doing, some small local changes in the form of domain
boundaries are made to specify more coarse partition of an approximated
graph and improve the quality of final partition of the original graph.

RAS Institute of Computer Simulation is developing a parallel library
implementing the above algorithms on the base of the widely used MPI
application interface. Besides, other special problems are solved that
kindle current interest with the number of system's processors up to
several hundreds.

When using systems consisting of several hundreds of processors and
computational grids of corresponding sizes, initial data loading from
their general array, as well as computation result recording take quite
much time. Besides, the data of sizes corresponding to the whole
computational domain may not be fitted into one system module's main
memory. In such case, there is a need in re-ordering original data
according to the partition obtained and in storing them separately. In
large parallel systems a situation is typical, where the number of
processors in use changes with each computation session and it is
necessary to re-distribute all the data for a short time.

The two-stage method of storing and processing large amounts of grid
data is proposed to overcome the above problems. First, the original
grid is subdivided into a large number of relatively small-size blocks,
the number of blocks has to significantly exceed a possible number of
processors used. Then data are re-ordered in compliance with such
division and will be stored separately in blocks. The second stage
consists in allocating to processors not all of the original set data,
but a sub-set of blocks of a considerably smaller size. Thus, a repeated
partition of the grid graph and its distribution over processors, as
well as access to data media during data reading/recording are performed
significantly faster.

PARALLELIZATION OVER DIRECTIONS FOR THE CALCULATION

OF THE 2D TRANSPORT EQUATION WITH SATURN-3 CODES

USING THE OpenMp INTERFACE

**A.I.Bochkov, V.A.Shumilin**

(RFNC-VNIIEF, Sarov)

The Saturn-3 package is designed for the numerical calculation of
stationary and non-stationary multidimensional transport problems on a
broad class of computers. Because of high dimension these problems take
much computational time. Saturn-3 implements various parallelization
algorithms to reduce the computational time.

The following takes place for the 2D non-stationary neutron transport
equation on a distributed memory computer using the MPI interface:

· parallelization over the energy groups

· parallelization over the computational domains

· fine grain parallelization for regular spatial grids.

For the 2D stationary neutron transport equation, the parallelization is
performed over the energy groups using both MPI and OpenMP.

OpenMP is a standard for programming on scalable symmetrical
multiprocessors (SMP systems) in the shared memory model. The OpenMP
standard includes the specifications of the compiler directives,
procedures and medium variables.

OpenMP is perfectly suitable for the developers who wish rapidly
parallelize their codes with large parallel cycles. The developer does
not generate a new parallel program, instead he simply adds OpenMP
directives to the sequential program text in a sequential manner.

The paper presents the method for parallelization over directions
implemented in Saturn-3 for the non-stationary neutron transport
equations using OpenMP interface.

The following functions were used as the parallelization efficiency
characteristics:

$S_{p} = \frac{t_{1}}{t_{p}}$ -- acceleration;
$E_{p} = \frac{t_{1}}{p \cdot t_{p}} \cdot \text{100}\text{\%}$ --
parallelization efficiency.

Depending on the problem formulation and the number of processors
activated, the parallelization efficiency of the 2D non-stationary
neutron transport equation varied from 50 to 100%.

# **TECHNOLOGIES FOR GENERATION**

# **OF 3D STRUCTURED GRIDS**

***T.N.Bronina, I.A.Gasilova, O.V.Ushakova***

(RAS UrB Institute of Mathematics and Mechanics, Ekaterinburg)

Technologies for generation of 3D structured grids in domains of complex
geometries are described. The technologies suggested are developed
within the approach from Ref./1/. The developed algorithms and codes for
grid computations are designed for numerically simulating 3D vortex
multicomponent media and solving other physical and engineering
problems.

The main idea of the approach has been presented in Refs./1/ and /2/.
Its main feature is the use of a special method for formalizing the
criterion of a grid similarity to uniform ones that provides the
smoothness of grids in combination with the orthogonality criterion,
realization of different boundary condition types for constructing grids
and a possibility to develop effective computational processes of grid
generation using discrete and variational formulations.

Within the approach /1/, we propose the algorithms for 3D structured
grid generation in domains resultant from the guide rotation about axis.
We propose the algorithms in which generation of a 3D grid is not
reduced to a 2D grid rotation about axis. When constructing 3D grids by
rotating 2Dgrids about axis (O-, C -type grids) some specific features
(degenerate cells) occur on the rotation axis, moreover these degenerate
cells become too small for small angles of rotation. Grids of such
quality are not desirable for the users, for this reason new algorithms
were suggested

The use of Lagrangian method of grid generation very often gives
severely distorted grids at some time points that are close to
degenerate ones. Further computations become possible only with the
improved quality of grids. For this purpose, we use the approach from
Ref./1/ for global reconstruction of a grid. The corresponding algorithm
is described.

Algorithms for grid generation in pipeline-like volumes with star-shaped
cross sections have been also developed. These algorithms and computer
codes automatically generate grids having no specific features inherent
in O-type grids.

To estimate the grid quality we use conditions from Ref./3,4/,
non-degeneracy criteria (Ref./5/) and criteria of optimality (Ref./1/),
as well as some other special quantitative characteristics of cell
sizes.

The corresponding codes have been developed to verify structured grids.
Examples of grids generated are presented.

This work has been financially supported by the Russian Foundation for
Fundamental Research, projects Nos.02-01-00236 and 00-15-96042.

References

1.  O.B. Khairullina, A.F. Sidorov, and O.V. Ushakova.

> *Variational methods of construction of optimal grids. Handbook of
> Grid* *Generation,* J.F.Thompson, B.K.Soni, and N.P.Weatherill, eds.,
> *CRC Press*, Boca Raton, FL, 1999, 36-1-36-25.

2.  A.F.Sidorov, O.B.Khairullina, A.F.Khairullin. *Parallel algorithms
    of generation of optimal multi-block-structed two-dimensional and
    three-dimensional grids of large size. Numerical Grid Generation in
    Computational Field Simulation, M.Cross, B.K.Soni and J.F.Thompson,
    J.Hauser, P.R.Eiseman eds*., ISGG, MS, 1998, 759-769.

3.  Ushakova O.V. *Conditions of nondegeneracy of three-dimensional
    cells. A formula of a volume of cells*. SIAM J. Sci. Comp, 23, 4,
    2001, 1273-1289.

4.  Ushakova O.V. *Nondegeneracy Conditions for Three-dimensional Cells
    and Formula for the Cell\'s Volume. Computational Mathematis and
    Mathematical Physics*. Vol.41, No. 6, 2001, pp.832-845.

5.  Ushakova O.V. *Nondegeneracy criteria for 3-D grid cells. Formulas
    for a cell volume. Grid Generation: New trends and applications in
    real-world simulations. Proceedings of the minisymposium at the
    International conference \"Optimization of finite-element
    approximations, splines and wavelets\"*. June 25-29, 2001.
    St.-Petersburg, Russia. Edited by S.A.Ivanenko, V.A.Garanzha,
    115-128.

**EFFICIENCY OF USING POWERFUL COMPUTERS**

**IN CONSTRUCTION OF COMPUTATIONAL METHODS**

**AND TREATMENT OF COMPUTED DATA FOR**

**ELECTROMAGNETIC WAVE IONOSPHERE REFLECTION WITH ACCOUNT FOR NONLINEAR
HEATING**

#### Yu.N. Brylev, N.F. Podderyugina, I.F. Podlivaev 

(Keldysh IAM)

### Setting up the problem of nonlinear interaction between electromagnetic wave and ionosphere plasma with inclusion of nonlinear heating. 

### Transformation of problem differential equations to the characteristic form convenient for development of computational algorithms.

Analysis of linearized local problem operator eigenelement properties
under different conditions using high-performance computers and
construction of an efficient computational method based on it.

Check for high-frequency components present in the reflected wave by
direct computation on a multiprocessor computer.

**COMBINED MODEL OpenMP & MPI FOR PARALLEL**

**COMPUTATIONS IN GAS DYNAMICS PROGRAM**

***A.N.Bykov, A.S.Zhdanov***

(FNC-VNIIEF, Sarov)

The combined model for parallel computations is studied by the example
of 3D gas dynamics problem of ellipsoid expansion. (The problem
dimension is 96 points in each of the directions.) A MPI interface base
program is used as a basis.

The testing of the parallel program in cluster computer system showed
that its efficiency dropped with increasing number of involved
processors. The efficiency drop is significantly slower, when one
computational node processor is used.

In the OpenMP&MPI program the computational work is distributed at the
procedure level. The algorithm for most procedures constitutes 3D
computational cycles.

The programs based on MPI and OpenMP&MPI have been tested on different
number of processors. Influence of hardware overheads on the parallel
program efficiency in the shared-memory operation has been determined.

The results of the effort are as follows:

1.  A technique for OpenMP interface base parallel computations is
    determined.

2 . Reasons for the overhead increase in operation of processes using
shared memory are studied.

3\. MPI and OpenMP& MPI base program speed of operation is comparatively
estimated.

4\. The practicability of using tools for automatic transformation of
sequential algorithms to parallel ones on base of OpenMP is
demonstrated. Possibilities of automatic parallelization tools Visual
KAP are studied.

5\. The mode OpenMP does not significantly change the program algorithm.

6\. The combined model OpenMP& MPI does not automatically save
significant time as compared to MPI.

7\. The principal reason for the decrease in actual program speed of
operation in the cluster multiprocessor system is hardware overheads in
the shared-memory operation.

**PROBLEMS OF LIQUID POLLUTIONS IN SOILS**

M.A.Vagin, A.V.Kosterin, K.A.Potashev

# (Kazan State University)

Problem 1. The task is to find the shape of a pollution jet propagating
under the force of gravity from a poor prolonged-action source on ground
surface to ground water table (GWT), as well as the field of flow
velocities in the jet.

Conventional models give a jet of a constant cross-section. However,
experiments show its significant enlargement. Consideration of this
effect is of great importance to the actual area of GWT pollution to be
determined and correct boundary conditions to be specified for the
problem of further migration of a groundwater pollutant.

A new macro-equation is proposed for the vertical jet pulse that
describes fluid leakage over local uncertainties of a porous array. A
self-similar problem of a jet from a point source has been formulated
and solved.

Problem 2. The process of absorbing a surface \"pool\" of liquid
hydrocarbons (HC) spilled during an accident by highly moistened soil is
considered. Experiments with soil columns carried out at the Kazan\'
University soil physics research laboratory show that some time later
(in two weeks) hydrocarbon penetration into soil completely ceases. To
give explanations to this surprising effect, it is necessary to extend
conventional flow models.

It is proposed to relate this \"protective effect\" to swelling of the
organic part of soil because of its interaction with hydrocarbons. The
corresponding computational model has been developed. A one-dimensional
problem of vertical HC migration through highly moistened soil has been
formulated and solved. Theoretic results are in a good agreement with
experimental data.

**FLOW FORMATION AND DEVELOPMENT**

**IN A SHOCK TUBE AT DIAPHRAGM OPENING**

**USING A LEAF**

***E.I.Vasilyev, D.V.Tsitser***

(Volgograd State University)

A method of the second order of accuracy was used for a numerical
simulation of non-stationary 2D flow of perfect gas, which originates in
a shock tube at diaphragm opening ensuring initial pressure
differential.

A study is made within three different models of diaphragm opening. The
first model implies opening of the diaphragm using a gate valve
perpendicular to the tube walls. The second model considers opening of
the diaphragm using a rigid leaf with a given law of motion. The third
model determines dynamics of inertial leaf diaphragm opening by pressure
difference on the inside and on the outside. The system of equations for
leaf motion derived herein is fixed. An implicit difference scheme of
the second order is used for its integration. Gas flow is calculated on
a curvilinear movable grid connected with the diaphragm. At the final
stage of opening the difference grid is reconstructed to adequately
describe the process of a leaf approaching to the wall.

A comparative analysis of computational data was made within above
models. It was found that the diaphragm opening method strongly affects
flow behavior both on the stage of opening and during long period after
opening. It is shown that within the model of leaf opening the
experimental effect of initial shock wave acceleration at large initial
pressure differences (\> 1000) is fully validated.

**AUTO-OSCILLATORY MODES OF SUPERSONIC FLOW ABOUT A DOUBLE WEDGE**

E.I.Vasilyev, A.V.Zenovich, A.V.Karpov

(Volgograd State University)

A plane supersonic flow about a double wedge (a wedge with an additional
inflection) by ideal perfect gas in the parameter range, where added
oblique shock waves emanate from the wedge vertex and from the
inflection, is considered. At high Mach numbers of the impinging flow
(M\>5) several interaction configurations of the above waves appear
which are highly dependent on the angle of the additional wedge
inflection.

The numerical simulation is performed with the method of establishment
using W modification of the Godunov method, with the modification being
of the second accuracy order in space and time. Two approaches were
adopted in the computations. The first used the two-grid
integrated-computation algorithm without shock wave separation and the
second employed curvilinear movable grids with separation of principal
shock fronts and triple points.

Most computational versions are in a good agreement with the results of
ref. /1/ obtained with another technique of the first order in time.
However, in contrast to /1/, stable auto-oscillatory modes of flow were
detected in a narrow angle range relating to a complex nature of the
interaction between several Mach fronts and subsonic flow zones behind
them. The auto-oscillatory modes are reliably detected by either
computational method: both with discontinuity separation and in the
integrated computation. A computational series revealed that the
auto-oscillatory modes were present in the entire range of impinging
flow Mach number range, 5\<M\<9, under study.

1.  Oleiniczak J., Wright M.J., Candler G.V. *Numerical study of
    inviscid shock interactions on double-wedge geometries*// J. Fluid
    Mech., 1997. Vol. 352, pp. 1-25.

**A CODE SIMULATING FRACTURE GROWTH**

**IN ELASTIC-PLASTIC MEDIUM**

# ***B.A.Voinov***

## (RFNC-VNIIEF, Sarov)

A code has been developed for fracture propagation computations in an
arbitrary 2D geometry for an elastic-plastic medium. Growth of fractures
is initiated by thermal and mechanical loads. The code has a built-in
generator of triangular grids of random complexity and a capability to
describe original geometry in a C-like metasyntax language. The code
basis is a module for solving a linear elasticity problem in finite
element approximation to elastic medium. Plastic deformation is
considered in plastic flow approximation. The code also has a built-in
module for solving heat conduction problems. Fractures are described
using two models: volume and surface models. For representing the
former, the method of enriched finite elements is used. For the latter,
there is a special model of mesh-aligned fractures. The number of any
fracture types in a problem is essentially unrestricted. The code allows
simulation of fracture interactions, conversion and branching.

The linear elasticity problem is solved using the direct method of
multi-front pseudo-elimination for sparse matrices. The method speedup
exceeds the speedup of the iterative method of conjugate gradients by a
factor of 10 to 30. To solve nonlinear equations with plastic
deformation consideration, the method of initial rigidity is used. The
fracture growth criterion is based on calculations of elastic energy
fluxes into a fracture top (Rice-Cherepanov integral).

The code also has capabilities for cyclic load simulation and evaluation
of fatigue failure probabilities basing on the statistical model.

**METHODS FOR STABILIZING PROGRAM MOTIONS**

# ***E.V.Voskresensky, T.F.Mamedova***

(N.P.Ogarev Mordovian State University, Saransk)

Let the equation of motion be of the form:

$\frac{\text{dy}}{\text{dt}} = f(t,y,v),\text{  }v \in K,\text{  }y \in R^{n},$
(1)

where $K$ is a class of admissible controls from the set of measurable
vector functions: $v:\lbrack T, + \infty) \rightarrow R^{m}$ and
$f\left( t,y,v(t) \right) \equiv F(t,y)$ is of Caratheodory type, so
that Cauchy problem (the initial value
problem)$\left( t_{0},y_{0} \right)$ has a single solution $y(t)$ in the
class of absolutely continuous vector functions with any *v*. Problems
of program motion control, construction and stabilization for a wide
range of equations (1) were solved by many authors, such as Zubov,
Pontryagin, Alekseyev, et al. The report presents a new concept of
stabilization with the program motion quality. In this case, Lyapunov
stability only (without asymptotic stability) is required. Conditions
sufficient to provide stabilization are presented for the equation of
motion, asymptotic features of solutions are examined.

**NUMERICAL ALGORITHM FOR THE CALCULATION OF 3D EQUATIONS DESCRIBING THE
PROPAGATION OF MAGNETO HYDRODYNAMIC WAVES IN IONOSPHERE**

***K.G.Gainullin, V.A.Zhmailo, Yu.F.Kirianov***

(RFNC-VNIIEF, Sarov)

Plasma dynamics equations are used for a wide spectrum of issues
involving the ion sphere dynamics such as the effects of the generation
and propagation of magneto hydrodynamic waves considering the impact of
the geomagnetic field. The derivation of the equations used and the
applicability conditions are described in reference /1/. The major role
for the analysis of the propagation of perturbations across the plasma
is played by Maxwell equations for the electromagnetic field, the
equation for the dynamics of the quasi-neutral ion sphere plasma and the
analog of Ohm law.

The paper presents the system of equations for the propagation of MHD
waves and the numerical algorithm for the calculations in 3D Cartesian
coordinates. In the numerical solution, the complicated nonlinear
system, except for the intensity vector H, conduction current density j
and plasma velocity v, is reduced to the second order equation for the
magnetic field intensity. To approximate the resulting 3D system with
respect to the projections of the vector E, a difference scheme was
developed based on the method of transversal and longitudinal directions
similar to the well known Douglas scheme /2/. The scheme is implicit in
each space variable with the weighting factor θ.

The numerical method was tested on a 1D problem given in /3/, on a 3D
problem having the analytical solution. The computational results for
the axisymmetrical problem obtained with the 3D method were compared
with the 2D data.

The test results are presented demonstrating the operability and the
efficiency of the numerical method.

References

1.  Sorokin V.M., Fedorovich G.V. *Physics of low speed MHD waves in
    ionospheric plasma*. Moscow, *Energoatomizdat*, 1982.

2.  Marchuk G.I. *Dissociation methods*. Moscow, *Nauka*, 1988.

3.  Carl Greifinger and Phyllys Greifinger. *Transmission of
    micropulsations through the lower ionosphere. Journal of Geophysical
    research*, 70, No. 9, 1965, 2217-2231.

**NUMERICAL SIMULATION OF GRAVITATIONAL MIXING USING MULTIPROCESSOR
COMPUTER SYSTEMS**

# ***S.M.Garina, N.V.Zmitrenko, M.E.Ladonkina, N.G.Proncheva, V.F.Tishkin***

## (RAS Institute for Computer Simulation, Moscow)

The paper presents the problem of Rayleigh-Taylor instability in its
traditional formulation: evolution of 3D, in general case, random
perturbations of a two-medium interface with constant gravitational
acceleration and various values of a parameter that determines the
medium compressibility effect is studied numerically. Computations have
been performed using both a single-processor computer and a
multiprocessor cluster of PCs.

Based on the computations performed, the effect of changes both in
physical (compressibility, initial data form) and computational
conditions (a code with computation parallelization, a difference scheme
of the second or the first order of time, etc.) has been analyzed. The
differences observed have been classified depending on the problem setup
and implementation form. It turns out that variations in numerical
implementation of the computation experiment cause no significant
changes in its result that still depends only on physical conditions
formulated in the problem.

**GENERALIZED ASSIGNMENT PROBLEM**

######### **A.E. Godin**

(MSU, Moscow)

####### The assignment problem consists in the following: in a quadratic matrix of *n* size one has to select *n* elements having a minimal sum provided that one element has been only selected from each column and each row. There are some variations of the assignment problem: a minimax assignment problem, an assignment problem on a rectangular matrix of *m* x *n* size (at a given *m\<n* one should first select *m* columns and then select elements from them), a two-level assignment problem etc. The problem has been widely used in practice and provoked interest since the 1920s. When solving the assignment problem, a technique of Hungarian method (a method of alternate chains) widely used in the problems of graph theory was developed. Algorithms with complexity on the order of *n^2.5^*ln *n* are available for solving this problem.

####### The paper suggests a natural generalization of the assignment problem: uniqueness of element selection from each row and column is not required, hence, the restriction on the matrix sizes is automatically removed. The notions of a cover (an admissible solution) and a deadlock cover (minimum in inclusion) are introduced. Theorems on classification of deadlock covers and reduction of a deadlock cover to a single canonical form are proved. Amount of deadlock covers is calculated, which demonstrates that a solution by item-by-item examination makes no sense because of immense computational complexity. 

The problem under consideration can be reduced to that of search for a
minimal cover of a complete bipartite weighted graph, for which exact
solution search algorithms with complexity of the order of *n^3^* are
available. However for problems of high dimensionality, when the number
*n* is comparable to the number of image points, such complexity can be
unacceptable. The paper presents an algorithm that includes the problem
structure and operates for quadratic time with an assured error of no
more than by a factor of two. Heuristic extensions of this algorithm
(working also for quadratic time) give more accurate results in many
cases.

References

1.  Dinits E., Konrod N. *An algorithm for solving the assignment
    problem* // *Dokl. AN SSSR*. Vol. 169, No. 1.

2.  Mainika E. *Optimization algorithms on nets and graphs*. M.: Mir,
    1981.

3.  Serdukov A. *On the problem of a cover // Upravlyaemye sistemy IM SO
    AN SSSR*. 1975. No. 14.

**CALCULATIONS OF LARGE ELASTIC-PLASTIC**

**DEFORMATIONS IN RATES WITH ADDITIVE AND**

######## MULTICATIVE DEFORMATION EXPANSIONS

######### **A.I. Golovanov, S.A. Kuznetsov**

(Kazan State University)

####### A problem of large continuum elastic-plastic deformations in rates is formulated (within Eulerian approach). A multiplicative decomposition of the deformation gradients to elastic and inelastic (plastic) components is used. Finger deformation and its rate are used as a measure for elastic deformations, while inelastic deformations are described using plastic deformation rate strain. A relation is derived that relates the Finger measure rate to total and plastic deformation rates. The determining relations are taken as:

-   a relation between the Finger measure and true Cauchy stresses for
    the elastic deformations;

-   as an expression for the deformation rate as a product of the
    quantity parameter and the tensor of their development "direction"
    (in the deformation space) for the plastic deformations.

Resolvent equations are derived from the equation of virtual powers by
its differentiation with respect to \"time\". The finite-element
procedure is used. The solution methodology consists in step-by-step
loading and construction of equations for node coordinate increments
(transfers from the current state to the next one) for the current
strained state. The residue from the current state disbalance is
introduced to the resolvent equation for the solution correction at the
loading step.

**A SCHEME OF SPLITTING FOR NUMERICAL SOLUTION OF**

**THE VLASOV KINETIC EQUATION**

# ***A.I.Golubev, T.G.Sysoyeva***

## (RFNC-VNIIEF, Sarov)

The paper presents a finite difference method for solving the Vlasov
equation using a new scheme of splitting, including Maxwell equations
splitting for electric and magnetic fields' components. In so doing, the
scheme of interpolation by cubic splines is used to solve
one-dimensional transport equations, when considering the effect of
spatial gradients on the distribution function. At the step of
considering the effect of electric field's components on the
distribution function, a monotone nonlinear scheme by A.E.Aloyan and
V.P.Dymnikov is used. For consideration of the magnetic field effects, a
special interpolation scheme is used that eliminates its operation
during plasma motion.

The scheme serviceability is demonstrated for the problem of modeling
Landau attenuation in heated collisionless plasma.

**THE INFORMATION AND REFERENCE SYSTEM**

**FOR NUCLEAR DATA NDX**

***A.N.Grebennikov, N.A.Krut\'ko, G.G.Farafontov***

(RFNC-VNIIEF, Sarov)

Currently, a fairly wide spectrum of sources is available and used
publishing in the electronic form various nuclear data in the area of
constant support of nuclear physics computations. In essence, the data
may be the results of direct measurements and estimates obtained from
various theoretical models. IAEA distributes in electronic form
individual files describing the atomic nucleus structure and the entire
libraries of estimated constants containing information about the
particle-nucleus interaction. Most existing libraries of estimated
constants represent data in the format of the US National Library ENDF.
This format has become international because of simplicity and
uniqueness.

Because of the importance of the estimated nuclear data as the most
significant part of the constant support a serious effort was made to
organize the structural storage of such data, develop access methods and
visualize the data in the graphics form. The above-mentioned goals have
been accomplished in the main in the archive of estimated data BEND
operated by manager BMU. Currently the BEND archive stores 13 libraries
of estimated constants in ENDF format - more than 5000 isotopes with the
total size exceeding 500 Mbytes.

The practices of using the BEND archive allowed formulation of two key
requirements for future development of service software for the constant
support:

-   integration of various data sources (not limited by estimated data)
    within a single software shell;

-   development of the shell such as the information and reference
    system with flexible contents and amount of supplied databases to be
    oriented to the needs of various user categories.

These requirements just serve the base for the NDX program concept.

The NDX system is implemented within Win32 environment (Windows 95, 98
and Windows NT 4.0) in the form of two executable modules - NDX.EXE and
NDX_NT.DLL. The requirements for the hardware resources are as follows:

-   The main memory size at least 16 Mbytes (for Windows NT this
    requirement is assured);

-   Video mode not lower than High Color (65536 colors).

The NDX program proper is implemented in Borland Delphi 4.0 environment.
The module for construction and calculation of resonance integrals,
NDX_NT.DLL runs in Watcom Fortran 10.6 environment.

**OPTIMAL COOLING OF HEATED THIN-WALL SYSTEMS**

**IN THERMAL SHIELD BREAKDOWN**

***Guk N.A*., *Obodan N.I.***

(Dnepropetrovsk National University, Dnepropetrovsk)

The paper discusses behavior of an insulated thin-wall system (plate)
located in thermal field and experiencing cooling on its thermal shield
breakdown. Upon the action, the system is in the stressed-strained
state. The forced cooling results in the residual stress field in the
system which "additionally loads" the system and can lead both to a
local clap of the damaged part and to the plate cracking in areas far
from the breakdown point. To avoid the above effects, a cooling law that
optimizes the process in terms of strength and resistance must be
determined. The condition
$Q(T) \Rightarrow \text{min}\left\lbrack \int_{\Omega}^{}{\sigma_{\text{res}}^{2}(T(q(t))d\Omega} \right\rbrack$,
where T(q(t)) is thermal field determined from solution to the problem
of heat conduction in cooling q(t), is used for the optimality
criterion. σ~res~ is assumed to be a parametric function of time.

The choice of the optimality criterion in its less rigid form does not
restrict the σ~res~ and depends on problem conditions: the plate should
be cooled down to ambient temperature for safety reasons. The heat
extraction function, q(t), is a control function and determines
temperature field change. The problem is formulated in its discrete form
to determine the control function. For the discrete model of a body its
finite element representation is used, the initial state upon heating is
described by stress vector, unknown function of the problem are
discretized. It is necessary to determine vector q(t) such, that the
optimality condition be met.

The optimization problem is solved using the principle of maximum as
applied to the heat conduction problem in its discrete form. The
optimality criterion, Q(T), is calculated in solving the problem of
plate stressed-strained state, which, reasoning from the theorem of
unloading, is linear for known stress field, Σ~init~(T~0~). Then
$\Sigma_{\text{res}}(T) = \Sigma_{\text{init}}(T_{0}) - \Sigma(T)$,
where Σ~init~(T~0~) is a stress vector at discretization points of the
plate heated to temperature Т~0~, Σ(Т) is a stress vector determined by
the linear model from action of cooling flow q(t) and calculated with
the finite element computation, Σ~res~(Т) is a vector of residual
stresses upon cooling.

The formulated problem leads to a linear dependence of the Hamiltonian
functional set up to obtain optimal solution on the cooling control
function. The suggested algorithm for determination of switch points of
q(t) consists of Cauchy problem solution to determine function T(t) and
conjugate system with respect to ψ(t) and an iterative procedure for the
switch point specification. Examples of optimal functions q(t) with
given geometrical parameters of the plate and heat spots, heat emission
coefficients and values of T~0~ are given. The function q(t) is
piecewise continuous and sequentially takes on values equal to those of
the heat emission coefficients.

**COMPUTATION OF PHOTON ABSORPTION CROSS SECTIONS**

**IN MIXTURES USING DATA FOR COMPONENTS**

***N.N.Degtyarenko G.M.Yeliseyev***

(RFNC-VNIIEF, Sarov)

VNIIEF has developed program MIXER and interpolation technology MIXER
for computing matrices of photon absorption cross sections and ranges in
mixtures /1,2/. Thomas-Fermi model for intra-atom potential,
hydrogen-like wave functions of discrete states and quasi-classical wave
functions of continuous ones are used. Compton scattering of photons,
free-free, bound-free and bound-bound absorptions are included. Spline
approximation is used to calculate the cross sections and line widths.

By considering the Thomas-Fermi equation set, a rule of relating states
can be formulated. Let atom Z~i~ be in state (T, ρ). Also, let the atom
Z~i~ be a component of some mixture at (T, ρ~i~). All absorption cross
sections in processes for the atom Z~i~ both of pure material and within
mixture will be the same, if mixture atom partial density, ρ~i~, equals
to ρ~.~ In fact, in these states the chemical potentials are the same
and the atom Z~i~ has an identical intra-atom potential. Hence, energy
levels and their occupation numbers for the atom Z~i~ will also be the
same. By virtue of additivity cross sections for the mixture will be
equal to the sum of the cross sections for the components. In MIXER
technology, the cross sections are calculated using data for the atom
Z~i~ that are stored in some archive through interpolation with splines
/4/ in density to ρ~i~ rather than by complex quantum mechanics
formulas. This approach reduces computation time by several orders of
magnitude.

Although our model is fairly simple, average Rosseland values obtained
with PERST-3 and MIXER are close to the average values from ref. /5/ and
agree with the measured data /5,6/ for almost all reference points. The
PERST-3 spectral and mean Rosseland and Plank opacities also agree with
the data of LEDCOP computations /5/ downloaded by us from InterNet. Ref.
/1/ demonstrates that full absorption cross sections are a fractal of
dimension 1.5. The PERST-3 and MIXER computations were used to simulate
target NIF /7/ and others, study ion beam energy conversion /8/, and
solve other problems.

References

1.  Degtyarenko N.N., Yeliseyev G.M. Fractals in applied physics. Ed. by
    A.E.Dubinov. Arzamas-16, VNIIEF, 1995. P.175.

2.  Nikiforov A.F., Novikov V.G., Solomyannaya A.D. Teplofizika Vysokikh
    Temperatur, 1996, V.34, No.2, P.220.

3.  Yeliseyev G., Degtyarenko N., Zhil\'nikova N., GSI-Rep.97-08 (1997)
    p.72

4.  Yeliseyev G. Physical mechanics of heterogeneous media. Ed. by
    Gadiyak G. Novosibirsk, USSR Academy of Sciences Siberian Branch
    ITAM, 1984, p.146.

5.  Third International Opacity Workshop & Code Comparison Study. Final
    Report, ed.A.Rickert et al, MPI fur Quantenoptik, Garching, Germany,
    1995.

6.  Yeliseyev G.M., Model I.Sh., Narozhny A.T. et al. Dokl. AN
    SSSR. 1986. V. 289. No.4. P. 850-851.

7.  Antonenko E.M., Degtyarenko N.N., Dolgoleva G.V., Yeliseyev G.M.,
    Yermolovich V.F. VANT. Ser. Matematicheskoye Modelirovaniye
    Fizicheskikh Protsessov. No.1, 1999, p.45.

8.  Degtyarenko N.N., Dolgoleva G.V., Yeliseyev G.M., Yermolovich V.F.
    VANT. Ser. Matematicheskoye Modelirovaniye Fizicheskikh Protsessov.
    No.3, 1999, p.51.

**NUMERICAL SOLUTION OF 3D GROUNDWATER IMPURITY TRANSPORT PROBLEMS**

***Yu.N.Deryugin, A.A.Bazin, O.M.Velichko, N.F.Gavrilov, V.V.Gorev,***

***I.V.Gorev, G.N.Gubkova, D.K.Zelensky, A.I.Panov, E.I.Ryabov,
V.I.Selin, L.I.Sizova, B.P.Tikhomirov, A.V.Kosterin, A.N.Chekalin,
V.M.Entov***

# (RFNC-VNIIEF, Sarov; RAS IMM, Kazan; RAS IMP, Moscow)

The paper considers 3D problems of groundwater migration and transport
of soluble impurities from type local sources.

Equations of the model include flow equations for a
saturated-unsaturated porous medium and mass transport equations. The
problem as a whole is therefore naturally subdivided into two problems.
The first of them -- flow problem -- consists in finding a
time-dependent velocity field for a carrier flow. The second one -- mass
transport problem -- is to determine the dynamics of impurity
concentration propagation against its background.

A single irregular grid consisting of parallelepipeds with oblique-angle
quadrangles in their bases is used to solve flow and mass transport
problems numerically. For digitizing differential equations, the method
of finite differences is used. When approximating the flow equation, the
implicit method of operators is used, where the difference analog to
grad operator refer to the grid nodes and the conjugate div operator is
specified at centers of the grid cells. To numerically solve the
equation of water-soluble impurity migration, the idea of splitting into
physical processes is used. The entire process is split into four
stages: impurity transport computations, diffusion computations with
regard to dispersion, consideration of sorption, consideration of
chemical reactions in a porous volume and at fluid/array interface. For
the transport equation, both the method of characteristics and
difference methods, such as the flow correction technique (FCT) and TVD
method of the third order of accuracy are used. For diffusion and
dispersion, an implicit method is used which is similar to the method of
solving the flow equation. For sorption and chemical kinetics, both
explicit and implicit difference schemes are used.

The code is implemented in C++ language.

Results of test and methodological computations for a number of problems
are presented.

The work is performed under ISTC Project 1565 funding.

**NUMERICAL SIMULATION OF THE FLEGMATIZED OCTOGEN**

**SAMPLE RESPONSE TO THE SHOCK CONSIDERING**

**THE DETONATION KINETICS IN THE D CODES**

***L.V.Dmitrieva, I.I.Karpenko, N.V.Korepova, V.G.Morozov***

(RFNC-VNIIEF, Sarov)

The paper is devoted to the numerical simulation of the experiments
studying the response of the flegmatized octogen samples to the shock
with the D codes considering the detonation kinetics /1/.

The paper goal is to carry out the numerical experiment based on the
experimentally obtained initial data and taking the experimental results
as reference points to reproduce the possible evolution pattern inside
the HE using the fact that the computations allow us to visualize the
density, pressure and burn-out distribution inside HE at any time.

The computations employed the detonation kinetics model /1/ successfully
used in the computations of the TATB-based HE behavior. The model is
phenomenological and contains the parameters that are calibrated with
some characteristic experiments.

A series of 1D computations gave a self-consistent collection of kinetic
constants allowing numerical reproduction of the experimentally known
depth of the detonation excitation in the flegmatized octogen samples
for various pressure values at the boundary. Then the resulting
constants were used to calculate the response of the flegmatized octogen
samples to the steel ball impact using the detonation kinetics model.

For 2D computations, the thickness of the protection aluminum plate at
the outer HE surface was varied. Depending on the plate thickness,
stationary detonation mode, damping mode or transient mode without
stationary detonation but with a considerable HE portion burn-out were
implemented.

A detailed modeling was carried out for the transient detonation growth
or damping process in the flegmatized octogen under impact depending on
the thickness of the protection shield. The time was found when the
detonation source formed in the HE at a considerable depth from the
outer surface impacted by the ball together with further evolution and
transition to the stationary detonation or damping. The computational
results are presented for acceleration of a steel indicator plate in
three detonation regimes - stationary, damping and transient.

The computational geometry and the experimental assembly parameters
correspond to the experiment data.

The acceleration rates were calculated for the steel indicator plates
when three detonation regimes were initiated by the impact of a
spherical fragment on the flegmatized octogen. Various regimes occur
when the thickness of the aluminum shield changes which corresponds to
the experimentally measured values. HE burn-out functions are given.

References

1.  Morozov V.G., Karpenko I.I., Kuratov S.E., Sokolov S.S., Shamrayev
    B.N., Dmitireva K.V. *Theoretical basis for the phenomenology model
    of shock sensitivity of TATB-based HE. Khimicheskaya fizika*,
    vol.14, No.2-3,1995.

2.  Dmitriev N.A., Dmitieva L.V., Malinovskaya E.V., Sofronov I.D.
    *Computational method for non-stationary 2D gas-dynamic calculations
    in Lagrangian coordinates. In: Teoreticheskie Osnovy i
    Konstruirovanie Chislennykh Algoritmov Zadach Mat. Fiziki*. Edited
    by K.I. Babenko. Moscow, *Nauka Publishers*. 1979, p.175-200.

3.  Belov G.V., ZimalinV.S., Oleinikov V.A., Yamolkin E.M. *Explosion
    parameters of the flegmatized octogen samples under shock. 12-th
    Deflagration and Detonation Symposium*, Pt.III, p.21-23,
    Chernogolovka, 2000.

**INCLUSION OF BREMSSTRAHLUNG IN GAMMA-RAY**

**TRANSPORT PROBLEMS**

***E.N.Donskoy, A.N.Zalyalov***

(RFNC-VNIIEF, Sarov)

A new model TTBIAS is suggested for inclusion of bremsstrahlung in
problems of gamma-ray transport in 'thick' regions. The model implies
local absorption of electrons and positrons and emission of all
bremsstrahlung gamma quanta generated by them. The quantity, spectral
angular and spatial distribution of the bremsstrahlung gamma quanta
energy per primary gamma quantum interaction event is pre-computed with
code ELIZA and tabulated. In terms of the transport equation, the
suggested model reduces to addition of one more integral nucleus to the
right-hand side of the gamma-ray transport equation.

The major difference of the model from model TTB involved in Los Alamos
code MCNP is the correct account for the angular distribution of
bremsstrahlung quanta escape and spatial distribution of the quanta
generation points: the model TTB assumes that all bremsstrahlung quanta
are generated at the electron or positron absorption point and fly in
the direction of the primary gamma quantum flight.

Comparative computations have been conducted for radiation transport in
thick lead plates. With the model TTBIAS, the problem is shown to be
computed several times faster than the problem of simultaneous transport
of gamma radiation and charged particles, but several times slower than
the gamma radiation transport problem without inclusion of
bremsstrahlung. By the example of the lead plate it is found that the
forward gamma quanta releases from thick plates computed with the model
TTBIAS are in a good agreement with those computed with inclusion of
charged particle transport.

The comparison between the ELIZA computations and MCNP computations
using the model TTB indicates that in this case the MCNP computations
yield about 2 times higher forward gamma quanta releases and do not
conserve average energy. The obtained difference is a result of two
competing factors: release underestimation because of neglecting the
spatial distribution of the bremsstrahlung quanta generation points and
release overestimation because of neglecting the angular distribution of
bremsstrahlung quanta escape.

The problem of imaging a gamma-rayed hollow lead sphere has been
computed with the model TTBIAS. It has been demonstrated that in the
image formation the principal role is played by the bremsstrahlung
component that significantly reduces the image contrast.

**A PARALLEL ROUTINE FOR THE GLOBAL SCALING**

**OF GRID VALUES IN \"LEGAK\" CODES**

***S.P. Egorshin, T.V. Rezvova***

(RFNC-VNIIEF, Sarov)

The Legak-MP codes based on the finite-difference LEGAK method is
designed for the calculation of gas dynamic flows of non-homogeneous
medium in lagrangian-eulerian variables on a regular quadrangular grid
in scalar and multiprocessor modes.

To maintain the grid during the computations in lagrangian-eulerian
variables, the Legak codes use a computational mechanism consisting of
two sequential routines - the grid interpolation routine and gas dynamic
value interpolation routine over the convection flows. This mechanism
operates as follows. First the grid interpolation routine generates a
new grid from the previous grid resulting from the gas-dynamic step
using an appropriate user-specified law. Then the value interpolation
routine scales the gas-dynamic values from the previous grid to the new
one over the convection flows. The use of gas-dynamic values over the
convection flows in this scaling mechanism restricts the displacement of
the computational nodes run by the grid interpolation routine. The
mechanism composed of these two routines runs at each step.

Usually, one grid interpolation iteration per step is sufficient to
retain a good quality of the grid. Otherwise, the computations stop and
the full scale iteration process is generated for the grid interpolation
with the scaling of values. In some cases (because of restrictions of
the displacement of the computational nodes) the number of iterations
reaches several hundreds. As a result, in the event of large production
applications the grid recovery takes much time.

The reduction of the grid recovery time necessitated the algorithm
development without restrictions to the displacement of the
computational nodes allowing the global grid restructuring and the
scaling of gas-dynamic values to the new grid within a single iteration.

In this connection, we developed a new mechanism to introduce into the
Legak codes. This allows both the generation of a new grid and the
global scaling of grid values to the new grid in the parallel mode on
distributed multiprocessors.

The new grid is generated in the grid interpolation routine with
iteration-free types that allows any displacement of the computational
nodes (required to restore the grid quality). The scaling of gas-dynamic
values to the new grid uses the so called \"raster method\" based on
splitting of the previous grid cells into small elementary parts of the
raster - pixels- followed by the count of the number of pixels in the
previous grid cell transferred to the new grid cells. The raster method
does not restrict the displacements of the grid nodes and allows the
scaling of gas-dynamic values to any new grid generated by the
interpolation routine.

The parallel routine version implementing the above mechanism of global
grid restructuring and subsequent scaling of the grid values is
currently introduced into the LEGAK codes.

**DESIGN OF MULTIPROCESSOR SYSTEMS**

***G.S Elizarov, V.V. Karatanov, V.V. Korneev, V.K Levin***

(Kvant Research Institute, Moscow)

Supercomputers as a separate class of specially designed, manufactured
in small series and correspondingly very expensive computers originated
from the development of Cray-1. The term \"supercomputer\" existed
earlier but clear distinguishing of high performance machines such as
CDC 6600 and others emerged since the development of Cray-1 exceeding
the existing systems more than by one order of magnitude in performance.

The electronics progress was initially reflected by ever growing
performance single processor systems. When this resource exhausted 2-,
4- and more shared memory processor systems emerged (Cray 2, Cray 3).

However, it was found several years later that single-chip
microprocessor industry developing in parallel exceeded so much the
specially developed unique multichip processors for supercomputers in
cost-performance ratio due to the enormous demand and extensive
investments that it became economically unreasonable to follow the
traditional way. The next step in the supercomputer development was as
follows. The highest performance commercially successful single-chip
microprocessor served the base for the development of the high computing
performance module. In parallel, a high performance communications
network was under development for the data transfer between such
computational modules together with the software. Cray-T3 is a classical
example of such approach. This approach still dominates.

However, approximately since 1995 the so called cluster approach to he
production of high performance multiprocessors started forming and is
extensively developing. This approach relies upon a very extensive
progress in the development technologies of local area networks.

Within a short period, the data transfer over local area networks formed
clear standards, the throughput of the channels increased from tens
kilobits per second to several Gigabits, multi-input non-blocking packet
switches appeared, efficient libraries were developed and extensively
used for the development of parallel applications. And this is with
cheap hardware resulting from massive production.

The up-to-date microprocessors also rapidly progressed to provide the
maximum efficient interaction with local area networks. The transfer
rate with local area network through the PCI bus reaches 500 Mbytes per
second.

The most important factor was a fast promotion of Linux operating
system. This system is initially designed for the network use and
reduces the designer expenses as a commonly used code allowing the
corrections, if necessary.

Thus the typical cluster multiprocessor includes the following:

-   several PCs or workstations as computational modules;

-   local area network (for example, Fastethernet) connecting these
    modules;

-   Linux OS for each module;

-   high level languages like C, C++ and MPI library for the development
    of parallel applications.

The clusters including several tens of cheap PCs connected by
Fastethernet were interconnected and are successfully used at many
scientific centers. The performance of these clusters for individual
applications may approach the total performance of all processors.

The global enthusiasm in developing the cluster systems resulted in the
emergence of some companies that developed and started the commercial
production of high performance LANs (Myrinet, Giganet, Quadrix elan, SCI
and others) specially designed for the development of high performance
cluster multiprocessors and oriented to the use at the basic platforms
such as DEC Alpha, Intel Pentium, Motorola Power PC.

The development of clusters including several hundreds and later
thousands processors showed, however, that despite existing background
and potential capabilities of some software and hardware tools the
creation of a cluster efficiently performing at a wide range of
scientific and engineering applications comparable with the highest
performance Teraflop class supercomputers is a difficult scientific and
engineering challenge. The choice of appropriate engineering design
requires a long-term experience in the development and implementation of
labor-intensive application algorithms on multiprocessors, in design and
operation of such computing systems.

The balanced choice of the type and specifications of the computational
module, the type and structure of the module interaction network
structure may improve by several times the ratio cost-performance of the
cluster system for the applications in the field of interest.

Instead, the lack of experience of computer system and applications
algorithm designers may dramatically reduce the performance and
reliability of the computer system and result in a difficult and
expensive maintenance.

All of the above problems were successfully resolved by developing the
cluster multiprocessor MVS-1000M with the peak performance of 10^12^
double precision floating point operations per second..

**THE METHOD OF LOCAL SPACE-TIME DECOMPOSITION IN PROBLEMS OF SHOT LASER
PULSE-PLASMA INTERACTION SIMULATIONS**

***N.V. Elkina, V.D. Levchenko***

(RAS M.V. Keldysh Instute of Applied Mathematics, Moscow)

An original algorithm, LSTD (Local Space-Time Decomposition), has been
developed for parallelization of problems of kinetic plasma simulations
with the macroparticle (particle-in-cell) method.

The kinetic approach to the plasma simulations uses the numerical
solution to Vlasov\'s equation set (the case of high-temperature
rarefied plasma) supplemented with Maxwellian equations for a
self-consistent field.

The problems of kinetic simulations are well parallelizable within the
conventional spatial decomposition method. However, this technique
possesses no required scalability degree because of nonuniform loading
of processors and the need in computation synchronization at each
timestep.

The suggested algorithm of local time-space decomposition resolves these
problems through inclusion of data and computation localization in the
space-time domain relations. The physical basis of the method is
finiteness of interaction extension, i.e. speed of light in vacuum.

In LSTD, each grid cell is assigned the current timestep and the current
state characteristic of a phase of the algorithm to. The condition of a
cell can change (with a new algorithm substep being executed) only when
adjacent cells have already reached the previous condition. The concept
of calculus window (CW) for localization of calculus is introduced. The
CW is an area, in which the calculus is carried out until the calculus
is independent on cells adjacent to the given area.

The code was tested using an actual problem of super-powerful short
laser pulse-plasma interaction as an example. The results are compared
to those of a model problem (with the total number of macroparticles
being N~p~^tot^ =10^6^ on the computational grid 32x32x32), whose data
can completely fitted into memory of a single computer. The code
scalability was studied in range N~p~^tot^ =32 10^3^ -0.13 10^9^ on a
single computer (Duron/600, 512M RAM, 2x30G HD). The efficiency made
70-90%. For the upper limit, the total problem data size is 5 GB, which
is larger by an order of magnitude than the computer RAM size. Moreover,
an increase in standard test problem efficiency by 30% was observed when
major computations were transferred to the processor cache through a
change in the calculus window size. In the code testing on
multiprocessor systems with weak connection (MVS-1000/16), the
computation speedup parameters were tested. The same problem was
computed on a different number of processor nodes (from 1 to 15). 1D
spatial decomposition was used (in addition to 1D LSTD). In so doing the
computational complexity per node ranged as N~p~^tot^ =(2-32) 10^6^,
which requires M~run~=80 MB - 1.2 GB for local memory size of 256 MB per
node. The problem computation speed increases by a factor of 12 as
compared to the one-processor system. Time diagrams of the algorithm
operation on 15 nodes of parallel system were studied. In this
computation, the load of different processors changes by a factor of 3
or more, which is due to the perturbation zone (laser pulse) extension
in the system (from left to right). Although no methods for dynamic
balance of processor load were used, the efficiency dropped as little as
by 30%. This is because the processors remain independent during the
synchronization. The conclusion was also made that the parallel
algorithm efficiency increases with long synchronization times.

**IMPLEMENTATION OF THE MIXED PARALLELIZATION**

**BASED ON MPI AND OPEN MP STANDARDS FOR**

**2D STATIONARY TRANSPORT EQUATION**

***A.S.Zapasskaya***

(RFNC-VNIIEF, Sarov)

Currently a major problem of the calculations for critical parameter λ
is to reduce the calendar computation time. This problem is particularly
pressing for the calculations on irregular grids. To resolve this
problem, a parallelization algorithm was developed for the 2D stationary
transport equation.

The parallelization over energy groups was used as the parallelization
method. This algorithm is general purpose, it works well for any number
of computational domains and on both regular and irregular grids which
is especially important for the calculations of different application
types.

Also, within the development of the parallel version of the code for the
calculation of the critical parameters the parallelization over the
anisotropy groups of secondary neutrons scattering taking 20-30 per cent
of the computation time has been performed.

The parallelization was performed both on shared memory using MPI
standard and on distributed main memory with the Open MP.

The testing was run on a multiprocessor in the Computational Division.

**\"ISOMIM\", A CODE FOR GRAPHICALLY PROCESSING**

**2D COMPUTATION RESULTS IN \"MIMOZA\" COMPLEX**

# ***V.V.Zmushko***

## (RFNC-VNIIEF, Sarov)

One of the important tasks during numerical simulations is to process
data obtained and represent simulation results. MIMOZA code complex uses
ISOMIM code for this purpose. The code is written in C++ using MFC and
ISR (Intel) graphic library under OS WINDOWS. The code is capable to
represent results of 2D computations on regular grids in the form of
isochart, both in Cartesian and polar coordinates in different modes.
There is also a capability to obtain various plots and to export data
into spreadsheets.

The paper gives the format of input data used in ISOMIM and functions of
the library used in MIMOZA complex for data recording in the required
form. Some capabilities of analyzing and representing computation
results are demonstrated. Supplementary capabilities developed to
support usable comparison of results for several problems are presented.
The paper is illustrated by results of solving a 2D magnetic
hydrodynamics problem.

**INTEGRABILITY OF (SU(n),S(U(k~1~)×U(k~2~)×U(k~3~))) PAIR**

**IN SMALL DIMENSIONS**

***K.M.Zuyev***

(M.V.Lomonosov MSU, Moscow)

An important differential geometry problem is the description of
Riemannian varieties with integrable geodesic flows. On the one hand,
there are topological obstacles to the integrability that forbid
existence of the flows on varieties of quite a complex topology. As for
positive results, only a few varieties are known, which it was possible
to construct integrable geodesic flows on. Topologically, almost all of
the varieties are homogeneous spaces.

Consider a simplectic variety, (M^2n^,ω), and a Hamiltonian system with
Hamiltonian H(x) on it. The system is called Liouville-integrable, if it
has pairwise commuting integrals f~1~,...,f~n~ that are functionally
independent almost everywhere. By Liouville theorem, in this case
nonsingular compact compatible surfaces of the first integral level are
n-dimensional tori, with the dynamics on the tori being very simple: the
system trajectories in appropriate coordinates are rectilinear tori
windings.

However, a situation is of frequent occurrence where the system has a
redundant set of integrals that are already not commuting with each
other. This kind of problems can be tackled with the non-commutative
integrability theory developed by A.S.Mishchenko and A.T.Fomenko /1/. If
a system is integrable in the non-commutative manner, then each
connected compact component of the compatible surface of the
f~1~,...,f~n~ ∈ F function level will be diffeomorphic to a torus of a
smaller dimension, however the dynamics on the tori will be like in the
classic case /1/.

A.B.Bolsinov and B.Jovanovic /2/ proved that the biinvariant metric
geodesic flow is non-commutatively integrable on any homogeneous space М
= G / H of compact Lie group G. In so doing, if ϑ is full commutative
algebra on Ad~G~ orbits of elements of general position v∈V and Ψ is
full commutative subalgebra of Ad~H~--invariant polynomials on V, then
α(ϑ) + β(Ψ) is the full commutative algebra of integrals of the
biinvariant metric geodesic flow on M. Here V is a special subalgebra of
the Lie algebra g and α and β are some maps. There is a well-known
construction referred to as an argument shift that constructs full
commutative polynomial families on an arbitrary orbit of the compact Lie
group. Pair (G, H) is called integrable, if a described algebra Ψ
exists. A.B.Bolsinov's hypothesis is that pair
*(SU(n),S(U(k~1~)×U(k~2~)×U(k~3~)))* can be integrated with the special
algebra Ψ constructed with the generalized method of chains. We were
able to prove that statement for some small dimensions.

***Theorem.** Pair (SU(n),S(U(k~1~)×U(k~2~)×U(k~3~))), where k~1~ +
k~2~ + k~3~ = n and k~i~ ≤ \[n/2\], is integrable with the generalized
method of chains for n=3,4,5,6.*

To provide a complete proof of the A.B.Bolsinov's hypothesis,
apparently, it is necessary to study basis reconstructions that appear
when small dimensions are changed.

References

1.  A.S. Mishchenko, A.F. Fomenko *Generalized Liouville method for
    Hamiltonian system integration// Funktsionalny Analiz i ego
    Prilozheniya. 1978. V.12 No. 2, pp.46-56*

2.  A.V.Bolsinov, B.Jovanovic *Integrable geodesic flows on homogeneous
    spaces// Matematichesky Sbornik. 2001. V.192, No.7, pp.21-40*

3.  A.V.Bolsinov, B.Jovanovic *Homogeneous spaces with completely
    integrable riemannian and sub-riemannian geodesic flows.*

**BRANCH COMPUTER CODES FOR PROBABILISTIC ANALYSIS**

**OF SAFETY AND UNCERTAINTY OF DETERMINISTIC MODELS**

***R.T.Islamov***

(RF MAE Departmental Center for Computation Codes, Moscow)

The paper presents the review of efforts carried out at Minatom to
develop the computer codes.

For probabilistic analysis of the first level safety, the RISK code was
developed. This product is the first branch domestic code for
probabilistic safety analysis and was developed in accordance with the
task order approved by the Department of Atomic Energy. The code is
compatible with the commonly used Swedish code RISK SPECTRUM in initial
data writing format which allows the use of previously developed safety
models for some units of Leningrad, Novovoronezhsk, Kolsk and other
nuclear power plants to analyze the safety of other units. The code uses
the trees of events and failures. The code permits to quantify the
minimum cross-sections, the sequence of accidents, the consequences of
accidents. The failures caused by general reasons are modeled explicitly
on the trees of failures.

The computer code was developed for the analysis of deterministic model
uncertainty. This product is the first domestic branch code of this kind
and is intended for the analysis of model uncertainty in the codes under
development. The program is based on the stochastic approximation model
and determines both the discrepancy between the experiment and theory
quantitatively and identifies the parameters and characteristics of the
model that make the greatest contribution to the discrepancy. One of the
first codes for which such analysis was carried out (comparison of the
experiment and theory) was thermal hydraulic code KORSAR.

**THE PARALLEL PROGRAMMING SYSTEM BASED**

**ON TYPE ALRORITHMIC STRUCTUERS**[^1]

# ***T.Istomin, O.Panin***

## (MSU, Moscow)

As it is known, research into the area of high-performance high-level
tools for parallel code development is continuing nowadays. Experimental
and commercial systems basing on different approaches but, however,
having the common goal that is to simplify and automate the process of
developing a parallel code are created.

Our development, the parallel programming system based on type
algorithmic structures is one of the experiments in this area. The
approach chosen is based on the ideas of re-using design solutions (type
algorithmic structures). A type algorithmic structure (TAS) is a pattern
of some algorithm intended for solving a certain class of problems. To
use such pattern in a code, it is necessary to perform its adjustment
(parametrization) for a specific subject area.

There are several projects in the world that use the same approach: for
example, SkIE (a joint development by Quarics Supercomputers World Ltd.
and the University of Pisa, Italy), CO~2~P~3~S project (the Alberta
University, Canada).

From the user's viewpoint, the system is a graphic constructor of a
parallel algorithm. Various tools for parametrizing patterns,
establishing links between patterns, adding new patterns are offered the
user. Each entity has a corresponding graphic representation in the
constructor. By operating graphic items, the user creates an algorithm's
skeleton in terms of TASs, then the skeleton obtained has to be filled
out, i.e. some serial code (which is the specific one for the problem
being solved) has to be introduced at certain places. Each pattern can
be allocated to a group of processors having one of the integrated
topologies. Currently, structure scaling is performed manually, however,
a possibility of automating the procedure is examined.

When all the required links have been established, patterns have been
adjusted correctly and all "empty places" have been filled, it becomes
possible to automatically generate a low-level representation of the
user's code in a target language with calls to some communication
library. Generation of C++ code with calls to MPI is currently
supported.

The system has a standard set of type algorithmic structures, for
example: apply to all, reduction, calculations using a grid, a pipeline,
a serial composition. The standard set of patterns will be further
augmented. Besides, there is a mechanism of adding user's structures.

Further examination of feasibility of TAS-based code automatic
optimization and a debugger development is planned.

# **ON LOWERING THE COMPLEXITY OF DIGITAL IMAGES**

# 

# ***I.M.Kanayev***

(MSU, Moscow)

During the current stage of computer technology development, processing
of images that comprise 70% of the total amount of data transmitted via
Internet is the problem of great importance. To prevent data losses,
users employ, as a rule, images in their full-color format, though in
the majority of cases it is sufficient to operate a less complex version
(an image having few colors, or even a monochromic one). This is caused
by the fact that they are offered either multipurpose methods of
reducing the number of image colors (for example, PhotoShop) which are
applicable to all images equally, or an opportunity to perform
operations manually.

The paper considers the input-adjusted reduction in the number of
colors. A more precise formulation of the problem is as follows: for
each image, a projector onto a fixed subset must be found from the
ensemble of all colors such, that the projection retain the local
intelligibility, i.e. the number of situations where the spacing between
the neighboring colors is greater than the difference threshold before
projecting and is smaller than the difference threshold after
projecting. This problem reduces to that the optimal painting vertex
finds the weighted graph.

The author considers an important specific case of projecting onto a
two-element set {0,1}that fully demonstrates the problem complexity.

The following theorems occur:

**Theorem 1.** The complexity of the algorithm for solving the problem
using the method of exhaustive search equals $(n^{2} - n)2^{n - 3}$
additions, where n is the number of the original image colors.

**Statement 1.** For any k\>0, there is an image, such that a greedy
algorithm constructs a projector for it which quality can be improved by
no less than k.

The author proposes the algorithm of complexity
$(n - 1)2^{n - 2} - 2$for finding the exact solution to the problem. The
complexity of the exhaustive search algorithm is worse than the proposed
algorithm complexity even with n\>3; asymptotically, with n→∞$$the ratio
of complexities of the exhaustive search and the proposed one algorithms
is
$\frac{2^{n - 3}(n^{2} - n)}{2^{n - 3}2(n - 1) - 2} \simeq \frac{n}{2}$.

**Theorem 2.** The problem under consideration belongs to the class of
NP-full problems.

The author developed two heuristic algorithms for finding approximate
solutions to the problem -- \"the highest contrast ratio\" algorithm and
the \"credent\" algorithm -- and wrote codes implementing them. The
following theorem occurs:

**Theorem 3.** The complexity of \"the highest contrast ratio\"
algorithm doesn\'t exceed

$2n^{2} + 6n + 3$and is no less than
$\frac{n^{2}}{2} + \frac{5n}{2} - 1$, the «credent» algorithm complexity
doesn\'t exceed
$\frac{9(2s - 1)n^{2}}{2} - \frac{n(\text{10}s - 1)}{2} - 4s + 4$
additions and $\frac{(2s - 1)n^{2}}{2} - \frac{n(2s - 1)}{2}$ divisions,
with this estimate being unimprovable, here *s* is a confidence
threshold.

Both algorithms have been tested. In the majority of cases, there are
evident advantages of the proposed adaptable methods against well-known
multipurpose methods (not depending on input material).

##### References

1.  N.Bakhvalov, N.Zhidkov, G.Kobel\'kov. Numerical Methods. -- Basis
    Knowledge Laboratory, 2000.

2.  M.Harry, D.Johnson. Computing Machines and Stubborn Problems. --
    Moscow, 1982.

3.  H.Papadimitrius, K.Staiglitz. Combinatory Optimization. -- Мoscow,
    1985.

**EMPLOYMENT OF DISTRIBUTED COMPUTERS**

**FOR LARGE NUMBER FACTORIZATION**

D.V.Kovalenko, D.P.Sidorov, S.A.Fedosin

# (N.P.Ogarev Mordovian State University, Saransk)

The swift development of information technologies also intensively
stimulates the studies and development of methods and hardware for
problems, a distinctive feature of which is a large number of simple
basic operations and feasibility of parallel execution of the
operations. As solution of this problem type on conventional computers
takes much computation time, a great many of options of parallel
high-performance computer systems are suggested: from universal
multiprocessor systems to specialized problem-specific speeding-up
boards.

Presently, the high-performance computer system development proceeds
along four major lines: vector-pipeline supercomputers, SMP systems, MPP
systems, and cluster systems /1/. Distributed computer systems
constitute a separate category.

The notions of parallel computations and distributed computations have
to be differentiated. In the parallel computations, operations are
parallelized at a micro-level, i.e. at the level of elementary
operations. In the distributed computations this approach in
unacceptable, as the data transfer time can be longer than the time of
the computations themselves, when, for example, individual computers
communicate through Internet. Therefore, the distributed computations
are applicable only to problems that can be parallelized to several
large computational blocks operated essentially independently.

A problem requiring huge computational costs is that of factorization
(decomposition to simple factors) of large integers /2/.

The common interest in the problem of decomposition to simple factors
abruptly increased in 1977, when R. L. Rivest, A. Shamir, and L. Adleman
developed the first public-key ciphering system presently known as RSA.
In the system, the ciphering/deciphering procedures are performed on
different keys. Knowledge of the ciphering key alone does not allow
message deciphering, thus it is not a private cipher component and is
typically published, so that any person could send a ciphered message to
the owner of the key. An individual who is aware of the public key only
can not deduce the private key /3/. The cryptoimmunity of the system is
based on the fact that large integers are hard to decompose to simple
factors. As no efficient algorithm for the number decomposition to
simple factors has been found thus far, the method almost absolutely
ensures secret data and message protection in computer network.

As an illustration of the distributed computation technique, the authors
of this paper have developed a distributed computer system for the
factorization of large numbers. The suggested program has been developed
in environment Borland Delphi 5.0 for platform Win32. The program has
been implemented by the "client-server" technology, and the computation
parallelization idea is concurrent search for divisors of the original
number with different methods on different clients. (The number of the
clients is unlimited.) Several methods have to be used because of the
following. There are a variety of factorization methods, and for a
random number it is impossible to predict what of them will find the
divisor sooner /2/. Besides, the bulk of the methods find divisors,
which, in their turn, are compound, so the original task breaks down
into several subtasks.

We took three methods (trivial division method, Pollard rho-method,
Fermat method) for the approach demonstration, which, however, does not
restrict the consideration generality. To extend possibilities of the
program, it is necessary to make only minor modifications to the code,
that is, to add procedures implementing the needed methods.

It should be noted that the standard tools of the Borland Delphi 5.0
language allow manipulation of integer numbers only in a very limited
range. So, a library of functions supporting large integer operations
has been developed. As this part of the system is very critical to the
speed of operation, the relevant functions are written in the embedded
assembler of Borland Delphi 5.0.

Consider the principle of operation of the program. To start the
computation, some number N has to be specified on the server. Then the
server begins to search for small simple divisors to the original number
with using the trivial division method. For this purpose a table of
simple numbers which is stored on the server is used. There is a
possibility to add numbers to the table, and in this version it contains
all simple numbers less than 10000 (a total of 1229 numbers). Despite
its simplicity, the method is optimal from the standpoint of search for
small divisors. (It should be noted that the Fermat method is most
suitable for search for large divisors, while for the Pollard rho-method
a situation is possible, where it finds no divisors whatsoever, and then
it makes sense to run it with other initial parameters.) The more so as
the divisors do exist, the other methods, as a rule, find just them only
with a much larger number of operations. When the table of simple
numbers has been exhausted, the server forms a list of client tasks
using other factorization methods (it makes no sense to use the trivial
division on clients) and numbers remaining upon the division by all
found simple divisors (if any). If some number of clients have already
connected to the server, then the tasks are sent to them. Given no
clients, the server will be waiting for their connection. If the number
of the clients is less than that of the methods used for the
factorization, than the task first in the queue of the tasks will be
sent to the newly connected clients. When the number of clients is
higher than that of the factorization methods used, the server will
issue tasks to calculate the divisor of the original number with
different methods to first ***m*** clients, where ***m*** is the number
of the methods, and the other clients will be entered into the waiting
queue.

If the client has found some divisor ***p***, information about the
latter will be transferred to the server, where numbers ***p*** and
***N***/***p*** are checked for simplicity (for example, by Fermat test
or one of probabilistic tests). If any one of the numbers is simple, it
will be entered into the list of simple divisors to the original number,
otherwise it will be used as a basis to form a task complex for search
for divisors to that number. In so doing, operation of the clients which
continue the search for divisors to the original number with other
methods will be stopped and newly formed tasks will be assigned to them.

The client's operation can be stopped at any moment, then the client
sends the state vector to the server and the server resends the task to
the first free client, which continues the computation from the stop
point.

It should be noted that the above divisor search technique implemented
in the program under discussion is efficient only when the quantity of
the divisors to the number remaining upon the operation of the trivial
division algorithm is larger than two. Should the divisor quantity be
two, it makes no sense to run several different methods in parallel, as
each of them, most probably, will find the same divisor. However,
efficiency of each algorithm is known, and what of them will be the
first to find the divisor can be predicted in advance in this case.

A merit of the discussed program is that its execution requires no
high-speed network link, in contrast to the cluster technologies, as the
client and the server communicate data packages insignificant in their
sizes, and the time it takes for the client to compute is significantly
longer than the data passing time. Hence, there is a possibility of
computations via Internet, which seems feasible thanks to the TCP/IP
protocol in the data communication between the client and the server.

1.  Datsyuk V.N., Bukatov A.A., Zhegulo A.I. *Methodological manual on
    the course: \"Multiprocessor systems and parallel programming\"*.
    Rostov-na-Donu, 2000.

2.  Knut D. *The art of programming*. Volume 2. Seminumerical
    algorithms. The 3^rd^ edition. Translated from English. Manual.
    Moscow, Williams Publishing House, 2000. -- 832 p.

3.  A. Menezes, P. van Oorschot, S. *Vanstone Handbook of applied
    cryptography, CRC Press*, Inc., 1997.

**APPLICATION OF MSPS AND MSH METHODS IS**

**AN EFFICIENT AND USABLE WAY T0 PARALLELIZE**

**NUCLEAR COMPUTATIONS**

# ***A.A. Kovalishin, N.I Laletin, N.V.Sultanov***

(RSC "Kurchatov Institute", Moscow)

Currently, two computational methods are in use, as applied to reactor
practice. The first one refers to so-called design computations. They
allow reactor characteristics for some limited class of problems to be
obtained with a sufficient precision and for a short enough time.
Benchmark computations refer to the second group of computations. Here,
the precision of computation results is restricted only by uncertainties
of our knowledge of micro-sections. These computations, however, require
much time and, for this reason, they are used only to estimate the
precision of design codes for a number of states. With such an approach,
it is not easy to analyze even intermediate states (interpolation) and,
in particular, states outside the limits of benchmark computations
(extrapolation). Application of a serial approach using methods with
initial approximations satisfying the design computation requirements
that would lead to benchmark results of computations is more preferable.
We think such approach is available, it has been developed by one of the
paper authors and is being improved by our research group. The first
method is a method of surface pseudo-sources (MSPS) which is used to
calculate parameters of cells. The second one is a method of surface
harmonics (MSH) which is used both for blanket computations and
coarse-grid computations of a reactor as a whole. The paper briefly
describes the methods above and the ways of applying them to
computations of VVER-440 and VVER-1000 reactors, and PWR benchmark. Both
methods have in common the following stages: selection of suitable
surfaces allowing convenient parametrization of the problem to be
solved; specification of a set of parameters; counting sequences of
solutions in regions bound by the surfaces selected; formulation of
finite-difference equations which coefficients are represented via
boundary values of the corresponding partial solutions (trial
functions). The both methods above are serial methods of solving the
transport equation, in sense that any desirable precision of solution
can be achieved by increasing the approximation number.

**PARALLELIZATION OF COMPUTATIONS**

**IN DVM SYSTEM**

N.A.Konovalov, V.A.Kryukov, N.V.Podderyugina, Yu.L.Sazanov

(RAS M.V.Keldysh Institute of Applied Mathematics, Moscow)

**The system purpose and composition.** DVM system developed at RAS
M.V.Keldysh IAM is intended for parallelizing serial codes written in C
and Fortran languages. These programming languages have been extended by
special comments (Fortran) and special macros ( C ) added which are the
serial code parallelism specification directives. Such directives are
"invisible" for standard compilers. So, there is a single copy of an
original language code, both for serial and parallel execution of it.
The two languages above have one and the same model of code parallelism.
That's why DVM has unified subsystems supporting parallel execution,
debugging, analysis and prediction of the parallel code performance for
both languages.

**Model of parallelism.** This is a high-level user-oriented code
parallelization model for those who program computation tasks (discrete
models of computational physics) in C and Fortran languages. The
parallelization model provides for the following two levels of
parallelism: one-dimensional arrays of problems and multidimensional
parallel cycles.

When parallelizing a serial code, the user should make the following
efforts: determination of an array of virtual processors (AVP);
specification of tasks and description of an array of tasks; description
of data and computations mapping for each of the tasks; description of
mapping tasks over AVP (each task is mapped on an AVP section);
specification of parallel cycles and description of mapping turns of
cycles over processors (a section of turns of a cycle is mapped over a
single processor).

**Parallel execution model.** Each processor is loaded by one and the
same code. The code execution begins with the first statement on each of
the processors. If a directive \"the beginning of the region of tasks\"
is encountered, the execution process will be subdivided into N branches
(where N is the number of tasks). Each task is executed in an AVP
section. If a parallel cycle directive is encountered, this cycle
execution is subdivided into M branches, where M is the number of
processors of the corresponding AVP section.

**Subsystem for supporting parallel code execution.** The subsystem
implements the following functions: implementation of the execution
model on distributed processor systems; distribution of data and
computations; control of remote data access; dynamic adjustment with
regard to application parameters, configuration of the system,
heterogeneity of computations, and heterogeneity of the system of
processors; data acquisition for debugging and studying them and for
making the performance predictions; verification of whether
specifications of parallelism are correct and consistent; modeling a
parallel code execution under OS WINDOWS and UNIX using one processor.

DVM system is installed and operated at the following institutions: RAS
Supercomputing Center (MSC), RAS M.V.Keldysh IAM, MSU Research Computer
Center, MSU' Department of Computational Mathematics & Cybernetics,
NNSU, RAS UrB IMM.

The complete set of documents can be found on Web-site
<http://www.keldysh.ru/dvm/>. The system can be downloaded from the site
above and installed in a PC or a workstation.

**PARALLELIZATION OF ALGORITHMS**

**EOR NUMERICAL GROUNDWATER DYNAMICS MODELS**

# ***P.K.Konosavsky, K.A.Soloveychik***

## (RAS Institute of Geoecology's Sanct-Petersburg Branch)

Development and improvement of numerical simulation of groundwater
dynamics consists not only in developing new numerical methods and
techniques, but also in providing hardware to support computational
processes that, it its turn, requires a new approach to the development
of algorithms and computation technologies. One of the most promising
lines here is to use parallel computations that could allow several
times increase of the computation performance even without replacement
(modernization) of an available computer park. Note that in order to
solve modern complex problems within the traditional serial approach,
there is an often need in unacceptably much simulation time (dozens of
hours and even more of continuos operation of a computer to solve
problems of different-density flow and migration of multicomponent
solutions).

Currently, the main approach to increasing the performance of
computations consists in algorithm parallelization and adaptation to
architectures of modern massively parallel supercomputers. Any parallel
algorithm development begins with decomposition of a problem to be
solved and determination of mechanisms of communications between
branches executable concurrently. In so doing, a supercomputer
architecture and implementation tools determine the final performance of
an algorithm proposed.

The paper proposes a parallel algorithm which is a new version of the
idea of parallelizing based on decomposing a problem into physical
processes. This algorithm's principal distinctive feature is that each
processor participating in computations performs computations for its
own physical process. To implement the algorithm, the modern MPI
parallel programming library has been chosen allowing the algorithm to
be made essentially independent on the hardware architecture of the
supercomputer in use. A parallel algorithm of such a kind is portable to
almost all platforms available (different UNIX, Windows 95, 98, NT
clones) and any hardware base (from personal computers and workstations
to large computer clusters and supercomputers of various architectures).

A problem of convective-dispersive impurity transport in sub-surface
waters accompanied by kinetic sorption/desorption processes has been
chosen as an example.

The model represents by itself two interconnected physical processes,
the latter (inter-phase interaction) also consists at least of two
processes (depending on the number of active centers in rocks under
consideration). The numerical solution performance for such problems is
mostly determined by computer equipment speedup and the use of
technological multiprocessor computations seems to be the most optimal
alternative to other approaches (say nothing of parallelization "inside"
each physical process).

The developed algorithm verification has been performed using
supercomputers of various architectures:

-   a local memory supercomputer Parsytec Power Mouse

-   shared-memory supercomputers Convex Exampler SPP-1600

-   a computer cluster built basing on 3 workstations SUN Ultra 10 Sparc
    Station under OS Solaris 2.7 integrated into 100 Mbit network

-   a computer cluster built basing on 3 personal computers Pentium
    III-300 under OS Windows98 integrated into 10Mbit network.

> The algorithm achieves its maximum speedup (2.1 times) on Convex
> SPP-1600 supercomputer. The result is representative enough with
> taking into consideration that 3 times speedup (using 3 processors) is
> maximum possible, in ideal case. At the same time, the more complex is
> the computation, the higher is the speedup. This is caused by a
> reduced percentage of initialization operations and different-kind
> data transmissions in computation process as a whole that allows
> speaking about further increase of the algorithm speedup with a
> computation complexity increasing. The speedup gained on the cluster
> of Sun Sparc stations shows the same trend, however, with noticeably
> lower absolute values. Here, of principal importance is the rate of
> data communications within the network. Computations performed using
> the cluster of PCs showed a yet lower speedup, this is caused by the
> low network throughput (10 Mbits).

**3D SIMULATION OF DIRECT AND INVERSE RALEIGH-BERNARD**

**AND RALEIGH-TAYLOR PROBLEMS**

# ***A.I.Korotkiy, I.A.Tsepelev***

(Institute for Mechanics and Mathematics of RAS UrB, Ekaterinburg)

The paper considers issues of numerical simulation of 3D direct and
inverse thermal convection problems for a heterogeneous viscid
incompressible fluid (Raleigh-Bernard problem), as well as direct and
inverse problems of gravitational instability of a heterogeneous viscid
incompressible fluid (Raleigh-Taylor problem). Such hydrodynamic
problems are the problems of active and prolonged use in geology and
geophysics for simulating various processes in earth interior (see, for
example, Ref./1/). Instabilities of geological structure types under
consideration are the result of lower density rocks bedding under higher
density rocks, or the result of temperature-dependent softening of lower
layers as compared to upper layers, or the result of material density
changes because of phase transformations. By an inverse problem of the
geological structure evolution is meant reconstruction of the history of
its motion resultant from evolution of perturbations in that medium
under gravitational and thermal fields. In other words, this is a
problem of reconstructing the medium state it had in the past using data
on its current state. Mathematically, the problem is formulated, as
follows: it is required to find solution to the corresponding system of
equations (Navier-Stokes equations, incompressibility equations,
equations of density and viscosity transport, the heat balance equation)
with corresponding boundary and initial conditions in backward time.

Numerical simulation of a heterogeneous viscous incompressible fluid
attracts undiverted attention on the side of geologists and
geophysicists, in particular, experts on distortions of sedimentary
rocks, mining engineers, earth interior explorers and engineers on
underground structures. This is because of such structures' capabilities
to cause landslips, create traps for hydrocarbons, to be storage
facilities for radioactive waste. For example, studying salt domes is of
great importance for national economy. Essentially all of the
oil-and-gas resources in the near-Kaspian basin are related to salt
structures. In order to understand the history of sediment accumulation,
erosion and deformation in sedimentary basins, backward reconstruction
of the basin evolution process is required. Improvement of the methods
of processing seismic and other-kind geological exploration data allows
more detailed studies of these processes. Studying such data in
combination with results of numerical simulations leads to a better
understanding of the nature of these phenomena.

Mathematically, the direct problem is stable in a sense that small
errors of physical value measurements and numerical errors of
computations provide small errors of finding a final status within an
essentially unlimited time period. However, the inverse problem (in
backward time) is unstable in a sense that small measurement errors or
numerical errors of computations may lead to large errors of
reconstructing the original status within relatively small time
intervals. In this case special methods of solving incorrect problems
are required.

Such problem solution is based on introduction of a vector potential for
the medium motion velocity (this allows the requirement of
incompressibility and pressure in motion equations to be eliminated).
The finite element method with the special basis of tricubic splines is
used for numerically approximating motion equations. The heat balance
equation is approximated using the difference method. Transport
equations are solved using the method of characteristics. Finding
solution to the heat conduction equation in backward time is reduced to
solving some special variational problem that provides the model
immunity to errors. Numerical simulation of such problems is
labor-consuming because of large dimensions of the corresponding
discrete approximations. Some progress has been already made due to the
use of special basic elements of the finite element method and special
representation of the medium motion velocity's vector potential. This
allowed obtainment of acceptable qualitative and quantitative results
with relatively small dimensions of discrete problems (Refs./2-4/).

The algorithms developed to solve the problems mentioned above have been
implemented for parallel computers MCS-100, MCS-1000, IBM SP2.
Computations demonstrate the algorithms as applicable to numerical
reconstruction of diapirs. Computation results are presented for some
model problems.

The work is supported by RFFR grants 01-07-90210 and 02-01-00354, as
well as the grant under the target program intended to support
inter-discipline projects in RAS UrB - RAS SB cooperation.

## References

1.  Chandrasekhar S. *Hydrodynamic and Hydromagnetic Stability*. Oxford:
    *Claredon Press*, 1961.

2.  Ismail-Zade A.T., Korotkiy A.I., Naimark B.M., et al.
    *Implementation of a 3D Hydrodynamic Model of Sedimentary Basin
    Evolution* // *ZhVM i MPh,* 1998, No.7, P.1190-1203.

3.  Ismail-Zade A.T., Korotkiy A.I., Naimark B.M., et al. *Numerical
    Simulation of 3D Viscid Flows Under Gravitational and Thermal
    Effects* // *ZhVM i MPh*, 2001, V.41, No.9, P.1399-1415.

4.  Korotkiy A.I., Tsepelev I.A., Ismail-Zade A.T., Naimark B.M. *3D
    Simulation of the Inverse Raleigh-Taylor Instability Problem*. //
    *Izvestia UrGU, Ser.: Mathematics and Mechanics*, 2002, Iss.4,
    No.22, P.81-90.

**COMPUTER SIMULATION OF NEAR-BREAST AREA DYNAMICS OF WATER WELL UNDER
ELECTRICAL DISCHARGE IN WATER**

***V.M. Kosenkov***

(Institute for Pulse Processes and Technologies, Ukrain)

The measurement of fast process parameters in two-phase medium within
the near-breast area of the water well makes the computer simulation the
basic method to obtain the complete data for the dynamic processes under
the electric discharge in water.

The computer model was formulated based on the following assumptions:
the porous medium is isotropic and saturated with a single-phase viscous
compressible fluid; the array material of the porous medium is
strengthened according to Coulomb law; the start of the separation
fissures formation corresponds to the criterion of the first strength
theory.

The computer model equations describe the mass, momentum and energy
conservation laws in the matrix and saturating fluid in the
heterogeneous approximation. The phase interaction of the media is taken
into account by Zhukovsky hypothesis. The coupling between the
deformations and stresses is described by the generalized Guk law using
the plastic flow theory. The surface of the plastic flow with respect to
Coulomb strengthening is specified by the modified Guber-Mieses
condition. The deformation and displacement rates are related by
geometrical ratios. The relation between the average stress and volume
deformation is governed by the equations of state for solid and liquid
phases. The resulting system of equations is described in reference /1/.

For more detailed description of fissure formation and growth processes,
the system /1/ additionally includes the relations between the fissure
lengths and number determining the occurrence of separation fissures and
their geometrical characteristics if the critical stress duration
exceeds the medium relaxation time /2/. This allows determination not
only of the magistral growth directions and length of fissures but their
thickness as well which permits to relate the dynamic loading parameters
to the variations of the porous medium permeability.

The tensile strength yield of the material in the area of shift
(compression-driven) matrix destruction is assumed to be zero. The
absence of contact between the fissure walls allows the assumption of
the normal and tangential stress components to be zero.

The system of equations includes additionally the appropriate initial
and boundary conditions. The liquid piston conditions were applied to
the interface between the well fluid and porous medium. Prior to the
electrical discharge the environment was in the non-perturbed state due
to the hydrostatic and rock pressure. The pressure at the well surface
was obtained from the hydrodynamic problem /3/.

The system of equations was solved with numerical methods /1/. The
solution revealed the effect of the pressure wave parameters impacting
the well surface, the structure of the near-breast area and physical
properties of the medium on the destruction processes of the porous
medium array and time-dependent fluid flow.

1.  Kosenkov V.M. *Computer model for the dynamics of porous
    elastic-plastic media saturated with fluid under electrical
    discharge in water. Problemy prochnosty*. 1997, N1.p.87-96.

2.  Bonushkin E. K., Zavada N.I., Novikov S.A., Uchaev A.Ya., *Kinetics
    of dynamic metal destruction under dynamic pulse heating*. Sarov,
    RFNC-VNIIEF, 1998,275 p.

3.  Barbashova G.A., Kosenkov V.M. *Determination of hydrodynamic
    loading on the oil well wall resulting from the electrical
    discharge*. PMTF. 2001, vol. 42, N 6, p.93-97.

**ON USING ERROR BACK PROPAGATION NEURON NETWORKS**

**IN SIMULATIONS OF MULTI-MATERIAL SYSTEMS**

# ***V.M.Kyashkin, S.V.Rod'kin, O.B.Tomilin, S.A.Fedosin***

# (N.P.Ogarev Mordovian State University, Saransk)

**Introduction**

Efficiency of a natural-science or technological process depends on its
ultimate final result achieved with given complete information about the
process mechanism. Such information is a set of effects on the process
proceeding as well as a set of quantitative relations among them, which
can be used as a basis to construct a completely determinate
computational model. Its availability allows us to specify necessary
conditions for achievement of the ultimate result.

**Setting up the problem**

When simulating multi-material systems, the investigator typically
resorts to the functional study of the process. Then the following sets
are considered: input measured and controlled process parameters
$\left\{ x_{i} \right\} = x_{1},x_{2},x_{3},\ldots,x_{k\ }$,
uncontrolled, randomly varying process parameters (so-called "noise")
$\left\{ Z_{i} \right\} = Z_{1},Z_{2},Z_{3},\ldots,Z_{m}$, output
process parameters
$\left\{ Y_{i} \right\} = Y_{1},Y_{2},Y_{3},\ldots,Y_{\ell}$.

With this standard problem formulation, the computational model is a
response function relating the input and output parameters:

$Y_{i} = \phi_{i}(x_{1},x_{2},x_{3},\ldots,x_{k})$ (1)

A specific form of the response function for one or another process is
determined by regressive analysis methods. As a rule, response function
(1) is represented as a polynomial of the following form:

$$Y = \beta_{0} + \sum_{j = 1}^{k}{\beta_{j}x} + \sum_{u = 1}^{k}{\sum_{j = 1}^{k}{\beta_{\text{uj}}x_{u}x_{j}}} + \sum_{j = 1}^{k}{\beta_{\text{jj}}x_{j}^{2}} + \ldots,$$

where

$$\beta_{j} = \frac{\partial\phi}{\partial x_{j}};\ \beta_{\text{uj}} = \frac{\partial^{2}\phi}{\partial x_{u}\partial x_{j}};\ \beta_{\text{uj}} = \frac{\partial^{2}\phi}{\partial x_{j}^{2}}\text{.}$$

Then the optimization is made with the method of steep ascendance across
the response surface, the near-extreme region is described, the response
surface is examined, and the desirability function is generated /1,2/.

As seen from all of the above, simulations of multi-material systems
with conventional methods result in quite cumbersome calculations to
obtain the final result representing actual simulated items. Moreover,
the computer simulations with such methods entail another kind of
difficulties: numerical methods are needed for approximate solution of
the differential equations (for example, transformation of a partial
differential equation to a finite difference equation).

A computer model hard to implement leads to long problem computation
times, which eventually results in a situation, where the simulation
results can prove not called for by the investigator or not satisfying
the posed problem /3/. This is especially true for fast processes /4/
characteristic, for example, of glass manufacturing. (The authors
encountered the problem when developing and introducing barium-lithium
glass manufacturing and molding for luminescent tubes of Lisma Stock
Company, Saransk.)

**Solution method**

It is suggested that the error back propagation neuron network
technology should be applied to this problem class.

The algorithm used is an extension of one of the simple perceptron
teaching procedures known as Widrow-Hoff rule (or delta rule) and
requires a teaching sample representation. The sample is composed of a
set of image pairs, between which a correspondence must be established,
and can be considered as an extensive setting of the vector function,
whose domain of definition is a set of input images and the set of
values is a set of outputs.

Before the teaching the connections are assigned small random values to.
Each iteration of the procedure consists of two phases. During the first
phase the input vector (image) is transferred to the network by
appropriately setting the input element state. Then the input signals
propagate through the network and generate some output vector. For the
algorithm operation it is necessary that the neuron-like element I/O
characteristic be non-decreasing and have a limited derivative. For this
purpose, a sigmoidal function of form (2) is used.

$y = \frac{1}{\text{1-e}^{- \left( \text{S-}\theta \right)}}$ (2)

The resultant output vector is compared to the required one. If they
coincide, no teaching occurs. Else, the difference between the actual
and required output values is calculated and transferred sequentially
from the output layer to the input one. This information about error is
used as a basis to modify the connections with the generalized delta
rule.

**Obtained results and their analysis**

The neuron network simulated by us to find functions of type (1) in
Delphi environment by experimental data has been implemented as a
separate component.

The data of experimental barium-lithium glass manufacturing was used for
the analysis. The neuron network was used to find a function, which is
described by the equation of linear multiple regression in conventional
methods. The neuron networks taught on the experimental data describe
the relation between the manufacturing properties of glasses and their
compositions. The data was analyzed using a three-layered neuron network
with the quantity of input layer neurons equal to that of components in
burden, the same neuron quantity is in the internal layer, which
characterizes the above-discussed neuron network functionality and is
hidden from the investigator. The output data characterizing the
barium-lithium glass manufacturing properties are located in the output
layer of the neuron network. The number of the epochs (teaching
iterations) for each of the five developed neuron networks is taken
equal to 10000. A major teaching parameter determining the teaching rate
is taken as 0.1 for all the above-mentioned networks. The accuracy of
the developed computational models can be judged by a change in the
mean-root square error during the teaching. If at the very beginning of
the teaching in the neuron network teaching for the first manufacturing
property (epoch 1) this is 0.602450631, then at the end of the teaching
(after the 10000-th epoch) this will be 0.000274846. Its fastest drop is
graphically displayed in Fig. 1.

![](media/image1.wmf)

**Conclusions**

As seen, the back propagation procedure resolves the problem of
simulation time and yields fairly good results in the implementation of
the computational models of multi-material systems /5,6/.

The error back propagation neuron network has been successfully used
when developing the software for analysis of the data of the experiments
conducted to study barium-lithium glass properties (Lisma Stock Company,
Saransk).

## References

1.  Yu.P.Adler, V.E.Markova, Yu.V.Granovsky. Planning the experiment in
    search for optimal conditions. Moscow, Nauka Publishers, 1975,279 p.

2.  V.V.Nalimov, N.A.Chernov. Statistical methods for planning of
    extreme experiments. Moscow, Nauka Publishers,1965,268 p.

3.  *A.A.Ligum, A.D.Malysheva.* Computer processing of experimental
    data*. Dneprodzerzhinsk,* DII Publishers*, 1992,47p.*

4.  *A.G.Ivakhnenko.* Long-term prediction and control of complex
    systems*. Kiev,* Tekhnika Publishers*, 1975,312p.*

5.  *F.Rozenblatt.* Principles of neurodynamics (perceptron and theory
    of brain mechanisms*). Moscow,* Mir Publishers*, 1965,480p.*

6.  R.J. McEliece, E.G Rosner, E.R. Rodemich, S.S. Venkatesh The
    capacity of Hopfield associative memory. //IEEE Transactions on
    Information Theory IT-33:461-82, 1987.

**PARALLEL COMPUTATION TECHNOLOGIES**

**FOR THE PROBLEM OF RADIATION TRANSPORT**

# ***E.F.Lelikova, L.I.Rubina, O.N.Ul'yanov, M.A.Chashchin***

## (Institute for Mechanics and Mathematics of RAS UrB, Ekaterinburg)

Radiation transport problems are of permanent interest for researchers,
because they occur in various areas, for example, when describing
processes in atmospheres of stars, in multiple-charge ion spectroscopy,
and laser designing.

The paper authors consider the radiation transport problems for a flat
immobile layer (or a layer that moves in some special manner) of a
finite thickness. For such problems, the method of numerically
simulating radiation transport processes has been developing for several
years. This method's specifics consists in that one and the same group
of investigators develops codes implementing two essentially different
approaches: MAPI that uses explicit formulas to solve transport
equations and MPLCh that uses Lagrangian method of interpolation
polynomials for the same purpose. For the period of investigations, the
method undertook significant changes to broaden the range of problems
(to make a physical model more complex, including consideration of a
large number of resonance lines, increase in the number of materials in
mixture, as well as the number of sub-layers in the region under
consideration in which the mixture differs in its composition, velocity
and some other parameters). Numerical simulation of radiation transport
processes is a difficult enough problem, because it is conjugated with
significant mathematical and, therefore, computational problems and
difficulties that require necessary large amounts of calculation
efforts. The method evolved depending on capabilities of computer
equipment and computation technologies (and was restricted by them).
Serial codes doesn't allow computation of informal problems, even with
the use of algorithms based on combination of numerical and analytical
methods and specially developed high-performance algorithms. Significant
advances have been achieved only with the use of parallel computation
technologies, when it has become possible to solve, for acceptable
computation times, not only model problems, but also those informal
enough that are of interest for engineers and physicists.

The experience of parallel algorithm and code development, as well as
the use of various multiprocessor computer systems (transputer systems,
computer clusters, multiprocessor parallel and massively parallel
systems, including MCS (from early models to MCS-1000M and MCS-1000/16))
allow some practical conclusions to be made and problems to be
formulated. The following can be mentioned among them, in particular:
selection of those preferable algorithms among alternative options that
can be parallelized more easily; acceptable computation times; selection
of algorithms that are independent on \"topology\" of the computer
system in use; the optimum number of processors used; minimizing the
amount of data to be transmitted during remote access; availability of
usable debugging tools. When developing a method for a long time by a
small team of investigators which is intended for solving a serious
problem, one of the aims is to provide a possibility of easy transition
from one computer system to another and concurrent use of various
computer systems. Another aim is to provide a possibility of running a
code using an arbitrary number of processors of one and the same
computer system. This is particularly the problem for programmers
operating communication environments MPI and Router. For this purpose,
special algorithms have been developed and transition from parallel
codes with fixed binding to the chosen number of processor in use, to
algorithms and codes intended for parallel computations on a relatively
random number of processors has been done. Later, a successful attempt
of transition to Fortran-DVM (DVM parallel programming system) has been
made. Numerical simulations are currently performed on parallel computer
systems MCS-1000/16 and MCS-1000/M.

The work has been performed under RFFR grants 00-15-9604 and
01-07-90210, as well as under the target program of supporting
inter-discipline projects being executed in cooperation of RAS SB and
RAS UrB scientists.

**A REMOTE DATA VISUALIZATION SYSTEM**

# ***N.V.Linich***

## (M.V.monosov MSU, Moscow)

Modern technologies for parallel code execution using multiprocessor
clusters can provide high computational capabilities to perform
large-scale computational experiments. Such clusters are not prevalent
and are now accessible only for a limited number of large institutions
because of their high cost. however, high priority users may have the
right of remote access to such multicomputer resources to use them.
Computational experiments, as a rule, generate numerical data required
for studying them visually. it makes no sense to transmit large amounts
of data through the network to a researcher's local computer. There is a
need in a remote visualization system allowing generated data images to
be constructed directly on the cluster by means of some graphic system
and transmitted to the user's computer. The paper presents such a
system, where data transmission is performed using web-technologies:
java-servlets, java-applets; xdraw system is used as a graphics code for
direct quick viewing of data during studies.

The remote visualization system presented is an interactive tool for
graphical data representation allowing data representation control using
ordinary web-tools.

**SOLVING A MANY-ELECTRON PROBLEM USING THE METHODS**

**OF MINIMIZING FUNCTIONS OF SEVERAL VARIABLES**

***Yu.B. Malykhanov***

(M.E. Evsev'yeva MSPI, Saransk)

The Hartree-Fock-Rutan method (HFR) is a most widely used quantum
mechanics method of calculating electron shells of atoms and electronic
structures of molecules. The method implementation in practice causes
severe difficulties in calculations and is reduced to solving two
interrelated problems: the problem of solving HFR equations directly and
the problem of finding an optimum basic set of atomic orbitals (AO)
allowing the achievement of Hartree-Fock limit of the system's
properties to be calculated. Optimization of nonlinear AO parameters
(orbital exponents) providing minimum energy is an especially difficult
problem. Traditionally, zero-order minimization methods (direct search)
of extremely slow convergence are used for this purpose, their
implementation causes tremendous calculations and is possible only with
the use of supercomputers in case of heavy atoms.

The author proposes a unified approach to solving both problems for
many-electron systems with both filled and open shells which is based on
advanced methods of the first and the second order for minimization of
functions of several variables possessing high convergence rates and
allowing the solution to be found with any a priory specified precision.
To implement minimization algorithms, formulas for calculation of the
first and the second system energy derivatives using the optimized
parameters which are elements of the density matrix and orbital
exponents of AO have been obtained. The proposed algorithms
serviceability has been verified by calculations of atoms of I&#8211;IV
periods with both filled and open shells; analytical Hartree-Fock
orbitals giving the Hartree-Fock limit for energy have been calculated
using Pentium II 400 MHz computer with a limited number of iterations
and with high precision (in terms of virial relation value).

**HARDWARE ARCHITECTURE IN THE INTERPROCESSOR COMMUNICATION SYSTEM OF
MP-X MULTIPROCESSOR SYSTEM**

# ***V.S.Popov, S.A.Stepanenko, A.A.Kholostov***

## (RFNC-VNIIEF, Sarov)

For the past ten years, RFNC-VNIIEF has been developing the line of
multiprocessor computer systems MP-X with \"hypercube\" architecture.
The practice of their development and operation allows the conclusion
that the most optimum option of such system structure is to apply widely
used processor platforms (further referred to as computational modules
(CM)) having the specialized interprocessor communication system (SIPC).
SIPC parameters (channel bandwidth and message delay time, first of all)
determine, to a significant extent, the performance of such complexes.
The paper describes major design principles and the structure of SIPC
system's hardware, shows main message transfer processes, gives SIPC
system parameters as compared to those of similar products.

Important features of MP-X SIPC are, as follows:

-   A capability to use it with computational modules built basing on
    various microprocessor types under various OS;

-   A capability to construct multiprocessor systems of an arbitrary
    topology; optimization for MP-X architecture (hypercube);

-   Introduction of a processor for communications that allows optimum
    control of interprocessor communications between CM and SIPC;

-   Hardware support for group data transmission operations;

-   Distinction between routes for transit data transfers and \"own\"
    data receipt/transmission, use of PCI 64/66 multiple-channel
    interface, availability of a large-capacity buffer memory in
    PCI-adapters allow the improvement of characteristics of the
    bottleneck, CM/SIPC interface which is a most vulnerable one, on our
    opinion, in such systems.

A SIPC channel throughput is 240 MB/s, a switch delay doesn't exceed 250
ns, throughput of a CM/SIPC communication channel is 2x512 MB/s. These
parameters correspond to the level of most recent products by a world
leader in manufacturing interprocessor communication means, Myricom,
Inc.(Ref./1/).

### References

1.  *Certification of a usable model* No.4266, Application
    No.95101975/20 of 02/09/95.

2.  Stepanenko S.A. *Estimation of temporal complexities of global
    communication algorithms for mpp systems* // VANT, Ser.:
    Math.model.phys.process., Iss.2, 1998, P.18-23.

3.  *Guide to myrinet-2000 switches and switch networks* / Myricom Inc.,
    2001 -- www.myri.com

**STUDYING TURBULENCE CHARACTERISTICS BY DIRECT**

**3D NUMERICAL SIMULATION OF GRAVITATIONAL MIXING**

## ***S.V.Rebrov, O.G.Sin'kova, A.L.Stadnik, V.P.Statsenko, Y.V.Yanilkin***

# (RFNC-VNIIEF, Sarov)

Turbulence evolution in gravitational field at a plane interface of two
incompressible fluids (gases) with a large density ratio 3 ≤ ρ~2~∕ρ~1~ ≤
40 is studied numerically using 3D hydrodynamic TREK code.

Computations were performed using a detailed computational grid
200x200x400 by parallelization to several dozens of processors.

Computation results were processed (averaged) to obtain moments of
hydrodynamic values: diagonal components of Reynolds tensor (turbulent
energy), turbulent flows, density profiles and its mean-square
pulsation. The results obtained are compared to predictions by
phenomenological turbulence models and known data of experiments.

**THE TECHNIQUES FOR NUMERICAL SOLUTION OF A 2D NONLINEAR**

**HEAT CONDUCTION EQUATION IN COMPLEX GEOMETRY PROBLEMS**

# ***N.M.Rud'ko, S.I.Skrypnik***

## (RFNC-VNIIEF, Sarov)

New techniques and a code for numerical solution of a 2D nonlinear heat
conduction equation have been developed in order to increase the
accuracy and producibility of 2D computations within RAMZES code complex
and make them more resistant to crashes.

A one-region version of the technique has been already described in
articles published in VANT journals in 1999 and 2000, and in papers
presented at the previous International Workshop on super-computations,
13-15 September 1999.

The technique development for solving complex problems which geometry
consists of several fragments, each of them being calculated using its
own more suitable computational grid, is described in the paper. The
principal equations used are given and the code that implements these
techniques is described.

The techniques presented allow a significant increase in the numerical
solution accuracy for 2D heat conduction problems, that is demonstrated
by numerical results for test problems having exact analytical
solutions.

**MODELING OF ISOTROPIC DESTRUCTION**

**OF ELASTIC-PLASTIC MEDIA WITH MIMOZA CODE**

***A.A.Sadovoi, N.V.Sokolova***

(RFNC-VNIIEF, Sarov)

The MIMOZA codes intended for the calculations of multidimensional solid
mechanics applications in hydrodynamic and elastic-plastic
approximations on Lagrangian and Eulerian-Lagrangian grids include the
isotropic kinetic model for the plastic destruction of the medium. It
relies upon the well known NAG model that experienced some modifications
and improvements basically associated with the modeling of deformation
processes in the destroyed material. The model treats the destruction as
a multiphase process including the formation, growth and merging of
microdefects and formation of fissures.

The pores are assumed to have the spherical form under plastic
destruction. The pores occur under the tensile stresses. The originally
integral material is unloaded following the adiabatic curve. However
since the negative threshold pressure P~g0~ the material exhibits the
tension and new pores form when the threshold pressure P~h0~ is reached.

When the critical damage, ω~cr~ is reached the material goes to the next
phase - merging of microdefects and formation of fissures. The material
no longer resists the expansion, that is the destruction occurs. The
material is assumed to dissociate into non-interacting fragments in this
case and the cell pressure is said to be zero. Further expansion of the
destroyed material does not change the pressure which remains zero.

If the damaged medium begins to compact the material exhibits
irreversible volume deformations that is compaction occurs.

When the compression proceeds, the full compaction may occur and then
the material compacts following the adiabatic curve of the solid
material.

This model was used to carry out the calculations of some experiments.
Generally we can outline a good agreement with the experimental data.

**A PROTOTYPE APPLICATION DEVELOPMENT SYSTEM**

**FOR PROCESSING EXPERIMENTAL DATA**

**USING HETEROGENEOUS MULTIPROCESSOR SYSTEMS**

# ***A. Sazonov, A.N. Salnikov***

(M.V.Lomonosov MSU, Moscow)

As practice shows, the use of a single-processor system of even the very
high performance to solve some problems may become difficult enough.
Such problems may require either large computational resources, or a
large storage capacity (or both of them at the same time). Very often,
such problems can be reduced to a set of smaller sub-problems whose data
dependence can be represented in the form of a graph. The problem of
processing images for which the processing process can be reduced to
application of a set of certain operations, such as boundary outlining,
changes of image brightness and contrast ratio, superposition, etc., is
an example.

After sub-problems of a given problem have been determined, it is
possible to construct a graph (being a data dependence graph) which is
equivalent to the process of data processing on the level of
sub-problems. A parallel code to be executed using a multiprocessor
system can be constructed using the dependence graph obtained. The
dependence graph nodes become functions that can be called by certain
processors comprising the multiprocessor system in use and its arcs
become data transmissions via communication channels, if necessary. The
vertexes and arcs of the dependence graph are marked with numerical
indexes characterizing the complexity of a sub-problem at a given vertex
and the amount of data transmitted via the graph edges, if it is the
case, when the graph's vertexes connected by the arc are called by
various processors.

The dependence graph nodes are dynamically selected for their execution
by a processor. When executing the code, decisions concerning a graph
node call by a processor are made basing on information about a current
process of the parallel code execution and data on the performance of
the multiprocessor system in use. The system of tools being developed to
support parallel programming will perform some part of OS functions. A
special set of tools responsible for the process of executing a code is
built into the parallel code. This set of tools will take care of proper
consideration of a problem information structure, will solve the problem
of composing a dynamic schedule for a parallel code execution by a
multiprocessor system, as well as the problem of processor load
balancing and communication channel loading.

The static schedule of executing the dependence graph nodes constructed
prior to the parallel code execution by a multiprocessor system should
be considered, if the built-in set of tools is unable to select a node
to be called by a processor basing on dynamics considerations. The
schedule is also considered, if a parallel code is run in static mode of
operation. The schedule is created using the genetic algorithm using
data acquired by testing the performances of processors, the rates of
data transmission via network, as well as a priori data on the
complexities of node operators in the algorithm's graph and on the
amount of data transmitted through an edge of the graph.

Utilities under development will provide flexible transition to a
multiprocessor system of a different architecture and help to lower, in
part, the reduction in the performance. Flexibility consists in that for
a new multiprocessor system there will be required to acquire data on
the network performance and the performances of processors, obtain a new
static schedule of executing a parallel code, and re-compile the code.

The work is being performed at Intel students' laboratory of
M.V.Lomonosov MSU.

**TECHNIQUES FOR SOLVING TIME-DEPENDENT ELASTIC-PLASTIC PROBLEMS USING
IRREGULAR POLYHEDRAL LAGRANGIAN GRIDS**

# ***S.S. Sokolov***

## (RFNC-VNIIEF, Sarov)

Numerical methods similar to Wilkins method are widely used to solve
multidimensional time-dependent elastic-plastic problems. Among most
popular multiple-purpose methods one can mention the method of finite
differences and the method of finite elements. The range of problems to
be solved is rather wide and it's impossible to cover it using a single
method, in a multidimensional case, in particular. For elastic-plastic
problems, among the most developed are finite-difference schemes mainly
using regular quadrangular (in 2D case) and hexahedral (in 3Dcase)
Lagrangian, Eulerian, and Eulerian-Lagrangian grids.

The paper presents a Lagrangian code for solving multidimensional
time-dependent elastic-plastic problems using irregular Lagrangian
grids. In a 2D case, the code, which is based on the Lagrangian
hydrodynamic DMK code, uses both an irregular difference grid consisting
of convex polygons with a random number of vertexes and a regular one
consisting of quadrangles. Each of these codes can be used for different
problem solution domains that remain convex during the computation
process. In a 3D case, the code for calculating elastic-plastic flows,
which is based on the Lagrangian hydrodynamic TMK code, uses an
irregular convex polyhedral Lagrangian grid of a random configuration.
Explicit finite-difference schemes are used to approximate equations in
3D and 2D codes.

The techniques described are implemented in TMK and DMK complexes that
allow continuum mechanics problems with severe distortions in complex
geometry regions to be solved.

**STUDYING THE EFFECT OF SOME ARTIFICIAL VISCOSITY FORMS**

**ON THE \"ENTROPY TRACE\" NATURE IN NUMERICAL SOLUTION**

**OF GAS DYNAMICS EQUATIONS**

***E.A.Solovyeva, A.M.Stenin***

(SarPTI, RFNC-VNIIEF, Sarov)

It is known that a serious demerit of the method for the numerical
computation of hydrodynamic jumps with artificial viscosity suggested by
Neumann and Richtmyer in 1950, which has recently become classical, is
appearance of a so-called \"entropy trace\" in the numerical solution.
The entropy trace appears as finite errors in densities and internal
energies in the vicinity of external boundaries of the computational
domains and in that of internal material interfaces on passage of shock
waves across the boundaries and interfaces as well as in shock wave
collision even in homogeneous material.

In 1987 W. F. Noh suggested a method of artificial heat conduction in
one of his papers to reduce the entropy trace in the numerical
computations. Computations show that the method facilitates a
significant reduction in the entropy trace with an appropriate
artificial heat conductivity factor in a number of problems.

This paper introduces artificial viscosities to all the three gas
dynamics equations. That is, the following is conservatively added to
the gas-dynamics equation system written in the form of the laws of
conservation in Eulerian coordinates: artificial mass diffusion to the
continuity equation, artificial momentum diffusion to the equation of
motion, and artificial total energy diffusion to the equation of energy.
The transition to the Lagrangian coordinates incidentally results in
relevant formulas for viscosities in each of the equations. Computations
for a number of problems suggest that in all instances the resultant
viscosity system is favorable to a significant reduction in the entropy
traces in the numerical solution.

Later, the viscosities were modified using shock Hugoniot conditions
written for the general case. Some typical computations for 1D ideal gas
flows possessing plane symmetry indicate that in most cases the maximum
deviation of density from its exact value in the entropy trace region is
no more than 3-4% when the suggested system of viscosities is used. In
some computations no entropy trace is observed whatsoever. The largest
deviation of density from its exact value occurs on the external region
boundary, when pressure is given for the boundary condition. But even in
this case the deviation is no more than 7-8%.

For 1D flows possessing cylindrical or spherical symmetry, in the
equation of motion it is appropriate to use tensor artificial viscosity
proportional to the deformation rate tensor deviator as well as its
modification obtained using the shock Hugoniot conditions.

The computations of the well-known Noh's problem for the cases of
cylindrical and spherical symmetry conducted with one of the
modifications to the new system of viscosities suggest that the maximum
deviation of density from its exact value is about 3%. In similar
computations with classic viscosity conducted with the same
dimensionless factors the deviation of density from its exact value
approaches 90% in the entropy trace zone. Note that all the computations
were conducted with the same dimensionless viscosity factors.

**TOPOLOGICAL REDUNDANCY**

**IN MULTIPROCESSOR ENVIRONMENTS**

S.A.Stepanenko

# (RFNC-VNIIEF, Sarov)

##  The paper proposes the redundancy method whose specific features are: first, full identity of backup components and those being reserved and, second, full preservation of the environment topology which means that when replacing a defective component by a backup one only logical numbering of components and links between them changes. No physical changes and degradation, in particular, take place. 

The method proposed is described for hypercube Γ^n^, where n=2^k-1^,
k=2,3,...,n is the hypercube dimension that coincides with Hamming code
dimension. This code properties are used in the method presented.

The method allows anyone failed component to be replaced. However, if
components failed are such that their vectors have been obtained by
changing one and the same component σ~i~ of the corresponding separated
vectors, then 2,3,...,2^k^ components can be reserved \"simultaneously\"
(by one and the same operation of logical numbering).

Probabilities of simultaneously reserving several processor elements
have been estimated. For m vectors we have $P_{m} = P_{m - 1}^{m}$.

It is known that 1D, 2D, 3D tori and other-type topologies can be
constructed in hypercube using Gray codes.

The method is usable to provide redundancy of such topologies with no
changes required.

The method can be used to provide redundancy of any structures
(automata, collectives, organisms) having the topologies above.

# References

1.  Von Neumann J. *Probabilistic logic*, California, Inst. of
    Technology, 1952,

> (in Russian: Automata, IL Publishers, 1956, Von Neumann J.
> Probabilistic logic and Synthesis of Reliable Entities Based on
> Non-Reliable Components).

2.  Moor E., Shannon C., *Reliable circuits using less reliable relays*,
    Journal of the Franklin Institute, N3 (1956), 191; N4, 281 (in
    Russian: C.Shannon. Papers on Information Theory and Cybernetics, IL
    Publishers, 1963).

3.  *Multiprocessor Computer Systems* /Edited by Khetagurov. - М.,
    Energiya Publishers, 1971.

4.  A.Avizhenis, J.-K.Lapris, *Dependable Computing*, TIIER, 1986, V.
    74, No. 5, P. 8-21.

5.  ASCI Red, http

6.  ASCI Blue Mountain, http

7.  V.Feller, *Introduction to the Probability Theory and Its
    Applications*.- М., Mir Publishers, 1964.

8.  U.Peterson, *Error-Correcting Codes*.- М., Mir Publishers, 1964.

9.  A.Kofmann, *Introduction to Applied Combinatorics*. -- M., Nauka
    Publishers, 1975.

**THE EFFICIENCY OF MULTIPROCESSOR ENVIRONMENTS**

**WITH GLOBAL COMMUNICATIONS**

# ***S.A.Stepanenko***

## (RFNC-VNIIEF, Sarov)

Multiprocessor efficiencies depend on time consumption to execute data
communications between processor elements.

MPI standard (Ref./1/) determines a set of global communication
operations consisting of most "popular" ones that are often encountered
during computations.

To implement this set of operations, various communication algorithms
are used. Such algorithms are developed with regard to

-   the topology of links between PEs of multiprocessors (ring, lattice,
    torus, hypercube, tree, ω-net and some others) (Refs./ 2,3/);

-   the data transmission mechanism implemented (channel switching,
    message buffering).

The paper presents analysis of the performance of various multiprocessor
environments with global communications. The estimates for algorithms'
temporal complexities have been obtained that allow the performance of
various multiprocessors to be determined for the set above. These
algorithms are shown to be close to perfect ones. The paper is
illustrated.

### References

1.  *MPI*: *A Message Passing Interface Standard. Message Passing
    Interface Forum, version1.1:* June 1995, ARPA&NSF grant ASC-9310330.

2.  G.T.Artamonov, V.D.Tyurin. *Computer Network and Multiprocessor
    System Topology*. -- M.: \"Radio i Svyaz\" Publishers, 1991, 248pp.

3.  S.A.Stepanenko. Temporal Complexity Estimates for Global
    Communications

> Algorithms on MPP systems. VANT, Ser.: Math.Model.Phys.Process., 1999,
> Iss.2, P.18-23.

**PARALLEL COMPUTATIONS WITH TOEPLITZ MATRICES**

***E.G. Sukhov***

(RAS Institute of Control Problems, Moscow)

Toeplitz matrices (T-matrices) are often encountered, when solving
problems of control, problems of item detection and their inner
reconstruction, problems of finding numerical solutions to integral
equations. Besides, such problems also have multiple-layer partitioned
matrices. Methods of solving Toeplitz systems of linear algebraic
equations of complexity O(n^2^) which are efficient ones for solving
small-size problems on serial computers, loose their efficiency with
parallelization \"straightforwardly\". The first cause of such a
phenomenon is that algorithms include the function of processing
variable-length vectors that complicates processor load balancing. The
second cause is related to the lack of acceptable matrix decompositions
in classic algorithms for T-matrices. This doesn't allow separation of
calculations with elements of matrices from calculations with their
right-hand sides during the process of solving systems of linear
algebraic equations.

Thus, it is necessary to develop parallel numerical methods of solving
Toeplitz systems that will be of high performances both in sense of the
number of floating point operations and in sense of parallelism. The
paper proposes recursive methods with partitioning to solve systems of
linear algebraic equations with T-matrices. It is shown that algorithms
basing on the composition of classic and \"fast\" methods are most
efficient ones.

**ON COMPUTATIONS OF TWO-COMPONENT FLOWS**

**IN COMPLEX-GEOMETRY CHANNELS**

***L.N.Sukhova***

(M.V.Keldysh IAM, Moscow)

Computation of two-component mixture flows is an important and urgent
problem. Simulation of processes in disperse gas flows through channels
of varying geometry can be used both for determination of optimal
conditions to purify disperse flows from a disperse phase (for example,
in aerosols) and at production facilities using such technologies for
transportation purposes to prevent sedimentation on walls of channels
(pipes). When solving problems of such a kind numerically, significant
difficulties arise, both with consideration of interactions between
phases and because of a computational domain's complex geometry. A
computational model of processes in channels of varying geometry for
disperse gas flows has been developed, as well as a numerical algorithm
for parallel computations.

**A MULTI-GRID METHOD IMPROVEMENT FOR HYDRODYNAMIC SIMULATIONS OF
OILFIELDS**

**WITH BOTTOMHOLE PRESSURES SPECIFIED**

# ***I.T.Usmanov***

## (Kazan State University)

At present time, multi-grid methods (Ref./1/) are among the most
\"fast\" finite-difference numerical schemes for solving elliptic
equations using large grids. However, computations using such methods
for regions having dotty sources and discharge outlets show that the
commonly used ways of their specification don't lead to satisfactory
results.

Ref./2/ describes the method of corrections using which one manages to
gain good agreement between the source power and the node pressure, the
latter can be interpreted as a \"bottomhole\" pressure (pressure on a
small-radius circle contour with a node being its center). For
subsurface hydraulic mechanics problems, of especial importance are the
two methods of determining the field of pool pressures: with specified
flow rates of wells (powers of sources) and in the mode of specified
bottomhole pressures. For the first method, a multi-grid scheme is
simply written with regard to corrections from Ref./2/. In the second
case, supplementary improvement of the linear equation system's factors
is needed after which a multi-grid technology becomes applicable to
solving the problem.

The algorithm developed has all of the required positive features of
classic multi-grid methods (fast convergence, low dependence on
gradients of factor fields and on the number of grid nodes). The
precision of solutions resulted from the algorithm operation has been
verified using test problems. Several computations have been performed
for real oilfields.

1.  K.Stuben, U.Trottenberg. *Multi-grid methods: fundamental
    algorithms, model problem analysis and applications*//*Multi-grid
    methods, Springer lecture notes in mathematica*. -- New York,
    *Springer Verlag,* No.960, 1982.

2.  A.N.Chekalin. *Numerical solutions of flow problems for water-oil
    pools*. -- Kazan, *KSU Publishers*, 1982.

**THE CHANNEL-TEMPORARY DIAGRAM**

**OF FAST CALCULATIONS OF INDIVIDUAL SUMS**

***Filippov N.A.***

(VyatSU, Kirov)

A system of **m** channels is discussed, in each of which there are
**k** operators, whose inputs one or two of **n** operands are moved to.
It is represented by the parallel-sequential time diagram of **k**
steps. At the first step, the problem, i.e. its function and operands,
is accepted and memorized. At the second step, the contents of the first
is analyzed and the structure of the fasted solution to the resultant
problem is constructed: all (k-2)\*m solvers and all connections between
them are determined (2 in (k-2) -- the first two blocks: the buffer of
memory and central logic device for analysis of the sum and synthesis of
the system of the solution to the problem with its structure of spatial
-- temporary interrelations of the operators).

Each operator solves a single or double elementary function: either a
specific trigonometric (single) function is calculated by one operand,
for example, angle, or one of arithmetic operations (double) is
performed on two arguments. After each operator a local logic key can be
put, which is also connected to the general logic input unit and
accounts for all preliminary results of operations on all channels and
has a managing connection with all local ones. Each of the latter
directs the result obtained before it to the input of an appropriate
operator located in an appropriate channel.

With each step the number of the channels becomes smaller. At the last
but one step, the final operation is performed (in a particular sum it
can occur earlier): **k** steps, like **m** channels, are taken the most
complex sum of a given type, which the structure under discussion is
just being constructed for.

A man, a one-processor digital computer, either arithmetic or logic
device, a functional converter based on permanent memory with all
expected results of intermediate computations can be used(as the
operator of each channel and step.

As there are many types of sums and, hence, the values of **n**, **m**,
**k** can significantly differ, for universality it is necessary to take
them large. On the other hand, this leads to unnecessary complexity of
the univezsal solution system for any sums from the standpoint of
specific sums requiring smaller **n**, **m**, **k** quantities. This
makes one to solve the optimization problem of the set selection. To
facilitate this, it would be reasonable if the application
mathematicians attribute similar problem types to separate categories
with their splitting into elementary, described ones, making problem
types from them with specification of their maximum sets **n**, **m**,
**k** for the solution of any of them.

The time analysis of the channel-time diagram allows us not only to
determine time of the solution to a specific sum, but also to secure the
mechanism for execution of intermediate elementary operations for
selection of shortest-time variants from them through possible updating.

**CONSTRUCTION OF OPTIMAL MESHES IN MULTIPLY**

**CONNECTED COMPLEX-TOPOLOGY DOMAINS**

**ON MULTIPROCESSOR COMPUTERS**

***A.F. Khairullin, O.B. Khairullina***

(RAS Urals Branch IMM, Ekaterinburg)

The suggested parallel algorithm for construction of 2D optimal
block-structured meshes with smooth grid lines at the block interfaces
in multiply connected complex-topology domains leans on block
overlapping method /1,2/. Criteria are specified for concurrent
computation of grids in the blocks (including with account for nonunique
mapping) on different processors, several modes are developed for
automatic block distribution over processors for uniform loading of the
processors used and reduction in data to be transferred. Independently
on the domain partition into blocks and the mapping topology, the grid
point coordinate array is closely packed. This allowed a significant
reduction in required memory to store it and making the overlap
arrangement most universal.

In complex-configuration domains the optimal meshes are constructed with
iterative method /1/ that implies minimization of the functional
responsible for the formalization of the criterion of closeness of the
grids to the ones uniform in distances between nodes and that of
closeness of the grids to the orthogonal ones at intersection points of
coordinate lines. To compute a grid of a needed quality, it is necessary
to give some initial approximation to it.

Refs. /3, 4/ develop a method for automatic partition of a
complex-geometry domain into blocks and a geometric method for
construction of an initial approximation to a curvilinear structured
grid in a block which is an single-connected domain of arbitrary type.
The methods and the above-described algorithm allowed us to resolve the
problem of automatic construction of a multi-block curvilinear grid
large in size (on the order of several hundreds of millions of nodes)
for an acceptable time in 2D domains of any connectivity and any
configuration using minimum input data, when the boundary in a physical
plane with nodes collocated in it and its image in a curvilinear
coordinate plane are given.

The algorithm extensively employs multiprocessor computer distributed
memory and is implemented in program complex MOPS-2A. The optimal grids
worked well in solving gas dynamics problems.

The work was carried out under financial support of Russian Fundamental
Research Foundation, Grants N 02-01-00236, N 00-15-96042.

References

1.  Khairullina O.B. *Method of constructing block regular optimal grids
    in two-dimensional multiply connected domains of complicated
    geometries* // *Russ. J. Numer. Anal. and Math. Modelling*. 1996. V.
    11, N 4. P. 343-358.

2.  Artemova N.A., Khairullin A.F., Khairullina O.B. *A parallel
    algorithm for computing optimal grids*// *Vychislitelnyje
    Tekhnologii*. 2001. V. 6, N 2. P. 3-13.

3.  Khairullin A.F., Khairullina O.B. *Automatic construction of initial
    approximation to curvilinear grid*// *Algoritmy i Programmnyje
    Sredstva Parallelnykh Vychislenij*. No. 4. Ekaterinburg, RAS Urals
    Branch IMM, 2000. P. 273-286.

4.  Khairullin A.F., Khairullina O.B. *Domain partition into blocks in
    construction of a multi-block curvilinear grid* // *Algoritmy i
    Programmnyje Sredstva Parallelnykh Vychislenij*. No. 5.
    Eekaterinburg, RAS Urals Branch IMM, 2000. P. 261-274.

**ON LOCAL REPRODUCIBILITY OF NONLINEAR**

**DIFFERENTIAL EQUATION SYSTEM**

***P.A. Shamanaev***

(N.P. Ogarev Mordovian State University, Saransk)

The paper considers a mechanical system, whose equation of motion is
described by nonlinear differential equation system

$\overset{˙}{x} = \text{Ax} + f(t,x)$ (1)

where $x \in R^{n}$,
$f \in C^{(p,q)}(\lbrack T, + \infty) \times V,R^{n})$,
$p \geq 0,q \geq 1$, $f(t,0) \equiv 0$, $V \subseteq R^{n}$ is a region
containing a vicinity of zero; $A$ is a constant ($n \times n$) matrix.
Also, let the following condition be held for vector-function $f$:

$\left\| f(t,u') - f(t,u\text{ \textbackslash)  rdline  <= ω \textbackslash( t,r \textbackslash)  ldline u' - u} \right\|$
$\forall u\text{',}u\text{ func ∈ S rSub \{ size 8\{r\} \} ,t func ∈  \textbackslash[ T,+ infinity  \textbackslash) \} \{}$

where
$\omega \in C(\lbrack T, + \infty) \times R_{+}^{1},R_{+}^{1}),\text{  }$$\omega(t,\alpha_{1}) \leq \omega(t,\alpha_{2})$
for $\alpha_{1} \leq \alpha_{2}$, $\omega(t,0) \equiv 0$; $S_{r}$ is a
sphere of radius $r$ with its center at the origin of coordinates.

Denote the set of all differential equation systems of type (1) by
$\Xi$. The solutions to the differential equation systems from the set
$\Xi$ are assumed to exist and be uniquely defined for any $t \geq T$
within some vicinity of zero.

The notion of reducibility for nonlinear differential equation systems
was introduced by E.V. Voskresensky /1/. Consideration of the
differential equation system solution behavior in the vicinity of zero
entails the problem of local reducibility of nonlinear differential
equation systems. By the local reducibility we mean the same as in ref.
/2/.

The problem of local reducibility of differential equation systems from
set $\Xi$ to linear system

$\overset{˙}{y} = \text{Ay}$

is posed. (2)

Denote different real parts of eigenvalues
$\lambda_{1},\lambda_{2},\ldots,\lambda_{n}$ of matrix $A$ by
$\Lambda_{1},\Lambda_{2},\ldots,\Lambda_{k}$, maximum orders of Jordan
cells of matrix $A$ that correspond to eigenvalues with real parts
$\Lambda_{1},\Lambda_{2},\ldots,\Lambda_{k}$ by
$m_{1} + 1,\text{  }m_{2} + 1,\ldots\text{, }m_{k} + 1$, respectively.
Without loss of generality, assume that
$\Lambda_{1} < \Lambda_{2} < \ldots < \Lambda_{k}$.

Theorem. Let the following conditions be met:

1\) matrix $A$ has eigenvalues of any multiplicity, whose real parts are
negative, and eigenvalues multiple of 1, whose real parts are zero;

2\) integral $\int_{1}^{+ \infty}{}s^{2m}\omega(s,r)\text{ds}$, where
$m = \text{max}\left\{ m_{i} \right\}$, is converging.

Then systems (1) and (2) are locally reducible.

#### 

References

1\. Voskresensky E.V*. On Lyapunov transformation groups*//
*Differentsialnyje Uravneniya* - 1996. - V. 32, No.11. - Pp. 1574-1575.

2\. Shamanaev P.A. *On Lyapunov transformations and motion stability*//
*Author's abstract of Ph.D. Thesis*. Saransk, 1997. -- 16p.

**STRUCTURE OF 3D STATIONARY ISOBARIC DOUBLE WAVES IN PERFECT
INCOMPRESSIBLE FLUID**

***V.E.Shemarulin***

(RFNC-VNIIEF, Sarov)

This paper pursues the structural studies of isobaric (inertial) flows
of perfect incompressible fluid (isobaric flows). The set of non-trivial
(non-constant) isobaric flows divides into two classes - rank 1 flows
(simple waves) and rank 2 flows (double waves) /1-4/. Reference /4/
completely describes the structure of simple waves that is 3D (rank 1)
and 2D non-stationary isobaric flows.

This paper contains an explicit geometrical structure description of 3D
isobaric stationary double waves. All such flows are proved to be the
combinations of domains of three basic flow types - shear, conic and
tangential flows (conic and tangential types with the general
description given in /4/). This result is completely similar to that for
the structure of isobaric simple waves obtained in /4/. The only
difference from rank 1 flows is that for rank 2 flows, the velocity has
a specific value on each line being the trajectory (streamline) of a
fluid particle in each plane (half-plane) naturally layering the flow
(see the definition of basic classes of isobaric flows in /4/) while for
rank 1 the flow is constant in each plane except for the flows along the
parallel lines where the velocity may also have an individual value on
each line.

Thus the results of this paper complete the local classification of
isobaric flows in the case of three independent variables - three space
or two space and one time variables.

The results obtained here were announced in /5/.

1.  1 Ovsyannikov L.V. *Lectures on gas dynamics fundamentals*. Moscow,
    *Nauka*, 1981. 368 p.

2.  Ovsyannikov L.V. *Isobaric gas motion. Differential equations*.
    1994Vol.30. N 10. P. 1792-1799.

3.  Bondarenko Yu.A. *Inertial 3D motion of inviscid incompressible
    fluid*. *VANT.Series: Mat.Mod.Fiz.Proc.* 1994. N3. P41-46.

4.  Shemarulin V.E. *Structure of isobaric simple waves in perfect
    incompressible fluid. VANT.Series: Mat.Mod.Fiz.Proc*. 2001. N3.
    P.26-36.

5.  Shemarulin V.E. *Structure of isobaric gas flows and perfect
    incompressible fluid. Doklady Akademii nauk*. 2002. Vol.383.N 2.
    P.206-210.

**NUMERICAL SIMULATIONS OF 2D GAS-DYNAMIC**

**AND ELASTIC-PLASTIC FLOWS ON ADAPTIVELY**

**EMBEDDING FRAGMENTED COMPUTATIONAL GRID**

***Yu.V.Yanilkin, A.V.Gorodnichev, S.P.Belyaev, L.I.Degtyarenko,***

***A.L.Stadnik, N.A.Khovrin, A.V.Volgin,***

(RFNC-VNIIEF, Sarov)

This paper describes program complex EGAK++ that is being developed for
simulations of 2D multi-material flows with severe strains, including
turbulent mixing. The complex employs a single Lagrangian-Eulerian
quadrangular regular computational grid allowing for multilevel
fragmentation in selected cells. The complex solves problems, in which
some small-scale flow features have to be separated, such as shock wave,
interface, etc.

For this purpose an approach is suggested that uses movable adaptively
embedding fragmented grids. Its idea is that movable Lagrangian-Eulerian
grid cells are dynamically cut and pasted to achieve a resolution degree
needed for description of one or another physical process or small flow
features.

The following requirements have to be met when cutting cells:

1.  The cell fragmentation method should be the same at all levels, new
    cells should be produced through cutting by straight lines drawn
    through middles of sides of larger cells.

2.  The number of the levels should not be larger than 5. (That is the
    minimum cell linear dimension can be 32 times as small as that of
    the main cell.)

3.  Neighboring cells can differ no more than by one level.

A novelty in the paper is that the approach is used for modeling
problems, in which a number of different-scale physical processes should
be taken into account at a time. This required development of programs
for automatic computational grid fragmentation depending on the
situation for simulations of a wide range of this kind of problems.

The programming in the complex involves an object-specific approach in
the C++ language. The paper discusses features both of the complex
arrangement and the difference schemes for simulations of gas-dynamic
and elastic-plastic flows on adaptively embedding grids, presents some
computed data.

The computations suggest that the complex EGAK++ is quite efficient and
promising.

[^1]: ^\*^