Модифицированная численная схема комплекса H3T

Modified Numerical Scheme for H3T Program Complex
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Гиззаткулов Н.М.
( N.M.Gizzatkulov)

ИПМ им. М.В.Келдыша РАН

Москва, 2003

Аннотация

Построена схема численного моделирования двумерного нестационарного течения теплопроводного газа в трехтемпературном приближении на криволинейных сетках. Дифференциальные уравнения модели расщеплены по процессам, газодинамическому и теплопроводному. Нелинейная система, полученная в теплопроводном этапе, линеаризуется методом Ньютона и решается пакетом SPARSKIT (решение систем алгебраических уравнений методом GMRES - обобщенный метод минимальных невязок). Выполнено сравнение данной схемы с явной схемой комплекса H3T, базирующейся на той же математической модели, получено полное совпадение результатов. Даны рекомендации по распараллеливанию схемы.

Abstract

Numerical scheme for numerical modelling 2D non stationary flow of the heat - conducting gas in three temperature approach on the curved grids is studied. Differential equations of the problem are splitted for the processes: gas dynamic process and heat - conductivity process. For linearizing nonlinear system obtained in the heat conductivity, it is used Nyuton approach. Algebraic equations are solved by GMRES approach using SPARSKIT software. Numerical scheme is verified by numerical tests.

Введение

В ИПМ им. М.В. Келдыша в течение двадцати лет разрабатывается многофункциональный программный комплекс H3T для расчета нестационарных двумерных уравнений газовой динамики с теплопроводностью на подвижных криволинейных сетках. Основное применение данного пакета - численное моделирование термоядерного синтеза в приближении теплопроводной газовой динамики. Расчетная область в H3T может быть достаточно сложной формы, до начала расчета она разбивается на элементарные области четырехугольной формы(ярусы). Элементарные области разбиваются двумя семействами координатных линий на криволинейные четырехугольники - ячейки сетки, стороны которых приближаются дугами окружностей (или отрезками прямых), сеточные функции задаются в центрах ячеек. Сетка строится на каждом шаге по времени и может быть достаточно сложной.
Разностная схема для системы дифференциальных уравнений газовой динамики с теплопроводностью строится на основе интегральных законов сохранения, записанных для каждой ячейки сетки с учетом ее движения. Алгоритм расчета основан на расщеплении сеточных уравнений на газодинамический, тепловой и релаксационный этапы. Газодинамический этап является обобщением схемы Годунова, предложенным в []. Расчет теплопереноса основан на вычислительной схеме, разработонной в []. Дискретизация оператора теплопроводности по пространству строится с помощью обобщения на случай произвольных сеток интегро - интерполяционного метода. Для получения схемы производится локальная аппроксимация решения многочленами, построенными по значениям решения на сеточных шаблонах, с последующей подстановкой многочленов в интегральное соотношение. Задавая определенным образом шаблон, вид многочлена и способ его построения, получаем схемы различного качества. По времени используется двухслойная дискретизации: переход на новый временной слой производится с помощью явно - итерационной схемы с чебышевским набором параметров, сохраняющей структуру связей, полученную при пространственной дискретизации.
Технические возможности ЭВМ предыдущих лет определили численные схемы, применяемые в H3T. При современных объемах памяти и тактовых частотах процессоров явная итерационная схема с Чебышевским набором параметров теплопроводного этапа является узким местом вычислительного комплекса, так как в областях с существенно различными коэффициентами теплопроводности счет приходится вести с обременительным ограничением на шаг по времени. Поэтому в данной работе предлагается модифицировать численную схему теплового этапа и сделать ее неявной. Далее будем называть эти схемы модифицированной и немодифицированной, соответственно. Также в работе приведено качественное сравнение результатов расчетов схем на тестовых задачах.

§ 1.  Система дифференциальных уравнений модели

Для трехтемпературной модели среда, в которой будут протекать рассчитываемые процессы, предполагается состоящей из трех видов частиц: ионов, электронов и фотонов
Среда описывается как газ с единой плотностью частиц r и общим вектором их скорости [(u)\vec]=(u,v). Температура T, давление p и удельная внутренняя энергия e (для единицы массы среды) ионов, электронов и фотонов предполагаются различающимися между собой в каждой точке пространства и времени. В связи с этим они будут отмечаться индексами i,e,f, соответственно.
Система дифференциальных уравнений для рассмативаемой модели среды в переменных Эйлера, записанная в форме законов сохранения, имеет следующий вид:

  r
t
 +div r
®
u
 
 =0,
(0.1)

 
t
 (r
®
u
 
 )+ div(r
®
u
 
   
®
u
 
 )+grad   p=0,
(0.2)

 
t
 (rei)+div(rei
®
u
 
 )+pidiv
®
u
 
 =
=div(Ki grad Ti)+cei(Te-Ti)+Qi,
(0.3)

 
t
 (ree)+div(ree
®
u
 
 )+pe div
®
u
 
 =
=div(Ke grad Te)+cei(Ti-Te)+cef(Tf-Te)+Qe,
(0.4)

 
t
 (ref)+div(ref
®
u
 
 )+pf div
®
u
 
 =
=div(Kf grad Tf)+cef(Te-Tf)+Qf,
(0.5)

Qi=Q1i+Q2i+Q3i,  Qe=Q1e+Q2e+Q3e,   Qf=Q1f.
Левые части этих уравнений описывают процессы, отвечающие за движение в газовой среде. Ki,  Ke,  Kf - коэффициенты теплопроводности. Соответствующие диффузионные слагаемые div(KgradT) описывают перенос тепла ионами и электронами и процесс диффузии энергии излучения в приближении "серой материи".
Следующие слагаемые, стоящие в правой части, имеют цель учесть релаксационные процессы обмена энергией между ионами и электронами и неравновесностью процесса переноса излучения фотонами. Участвующие в их описании коэффициенты cei,  cef также как и коэффициенты теплопроводности Ki, Ke, Kf могут зависеть от искомых функций, а также пространственных переменных и времени.
Q={Qi,  Qe,  Qf} - эти члены отвечают за энерговложение от различных источников в ионы, электроны и фотоны. Источники могут быть следующими:

Q1e,i,f=f(r,{e}i,e,f,x,r,t);
Q2i,e - энерговыделение от термоядерных реакций;
Q3i,f - энерговыделение от нейтронов, которое расчитывается в комплексе NPT[].
Система (1) - (5) замыкается заданием уравнений состояния для всех трех компонент среды, например, в виде:

ek
=
ek(r,Tk),  pk=pk(r,Tk),  k=i,e,
(0.6)
ef
=
sTf4,      pf=  sTf4
3
 .
Уравнения состояний для электронов и ионов должны удовлетворять известным термодинамическим тождествам:
r2  ek
r
 -pk+Tk  pk
Tk
 =0,    k=i,e.
а также начальным и граничным условиям на границе области течения газа.
В уравнении (2) p - суммарное давление частиц:
p = pi+pe+pf
(0.7)
Полное описание системы дифференциальных уравнений дано в [].

§ 2.  Разностная схема

Схема счета основана на декомпозиции(при необходимости) области решения на подобласти. Эта операция включает три этапа: геометрическое разбиение исходной расчетной области на подобласти(введением при необходимости дополнительных границ), движения границ, построение сетки в каждой подобласти. Разбиение исходной области на подобласти определяется, в сновном необходимостью ее представления набором "топологических четырехугольников". Эти "четырехугольники" должны в ходе расчета покрывать всю область искомого решения. Несколько подобластей могут объединяться в структуру, называемую ярусом, если при этом для объединения сохраняется топология четырехугольника, а сами подобласти составляют прямоугольную решетку. Ярус является основной вычислительной структурой: на каждом шаге по времени уравнения (1) - (5) решаются в ярусе с определенными краевыми условиями. На этапе счета теплопроводности несколько ярусов могут объединяться в счетные объекты называемые кластерами, в которых сквозным образом осуществляется расчет теплопереноса.
Каждая из подобластей(счетная область) является либо физической областью, либо ее частью. Принадлежностями области являются ячейки, границы, физические параметры(параметры уравнения состояния, коэффициента теплопроводности и т.д.). В каждом ярусе мы строим четырехугольную криволинейную сетку. Особенностью метода является процедура учета взаимосвязи между ярусами, не предполагащая гладкости сетки через границы ярусов. В данном параграфе будет описана модификация теплового этапа расчета.
Система дифференциальных уравнений
r  Ei
t
=
div(Ki grad Ti)+cei(Te-Ti)+Qi,
r  Ee
t
=
div(Ke grad Te)+cei(Ti-Te)+cef(Tf4-Te4)+Qe,
(0.8)
r  Ef
t
=
div(Kf grad Tf)+cef(Te4-Tf4)+Qf,
описывает теплопроводные и релаксационные процессы. Разностная аппроксимация имеет вид
^
E
 
i 
=
Ei+  t
r
 D(
^
Ti
 
 )+  t
r
 cei(
^
T
 
e 
 -
^
T
 
i 
 )+  t
r
 Qi,
^
E
 
e 
=
Ee+  t
r
 D(
^
Te
 
 )+  t
r
 cei(
^
T
 
i 
 -
^
T
 
e 
 )+  t
r
 cef(
^
T
 
 4
f 
 -
^
T
 
 4
e 
 )+  t
r
 Qe,
(0.9)
^
E
 
f 
=
Ef+  t
r
 D(
^
Tf4
 
 )+  t
r
 cef(
^
T
 
 4
e 
 -
^
T
 
 4
f 
 )+  t
r
 Qf.
Энергии на верхнем временном слое линеаризуются по температуре:
^
E
 
i 
 =ai
^
T
 
i 
 +bi,  
^
E
 
e 
 =ae
^
T
 
e 
 +be,  
^
E
 
f 
 =af
^
T
 
f 
 +bf;
соответствующие формулы для коэффициентов ak и bk (k=i,e,f) выписаны в []. Температуру [^(T)]e4 можно линеаризовать двумя способами: либо [^(T)]e4=z1m[^(T)]f4+z1n, где
a =  4cefTe3
ae
 ,  z1m =  a
1+a
 ,  z1n =  Te4
1+a
(0.10)
либо [^(T)]e4=z2m[^(T)]e+z2n, где
z2m = 4 Te3,  z2n = -3 Te4
(0.11)
Линеаризация (10) была получена А.В. Забродиным [] упрощением исходной дифференциальной системы.
Cистему (9) с учетом (10),(10) можно переписать в виде

ai
^
T
 
i 
 -
-
t
 
 D(
^
T
 
i 
 )+
-
cei
 
 (
^
T
 
i 
 -
^
T
 
e 
 )
=
-bi+Ei+
-
Q
 
i 
 ,
ae
^
T
 
e 
 -
-
t
 
 D(
^
T
 
e 
 )+
-
cei
 
 (
^
T
 
e 
 -
^
T
 
i 
 )-(1-z1m)
-
c
 
ef 
 Tf4
=
(0.12)
= -be+Ee-z1n
-
c
 
ef 
 +
-
Q
 
e 
 ,
af
^
T
 
 4
f 
 -
-
t
 
 D(
^
T
 
 4
f 
 )+(1-z1m)
-
c
 
ef 
 Tf4
=
-bf+Ef+z1n
-
c
 
ef 
 +
-
Q
 
f 
 ,
где [`(ck)]= [(t)/(r)]ck,  [`(Qk)]= [(t)/(r)]Qk   (k=ei,ef) и [`(t)] = [(t)/(r)].
Перепишем систему (12) в матричном виде. Теплопроводные члены
D([^(T)]k),  k=i,e,f, соберем в матрицах A1,  b1 вида
A1= ж
з
з
и
C11
0
0
0
C22
0
0
0
C33
 ц
ч
ч
ш 		,
(0.13)

b1= ж
з
з
и
D1
D2
D3
 ц
ч
ч
ш 		.
(0.14)
Остальные элементы образуют матрицу и вектор правых частей:

A2= ж
з
з
з
з
и
ai+
-
c
 
ei 
-
-
c
 
ei 
0
-
-
c
 
ei 
ae+
-
c
 
ei 
-(1-z1m)
-
c
 
ef 
0
0
af+(1-z1m)
-
c
 
ef 
ц
ч
ч
ч
ч
ш 		,
(0.15)

b2= ж
з
з
з
з
и
-bi+Ei+
-
Q
 
i 
-be+Ee-z1n
-
c
 
ef 
 +
-
Q
 
e 
-bf+Ef+z1n
-
c
 
ef 
 +
-
Q
 
f 
ц
ч
ч
ч
ч
ш 		,
(0.16)

^
T
 
 = ж
з
з
з
з
и
^
T
 
i 
^
T
 
e 
^
T
 
f 
ц
ч
ч
ч
ч
ш 		.
(0.17)
Таким образом, с учетом (13)-(16) система (12) будет иметь вид
(A1+A2)
^
T
 
 =b1+b2.
(0.18)
Решая систему (18) получим вектор температуры [^(T)] на момент времени t+t. Линейная система является плохим приближением исходной нелинейной системы, поэтому для того, чтобы гарантировать хорошую аппроксимацию, необходимо организовать итерационный процесс по методу Ньютона []. А именно, температуру ([^(T)]e)4 линеаризуем согласно (11) в окрестности [^(T)]ek-1. Тогда система (9) с учетом (10) и линеаризации переписывается в виде

ai
^
T
 
 k
i 
 -
-
t
 
 D(
^
T
 
 k
i 
 )+
-
cei
 
 (
^
T
 
 k
i 
 -
^
T
 
 k
e 
 )
=
-bi+Ei+
-
Q
 
i 
 ,
ae
^
T
 
 k
e 
 -
-
t
 
 D(
^
T
 
 k
e 
 )+
-
cei
 
 (
^
T
 
 k
e 
 -
^
T
 
 k
i 
 )-
-
c
 
ef 
^
Tk
 
 4
f 
 +z2m
-
c
 
ef 
-
T
 
 k
e 
=
(0.19)
= -be+Ee-z2n
-
c
 
ef 
 +
-
Q
 
e 
 ,
af
^
Tk
 
 4
f 
 -
-
t
 
 D(
^
Tk
 
 4
f 
 )-z2m
-
c
 
ef 
^
T
 
 k
e 
 +
-
c
 
ef 
^
Tk
 
 4
f 
=
-bf+Ef+z2n
-
c
 
ef 
 +
-
Q
 
f 
 .
Запишем систему (19) матричном виде. Для этого нужно внести поправки в матрицу A2 и вектор b2:

A2k= ж
з
з
з
з
и
ai+
-
c
 
ei 
-
-
c
 
ei 
0
-
-
c
 
ei 
ae+
-
c
 
ei 
 +z2m
-
c
 
ef 
-
-
c
 
ef 
0
-z2m
-
c
 
ef 
af+
-
c
 
ef 
ц
ч
ч
ч
ч
ш 		,
(0.20)

b2k= ж
з
з
з
з
и
-bi+Ei+
-
Q
 
i 
-be+Ee-z2n
-
c
 
ef 
 +
-
Q
 
e 
-bf+Ef+z2n
-
c
 
ef 
 +
-
Q
 
f 
ц
ч
ч
ч
ч
ш 		,
(0.21)

^
Tk
 
 = ж
з
з
з
з
и
^
Tk
 
i 
^
Tk
 
e 
^
Tk
 
f 
ц
ч
ч
ч
ч
ш 		.
(0.22)
Таким образом, учитывая (13),(14),(20)-(22) система (19) будет иметь вид
(A1+A2k)
^
Tk
 
 =b1k+b2k.
(0.23)
Выполняя k итераций, получим вектор температур [^(Tk)], где [^(T0)]=[^(T)].
На первом этапе казалось, что контроля количества итерации по невязке достаточно для решения нелинейной системы. Практика показала, что в задачах физики плазмы при прохождении фронта тепловой волны границы раздела двух существенно разных веществ абсолютные невязки после Ньютоновских итераций порядка(10-3 -101), при начальной абсолютной невязке порядка (108 -1010). Возникает характерная картина - прогрев газа в режиме тепловой волны, фронт которой движется с конечной скоростью. Расчет такого режима затруднен тем, что температура существенно негладкая около точки фронта. Поэтому необходимо было ввести автомат, который вел бы контроль шага по времени от величины невязки. Известно, что метод Ньютона сходится квадратично []:
rk Ј C-1(C  r0)2k; ||fxx(x)||   Ј   C2;  ||fx-1(x)||   Ј   C1;  C=C2C12;
Теперь можно более точно указать, насколько хорошим должно быть начальное приближение x0, чтобы процесс заведомо сходился. Очевидно, для этого достаточно выполнения неравенстваю
C||f(x0)|| Ј   q <   1.
Исходя из вышесказанного был предложен следующий алгоритм. На каждой Ньютоновской итерации рассчитываются невязки ri,       re,      rf нелинейной системы (8), в норме L2 по счетному пространству. Итерации по нелинейности выполняются до тех пор, пока одна из невязок убывает в заданное число раз (один, два порядка или более) или не достигается величина nmax=10(максимальное количество Ньютоновских итераций). В случае, когда одна из невязок после итерации отличается от средней величины невязки, выработанной за предыдущие сто шагов, более чем в заданное число раз(103 - 104) то шаг по времени уменьшается в двое и процедура измельчения повторяется до тех пор, пока невязка не достигает средней велечины. В [] приведено много модифицированных алгоритмов Ньютона, которые оптимальнее вырабатывают шаг приращения временной переменной. На данном этапе реализован самый простой из всех вариантов для проведения методических расчетов и верификации численной схемы.
Для решения системы линейных алгебраических уравнений с разряженной матрицей использовался программный комплекс Sparskit [], использующий метод обобщенной минимальной невязки(GMRES). GMRES считается одним из самых эффективных методов решения несимметричных линейных систем. В пакете использовались следующие настройки: e = 2.22D-20 (локальный толеранс),       maxits=200,      im=15, где im - размерность подпространства Крылова, maxits - максимальное количество итераций в методе обобщенных невязок. В процессе расчетов Sparskit зарекомендовал себя как надежный алгебраический решатель. К тому же, данный пакет имеет параллельный вариант Sparslib, что является существенным заделом на будущее. О пакете Sparslib и его возможностях будет написано ниже.

§ 3.  Тестовый расчет 1

Целью данного расчета является сравнение модифицированной и не модифицированной схем счета на "реальной" задаче. Рассматривается цилиндрическая мишень, состоящая из нескольких цилиндрических слоев разных веществ бесконечной длины с открытыми торцами Рис. . В начальной геометрии конструкция обладает осевой симметрией с сечением, перпендикулярным оси вращения.
Figure 0.1: Геометрия задачи 1

Внешнее энерговложение осуществляется в область M2 по следующей формуле:
Q( t )
=
 2gc03 Gx - [(2g)/(g+1)]
( g- 1 )2r0 m2
{ ( 1 -x - [(g- 1)/(g+ 1)] )
  
+
 x - 1
g+ 1
 й
л 		 V2 ( 0 )
aV0 ( 0 )
 +( g- 1 ) + 2x- ( x+ 1)x[ 2/(g+ 1)] щ
ы 		ь
э
ю 		·zn
(0.24)

a
=
 2( m1 + m2 + m3) + ( g+ 1 )m0 /g
(g- 1 )( 2m3 + m2 )
 ,   
G
=
 m0 ( 1 - 2k2 )
g
+  2k1
g+ 1
 +  g+ 1
g2
 k3 m02     ,
(0.25)

k1=m1+  m2
3
 +
m1+  m2
2
(m3+  m2
2
 )2
 й
л 		 m2
3
 ж
и 		m3+  m2
4
 ц
ш 		+m1 ж
и 		m3+  m2
3
 ц
ш 		щ
ы 		,
k2 = -  1
( m3 + m2 /2 )2
 й
л 		 m2
3
 ж
и 		m3 +  m2
4
 ц
ш 		+ m1 ж
и 		m3 +  m2
3
 ц
ш 		щ
ы 		,   
k3 =  m3 + m2 /3
2( m3 + m2 /2)2
 ,
(0.26)
где mi=(ri-ri-1)ri плоские массы областей Mi,      i=0,1,2,3;       r0=20,      r1=21;     r3=39,      r4=44,4; V2(0)=[(m2)/(r2)]; x = 1-[(c0 t)/(r0)];        0 < x Ј 1, r0 - внешний начальный радиус области DT, c0 - начальная скорость звука в области DT, t - текущий момент времени,g - показатель адиабаты области с массой m2.
В начальный момент в первой области M0 задается давление P|(t=0)=[(r0 c02)/(g1)](g1 - показатель адиабаты в области с массой m0). Такое же давление задается и в остальных областях. Рассчитывается движение системы (уравнения газодинамики) и кинетика термоядерных реакций. Уравнения состояния - идеальный газ с показателем g = [ 5/3] . В области Pb или Be задается энерговложение Q(t) на основе вышеприведенной формулы. На поверхности области поток тепла полагается равный 0.
Ниже будет выполнено качественное сравнение немодифицированной схемы с модифицированной схемой, в которой количество итераций по нелинейности на временном шаге зависит от поведения невязки в L2 норме и модифицированной схемой с постоянным количеством итераций = 7, отметим, что семь итераций взяты с запасом для подтверждения алгоритма выбора количества итераций на временном шаге. Далее будем называть эти три метода как метод 1, 2 и 3, соответственно.
Сравним качественно три этих метода. На следующих рисунках приведено распределение газодинамических велечин по области на момент времени t=103 для трех различных методов счета. Непрерывной линией обозначена схема 1, пунктирной схема 2, штрих пунктирной схема 3.
Figure 0.2: На левом рисунке: распределение плотности вдоль радиуса;
На правом рисунке: изображен фрагмент левого рисунка.

Ниже в таблице приведено распределение полной энергий по счетным областям для трех различных вариантов счета на момент времени t=103.5 на сетке 3 на 160. Сравнительная таблица расчетов:
Вариант расчета 1 2 3
E полнаяE полнаяE полная
M00.11035E+060.11388E+060.11391E+06
M10.17615E+070.14464E+070.14471E+07
M20.40556E+070.20010E+080.20001E+08
M30.32870E+070.40551E+070.40541E+07
общая по областям0.26029E+080.25625E+080.25616E+08
общее время счета в минутах 1105066
Figure 0.3: На левом рисунке: распределение полной энергии вдоль радиуса;
На правом рисунке: изображен фрагмент левого рисунка.
Figure 0.4: На левом рисунке: распределение полного давления вдоль радиуса;
На правом рисунке: изображен фрагмент левого рисунка.
Как видно из рисунка 0.2 левого плотности совпадают. На правом рисунке рисунка 0.2 в масштабе приведена окрестность максимума по плотности, где можно разглядеть небольшую погрешность методов. На рисунках: рис. 0.3 приводится распределения полной энергии, рис. 0.4 распределения давления вдоль радиуса. На рис. приводится распределения энергии электронов вдоль радиуса.
Figure 0.5: Распределение энергии электронов вдоль радиуса
Figure 0.6: На левом рисунке: распределение температуры электронов вдоль радиуса;
На правом рисунке: изображен фрагмент левого рисунка.
Figure 0.7: Зависимость процессорного времени от переменной t
На рис. 0.6 приводится распределения температуры электронов вдоль радиуса. На рис. 0.7 приведена зависимость времени счета от временной переменной: такое поведение кривой можно объяснить тем, что в начальный момент времени градиенты теплопроводных членов не велики и шаг по пространству достаточно большой, соответственно количество итераций требуемых для решения задачи явным методом не велико. Далее, по мере сжатия мишени, ситуация меняется, схема 1 начинает отставать и даже проигрывает методу, где итерации по нелинейности постоянны и взяты с запасом. Как видно из графиков время затрачиваемое алгебраическим алгебраическим солвером на разрешение одного шага постоянно.
Figure 0.8: Зависимость количества Ньютоновских итераций от шага по времени
В ходе расчетов было обнаружено, что существуют моменты времени, когда невязки ri,  re,  rf в точках где располагается фронт тепловой волны имеют высокий порядок (102 - 104), как отмечалось ранее, в данных точках происходит измельчение шага по времени. В данном расчете из 9000 шагов схемы данная проблема возникла в 12 из них, причем шаг по времени уменьшался в среднем в 32 раза.
Ниже приведены графики невязок для ионов, электронов и фотонов в L2 норме по счетному пространству рис. , рис. .
Figure 0.9: Зависимость величины невязки от шага по времени
Figure 0.10: Зависимость величины невязки от шага по времени
Так же был проведен расчет на вдвое измельченной сетке. Результаты расчета приведены на момент времени t=103.
Явная схема:
Сетка 4x160 4x320
E полнаяE полная
M00.14785E+060.13150E+06
M10.17937E+070.18884E+07
M20.19309E+080.19243E+08
M30.44719E+070.45281E+07
общая по областям0.25722E+080.25791E+08
Неявная схема:
Сетка 4x160 4x320
E полнаяE полная
M00.11462E+060.10485E+06
M10.17195E+070.17987E+07
M20.19326E+080.19280E+08
M30.45428E+070.45972E+07
общая по областям0.25703E+080.25781E+08

§ 4.  Параллельные аспекты схемы

При численном моделировании задач "Тяжелоионного термоядерного синтеза" время счета является критическим параметром. Поэтому вопрос написания Параллельного варианта модифицированной схемы возникает сам по себе. Многопроцессорная программа для не модифицированной схемы[] была реализована и подробно описана в работе []. В данной работе построена схема, которая приводит к необходимости решения системы алгебраических уравнений большой размерности. Трудность заключается в том, что написание параллельного солвера для системы уравнений достаточно не тривиальная задача. Один из методов выхода из тупиковой ситуации - это использование параллельных пакетов, которые представлены на сервере www.parallel.ru. Изучив по отдельности возможности каждого из пакетов, выбор пал на Sparslib. Ниже приводится таблица зависимости времени счета для модельной задачи (2 - х мерное уравнение теплопроводности на прямоугольной сетке размером 200 на 200, первая краевая задача) от количества процессоров для данного программного продукта.
Параметры/Кол-во проц. 1 4 9
Время счета в сек.3.29091.20840.8093
Коэфф. ускорения12.74
Как видно из таблицы, пакет вполне пригоден для использования в качестве параллельного солвера. Таким образом в дальнейшем предполагается реализовать параллельный вариант комплекса H3T с модифицированной схемой.

Заключение

Построена неявная схема численного моделирования двумерного нестационарного течения теплопроводного газа в трехтемпературном приближении на криволинейных сетках. Линеаризация уравнений газовой динамики с теплопроводностью выполнена по методу Ньютона. Для решения системы линейных алгебраических уравнений с разряженной матрицей использовался программный комплекс Sparskit [], использующий метод обобщенной минимальной невязки(GMRES). Ньютоновский итерационный алгоритм построен с автоматическим выбором шага по временной переменной. На основе тестового расчета выполнено сравнение данной схемы с явной схемой, базирующейся на той же математической модели. Получено совпадение результатов не модифицированной и модифицированной схем.

Bibliography

[]
В.Т. Жуков, А.В. Забродин, О.Б. Феодоритова Схема решения нестационарных двумерных уравнений газовой динамики с теплопроводностью на подвижных криволинейных сетках. Препринт N 18 М.: ИПМ им. М.В. Келдыша РАН, 1991.
[]
В.Т. Жуков, О.Б. Феодоритова Разностные схемы для уравнений теплопроводности на основе локальных среднеквадратичных приближений. Препринт N 97 М.: ИПМ им. М.В. Келдыша РАН, 1989.
[]
В.Т. Жуков, А.В. Забродин, В.С. Имшенник ,О.Б. Феодоритова Численное моделирование мишени тяжелоионного термоядерного синтеза в приближении теплопроводной газодинамики. Препринт N 41 М.: ИПМ им. М.В. Келдыша РАН, 1993.
[]
А.В. Забродин, Г.П. Прокопов. Методика численного моделирования двумерных нестационарных течений теплопроводного газа в трехтемпературном приближении. ВАНТ, сер. Математическое моделирование физических процессов. 1998, вып.3.
[]
Г.В. Долголева, А.В. Забродин. Разработка термоядерных мишеней на основе реализации концепции безударного сжатия. Аэромеханика и газовая динамика, 2002 г., N2, стр. 48 - 54.
[]
В.С. Рябенький. Введение в вычислительную математику: Учеб. пособие: Для вузов. - М.: Физматлит, 1994.
[]
Р.П. Федоренко Введение в вычислительную физику: Учеб. пособие: Для вузов. - М.: Изд - во Моск. физ. - техн. ин - та, 1994.
[]
SPARSKIT: a basic tool-kit for sparse matrix computations (Version 2). - University of Minnesota Department of Computer Science and Engineering. www.cs.umn.edu/Research/arpa/SPARSKIT

File translated from TEX by TTH, version 3.40.
On 14 May 2004, 17:10.