1. Введение *

2. Математическая постановка задачи об отражении наклонной ударной волны и метод ее решения. *

3. Основные конструкции языка Норма *

4. Программа для распределенной системы: общие проблемы и их автоматическое решение в системе Норма *

5. Разработка программы: от Фортрана к Норме *

6. Трансляция и запуск Норма-программы. *

7. Литература *

8. Приложение 1. Текст последовательной Фортран-программы *

9. Приложение 2. Текст Норма-программы *

10. Приложение 3. Текст Фортран-программ, используемых в Норма-программе *

1. Введение

Цель данного материала – достаточно подробно описать основные приемы разработки параллельной программы на языке Норма [1-3]. Ориентирован этот материал в первую очередь на прикладного специалиста, который сталкивается с необходимостью разработки параллельной программы для решения своей расчетной задачи.

Надо отметить, что в настоящее время имеется много различных методов и средств разработки параллельных программ (с ними можно познакомиться, например, в [4]). Эти способы требуют, в различной форме, достаточно сложной и трудоемкой работы программиста по выявлению параллелизма и представлению этого параллелизма в программе. Надежды на автоматическое распараллеливание последовательной программы пока не оправдываются: для параллельных вычислительных систем с общей памятью распараллеливающие компиляторы не обеспечивают эффективное распараллеливание, а для параллельных вычислительных систем с распределенной памятью проводятся лишь исследовательские работы в данном направлении.

Все это приводит к тому, что у прикладного специалиста возникают трудности при освоении области параллельных вычислений, которые часто осложняются наличием у него опыта "простой" разработки последовательных программ и попытками использовать этот опыт при разработке параллельных программ.

Программирование на языке Норма также требует усилий – не должно создаваться иллюзии, что применение этого языка снимает все проблемы параллельного программирования, конечно, это не так. Самое очевидное соображение против использования этого языка – необходимость изучать этот новый язык. Однако усилия по написанию Норма-программы, по нашему мнению, оправданы - вся работа с программой на многопроцессорных вычислительных системах (изменение размеров сетки, изменение числа процессоров, изменение выходного языка) поддерживается системой, а проблемы по синтезу этой параллельной программы остаются только у Норма-компилятора, а не у программиста. Кроме того, имеется положительный опыт применения этого языка для решения практических задач математической физики (например, [5-9]).

В качестве исходной информации для разработки Норма-программы взята последовательная программа на Фортране 77 для численного решения задачи об отражении сильной ударной волны от стенки, предоставленная С.Д.Устюговым (ИПМ РАН). Следует отметить, что такой подход является не вполне правильным - разработку Норма-программы следует вести по расчетным формулам, а не по готовой последовательной программе, при этом система Норма позволяет получать по одной и той же Норма-программе как параллельную, так и последовательную реализацию. Перевод последовательной программы на язык Норма является трудоемким процессом — требуется, по крайней мере, детальный анализ исходной последовательной программы, чтобы понять, а что же необходимо вычислить. При этом приходится "распутывать" информационные зависимости, часто не связанные с существом расчета, а возникшие как результат определенного стиля программирования и богатства возможностей Фортрана (экономия памяти, косвенная адресация, массовое использование COMMON-блоков и т.п.). Тем не менее, мы опишем процесс перехода от последовательной Фортран-программы к Норма-программе, поскольку эта ситуация является типичной – прикладной специалист чаще всего предпочитает именно такую форму представления необходимых расчетов.

Сначала будет описана математическая постановка задачи об отражении волны и метод ее решения (п.2), затем будут кратко определены основные конструкции языка Норма (п.3), далее рассмотрены общие вопросы, возникающие при программировании для распределенных вычислительных систем и способы их решения в системе Норма (п.4), и после этого описана сама методики перехода от рассматриваемой Фортан-программы к Норма-программе (п.5).

2. Математическая постановка задачи об отражении наклонной ударной волны и метод ее решения

Рассматривалось численное решение тестовой задачи: отражение наклонной ударной волны от стенки. Расчетная область представляла собой двумерную прямоугольную область размером . Плоскость ударной волны наклонена под углом к оси . Нижний край ударной волны находится в точке , верхний край в точке . Значения всех физических величин на левой границе и на нижней границе от до всегда равны значениям величин за фронтом ударной волны, в области, определенной соотношением . В области определены отражающие граничные условия. На правой границе удерживались значения величин, определённые в начальный момент, перед фронтом ударной волны. На верхней границе, используя определение положения фронта ударной волны в момент времени по формуле - число Маха, - скорость звука), для брались значения величини за фронтом ударной волны, для - перед фронтом ударной волны.

Начальные значения физических величин в расчётной области определялись с помощью соотношений Ренкина-Гюгонио для ударной волны с числом Маха при показателе адиабаты и имели вид



перед фронтом ударной волны и



за фронтом ударной волны.

Моделирование во времени процесса отражения ударной волны проводилось численным решением системы уравнений идеальной гидродинамики, записанной в дивергентной форме



где , , и - вектор консервативных переменных, , - вектора потоков, - плотность, - компоненты скорости, - давление, - полная энергия, - удельная внутренняя энергия.

Эта система решалась с помощью явной консервативной конечноразностной TVD схемы Годуновского типа второго порядка по пространству и времени [10]









где - правый собственный вектор матрицы , - правый собственный вектор матрицы . Элементы вектора находятся по формуле



.

Значения есть собственные значения матрицы , соответствующие характеристичеким скоростям и определялись для некоторого среднего состояния по известным значениям и по методу Рое [11]. Функция есть энтропийная коррекция к и имеет вид [12]:

,

где - константа, равная в расчёте .

В вычислениях использовалась функция-лимитер

.

Численный поток определялся по аналогичным формулам.

Продвижение во времени от к осуществлялось методом Рунге-Кутта третьего порядка сохраняющего TVD свойство исходной конечноразностной схемы [13].





Шаг по времени определялся из условия Куранта



В вычислениях число Куранта бралось равным 0.5. Расчёты проводились на равномерной сетке размерностью .

3. Основные конструкции языка Норма

Специализированный декларативный язык нового поколения Норма [1-3] (разработка ИПМ им.М.В.Келдыша РАН) позволяет описывать решение широкого класса задач математической физики и при этом использовать привычные для прикладного специалиста понятия. Такое описание не ориентировано на конкретную архитектуру компьютера, поэтому оно предоставляет большие возможности для выявления естественного параллелизма и организации вычислений. Программа на языке Норма фактически является представлением математических расчетных формул, полученных прикладным специалистом.

1. Способ задания областей

Понятие области введено в языке Норма для представления понятия индексного пространства. Неформально, область – это множество значений, которые могут принимать (или принимают) индексы величин, используемых в расчетных формулах. Точнее, область - это совокупность целочисленных наборов {i₁,....,i_n }, n>0, i_j>0, j=1,...,n, каждый из которых задает координаты точки n-мерного индексного пространства. С каждым направлением (осью координат) n-мерного пространства задачи связывается уникальное имя - имя индекса (имя оси координат индексного пространства).

Следует отметить, что область определяет значения координат точек индексного пространства, а не значения расчетных величин в этих точках. Например, если требуется вычислить значения величины Y_i,j, i, j =1,...,n на некоторой сетке X_i,j, i,j =1,...n, заданной, например, формулой X_i,j =F(h,i,j), (F - заданная функция, h - заданный параметр), то следует:

1. описать область, состоящую из точек (i,j), i,j = 1,...n;

(2) описать на этой области величины X и Y;

(3) задать на этой области правило вычисления значений сетки X_i,j : X_i,j = F(h,i,j) и правило вычисления значений Y_i,j : Y_i,j = G(X_i,j ) (функция G также некоторым способом задана).

В языке Норма область может иметь имя. Индексы областей специально не описываются - они вводятся при определении областей. Области используются при описании величин, определенных на области, при задании вычислений в операторах ASSUME, при описании входных или выходных величин, при задании областей фактических параметров в вызовах разделов или функций, в функциях редукции.

Ключевым понятием при описании прямоугольных областей является понятие одномерной области. Одномерная область служит для задания диапазона точек на некоторой оси координат индексного пространства. В простейшем случае при описании одномерной области указывается имя одномерной области, имя индекса и границы изменения значений индекса, например ниже описаны одномерные области с именами oi, oj, om:

Имя-одномерной-области может использоваться для ссылки на эту область. Имя-индекса есть индексная переменная, множество значений которой определяется диапазоном. Границами диапазона являются целые положительные константные выражения, построенные из целых констант, параметров области и арифметических операций. Требуется, чтобы в программе всем параметрам области были присвоены конкретные значения в описании параметров области. Для этого в языке предусмотрена конструкция DOMAIN PARAMETER, например, DOMAIN PARAMETER Nx=388, MyJ=4, Ny=122. Таким образом, области oi и oj представляют собой наборы точек 1..388 и 4..123 соответственно.

Многомерная область строится при помощи операции “;” произведения прямоугольных областей. Например, описания трехмерной области oijm и двумерной oijstep, заданных с помощью операции произведения одномерных областей, могут быть описаны следующим образом:

1. Способ задания величин

Скалярные величины (скаляры) и величины на области относятся к арифметическим величинам. Описание ставит в соответствие каждой арифметической величине уникальное в текущей программной единице имя величины, а также задает тип величины: REAL, INTEGER или DOUBLE (по умолчанию - тип REAL). Пример описания скаляров:

VARIABLE dt0,dk REAL.

Каждая величина на области при описании связывается с указанной в описании областью. Эта область определяет имена индексов, которые могут использоваться в индексных выражениях при обращении к данной величине.

Например, описание

VARIABLE qu DEFINED ON oijm

говорит о том, что определена величина qu на области oijm, то есть эта величина имеет индексы i,j,m и может принимать значения во всех точках области oijm.

Если индексное выражение у величины не содержит смещения, то оно может быть опущено, например, записи qu[i, j, m-1] и qu[m-1] эквивалентны.

Если значение некоторого индекса задается константой, необходимо явно указать, к какому направлению относится константа, например, обращение к величине qu_i,j,m с индексами i=5,j=3,m=2 выглядит следующим образом: qu[i=5,j=3,m=5].

Если необходимо значения одного индекса связать со значениями другого, необходимо сделать явное указание, например, диагональные элементы в плоскости m=2 можно определить следующим образом: qu[i,j=i,m=2].

Существенно, что в языке Норма введено ограничение на вид индексных выражений - допускаются выражения только вида , где - имя индекса, а - константное выражение целого типа.

2. Задание вычисления переменной

В Норме определены три вида операторов: скалярный оператор, оператор ASSUME и оператор COMPUTE. Операторы предназначаются для описания вычислительных действий, необходимых для решения задачи.

Скалярный оператор предназначен для вычисления арифметических значений скаляров. По существу, это аналог оператора присваивания традиционных языков программирования.

1. Оператор ASSUME

Оператор ASSUME используется для описания вычисления арифметических значений величин, определенных на областях. Оператор ASSUME записывается следующим образом:

FOR <область> ASSUME <соотношение> .

и является запросом на вычисление соотношения во всех точках области. В результате выполнения оператора ASSUME величина, указанная в соотношения, получает значения во всех точках области, указанной в заголовке оператора, однако каким образом это вычисление реализуется – последовательно, параллельно, в каком порядке – в языке не определяется. Например, в результате выполнения оператора

FOR ((i=3..Nx-2);(j=3..Ny-2);(m=1..4)) ASSUME qu=1.0.

во всех точках области, указанной после слова FOR , величине qu будет присвоено значение1.0 (индексы у величины qu в записи соотношения опущены ).

Если на одной и той же области надо вычислить несколько величин, можно задать несколько соотношений, например,

FOR onedomain ASSUME qu=1.0; e=p+v; u=w[i-1]/2.

После ключевого слова FOR в операторе ASSUME могут быть указаны несколько областей, что позволяет сократить запись в случае вычисления одних и тех же формул на различных областях. Например, пусть определены две области:

og1:(i=3..Nx-2);(j=Ny-1..Ny);(m=1..4)).

og2:(i=3..Nx-2);(j=1..2);(m=1..4)).

Тогда оператор

FOR og1,og2 ASSUME qu=qu1.

описывает то же вычисление, что и два оператора

FOR og1 ASSUME qu=qu1.

FOR og2 ASSUME qu=qu1.

Язык Норма является языком с однократным присваиванием (single assignment language), величины могут принимать значения только один раз, переприсваивание значений невозможно по определению. Таким образом, некоторое внешнее сходство оператора ASSUME с оператором присваивания не должно вводить в заблуждение.

2. Оператор COMPUTE

Оператор COMPUTE вызова раздела в определенном смысле является аналогом понятия вызова подпрограммы в традиционных языках программирования, и, по существу, является обобщением понятия соотношения, используемого в операторе ASSUME, так как дает возможность получать несколько величин-результатов.

При задании оператора COMPUTE надо задать значения величин, которые являются исходными данными для вычисления (фактические in-параметры) и имена величин, которые являются результатами вычисления (фактические out-параметры). Например, оператор

COMPUTE step(gam1,gamma,c0,dx,qu ON oijm RESULT umin ON oijstep).

является запросом на вычисление с именем step (подробности которого описаны в разделе с именем step), с in-параметрами gam1,gamma,c0,dx,qu (последний параметр – все значения величины qu на области oijm), по которым будут вычислены значения величины-результата umin (величина umin будет вычислена во всех точках области oijstep).

Фактическими in-параметрами раздела могут быть арифметические выражения, имена внешних разделов и внешних функций, либо величины на областях.

Фактическими out-параметрами могут быть скаляры, величины с индексами (быть может, заданными правилом умолчания индексов), величины на областях.

Вычисление раздела возможно, если все фактические in-параметры уже вычислены (приняли некоторые значения), побочный эффект в языке Норма невозможен.

3. Оператор COMPUTE в операторе ASSUME

Оператор COMPUTE может использоваться для задания сложного соотношения в операторе ASSUME, например

FOR oijstep ASSUME

COMPUTE step1(gam1,gamma,c0,dx,qu ON oijm/(i=i,j=j) RESULT umin).

Эта запись означает, что в каждой точке области oijstep надо провести вычисление step1, (подробности которого описаны в разделе с именем step1), с in-параметрами gam1,gamma,c0,dx,qu (последний параметр – значения величины qu на динамической области), по которым будет вычислено значение единственного результата - величины umin.

Фактический in-параметр qu ON oijm/(i=i,j=j) задает значения величины qu на динамической области, то есть области, которая меняется при различных значениях индексов из заголовка oijstep оператора ASSUME. Неформально, для каждого фиксированного значения индексов (i,j) области oijstep, например, i=i0, j=j0, при вызове раздела step1 в качестве фактического параметра передается рабочая величина R(1:4), такая, что R(1)=qu(i0,j0,1); R(2)=qu(i0,j0,2); R(3)=qu(i0,j0,3); R(4)=qu(i0,j0,4).

В целом, осуществление вычисления step1 в каждой точке области oijstep позволяет вычислить значения величины-результата umin во всех точках области oijstep, то есть вычисление step из предыдущего пункта и вычисление step1 в операторе ASSUME дадут одинаковый результат - величина umin будет вычислена во всех точках области oijstep. Разделы step и step1 при этом, конечно, будут отличаться.

Таким образом, в п.3.3.2 и п.3.3.3 приведены два способа описания вычислений с помощью оператора вызова раздела.

Отметим, что если вызов раздела находится в теле оператора ASSUME, то параметрами-результатами не могут быть скаляры и величины на статических областях (которые не меняются при различных значениях индексов из заголовка оператора ASSUME) - это заведомо приводит к переприсваиванию, что недопустимо в языке Норма.

4. Интерфейс с программами, написанными на Фортране

Подпрограммы и функции, написанные на Фортране, могут быть вызваны из программы, написанной на языке Норма, при помощи обычных средств вызова разделов и обращения к функциям. При этом осуществляется контроль переприсваивания при анализе фактических параметров.

В системе Норма также реализована подсистема конфигурации программ, которая, в частности, дает возможность проверки соответствия фактических и формальных параметров как для Норма-разделов, так и в рамках интерфейса Нормы и Фортрана.

 

3. Итерационный цикл

В языке Норма для описания итерационных вычислений используется специальная конструкция ITERATION, которая имеет следующий вид:

ITERATION <список итерируемых величин> ON <индекс итерации>.

INITIAL <индекс итерации>=0:

<операторы присваивания начальных значений итерируемым величинам>

EXIT WHEN (<условие выхода из итерации>).

END ITERATION <индекс итерации>.

Граничные значения итерируемой величины задаются при помощи обычных операторов языка Норма внутри блока BOUNDARY...END BOUNDARY. Эти значения считаются неизменными на протяжении всей итерации и определены на каждом шаге итерации.

Начальные значения для итерируемых величин задаются операторами языка Норма при помощи блока INITIAL...END INITIAL.

Индекс итерации может использоваться в списке индексов, задаваемом для итерируемой переменной, то есть фактически можно считать, что итерируемая величина в пределах итерации имеет дополнительный итерационный индекс по фиктивному направлению. Правило умолчания индекса справедливо и для итерационного индекса. Например, запись qu[itime-1], где qu – итерируемая величина, itime – индекс итерации, задает значения величины qu с предыдущего шага итерации. Обычно в пределах итерации сохраняются значения итерируемых переменных для двух шагов итерационного цикла: текущего и предыдущего.

1. Программа для распределенной системы: общие проблемы и их автоматическое решение в системе Норма

Будем использовать следующую простую модель задачи. Заданы статические регулярные области (сетки), в точках сеток определены расчетные величины (скорость, давление, температура и т.п.) и вспомогательные величины, которые вводятся по мере необходимости при определении расчетных формул. Для всех этих величин определены расчетные формулы, определяющие необходимые в задаче вычисления. С точки зрения объема вычислений, точки, лежащие внутри и на границах сеток, обычно различаются, однако для простоты будем считать, что в каждой точке сетки проводятся одинаковые (по времени) вычисления. Задачу распараллеливания для распределенных систем будем рассматривать как задачу статического распределения в модели параллелизма по данным: вычисления выполняются на том процессоре, где находятся вычисляемые данные. Таким образом, необходимо распределить точки сеток (расчетные и вспомогательные величины) между процессорами вычислительной системы и обеспечить доступ к значениям величин, используемым в расчетных формулах. Будем предполагать также, что каждый процессор должен работать по одной и той же программе (типичная ситуация при использовании модели передачи сообщений и библиотеки MPI, поддерживающей программирование в этой модели).

При построении фрагмента параллельной распределенной программы, реализующей параллельное вычисление на области O некоторой переменной X, определенной на области D, необходимо определить:

1) номера процессоров, которые должны выполнять данный фрагмент программы,

2) границы изменения параметров циклов, обеспечивающих обход области O,

3) описание массива X,

4) вид и расположение операторов передачи сообщений (обменов), которые обеспечивают доступ к значениям величин, используемым в расчетных формулах для вычисления X.

Отметим, что задачи 1)-4) решаются автоматически в процессе компиляции Норма-программы, что сокращает объем технической, а часто и интеллектуальной работы по программированию.

Проиллюстрируем на примере, в чем заключается решение задач 1)-4). Пусть надо вычислить следующие соотношения:

1. U(i,j)=i+j, для i=1..N; j=1..N.

2. W(i,j)=U(i,j)-1, для i=1..N; j=1..N.

3. V(i,j)= U(i-1,j)+1, для i=25..N-25; j=1..N.

Пусть N=100 и распараллеливание следует проводить на линейку из 10 процессоров по направлению i. На языке Норма эти вычисления описываются следующим образом:

! Описание областей Oi, Oi1, Oj

Oi:(i=1..N). Oi1:(i=25..N-25). Oj:(j=1..N).

! Область Oij есть декартово произведение (i=1..N)X(j=1..N)

Oij:(Oi;Oj).

! Область Oij1 есть декартово произведение (i=25..N-25)X(j=1..N)

Oij1:(Oi1;Oj).

! Определение параметра областей

DOMAIN PARAMETERS N=100.

! Описание переменных на областях

VARIABLE U,V,W DEFINED ON Oij.

! Число процессоров по направлению i равно 10.

DISTRIBUTION INDEX i=1..10, j=1.

! Расчетные формулы (индексы без смещений можно опускать)

FOR Oij ASSUME U=i+j; W=U-1.

FOR Oij1 ASSUME V=U[i-1]+1.

Фрагмент выходной программы, который будет построен компилятором с языка Норма при трансляции на язык Фортран-77 с использованием библиотеки MPI, представлен ниже.

C Вычисление параметров циклов в зависимости от номера процессора IN01 = 1 IK01 = 10 IN02 = 1 IK02 = 10 IF(IPRWWW.EQ.3) THEN IN02=5 IF(IPRWWW.EQ.8) THEN IK02=5 DO j=1,100 DO i=IN01,IK01 C Выражение i+j в процессоре с номером IPR U(i,j)=((IPRWWW-1)*10+i)+j ENDOO ENDDO C Требуются вычисленные значения из процессоров с номерами 3..7 IF(IPRWWW.LT.3.OR.IPRWWW.GT.7.OR.JPRWWW.NE.1) GOTO 4 DO j=1,100 RM1(j)=U(10,j) ENDDO CALL MPI_SEND(RM1,100,MPI_REAL, ITASKJ(JPRWWW+(IPRWWW+1-1)), >2,MPI_COMM_WORLD,IER1) 4 CONTINUE DO j=1,100 DO i=IN01,IK01 W(i,j)=U(i,j)-1 ENDOO ENDDO C Используются вычисленные значения в процессорах с номерами 4..8 IF(IPRWWW.LT.4.OR.IPRWWW.GT.8.OR.JPRWWW.NE.1) GOTO 7 CALL MPI_RECV(RM2,100,MPI_REAL,ITASKJ(JPRWWW+(IPRWWW-2)), >2,MPI_COMM_WORLD,STTS1,IER1) DO j=1,100 U(0,j)=RM2(j) ENDDO 7 CONTINUE IF(IPRWWW.LT.3.OR.IPRWWW.GT.8) GOTO 8 DO j=1,100 DO i=IN02,IK02 V(i,j)= U(i-1,j)+1 ENDOO ENDDO CONTINUE

Видно, что эта программа менее проста и понятна, чем программа на Норме, однако это скорее всего не должно беспокоить пользователя – эта программа построена автоматически компилятором. Поясним, используя текст этой программы, как решены задачи 1)-4). Идентификатор IPRWWW обозначает номер процессора из линейки процессоров и принимает значения от 1 до 10, идентификатор JPRWWW=1 (представление линейки процессоров).

Определение номеров процессоров. Пусть некоторое вычисление необходимо выполнить при значениях индекса i, принимающего значения в диапазоне от до . В языке Норма значения и известны на этапе трансляции, поэтому компилятор по значению () и числу точек, распределенных на процессор, определяет номер процессора () в котором находится точка (). Таким образом, диапазон процессоров, участвующих в этом вычислении, задается сегментом [, ]. Процессоры, номера которых находятся за пределами этого сегмента, не участвуют в реализации данного вычисления. В нашем примере операторы, реализующие соотношения 1-2, выполняются во всех 10 процессорах. Соотношение 3 вычисляется при значениях индекса i в диапазоне от 25 до 75. Точка i=25 принадлежит процессору с номером 3, а точка i=75 принадлежит процессору с номером 8. Поэтому оператор, реализующий соотношение 3, выполняется только в процессорах, номера которых лежат в диапазоне [3,8] – это условие реализуется последним условным оператором фрагмента программы.

Определение границ изменения параметров циклов. Аналогично тому, как определяется диапазон процессоров, реализующих заданное соотношение (вычисление), определяются начальные и конечные значения параметра цикла i, при которых надо выполнить оператор, реализующий соотношение. В нашем случае, операторы 1 и 2 выполняются (по индексу i) от IN01 до IK01, а оператор 3 - от IN02 до IK02. При этом, во всех 10 процессорах, значения IN01=1, IK01 = 10. Более сложный вид имеют значения IN02 и IK02. В процессоре с номером 3 значение IN02 = 5, во всех остальных IN02 = 1. В процессоре с номером 8 значение IK02 = 5, а в остальных IK02 = 10.

Описания массивов. Построение описаний для величин, определенных на областях, проводится достаточно простым образом. В последовательной программе переменным U,V,W соответствуют описания

DIMENSION U(100,100), V(100,100), W(100,100)

При построении распределенной программы учитывается, что эти переменные должны быть распределены в соответствии с конструкцией

DISTRIBUTION INDEX i=1..10, j=1

то есть "разрезаны" на 10 процессоров по частям, состоящим из 10´ 100 элементов. Кроме этого, поскольку в формуле для вычисления V требуются значения U[i-1], необходимо выделить память ("теневую" грань) для приема значений U[i-1] из соседнего левого процессора. Это делается за счет расширения влево диапазона при описании массива U. Таким образом, описания для распределенной программы переменных U,V,W выглядят так:

DIMENSION U(0:10,100), V(10,100), W(10,100)

Генерация операторов обмена. Интерфейс передачи сообщений MPI предоставляет разнообразные средства реализации взаимодействия в модели передачи сообщений. В рассматриваемой программе взаимодействие между процессорами реализовано при помощи подпрограмм MPI_SEND (отправить сообщение) и MPI_RECV (принять сообщение), которые имеют следующие параметры:

CALL MPI_SEND(<адрес посылки>,<длина посылки>,<тип посылаемых элементов>,

<имя процесса-получателя>,<тэг>,<коммуникатор >,<код возврата>)

CALL MPI_RECV (<адрес приема>,<длина приема>, <тип принимаемых элементов>,

<имя процесса-отправителя>,<тэг>, <коммуникатор>,<статус>,<код возврата>.

В других системах передачи сообщений (PVM, GNS и т.п.) формат подпрограмм обмена сообщениями имеют немного другой вид. Однако, во всех этих системах для процесса обмена определены параметры, указывающие процесс-отправитель (процесс-получатель), адрес посылаемого объекта, длину сообщения и способ идентификации сообщения (тэг).

Для обеспечения взаимодействия в каждой распределенной программе компилятор с языка Норма автоматически определяет и использует следующие структуры данных:

ITASKJ(n,m) - матрица имен параллельных процессов, участвующих в вычислении (ее размеры n и m фактически совпадают с параметрами распределения, заданными в конструкции DISTRIBUTION INDEX, поэтому в данном примере описание матрицы имен имеет вид INTEGER ITASKJ(10,1);

IPRWWW, JPRWWW - значения индексов матрицы ITASKJ, определяющие имя каждого из процессов (другими словами, каждый процесс "знает", что его имя - ITASKJ(IPRWWW,JPRWWW)).

Таким образом, вызов

CALL MPI_SEND (RM1,100,MPI_REAL, ITASKJ(JPRWWW+IPRWWW),

2,MPI_COMM_WORLD,IER1)

приводит к пересылке 100 элементов массива RM1 типа REAL соседу справа. Посылку осуществляет каждый из процессов с номерами JPRWWW=1, IPRWWW=3..7.

Вызов

CALL MPI_RECV (RM2,100,MPI_REAL,ITASKJ(JPRWWW+IPRWWW-2),

2,MPI_COMM_WORLD,STTS1,IER1)

приводит к приему 100 элементов типа REAL в массив RM2 от соседа слева. Прием осуществляет каждый из процессов с номерами JPRWWW=1, IPRWWW=4..8. Тэг со значением 2 позволяет “связать” прием и отправку сообщения (это важно , когда один и тот же процессор посылает одному и тому же процессору несколько сообщений – разные тэги позволяют однозначно идентифицировать эти сообщения).

Определение параметров операторов MPI_SEND и MPI_RECV проводится автоматически на этапе компиляции. Место в выходной программе, где должны стоять эти операторы, определяется также компилятором из следующих соображений. Оператор MPI_SEND ставится непосредственно за оператором, в котором вычисляется требуемое для передачи значение. Оператор MPI_RECV ставится непосредственно перед оператором, в котором используется принимаемое (требуемое) значение. Это позволяет надеяться на перекрытие времени счета и обмена данными. В нашем примере, сразу после вычисления значений переменной U стоит оператор записи значений величины U в рабочий массив RM1 и оператор отправки этого массива MPI_SEND. При этом, как отмечено выше, отправляют значения не все процессоры. Соответствующий оператор приема MPI_RECV расположен перед оператором, где используется значение U[i-1,j]. Соответствие обеспечивается наличием уникального тэга 2, причем генерацию тэгов также поддерживает компилятор. Отметим, что компилятор с языка Норма поддерживает несколько методов расстановки операторов обмена сообщениями.

В заключение несколько слов о схеме распараллеливания вычислений при трансляции Норма-программы.

Пусть распределенная вычислительная система представлена как матрица процессорных элементов. Как и при задании матриц можно считать, что каждый элемент такой вычислительной системы характеризуется своими координатами: номером строки (i) и номером столбца (j).

В программе на языке Норма программист сам определяет, точки какого индексного направления будут распределяться по процессорам системы, и задает число процессоров, на которых будет решаться исходная задача. Направление "разрезания" и число процессоров задается, например, описанием следующего вида: DISTRIBUTION INDEX i=1..6;j=1..5. Такая запись означает, что все области, в описании которых используются направления i или j, будут распределяться на 6 процессоров по направлению i и на 5 процессоров по направлению j.

Для эффективного распараллеливания при задании индексных направлений i, j в конструкции DISTRIBUTION INDEX обычно используются следующие два критерия:

1) максимальное число точек по индексному направлению (что приводит к распределению между процессорами наиболее "тяжелой работы"),

2) минимальное число вхождений индекса, заданного со смещением, в индексные выражения (это обычно минимизирует число обменов между процессорами, так как если точка с координатой i находится в процессоре с номером , то точка с координатой i-1, в зависимости от значения i, может находиться либо в процессоре, либо в процессоре .

Разрезание исходных данных проводится с учетом выбранных направлений. Обозначим верхнюю границу по направлению i для величины через . Пусть . Используемая процедура разрезания основывается на размещении постоянного числа компонентов всех массивов в каждом процессоре (по фиксированному направлению). Пусть по направлению i используется m процессоров, а по направлению j - n процессоров. Тогда число компонентов массива по направлению i (соответственно j) задается значением ().

Пусть координата процессора по направлению i есть iprwww, по направлению j - jprwww. Тогда в процессоре с номером (iprwww, jprwww) вычисляются значения по направлению i от (iprwww-1)ipoints+1 до iprwwwipoints, а по направлению j от (jprwww-1)jpoints+1 до jprwwwjpoints. В нашем примере процессоры образуют матрицу размерностью 6 X 5. Исходная матрица значений имеет размерность . Таким образом, каждый процессор проводит расчет величины qu размерностью .

3. Разработка программы: от Фортрана к Норме
1. Общая схема вычислений исходной программы

Текст исходной программы на Фортране приведен в Приложении 1. Вычисления в этой программе организованы по следующей схеме:

Присваивание начальных значений расчетным переменным Итерационный цикл по времени Вычисление временного шага для итерации по времени Итерационный цикл Рунге-Кутта (3 шага) Вычисление потока F Вычисление потока G Вычисление величины qu Конец итерации Рунге-Кутта Вычисление значений величины qu на границах Конец итерации по времени Запись вычисленных значений величины qu в файл

Определим общую структуру программы, а затем уточним способы вычисления необходимых переменных, а также способы описания и использования необходимых разделов (подпрограмм).

Программу для вычислительной системы с распределенной архитектурой представим в виде двух программ: главной программы и распределенной программы. Главная программа выполняется на одном процессоре системы и в основном используется для ввода/вывода данных. Основной расчет осуществляет распределенная программа, которая выполняется параллельно всеми остальными процессорами, затребованными для счета.

2. Главная программа

В исходной программе значения переменной qu многократно перевычисляются, так как Фортран допускает переприсваивание значений переменных. В программе на Норме в таких случаях будут вводиться дополнительные “копии” величин для исключения переприсваивания (в выходной программе компилятор оптимизирует использование требуемой из-за введения лишних переменных лишней памяти). Для определения начальных значений величины qu на первом этапе вводится величина qu0. Главная программа на Норме выглядит следующим образом:

MAIN PART GDST.

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

oijm:((i=1..Nx);(j=1..Ny);(m=1..4)). ok1:(k=1..Nx). ok2:(k=1..Ny).

VARIABLE QU0,QU DEFINED ON oijm.

VARIABLE x1 DEFINED ON ok1. VARIABLE x2 DEFINED ON ok2.

VARIABLE gamma,gam1,c0,dt0,an,fi,pi,dx REAL.

pi=3.1415926.

gamma=1.4.

gam1=gamma-1.0.

c0=0.5.

an=1.0/6.0.

fi=pi/3.0.

COMPUTE init(pi,gam1,fi,an

RESULT dx,x1 ON ok1, x2 ON ok2, QU0 ON oijm,dt0).

COMPUTE GD(pi,gamma,gam1,c0,fi,an,dx,dt0,x1 ON ok1,x2 ON ok2,QU0 ON oijm

RESULT QU ON oijm).

COMPUTE writer(x1 ON ok1, x2 ON ok2, QU ON oijm).

END PART.

В этой программе, кроме присваивания значений скалярным переменным, заданы вызовы трех вычислений - init, GD, и writer. Вычисление init является подпрограммой, написанной на Фортране и используемой в рамках интерфейса Нормы и Фортрана. Вычисление GD является разделом языка Норма и будет описано на Норме. Вычисление writer также является подпрограммой, написанной на Фортране. Рассмотрим эти три части программы.

Вызов раздела

COMPUTE init(pi,gam1,fi,an

RESULT dx,x1 ON ok1, x2 ON ok2, QU0 ON oijm,dt0).

с входными параметрами pi, gam1, fi приводит к вычислению значений величины qu0 во всех точках области oijm (i=1..Nx;j=1..Ny;m=1..4), значений векторов x1 и x2; значения шага по пространству dx и начального значения dt0. Так как расчет начальных данных будет выполняться в главной программе и будет выполняться на одном процессоре, то можно просто использовать соответствующий код последовательной Фортран-программы (см. Приложение 1), оформив его в виде подпрограммы init, который практически совпадает с начальным фрагментом текста исходной программы SHATL на Фортране.

SUBROUTINE init(pi,gam1,fi,an,dx,x1,x2,QU,dt0) INCLUDE 'TASK.PAR' real QU(Nx,Ny,4) real x1(Nx),x2(Ny) db=1.0 dx=db/(Ny-1) x1(1)=-2.*dx DO I=2,Nx x1(I)=x1(I-1)+dx ENDDO x2(1)=-2.*dx DO I=2,Ny x2(I)=x2(I-1)+dx ENDDO DO J=1,Ny DO I=1,Nx IF(x1(i).LT.AN+x2(J)/TAN(fi)) THEN QU(I,J,1)=8. QU(I,J,2)=8.25*SIN(PI/4.)*QU(I,J,1) QU(I,J,3)=-8.25*COS(PI/4.)*QU(I,J,1) QU(I,J,4)=116.5/gam1+(QU(I,J,2)**2+QU(I,J,3)**2)/2./QU(I,J,1) ELSE QU(I,J,1)=1.4 QU(I,J,2)=0. QU(I,J,3)=0. QU(I,J,4)=1./gam1 ENDIF ENDDO ENDDO dt0=0.0 RETURN END

В данной подпрограмме, а также в приводимых ниже подпрограммах, используется ссылка на файл 'TASK.PAR'. Этот файл содержит следующие описания:

PARAMETER(NX=388,NY=122,ipoints=65,jpoints=25) C ipoints=[Nx/iproc], jpoints=[NY/jproc] PARAMETER(iproc=6,jproc=5) INTEGER ITASKJ(iproc*jproc) COMMON /NORMAMESSPASS/ITASKJ, IPRWWW, JPRWWW, IVVVVJ,IWWWWJ

В описаниях задаются значения основных параметров задачи: Nx, Ny - размерности сеток по направлению i и по направлению j; iproc (jproc) - число процессоров по направлению i(j); ipoints(jpoints) - число точек сетки в одном процессоре по направлению i(j). Последняя строка содержит описание стандартного (для Норма-системы) COMMON-блока. В частности, iprwww указывает номер процессора по направлению i, а jprwww указывает номер процессора по направлению j, а массив ITASKJ задает все процессы, участвующие в параллельной работе. Неформально, этот файл обеспечивает интерфейс программ на Фортране, построенных компилятором с языка Норма, и программ на Фортране, написанных вручную.

Оператор

COMPUTE GD(pi,gamma,gam1,c0,fi,an,dx,dt0,x1 ON ok1,x2 ON ok2,QU0 ON oijm

RESULT QU ON oijm).

задает вызов раздела GD, который будет определен как распределенный (предназначенный для параллельного выполнения) и в котором будут описаны основные вычисления задачи. Исходными параметрами для раздела GD являются значения величины qu0, вектора x1 и x2 и ряд констант. Результатом вычисления раздела являются значения расчетной величины qu во всех точках области определения.

Оператор

COMPUTE writer(x1 ON ok1, x2 ON ok2, QU ON oijm).

задает вызов раздела writer для записи вычисленных значений величины qu в файл. Раздел writer, как и раздел init, практически совпадает с заключительным фрагментом текста исходной программы SHATL на Фортране и оформлен в виде подпрограммы.

SUBROUTINE writer(x1,x2,QU) INCLUDE 'TASK.PAR' real QU(Nx,Ny,4) real x1(Nx) real x2(Ny) open(40,file=’REZ’) WRITE(40,*) 'VARIABLES="X","Y","D","U","V"' WRITE(40,*) 'ZONE T="XX",I=118,J=384,F=POINT' DO I=3,NX-2 DO J=3,NY-2 TT=QU(I,J,1) UU=QU(I,J,2)/TT VV=QU(I,J,3)/TT WRITE(40,*) X1(I),X2(J),TT,UU,VV ENDDO ENDDO close(40) RETURN END

3. Распределенная программа

Теперь перейдем к определению распределенного раздела GD. Величина qu вычисляется в результате итерационного процесса (цикла) по времени itime, причем на каждом шаге итерации вычисление этой величины разбивается на 2 этапа. На первом этапе вычисление проводится с помощью 3-х шагов итерации по методу Рунге-Кутта (по переменной irk). На втором этапе перевычисляются значения на нижней и верхней границах исходной области.

Область вычисления величины qu схематично (не изображен индекс m) представлена на Рис.1.

Величина qu остается неизменной в областях (i=1..2);(j=1..Ny);(m=1..4) и (i=Nx-1..Nx);(j=1..Ny);(m=1..4) - эти области в программе на языке Норма обозначены соответственно og1 и og2. Поэтому значения величины qu в этих областях необходимо определять в конструкции граничных значений (BOUNDARY...ENDBOUNDARY) итерации по itime. В областях og3: ((i=3..Nx-3);(j=1..2);(m=1..4)) и og4: ((i=3..Nx-3);(j=Ny-1..Ny);(m=1..4)) величина qu перевычисляется на каждом шаге итерации по времени itime. При этом на каждой области вычисления проводятся по различным формулам и реализуются в разделах downbound и upbound.

	Ny					og4
	Ny-1
		og1				oin			og2
J ­
	2					og3
	1
		1	2						Nx-1	Nx
					I ®

Рис.1 Схема расчетных областей задачи

В области oin: ((i=3..Nx-3);(j=3..Ny-2);(m=1..4)) вычисление величины qu проводиться с помощью 3-х шагов итерации методом Рунге-Кутта. В остальных точках области при проведении итерации методом Рунге-Кутта значения величины qu остаются неизменными, поэтому в программе на Норме эти значения также задаются в конструкции BOUNDARY...ENDBOUNDARY итерации по irk. Схема расчетных областей с точки зрения организации итерационных вычислений по методу Рунге-Кутта приведены на Рис.2.

	Ny					og4
	Ny-1	Граница для итерации Рунге-Кутта
				Область итеративных вычислений Рунге-Кутта
		og1				oin				og2
J ­
	2					og3
	1
		1	2						Nx-1	Nx
					I ®

Рис.2 Схема расчетных областей для итерации Рунге-Кутта

Программа распределенного раздела GD приведена ниже. Строки, начинающиеся с символа "!", в программе на языке Норма являются комментарием.

PART GD.

! входные параметры для распределенной программы

pi,gamma,gam1,c0,fi,an,dx,dt0,x1,x2,QU0

! результаты вычисления раздела

RESULT QU

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

DISTRIBUTION INDEX i=1..6,j=1..5.

oi:(i =1..Nx). oi3nm2:(i=3..Nx-2). oi2nm2:(i=2..Nx-2).

oj:(j =1..Ny). oj3nm2:(j=3..Ny-2). oj2nm2: (j =2..Ny-2).

m4:(m =1..4). m5: (m =1..5).

! i=1,Nx j=1,Ny,m=1,4

oijm:(oi;oj;m4).

!i=3,Nx-2 j=3,Ny-2,m=1,4 i=3,Nx-2 j=3,Ny-2

oin:(oi3nm2;oj3nm2;m4). oijstep:(oi3nm2;oj3nm2).

!i=3,Nx-2 j=1,Ny,m=1,4

oin1:(oi3nm2;oj;m4).

!i=1,2 j=1,Ny,m=1,4 i=Nx-1,Nx j=1,Ny; m=1,4

og1:((i=1..2);oj;m4). og2:((i=Nx-1..Nx);oj;m4).

!i=3,Nx-2 j=1,2;m=1,4 i=3,Nx-2, j=Ny-1,Ny; m=1,4

og3:(oi3nm2;(j=1..2);m4). og4:(oi3nm2;(j=Ny-1..Ny);m4).

og33:(oi3nm2;(j=1..2)). og44:(oi3nm2;(j=Ny-1..Ny)).

!i=3,Nx-2 j=1,2;m=1,4 i=3,Nx-2, j=Ny-1,Ny; m=1,4

og5:(oi3nm2;(j=1..2);m4). og6:(oi3nm2;(j=Ny-1..Ny);m4).

!i=1,Nx j=1,Ny i=1,Nx j=1,Ny m=1,5

oprx:(oi;oj). oprw:(oprx;m5).

!i=2,Nx-2 j=3,Ny-2 m=1,4 i=3,Nx-2 j=2,Ny-2 m=1,4

oif:(oi2nm2;oj3nm2;m4). oig:(oi3nm2;oj2nm2;m4).

ok1:(k=1..Nx). ok2:(k=1..Ny).

VARIABLE QU0,QU,QU2 DEFINED ON oijm.

VARIABLE umin DEFINED ON oijstep.

VARIABLE x1 DEFINED ON ok1.

VARIABLE x2 DEFINED ON ok2.

VARIABLE f DEFINED ON oif.

VARIABLE g DEFINED ON oig.

VARIABLE prx DEFINED ON oprx.

VARIABLE prw DEFINED ON oprw.

VARIABLE pi,gam1,gamma,c0,dt,dt0,dt0it,dk,dx,dtr,an,fi REAL.

VARIABLE A1,A2,A9,B9,C9,D9 REAL.

A9=0.25. B9=0.75. C9= 0.33333333. D9=0.66666666. dk=0.15.

ITERATION QU,dt0it ON itime.

BOUNDARY

!i=1..2; j=1..Ny;m=1..4 & i=Nx-1..Nx; j=1..Ny;m=1..4

FOR og1;og2 ASSUME QU=QU0.

END BOUNDARY

INITIAL itime=0:

dt0it=dt0.

! i=3..Nx-2; j=1..Ny;m=1..4

FOR oin1 ASSUME QU=QU0.

END INITIAL

!i=3..Nx-2; j=3..Ny-2

FOR oijstep ASSUME

COMPUTE step1(gam1,gamma,c0,dx,QU[itime-1] ON oijm/(i=i,j=j) RESULT umin).

dt =MIN((oijstep)umin).

dt0it=dt0it[itime-1]+dt.

OUTPUT dt(FILE='tst').

dtr=dt/dx.

ITERATION QU2 ON irk.

BOUNDARY

! i=1..2; j=1..Ny;m=1..4 & i=Nx-1..Nx; j=1..Ny;m=1..4

! i=3..Nx-2; j=1..2;m=1..4 & i=3..Nx-2; j=Ny-1..Ny;m=1..4

FOR og1;og2;og3;og4 ASSUME QU2=QU[itime-1].

END BOUNDARY

INITIAL irk=0:

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME QU2=QU[itime-1].

END INITIAL

! i=1..Nx;j=1..Ny

FOR oprx ASSUME prx=gam1*(QU2[irk-1,m=4]

-(QU2[irk-1,m=2]**2+QU2[irk-1,m=3]**2)/2.0/QU2[irk-1,m=1]).

! i=1,Nx;j=1,Ny; m=1,3

FOR oprw/m=1..3 ASSUME prw=QU2[irk-1].

! i=1,Nx;j=1,Ny; m=4

FOR oprw/m=4 ASSUME prw=QU2[irk-1]+prx.

! i=1,Nx;j=1,Ny; m=5

FOR oprw/m=5 ASSUME prw=gamma*prx/QU2[irk-1,m=1].

! i=2..Nx-2;j=3..Ny-2;m=1..4

COMPUTE mhdx (gam1,prw ON oprw, prx ON oprx, QU2[irk-1] ON oijm

RESULT f ON oif).

! i=3..Nx-2;j=2..Ny-2;m=1..4

COMPUTE mhdy (gam1,prw ON oprw, prx ON oprx, QU2[irk-1] ON oijm

RESULT g ON oig).

COMPUTE coeff (irk, A9,B9,C9,D9 RESULT A1,A2).

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME

QU2=A1*QU[itime-1]+A2*(QU2[irk-1]-dtr*(f-f[i-1]+g-g[j-1])).

EXIT WHEN (irk=3).

END ITERATION irk.

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME

COMPUTE absqu (QU2 ON oin/(i=i,j=j),i,j,an,x1 ON ok1

RESULT QU ON oin/(i=i,j=j)).

! i=3..Nx-2; j=1..2;m=1..4

FOR og33 ASSUME

COMPUTE downbound (QU ON og3/(i=i,j=3),i,x1 ON ok1,an,pi,gam1

RESULT QU ON og3/(i=i,j=j)).

! i=3..Nx-2; j= Ny-1,Ny;m=1..4

FOR og44 ASSUME

COMPUTE upbound (i,j,x1 ON ok1, x2 ON ok2,dt0it,fi,an,pi,gam1

RESULT QU ON og4/(i=i,j=j)).

EXIT WHEN (dt0it[itime]> dk).

END ITERATION itime.

END PART.

Все величины, определенные на областях, в которых индексы совпадают с индексами из описания DISTRIBUTION INDEX, распределяются равномерно по процессорам. В нашей программе величины x1 и x2 мы предпочитаем не распределять по процессорам. Таким образом, в каждом процессоре будут находиться свои копии этих векторов.

Величины prx и prw в исходной Фортан-программе определялись по одним и тем же формулам, как при вычислении потока f, так и при вычислении потока g. Разница состояла только в том, что вычисление проводилось при различных границах операторов цикла. В первом случае диапазон имел вид i=1,...,Nx; j=3,...,Ny-2 а во втором случае - i=3,...,Nx-2; j=1,Ny. В нашей программе мы вычисляем их один раз на области i=1,Nx;j=1,Ny.

Рассмотрим некоторые части этой программы более подробно. В начале программы приведены описания необходимых областей (эти описания определяются областями значений, которые принимают индексы в исходной программе на Фортране) и описания расчетных величин. Далее описываются необходимые вычисления.

1. Вычисление временного шага dt для итерации по времени

Значение временного шага dt зависит от значений величины qu в точках с координатами по оси i от 3 до Nx-2; а по оси j - от 3 до Ny-2.

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

В нашем случае вычисление величины umin выполняется в подпрограмме, в которой вычисляется один элемент величины umin. Вызов подпрограммы производится в операторе ASSUME на области oijstep:

FOR oijstep ASSUME

COMPUTE step1(gam1,gamma,c0,dx,QU[itime-1] ON oijm/(i=i,j=j) (*)

RESULT umin).

Выражение qu[itime-1] ON oijm/(i=i,j=j) означает, что в качестве параметра при вызове раздела step1, при фиксированных значениях параметров цикла - i=i0,j=j0, передается рабочий вектор R, состоящий из 4-х компонентов. При этом, значения соответствующих компонентов вектора R вычисляются по следующему правилу: R(1)=qu(i0,j0,1); R(2)=qu(i0,j0,2); R(3)=qu(i0,j0,3); R(4)=qu(i0,j0,4).

Возможен и другой подход к организации вычислений значений величины umin. При этом подходе в подпрограмме вычисляются все значения величины umin во всех точках области oijstep. Такое вычисление задается оператором:

COMPUTE step2(gam1,gamma,c0,dx,QU ON oijm RESULT umin ON oijstep). (**)

Ниже в таблице представлены описания двух подпрограмм, реализующих рассматриваемое вычисление величины umin.

SUBROUTINE step1(gam1,gamma,c0,dx,QU,umin) real QU(4) DTM=1.E+10 PRR=GAM1*(QU(4)-(QU(2)**2+ QU(3)**2)/2./QU(1)) CL=SQRT(GAMMA*PRR/QU(1)) UCSX=QU(2)/QU(1) UCSY=QU(3)/QU(1) DTMX=1./(ABS(UCSX)+CL) DTMY=1./(ABS(UCSY)+CL) umin=C0*DX*AMIN1(DTM,DTMX,DTMY) RETURN END

SUBROUTINE step2(gam1,gamma,c0,dx,QU,umin) INCLUDE 'TASK.PAR' real QU(ipoints,jpoints,4) real umin(ipoints,jpoints) in=1 ik=ipoints jn=1 jk=jpoints if(iprwww.eq.1)in=3 if(iprwww.eq.iproc)ik=NX-(iproc-1)*ipoints if(jprwww.eq.1)jn=3 if(jprwww.eq.jproc)jk=NY-(jproc-1)*jpoints DTM=1.E+10 do I=in,ik do J=jn,jk PRR=GAM1*(QU(I,J,4)-(QU(I,J,2)**2+ QU(I,J,3)**2)/2./QU(I,J,1)) CL=SQRT(GAMMA*PRR/QU(I,J,1)) UCSX=QU(I,J,2)/QU(I,J,1) UCSY=QU(I,J,3)/QU(I,J,1) DTMX=1./(ABS(UCSX)+CL) DTMY=1./(ABS(UCSY)+CL) umin(i,j)=C0*DX*AMIN1(DTM,DTMX,DTMY) enddo enddo RETURN END

Разница между подпрограммами step1 и step2 определяется следующими факторами.

Подпрограмма step1 вычисляет одно скалярное значение umin, не зависит от значений переменных по индексам распределения i и j, всю организацию распределенных вычислений берет на себя компилятор при трансляции оператора ASSUME (*).

Подпрограмма step2 имеет более сложный вид, так как она вычисляет значения umin на области, зависящей от индексов распределения i и j, зависит от значений переменной qu с индексами распределения i и j, а организация распределенных вычислений, сделанная компилятором при трансляции оператора (**) и отраженная в файле 'TASK.PAR', должна быть учтена при написании тела подпрограммы. Действительно, эта подпрограмма выполняется на каждом процессоре и поэтому границы используемых массивов по индексам распределения i и j зависят от размерности используемой сетки и числа процессоров, используемых для решения задачи (эти параметры содержатся в файле 'TASK.PAR'). Кроме того, существенным является вопрос определения границ для циклов, в которых вычисляются значения величины umin. В исходной программе диапазоны изменения параметров циклов имеет вид: i=3..Nx-2, j=3..Ny-2. При распределении на процессор по ipoints´ jpoints точек диапазоны изменения параметров циклов имеют вид i=1..ipoints, j= 1..jpoints для всех процессоров, кроме 4-х угловых (с координатами (1,1),(iproc,1),(1,jproc),(iproc,jproc) в матрице процессоров). Первые операторы подпрограммы step2 как раз и вычисляют эти значения – здесь фактически вручную программируются действия, которые для подпрограммы step1 будет сделана компилятором.

Два варианта организации вычислений - подпрограмма step1 и подпрограмма step2 – приведены для того, чтобы показать два способа использования интерфейса Нормы и Фортрана. Хотя второй способ сложнее, он также может оказаться полезным. Для разрабатываемой нами программы на Норме мы в дальнейшем будем использовать первый вариант вызова из распределенного раздела на Норме подпрограммы на Фортране.

Отметим также, что вызов подпрограмм на Фортране из программы на Норме используется для того, чтобы показать, как можно использовать уже существующий код на Фортране. Вычисления, задаваемые в подпрограмме на Фортране, могут быть также описаны и на Норме – ниже мы приведем примеры вызова из распределенного раздела на Норме другого распределенного раздела, описанного на Норме.

Вторая фаза расчета временного шага dt - вычисление минимального значения в массиве umin по всем точкам области oijstep - производится с помощью вызова стандартной функции MIN нахождения минимума по области, определенной в языке Норма:

dt =MIN((oijstep)umin).

Здесь для пользователя никаких проблем больше нет – реализует это вычисление компилятор.

Отметим в заключение, что в программе кроме величины временного шага dt используется значение текущего времени dt0it, которое перевычисляется на каждом шаге итерации по itime. Это вычисление описано оператором

dt0it=dt0it[itime-1]+dt.

в котором у итерируемой величины dt0it в правой части задано индексное выражение [itime-1], что означает значение dt0it с предыдущего итерационного шага по itime.

Значение dt, вычисляемое на каждом шаге итерации по itime, выводится в файл 'tst' при помощи описания

OUTPUT dt(FILE='tst').

2. Вычисление потоков F и G

Вычисление потоков F и G осуществляется в результате вложенной итерации по методу Рунге-Кутта, выполняемой по индексу итерации irk. В процессе вычисления потоков перевычисляется величина qu, поэтому в программе на Норме вводится дополнительная переменная-“копия” qu2, чтобы избежать переприсваивания. В программе на Норме величины для представления потоков обозначены переменными f и g.

Вычисление потока F оформим в виде раздела mhdx, написанного на языке Норма, а потока G - в виде раздела mhdy, также написанного на языке Норма.

Сначала рассмотрим, какие величины являются исходными для вычисления раздела mhdx, а какие величины являются результатом вычисления раздела.

Исходными данными для вычисления потока F являются значения величины qu2 с предыдущего итерационного шага, prx и prw, вычисление которых описано сразу за блоком INITIAL...ENDINITIAL итерации по irk, и константы gam1. Результатом выполнения раздела mhdx являются значения величины f в точках области oif: ((i=2,Nx-2);(j=3,Ny-2)). Запрос на вычисление величины f в нашем случае записан в следующем виде:

COMPUTE mhdx (gam1,prw ON oprw,prx ON oprx,QU2[irk-1] ON oijm

RESULT f ON oif).

Теперь рассмотрим подробнее процесс вычисления потока f, заданный в исходной программе на Фортране. Рассмотрим только операторы циклов, обращая особое внимание на границы изменения индексов циклов, и выделим те переменные, которые автор программы использует с индексами. Схематично определяется следующая последовательность вычислений.

J=3,Ny-2

I=1,Nx; PRW=

I=1,Nx-1 QSR=

I=1,Nx; EQ=

I=1,Nx-1 SLP=

I=1,Nx-2; GGF=

I=2,Nx-2 AS, GMX, XFI, ED, F=

Вся программа представляет цикл с параметром j, изменяющимся от 3 до Ny-2. Далее идет последовательность вычислений ряда вспомогательных величин. При этом каждая величина вычисляется в цикле по i со своими границами. Автор Фортран-программы с целью экономии памяти описывает используемые переменные как вектора (конечно, за исключением величины f, которая является двумерной). Фактически же все величины являются двумерными - это определяется использованием внешнего цикла по направлению j и вложенных циклов по направлению i.

Конечно, у этих величин есть еще размерность по m от 1 до 5 (или от 1 до 4), но эта размерность используется только для обозначения различных физических величин, которые задаются в одной переменной.

Введем описания областей для определения этих вспомогательных величин. Сначала опишем одномерные области:

oi:(i=1..Nx). oi1:(i=1..Nx-1). oi2nm2:(i=2..Nx-2). oi2:(i=1..Nx-2).

oj:(j=1..Ny). oj3nm2:(j=3..Ny-2).

m5:(m=1..5). m4:(m=1..4).

из которых построим необходимые многомерные:

! i=1,Nx;j=1,Ny;m=1,4 i=2,Nx-2;j=3,Ny-2;m=1,4

oijm: (oi;oj;m4). of: (oi2nm2;oj3nm2;m4).

!i=1,Nx;j=1,Ny i=1,Nx;j=1,Ny;m=1,5

oprx:(oi;oj). oprw:(oprx;m5).

!i=1,Nx;j=3,Ny-2;m=1,4 i=1,Nx-1;j=3,Ny-2;m=1,5

oeq: (oi;oj3nm2;m4). oslp: (oi1;oj3nm2;m5).

oslp1: (oi1;oj3nm2).

oslp2: (oi1;oj3nm2;m4).

!i=1,Nx-2;j=3,Ny-2;m=1,4; i=1,Nx-1;j=3;Ny-2

oggf: (oi2;oj3nm2;m4). orsl: (oi1;oj3nm2).

и две условные области oft и off:

oft,off:of/slp<>0.

Опишем на введенных областях необходимые величины:

VARIABLE QU DEFINED ON oijm.

VARIABLE eq DEFINED ON oeq.

VARIABLE prw DEFINED ON oprw.

VARIABLE prx DEFINED ON oprx.

VARIABLE qsr,rsr,rsl,alf,alg DEFINED ON oslp.

VARIABLE slp,xxx DEFINED ON oslp2.

VARIABLE csr DEFINED ON oslp1.

VARIABLE rsr,rsl,alf,alg DEFINED ON oslp.

VARIABLE ggf DEFINED ON oggf.

VARIABLE as,gmx,xfi,ed,f DEFINED ON of.

VARIABLE gam1 REAL.

VARIABLE sl, stx DEFINED ON oggf.

Неформально схему вычисления можно представить в виде следующей последовательности операторов :

FOR oeq ASSUME eq= .

FOR oslp ASSUME rsl= ;

rsr= ;

alf= ;

alg= .

FOR oslp ASSUME qsr= .

FOR oslp2 ASSUME slp= .

FOR oggf ASSUME sl= ;

stx= ;

ggf= .

FOR oft ASSUME gmx= .

FOR of ASSUME as= ;

xfi= ;

ed= ;

f= .

Отметим, что порядок записи операторов ASSUME несущественен, так как транслятор, по зависимостям переменных, определяет правильную последовательность вычислений.

Ниже вид некоторых из этих операторов конкретизируется, причем приводится как фрагмент текста на Фортране их исходной программы, так и соответствующий ему фрагмент на Норме. После этого приводится полный текст программы раздела mhdx на Норме.

1. Вычисление величины eq

В приведенной ниже таблице, в левом столбце приведен фрагмент исходной программы. В правом столбце приведены операторы языка Норма, реализующие те же вычисления. В операторах языка Норма мы отказались от использования вспомогательных переменных RR0 и RRA.

DO I=1,NX RR0=QU(I,J,1) RRA=QU(I,J,2) EQ(I,1)=RRA EQ(I,2)=RRA**2/RR0+PRX(I) EQ(I,3)=RRA*QU(I,J,3)/RR0 EQ(I,4)=(QU(I,J,4)+PRX(I))*RRA/RR0 ENDDO

FOR oeq/m=1 ASSUME eq=QU[m=2]. FOR oeq/m=2 ASSUME eq=QU[m=2]**2/QU[m=1]+prx. FOR oeq/m=3 ASSUME eq=QU[m=2]*QU/QU[m=1]. FOR oeq/m=4 ASSUME eq=(QU+prx)*QU[m=2]/QU[m=1].

Конечно, как и прежде, можно оформить вычисление величины eq в виде подпрограммы и и описать это вычисление оператором COMPUTE:

FOR oeq ASSUME

COMPUTE eq1(QU ON oijm/(i=i,j=j),prx RESULT eq/(i=i,j=j)).

или

COMPUTE eq2(QU ON oijm,prx ON oprx RESULT eq ON oeq).

В первом случае при каждом вызове подпрограммы eq1 вычисляется значение величины eq в одной точке области. Во втором случае при вызове подпрограммы eq2 значение величины eq вычисляется во всех точках области oeq. Ниже (в качестве возможных вариантов) приведены тексты этих подпрограмм (см. также обсуждение вариантов интерфейса с Фортраном из п.5.3.1).

SUBROUTINE eq1(QU,prx,eq) real QU(4) real eq(4) real prx RR0=QU(1) RRA=QU(2) EQ(1)=RRA EQ(2)=RRA**2/RR0+PRX EQ(3)=RRA*QU(3)/RR0 EQ(4)=(QU(4)+PRX)*RRA/RR0 RETURN END

SUBROUTINE eq2(QU,prx,eq) include 'task.par' real QU(ipoints,jpoints,4) real eq(ipoints,jpoints,4) real prx(ipoints,jpoints) in=1 ik=ipoints jn=1 jk=jpoints if(iprwww.eq.1)ik=NX-(iprwww-1)*ipoints if(jprwww.eq.iproc)jn=3 if(jprwww.eq.jproc)jk=NY-(jprwww-1)*jpoints do I=in,ik do J=jn,jk RR0=QU(I,J,1) RRA=QU(I,J,2) EQ(I,J,1)=RRA EQ(I,J,2)=RRA**2/RR0+PRX(I,J) EQ(I,J,3)=RRA*QU(I,J,3)/RR0 EQ(I,J,4)=(QU(I,J,4)+PRX(I,J))*RRA/RR0 enddo enddo RETURN END

 

В программе eq2 циклы выполняются от in до ik по направлению i и от jn до jk по направлению j. Конечно, значения этих параметров являются различными в каждом процессоре. В исходной программе цикл по направлению i выполняется от 1 до 388. В процессоре с номером iprwww=6 искомые значения вычисляются по направлению i в диапазоне от (iprwww-1) ipoints+1 до iprwwwipoints или от 565+1=326 до 665 =390. Но нам надо проводить счет до i=388, поэтому для процессора, у которого iprwww =равен 6, конечное значение цикла по i (ik) должно быть равно 63 или ipoints-2. Аналогично должны быть определены значения величин jn и jk для процессоров у которых jprwww=1 и jprwww=5.

Замечание 1. Автоматическое копирование операторов исходной программы невозможно в обоих случаях. В программе eq1 параметр qu является вектором, состоящим из 4-х компонентов. Поэтому необходимо в исходном тексте в каждом вхождении величины qu удалить индексные переменные i и j. Аналогичное замечание справедливо и для величины eq.

Замечание 2. Вопрос о том, надо ли описывать некоторое вычисление в виде операторов ASSUME, или в виде раздела, зависит от принятой технологии разработки программы и аналогичен вопросу о том, когда в программе следует использовать подпрограмму. В данном материале мы эти вопросы не рассматриваем и иллюстрируем лишь различные способы организации программы на Норме.

В итоге, отдадим предпочтение первому способу описания вычисления величины eq – при помощи операторов ASSUME.

2. Вычисление величины qsr

При вычислении величины qsr в исходной Фортран программе используется подпрограмма CSXK, которая при вычислении потока F вызывается с параметром 1: CALL CSXK(1), а при вычислении потока G - с параметром 2: CALL CSXK(2). При этом параметр-индикатор, равный 1 или 2, фактически определяет, по какому направлению индексного пространства – i или j – надо проводить вычисление, а требуемые для расчетов данные – массив prw – для каждого из потоков предварительно вычисляется и передается через COMMON-блок.

Такой способ организации вычислений является примером введения избыточных информационных зависимостей и ограничений при разработке последовательной программы, которые, по-видимому, связаны с экономией памяти (для программы и для данных) и которые требуют дополнительного сложного анализа состояния памяти данных при попытке автоматического распараллеливания такой программы. Такой анализ не всегда возможен – определить, есть зависимость по данным или нет, в общем случае часто не удается. Этот пример иллюстрирует одну из проблем преобразования последовательной программы в параллельную – критерии и методы разработки последовательной программы и параллельной программы различаются, и часто конфликтуют.

При трансляции программы на языке Норма проблема анализа состояния памяти не возникает – переприсваивание по определению языка невозможно. И уже в условиях разрешимости задачи распараллеливания могут решаться вопросы экономии памяти (автоматически компилятором или пользователем).

При вычислении величины qsr используются вспомогательные величины rsl и rsr. В исходной программе эти величины описаны как скаляры. В нашей программе мы определим их как величины, определенные на области oslp: (i=1,Nx-1;j=3,Ny-2). Фрагменты соответствующих вычислений на Фортране (для случая вызова CSXK(1), то есть NK=1, NN=Nx) и Норме приведены ниже:

 

DO I=1,Nx-1 RR0=PRW(I,1) RR1=PRW(I+1,1) RSL=SQRT(RR0) RSR=SQRT(RR1) ALF=RSL/(RSL+RSR) ALG=1.-ALF QSR(I,1)=RSL*RSR QSR(I,2)=ALF*PRW(I,2)/RR0+ALG*PRW(I+1,2)/RR1 QSR(I,3)=ALF*PRW(I,3)/RR0+ALG*PRW(I+1,3)/RR1 QSR(I,4)=ALF*PRW(I,4)/RR0+ALG*PRW(I+1,4)/RR1 QSR(I,5)=ALF*PRW(I,5)+ALG*PRW(I+1,5)+GAM1/2. 1 *QSR(I,1)/(RSL+RSR)**2* 1 ((PRW(I+1,2)/RR0-PRW(I,2)/RR0)**2 1 +(PRW(I+1,3)/RR0-PRW(I,3)/RR0)**2) ENDDO

FOR oslp ASSUME rsl=SQRT(prw[m=1]); rsr = SQRT(prw[i+1,m=1]); alf=rsl/(rsl+rsr); alg=1.0-alf. FOR oslp/m=1 ASSUME qsr=rsl*rsr. FOR oslp/m=2..4 ASSUME qsr=alf*prw/prw[m=1] +alg*prw[i+1]/prw[i+1,m=1]. FOR oslp/m=5 ASSUME qsr=alf*prw+alg*prw[i+1] +gam1/2.0*qsr[m=1]/(rsl+rsr)**2* ((prw[i+1,m=2]/prw[m=1]-prw[m=2]/prw[m=1])**2 +(prw[i+1,m=3]/prw[m=1] -prw[m=3]/prw[m=1])**2).

При вычисление величины qsr, как и выше для величины eq, мы используем только операторы языка Норма, а не подпрограмму, к которой производится обращение из Норма-программы.

Отметим, что при расчете величин qsr и rsr используются значения величины prw, заданные с индексным смещением i+1. При распределении величин по процессорам в качестве направлений, по которым проводится “разрезание”, в конструкции DISTRIBUTED INDEX выбраны направления i и j. Поэтому использование в правой части оператора ASSUME величины prw c индексом i+1 означает, что для вычисления требуются значения величины prw, находящиеся в соседнем справа (по направлению i) процессоре и они должны быть из этого процессора переданы . Компилятор с языка Норма автоматически определяет необходимость таких обменов и строит в выходной программе соответствующие операторы.

Конструкция DISTRIBUTED INDEX, заданная в разделе mhdx, с точки зрения пользователя эквивалентна конструкции DISTRIBUTION INDEX. При помощи конструкции DISTRIBUTED INDEX должны быть заданы индексы распределения в распределенных разделах, вызываемых из распределенных разделов – эта информация используется компилятором для необходимой организации параллельных вычислений.

3. Вычисление величин slp и csr

Вычисление slp и csr оформим в виде подпрограммы на Фортране. Как обычно, определим, какие величины необходимы для вычисления (in-параметры), и для каких значений индексов (на каких областях) они требуются.

Цикл из исходной программы приведен в левом столбце таблицы. Тест подпрограммы cslx приведен в правом столбце. В подпрограмме вычисляются значения величин slp и csr при фиксированных значениях индексных переменных (i,j). В исходной программе вычисление величин csr и slp проводится в цикле по i от 1 до Nx-1 и по j от 1 до Ny-2. Поэтому обращение к подпрограмме cslx выполняется с помощью оператора COMPUTE, помещенного в оператор ASSUME:

FOR oslp1 ASSUME

COMPUTE cslx(gam1, QU ON oijm/(i=i+1,j=j),

QU ON oijm/(i=i,j=j), qsr ON oslp/(i=i,j=j)

RESULT slp ON oslp2/(i=i,j=j),csr).

DO J=1,NY-2 ... DO I=1,NX-1 QR1=QU(I+1,J,1)-QU(I,J,1) QR2=QU(I+1,J,2)-QU(I,J,2) QR3=QU(I+1,J,3)-QU(I,J,3) QR4=QU(I+1,J,4)-QU(I,J,4) CSR(I)=SQRT(QSR(I,5)) TMX=QSR(I,2)/CSR(I) TMY=QSR(I,3)/CSR(I) TM2=GAM1*(TMX**2+TMY**2)/2. TMC=GAM1*TMX/CSR(I) TMD=GAM1*TMY/CSR(I) GMC=GAM1/CSR(I)**2 CS2=1./CSR(I) SLP(I,1)=((TMX+TM2)*QR1-(CS2+TMC)*QR2- > TMD*QR3+GMC*QR4)/2. SLP(I,2)=((1.-TMY- > TM2)*QR1+TMC*QR2+(CS2+TMD)*QR3-GMC*QR4)/2. SLP(I,3)=((1.+TMY-TM2)*QR1+TMC*QR2-(CS2-> TMD)*QR3-GMC*QR4)/2. SLP(I,4)=((-TMX+TM2)*QR1+(CS2-TMC)*QR2->TMD*QR3+GMC*QR4)/2. ENDDO ... ENDDO

SUBROUTINE cslx(gam1,QU1,QU2,qsr,slp,csr) real QU1(4),QU2(4),qsr(5),slp(4),csr QR1=QU(1)-QU2(1) QR2=QU1(2)-QU2(2) QR3=QU1(3)-QU2(3) QR4=QU1(4)-QU2(4) csr=SQRT(qsr(5)) TMX=QSR(2)/csr TMY=QSR(3)/csr TM2=gam1*(TMX**2+TMY**2)/2. TMC=gam1*TMX/csr TMD=gam1*TMY/csr GMC=gam1/csr**2 CS2=1./csr slp(1)=((TMX+TM2)*QR1-(CS2+TMC)*QR2- > TMD*QR3+GMC*QR4)/2. slp(2)=((1.-TMY-TM2)*QR1 > +TMC*QR2+(CS2+TMD)*QR3-GMC*QR4)/2. slp(3)=((1.+TMY-TM2)*QR1+TMC*QR2- > (CS2-TMD)*QR3-GMC*QR4)/2. slp(4)=((-TMX+TM2)*QR1+(CS2-TMC)*QR2- > TMD*QR3+GMC*QR4)/2. RETURN END

4. Вычисление величины ggf

При вычислении величины ggf в программе на Норме вместо переменных SL1 и SL2 мы подставили в формулу их выражения для сокращения числа используемых промежуточных величин. Вычисление описывается одним оператором ASSUME, а при записи соотношений использованы стандартные функции, доступные как в Фортране, так и в Норме:

DO I=1,NX-2 DO M=1,4 SL=2.*SLP(I,M) SL1=2.*SLP(I+1,M) SL2=(SLP(I,M)+SLP(I+1,M))/2. STX=SIGN(1.,SL) GGF(I,M)=STX*AMAX1(0., AMIN1(ABS(SL),STX*SL1,STX*SL2)) ENDDO ENDDO

FOR oggf ASSUME sl=2.0*slp;   stx=SIGN(1,sl); ggf=stx*AMAX1(0.0, AMIN1(ABS(sl),stx*2.0*slp[i+1], stx*(slp+slp[i+1])/2.0)).

5. Вычисление величины as, gmx, ed, f

Оставшиеся вычисления, необходимые для определения потока f, приведем без комментариев – они соответствуют операторам исходной Фортран-программы, приведенным в левом столбце

DO I=2,NX-2 AS(1)=QSR(I,2)-CSR(I) AS(2)=QSR(I,2) AS(3)=QSR(I,2) AS(4)=QSR(I,2)+CSR(I) DO M=1,4 IF(SLP(I,M).NE.0.) THEN MH=M GMX(M)=FF(AS(M)) >*(GGF(I,M)-GGF(I-1,M))/SLP(I,M)/2. ELSE GMX(M)=0. ENDIF ENDDO DO M=1,4 MH=M AMX=AS(M)+GMX(M) XFI(M)=FF(AS(M))*(GGF(I,M)+GGF(I-1,M))/2. >-FF(AMX)*SLP(I,M) ENDDO ED(1)=XFI(1)+XFI(2)+XFI(3)+XFI(4) ED(2)=(QSR(I,2)-CSR(I))*XFI(1) >+QSR(I,2)*XFI(2)+QSR(I,2)*XFI(3) >+(QSR(I,2)+CSR(I))*XFI(4) ED(3)=QSR(I,3)*XFI(1)+(QSR(I,3)+CSR(I)) >*XFI(2) +(QSR(I,3)-CSR(I))*XFI(3) >+QSR(I,3)*XFI(4) ED(4)=(QSR(I,4)-QSR(I,2)*CSR(I))*XFI(1) >+((QSR(I,2)**2+QSR(I,3)**2)/2. >+QSR(I,3)*CSR(I))*XFI(2)+((QSR(I,2)**2 >+QSR(I,3)**2)/2.->QSR(I,3)*CSR(I))*XFI(3) >+(QSR(I,4)+QSR(I,2)*CSR(I))*XFI(4) DO M=1,4 F(I,J,M)=(EQ(I,M)+EQ(I+1,M)+ED(M))/2. ENDDO ENDDO

FOR of/m=1 ASSUME as=qsr[m=2]-csr. FOR of/m=2 ASSUME as=qsr. FOR of/m=3 ASSUME as=qsr[m=2]. FOR of/m=4 ASSUME as=qsr[m=2]+csr.   FOR oft ASSUME gmx=FF(as,m) *(ggf-ggf[i-1])/slp/2.0.   FOR off ASSUME gmx=0.0.         FOR of ASSUME xfi=FF(as,m)*(ggf+ggf[i-1])/2.0 -FF(as+gmx,m)*slp. FOR of/m=1 ASSUME ed=xfi+xfi[m=2]+xfi[m=3]+xfi[m=4]. FOR of/m=2 ASSUME ed=(qsr-csr)*xfi[m=1] +qsr*xfi+qsr*xfi[m=3]+(qsr+csr)*xfi[m=4]. FOR of/m=3 ASSUME ed=qsr*xfi[m=1]+(qsr+csr)*xfi[m=2] +(qsr-csr)*xfi+qsr*xfi[m=4]. FOR of/m=4 ASSUME ed=(qsr-qsr[m=2]*csr) *xfi[m=1]+((qsr[m=2]**2+qsr[m=3]**2)/2.0 +qsr[m=3]*csr)*xfi[m=2]+((qsr[m=2]**2 +qsr[m=3]**2)/2.0-qsr[m=3]*csr)*xfi[m=3] +(qsr+qsr[m=2]*csr)*xfi. FOR of ASSUME f=(eq+eq[i+1]+ed)/2.0.

3. Программа раздела mhdx для вычисления потока F

Таким образом, полная программа раздела с необходимыми описаниями областей и величин имеет следующий вид.

PART mhdx.

!входные параметры для распределенной программы

gam1,prw,prx,QU

!результаты вычисления раздела

RESULT f

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

DISTRIBUTED INDEX i=1..6,j=1..5.

oi:(i=1..Nx). oi1:(i=1..Nx-1). oi2nm2:(i=2..Nx-2). oi2:(i=1..Nx-2).

oj: (j =1..Ny). oj3nm2: (j =3..Ny-2).

m5: (m=1..5). m4: (m=1..4).

! i=1,Nx;j=1,Ny;m=1,4 i=2,Nx-2;j=3,Ny-2;m=1,4

oijm: (oi;oj;m4). of: (oi2nm2;oj3nm2;m4).

!i=1,Nx;j=1,Ny i=1,Nx;j=1,Ny;m=1,5

oprx:(oi;oj). oprw:(oprx;m5).

!i=1,Nx;j=3,Ny-2;m=1,4 i=1,Nx-1;j=3,Ny-2;m=1,5

oeq: (oi;oj3nm2;m4). oslp: (oi1;oj3nm2;m5).

oslp1: (oi1;oj3nm2).

oslp2: (oi1;oj3nm2;m4).

!i=1,Nx-2;j=3,Ny-2;m=1,4; i=1,Nx-1;j=3;Ny-2

oggf: (oi2;oj3nm2;m4). orsl: (oi1;oj3nm2).

oft,off:of/slp<>0.

EXTERNAL FUNCTION FF.

VARIABLE QU DEFINED ON oijm.

VARIABLE eq DEFINED ON oeq.

VARIABLE prw DEFINED ON oprw.

VARIABLE prx DEFINED ON oprx.

VARIABLE qsr,rsr,rsl,alf,alg DEFINED ON oslp.

VARIABLE slp DEFINED ON oslp2.

VARIABLE csr DEFINED ON oslp1.

VARIABLE rsr,rsl,alf,alg DEFINED ON oslp.

VARIABLE ggf DEFINED ON oggf.

VARIABLE as,gmx,xfi,ed,f DEFINED ON of.

VARIABLE gam1 REAL.

VARIABLE sl, stx DEFINED ON oggf.

!i=1,Nx;j=3,Ny-2;m=1,4

FOR oeq/m=1 ASSUME eq=QU[m=2].

FOR oeq/m=2 ASSUME eq=QU[m=2]**2/QU[m=1]+prx.

FOR oeq/m=3 ASSUME eq=QU[m=2]*QU/QU[m=1].

FOR oeq/m=4 ASSUME eq=(QU+prx)*QU[m=2]/QU[m=1].

!i=1,Nx-1;j=3,Ny-2;m=1,5

FOR oslp ASSUME rsl=SQRT(prw[m=1]);

rsr = SQRT(prw[i+1,m=1]);

alf=rsl/(rsl+rsr);

alg=1.0-alf.

FOR oslp/m=1 ASSUME qsr=rsl*rsr.

FOR oslp/m=2..4 ASSUME qsr=alf*prw/prw[m=1]+alg*prw[i+1]/prw[i+1,m=1].

FOR oslp/m=5 ASSUME qsr=alf*prw+alg*prw[i+1]

+gam1/2.0*qsr[m=1]/(rsl+rsr)**2*((prw[i+1,m=2]/prw[m=1]

-prw[m=2]/prw[m=1])**2+(prw[i+1,m=3]/prw[m=1]

-prw[m=3]/prw[m=1])**2).

FOR oslp1 ASSUME

COMPUTE cslx(gam1, QU ON oijm/(i=i+1,j=j),

QU ON oijm/(i=i,j=j), qsr ON oslp/(i=i,j=j)

RESULT slp ON oslp2/(i=i,j=j),csr).

!i=1,Nx-2;j=3,Ny-2;m=1,4

FOR oggf ASSUME sl=2.0*slp;

stx=SIGN(1.0,sl);

ggf=stx*AMAX1(0.0,AMIN1(ABS(sl),stx*2.0*slp[i+1],

stx*(slp+slp[i+1])/2.0)).

!i=2,Nx-2;j=3,Ny-2;m=1,4

FOR of/m=1 ASSUME as=qsr[m=2]-csr.

FOR of/m=2 ASSUME as=qsr.

FOR of/m=3 ASSUME as=qsr[m=2].

FOR of/m=4 ASSUME as=qsr[m=2]+csr.

FOR oft ASSUME gmx=FF(as,m)*(ggf-ggf[i-1])/slp/2.0.

FOR off ASSUME gmx=0.0.

FOR of ASSUME xfi=FF(as,m)*(ggf+ggf[i-1])/2.0-FF(as+gmx,m)*slp.

FOR of/m=1 ASSUME ed=xfi+xfi[m=2]+xfi[m=3]+xfi[m=4].

FOR of/m=2 ASSUME ed=(qsr-csr)*xfi[m=1]+qsr*xfi

+qsr*xfi[m=3]+(qsr+csr)*xfi[m=4].

FOR of/m=3 ASSUME ed=qsr*xfi[m=1]+(qsr+csr)*xfi[m=2]

+(qsr-csr)*xfi+qsr*xfi[m=4].

FOR of/m=4 ASSUME ed=(qsr-qsr[m=2]*csr)*xfi[m=1]

+((qsr[m=2]**2+qsr[m=3]**2)/2.0+qsr[m=3]*csr)*xfi[m=2]

+((qsr[m=2]**2+qsr[m=3]**2)/2.0-qsr[m=3]*csr)*xfi[m=3]

+(qsr+qsr[m=2]*csr)*xfi.

FOR of ASSUME f=(eq+eq[i+1]+ed)/2.0.

END PART.

Программа для вычисления потока G – раздел mhdy – приведена в Приложении 2. Обсуждать ее мы не будем, так как она аналогична разделу mhdx с той лишь разницей, что вычисления проводятся по другим (симметричным) областям, в которых индексы i и j принимают другие (симметричные) значения, а величины в расчетных формулах имеют смещения j+1по индексу j (вместо смещения i+1по индексу i). Эти модификации вызваны изменением индексов в исходной программе в циклах для вычисления G по закону

DO I=3,NX-2

DO J=...,...

в то время как для потока F циклы имели вид

DO J=3,NY-2

DO I=...,...

4. Завершение итерации по методу Рунге-Кутта

Вычисленные в итерации по методу Рунге-Кутта, выполняемой по индексу итерации irk, значения потоков f и g используются для вычисления значений величины qu2 (во внутренней области расчета, см. Рис. 2 в п.5.3). Напомним, что qu2 является дополнительной переменной-“копией”, которая введена, чтобы избежать переприсваивания при вычислении величины qu. Эти вычисления описаны двумя операторами языка Норма:

COMPUTE coeff (irk, A9,B9,C9,D9 RESULT A1,A2).

FOR oin ASSUME

QU2=A1*QU[itime-1]+A2*(QU2[irk-1]-dtr*(f-f[i-1]+g-g[j-1])).

первый из которых вызывает подпрограмму на Фортране, реализующую определение констант в схеме Рунге-Кутта 3-го порядка. Тело этой подпрограммы совпадает с соответствующими операторами исходной Фортран-программы (приведенными в правй колонке таблицы):

SUBROUTINE coeff(ijk,A9,B9,C9,D9,A1,A2) IF(IJK.EQ.1) THEN A1=0. A2=1. ELSE IF(IJK.EQ.2) THEN A1=B9 A2=A9 ELSE A1=C9 A2=D9 ENDIF ENDIF RETURN END

IF(IJK.EQ.1) THEN A1=0. A2=1. ELSE IF(IJK.EQ.2) THEN A1=B9 A2=A9 ELSE A1=C9 A2=D9 ENDIF ENDIF

5. Вычисление величины qu на границах

Далее, на очередном шаге итерации по времени itime осуществляется корректировка значений величины qu во внутренней области oin расчета (см. Рис. 1 из п.5.3), вычисление значений величины qu на верхней границе og4 области расчета и вычисление значений величины qu на верхней границе og3 области расчета. Описание этих вычислений на Норме приведено ниже.

! Корректировка значений QU

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME

COMPUTE absqu (QU2 ON oin/(i=i,j=j),i,j,an,x1 ON ok1

RESULT QU ON oin/(i=i,j=j)).

! Расчет на нижней границе

! i=3..Nx-2; j=1..2;m=1..4

FOR og33 ASSUME

COMPUTE downbound (QU ON og3/(i=i,j=3),i,x1 ON ok1,an,pi,gam1

RESULT QU ON og3/(i=i,j=j)).

! Расчет на верхней границе

! i=3..Nx-2; j= Ny-1,Ny;m=1..4

FOR og44 ASSUME

COMPUTE upbound (i,j,x1 ON ok1, x2 ON ok2,dt0it,fi,an,pi,gam1

RESULT QU ON og4/(i=i,j=j)).

Подпрограммы absqu, downbound и upbound фактически совпадают с соответствующими операторами исходной Фортран-программы и приведены в Приложении 3.

4. Трансляция и запуск Норма-программы.

В итоге, на основании исходной последовательной программы на Фортране (Приложение 1), написана программа на Норме (Приложение 2), вызывающая подпрограммы на Фортране (Приложение 3 ). Для того, чтобы получить параллельную программу для счета на 30 процессорах (это число процессоров определено в Норма-программе) распределенной системы, надо странслировать Норма-программу (пусть она находится в файле example.hop), задав команду

norma example.hop mpi

затем собрать параллельную программу, задав команду

normacnf example.fmp exampleadd.for /mpi

( здесь example.fmp – результат трансляции Норма-программы, exampleadd.for – файл с дополнительными Фортран-подпрограммами). В результате будет получена параллельная программа example.f, которая может быть странслирована Фортран-компилятором, например, так:

mpif77 -o norma_go example.f

после чего запущена на счет:

mpirun –np 31 norma_go

Подробнее о компиляции и запуске программ на языке Норма можно прочитать в [14].

 

5. Литература

1. А.Н.Андрианов, А.Б.Бугеря, К.Н.Ефимкин, И.Б.Задыхайло. Норма. Описание языка. Рабочий стандарт. Препринт ИПМ им. М.В.Келдыша РАН, N 120, 1995, 50 с.
2. И.Б.Задыхайло, К.Н.Ефимкин. Содержательные обозначения и языки нового поколения. Информационные технологии и вычислительные системы, 1996, N 2, с.46-58.
3. http: //www.keldysh.ru/norma
4. Информационно-аналитический центр PARALLEL.RU. http://parallel.ru
5. А.Н.Андрианов, Г.Н.Гусева, И.Б.Задыхайло. Применение языка Норма для расчета дозвукового течения вязкого газа. Математическое моделирование, т.11, N 9, 1999, с. 45-53.
6. A.N.Andrianov, S.B.Bazarov, A.B.Bugerya, K.N.Efimkin. Solution of three-dimensional problems of gas dynamics on multiprocessors computers. Computational Mathematics and Modeling, v.10, N2, 1999, p.140-150.
7. A.N.Andrianov, K.N.Efimkin, V.Y.Levashov, I.N.Shishkova. The Norma Language Application to Solution of Strong Nonequilibrium Transfer Processes Problem with Condensation of Mixtures on the Multiprocessor System. Computational Science - ICCS 2001, Proceedings of International Conference, San Francisco, May 2001, Lecture Notes in Computer Science, v.2073, p.502-510.
8. А.Н.Андрианов. Использование языка Норма для решения вычислительных задач на нерегулярных сетках. Тез. Всеросс. научн. конф. "Фундаментальные и прикладные аспекты разработки больших распределенных программных комплексов", Новороссийск, 1998, с.120-123.
9. А.Н.Андрианов, А.В.Жохова, Б.Н.Четверушкин. Использование параллельных алгоритмов для расчетов газодинамических течений на нерегулярных сетках. В сб. Прикладная математика и информатика, М., Диалог-МГУ, 2000, с. 68-76.
10. H.C.Yee. Construction of explicit and implicit symmetric TVD schemes and their applications, J. Comput. Physics, v.68, 1987, p.151.
11. P.L.Roe. Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Physics, v. 43, 1981, p. 357.
12. A.Harten, J.M.Hyman. A self-adjusting grid for the computation of weak solutions of hyperbolic conservations laws, , J. Comput. Physics, v. 50, 1983, p. 235.
13. S.-W.Shu. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II, J. Comput. Physics, v. 83, 1989, p. 32.
14. Бугеря А.Б. Компиляция и запуск программ на языке Норма. Препринт ИПМ им. М.В.Келдыша РАН, N 34, 2001, 12 с.

2. Приложение 1. Текст последовательной Фортран-программы

PROGRAM SHATL

PARAMETER (NX=388,NY=122)

EXTERNAL FF

REAL AS(4),GMX(4),XFI(4),ED(4)

REAL CSR(NX),PRX(NX)

COMMON/WQU/QU(NX,NY,4)

COMMON/WQU0/QU0(NX,NY,4)

COMMON/WFF/F(NX,NY,4)

COMMON/WFG/G(NX,NY,4)

COMMON/WRMX/QSR(NX,5)

COMMON/WRMS/SLP(NX,4)

COMMON/WRMY/EQ(NX,4),GGF(NX,4)

COMMON/SWR/PRW(NX,5)

COMMON/WQ/X1(NX),X2(NY)

COMMON/DEP/EPS,MH

COMMON/GWGW/GAM1

C Задание начальных параметров

C DB - Размер области по Y

C C0 - Постоянная Куранта

C FI - Угол наклона фронта ударной волны к направлению X

C DK - Полное время расчета

C AN - Точка пересечения фронта ударной волны с осью X

C GAMMA - Постоянная адиабаты

DATA DB,C0,EPS/1.,0.5,0.05/

DATA A9,B9,C9,D9/0.25,0.75,0.33333333,0.66666666/

PI=3.1415926

FI=PI/3.

DK=0.015

AN=1./6.

GAMMA=1.4

GAM1=GAMMA-1.

C Вычисление шага по пространству

DX=DB/(NY-1)

C Задание массивов координат центров ячеек сетки по X

X1(1)=-2.*DX

DO I=2,NX

X1(I)=X1(I-1)+DX

ENDDO

C Задание массивов координат центров ячеек сетки по Y

X2(1)=-2.*DX

DO I=2,NY

X2(I)=X2(I-1)+DX

ENDDO

C Начальные данные для вектора консервативных переменных

DO J=1,NY

DO I=1,NX

IF(X1(I).LT.AN+X2(J)/TAN(FI)) THEN

QU(I,J,1)=8.

QU(I,J,2)=8.25*SIN(PI/4.)*QU(I,J,1)

QU(I,J,3)=-8.25*COS(PI/4.)*QU(I,J,1)

QU(I,J,4)=116.5/GAM1+(QU(I,J,2)**2+QU(I,J,3)**2)/2./QU(I,J,1) ELSE

QU(I,J,1)=1.4

QU(I,J,2)=0.

QU(I,J,3)=0.

QU(I,J,4)=1./GAM1

ENDIF

ENDDO

ENDDO

DT0=0.

C Задание значений вектора консервативных переменных с 1 шага

DO I=1,NX

DO J=1,NY

DO L=1,4

QU0(I,J,L)=QU(I,J,L)

ENDDO

ENDDO

ENDDO

C Начало тела основного цикла

777 CONTINUE

IF(DT0.LT.DK) THEN

C Вычисление шага по времени

DTM=1.E+10

DO I=3,NX-2

DO J=3,NY-2

PRR=GAM1*(QU(I,J,4)-(QU(I,J,2)**2+QU(I,J,3)**2)

1 /2./QU(I,J,1))

CL=SQRT(GAMMA*PRR/QU(I,J,1))

UCSX=QU(I,J,2)/QU(I,J,1)

UCSY=QU(I,J,3)/QU(I,J,1)

DTMX=1./(ABS(UCSX)+CL)

DTMY=1./(ABS(UCSY)+CL)

DTM=AMIN1(DTM,DTMX,DTMY)

ENDDO

ENDDO

DT=C0*DX*DTM

DT0=DT0+DT

DTR=DT/DX

DO IJK=1,3

C Вычисление потока в направлении X

DO J=3,NY-2

DO I=1,NX

DO L=1,3

PRW(I,L)=QU(I,J,L)

ENDDO

PRX(I)=GAM1*(QU(I,J,4)-(QU(I,J,2)**2+QU(I,J,3)**2)

1 /2./QU(I,J,1))

PRW(I,4)=QU(I,J,4)+PRX(I)

PRW(I,5)=GAMMA*PRX(I)/QU(I,J,1)

ENDDO

CALL CSXK(1)

DO I=1,NX

RR0=QU(I,J,1)

RRA=QU(I,J,2)

EQ(I,1)=RRA

EQ(I,2)=RRA**2/RR0+PRX(I)

EQ(I,3)=RRA*QU(I,J,3)/RR0

EQ(I,4)=(QU(I,J,4)+PRX(I))*RRA/RR0

ENDDO

DO I=1,NX-1

QR1=QU(I+1,J,1)-QU(I,J,1)

QR2=QU(I+1,J,2)-QU(I,J,2)

QR3=QU(I+1,J,3)-QU(I,J,3)

QR4=QU(I+1,J,4)-QU(I,J,4)

CSR(I)=SQRT(QSR(I,5))

TMX=QSR(I,2)/CSR(I)

TMY=QSR(I,3)/CSR(I)

TM2=GAM1*(TMX**2+TMY**2)/2.

TMC=GAM1*TMX/CSR(I)

TMD=GAM1*TMY/CSR(I)

GMC=GAM1/CSR(I)**2

CS2=1./CSR(I)

SLP(I,1)=((TMX+TM2)*QR1-(CS2+TMC)*QR2-TMD*QR3+GMC*QR4)/2.

SLP(I,2)=((1.-TMY-TM2)*QR1+TMC*QR2+(CS2+TMD)*QR3-GMC*QR4)/2.

SLP(I,3)=((1.+TMY-TM2)*QR1+TMC*QR2-(CS2-TMD)*QR3-GMC*QR4)/2.

SLP(I,4)=((-TMX+TM2)*QR1+(CS2-TMC)*QR2-TMD*QR3+GMC*QR4)/2.

ENDDO

DO I=1,NX-2

DO M=1,4

C WOODWARD

SL=2.*SLP(I,M)

SL1=2.*SLP(I+1,M)

SL2=(SLP(I,M)+SLP(I+1,M))/2.

STX=SIGN(1.,SL)

GGF(I,M)=STX*AMAX1(0.,AMIN1(ABS(SL),STX*SL1,STX*SL2))

ENDDO

ENDDO

DO I=2,NX-2

AS(1)=QSR(I,2)-CSR(I)

AS(2)=QSR(I,2)

AS(3)=QSR(I,2)

AS(4)=QSR(I,2)+CSR(I)

DO M=1,4

IF(SLP(I,M).NE.0.) THEN

MH=M

GMX(M)=FF(AS(M))*(GGF(I,M)-GGF(I-1,M))/SLP(I,M)/2.

ELSE

GMX(M)=0.

ENDIF

ENDDO

DO M=1,4

MH=M

AMX=AS(M)+GMX(M)

XFI(M)=FF(AS(M))*(GGF(I,M)+GGF(I-1,M))/2.-FF(AMX)*SLP(I,M)

ENDDO

ED(1)=XFI(1)+XFI(2)+XFI(3)+XFI(4)

ED(2)=(QSR(I,2)-CSR(I))*XFI(1)+QSR(I,2)*XFI(2)+QSR(I,2)*XFI(3)

1 +(QSR(I,2)+CSR(I))*XFI(4)

ED(3)=QSR(I,3)*XFI(1)+(QSR(I,3)+CSR(I))*XFI(2)

1 +(QSR(I,3)-CSR(I))*XFI(3)+QSR(I,3)*XFI(4)

ED(4)=(QSR(I,4)-QSR(I,2)*CSR(I))*XFI(1)

1 +((QSR(I,2)**2+QSR(I,3)**2)/2.+QSR(I,3)*CSR(I))*XFI(2)

1 +((QSR(I,2)**2+QSR(I,3)**2)/2.-QSR(I,3)*CSR(I))*XFI(3)

1 +(QSR(I,4)+QSR(I,2)*CSR(I))*XFI(4)

DO M=1,4

F(I,J,M)=(EQ(I,M)+EQ(I+1,M)+ED(M))/2.

ENDDO

ENDDO

ENDDO

C Вычисление потока в направлении Y

DO I=3,NX-2

DO J=1,NY

DO L=1,3

PRW(J,L)=QU(I,J,L)

ENDDO

PRX(J)=GAM1*(QU(I,J,4)-(QU(I,J,2)**2+QU(I,J,3)**2)

1 /2./QU(I,J,1))

PRW(J,4)=QU(I,J,4)+PRX(J)

PRW(J,5)=GAMMA*PRX(J)/QU(I,J,1)

ENDDO

CALL CSXK(2)

DO J=1,NY

RR0=QU(I,J,1)

RRA=QU(I,J,3)

EQ(J,1)=RRA

EQ(J,2)=RRA*QU(I,J,2)/RR0

EQ(J,3)=RRA**2/RR0+PRX(J)

EQ(J,4)=(QU(I,J,4)+PRX(J))*RRA/RR0

ENDDO

DO J=1,NY-1

QR1=QU(I,J+1,1)-QU(I,J,1)

QR2=QU(I,J+1,2)-QU(I,J,2)

QR3=QU(I,J+1,3)-QU(I,J,3)

QR4=QU(I,J+1,4)-QU(I,J,4)

CSR(J)=SQRT(QSR(J,5))

TMX=QSR(J,2)/CSR(J)

TMY=QSR(J,3)/CSR(J)

TM2=GAM1*(TMX**2+TMY**2)/2.

TMC=GAM1*TMX/CSR(J)

TMD=GAM1*TMY/CSR(J)

GMC=GAM1/CSR(J)**2

CS2=1./CSR(J)

SLP(J,1)=((TMY+TM2)*QR1-TMC*QR2-(CS2+TMD)*QR3+GMC*QR4)/2.

SLP(J,2)=((1.-TMX-TM2)*QR1+(CS2+TMC)*QR2+TMD*QR3-GMC*QR4)/2.

SLP(J,3)=((1.+TMX-TM2)*QR1-(CS2-TMC)*QR2+TMD*QR3-GMC*QR4)/2.

SLP(J,4)=((-TMY+TM2)*QR1-TMC*QR2+(CS2-TMD)*QR3+GMC*QR4)/2.

ENDDO

DO J=1,NY-2

DO M=1,4

C WOODWARD

SL=2.*SLP(J,M)

SL1=2.*SLP(J+1,M)

SL2=(SLP(J,M)+SLP(J+1,M))/2.

STX=SIGN(1.,SL)

GGF(J,M)=STX*AMAX1(0.,AMIN1(ABS(SL),STX*SL1,STX*SL2))

ENDDO

ENDDO

DO J=2,NY-2

AS(1)=QSR(J,3)-CSR(J)

AS(2)=QSR(J,3)

AS(3)=QSR(J,3)

AS(4)=QSR(J,3)+CSR(J)

DO M=1,4

IF(SLP(J,M).NE.0.) THEN

MH=M

GMX(M)=FF(AS(M))*(GGF(J,M)-GGF(J-1,M))/SLP(J,M)/2.

ELSE

GMX(M)=0.

ENDIF

ENDDO

DO M=1,4

MH=M

AMX=AS(M)+GMX(M)

XFI(M)=FF(AS(M))*(GGF(J,M)+GGF(J-1,M))/2.-FF(AMX)*SLP(J,M)

ENDDO

ED(1)=XFI(1)+XFI(2)+XFI(3)+XFI(4)

ED(2)=QSR(J,2)*XFI(1)+(QSR(J,2)+CSR(J))*XFI(2)

1 +(QSR(J,2)-CSR(J))*XFI(3)

1 +QSR(J,2)*XFI(4)

ED(3)=(QSR(J,3)-CSR(J))*XFI(1)+QSR(J,3)*XFI(2)

1 +QSR(J,3)*XFI(3)+(QSR(J,3)+CSR(J))*XFI(4)

ED(4)=(QSR(J,4)-QSR(J,3)*CSR(J))*XFI(1)

1 +((QSR(J,2)**2+QSR(J,3)**2)/2.+QSR(J,2)*CSR(J))*XFI(2)

1 +((QSR(J,2)**2+QSR(J,3)**2)/2.-QSR(J,2)*CSR(J))*XFI(3)

1 +(QSR(J,4)+QSR(J,3)*CSR(J))*XFI(4)

DO M=1,4

G(I,J,M)=(EQ(J,M)+EQ(J+1,M)+ED(M))/2.

ENDDO

ENDDO

ENDDO

C Выбор констант в схеме Рунге-Кутта 3-го порядка

IF(IJK.EQ.1) THEN

A1=0.

A2=1.

ELSE

IF(IJK.EQ.2) THEN

A1=B9

A2=A9

ELSE

A1=C9

A2=D9

ENDIF

ENDIF

C Вычисление значений вектора переменных на n+1 шаге по времени

DO I=3,NX-2

DO J=3,NY-2

DO L=1,4

FX=F(I,J,L)-F(I-1,J,L)

GY=G(I,J,L)-G(I,J-1,L)

QU(I,J,L)=A1*QU0(I,J,L)+A2*(QU(I,J,L)-DTR*(FX+GY))

ENDDO

ENDDO

ENDDO

ENDDO

C Корректировка значений вектора переменных на n+1 шаге по времени

DO I=3,NX-2

IF(X1(I).GT.AN) THEN

QU(I,3,3)=ABS(QU(I,3,3))

ENDIF

ENDDO

C************** Верхняя граница ***************************

DO J=NY-1,NY

XT1=10.*DT0/SIN(FI)+AN+X2(J)/TAN(FI)

DO I=3,NX-2

IF(X1(I).LT.XT1) THEN

QU(I,J,1)=8.

QU(I,J,2)=8.25*SIN(PI/4.)*QU(I,J,1)

QU(I,J,3)=-8.25*COS(PI/4.)*QU(I,J,1)

QU(I,J,4)=116.5/GAM1+(QU(I,J,2)**2+QU(I,J,3)**2)/2./QU(I,J,1) ELSE

QU(I,J,1)=1.4

QU(I,J,2)=0.

QU(I,J,3)=0.

QU(I,J,4)=1./GAM1

ENDIF

ENDDO

ENDDO

C************** Нижняя граница ***************************

DO I=3,NX-2

IF(X1(I).GT.AN) THEN

QU(I,2,1)=QU(I,3,1)

QU(I,1,1)=QU(I,3,1)

QU(I,2,2)=QU(I,3,2)

QU(I,1,2)=QU(I,3,2)

QU(I,2,3)=0.

QU(I,1,3)=0.

QU(I,2,4)=QU(I,3,4)

QU(I,1,4)=QU(I,3,4)

ELSE

DO J=1,2

QU(I,J,1)=8.

QU(I,J,2)=8.25*SIN(PI/4.)*QU(I,J,1)

QU(I,J,3)=-8.25*COS(PI/4.)*QU(I,J,1)

QU(I,J,4)=116.5/GAM1+(QU(I,J,2)**2+QU(I,J,3)**2)/2./QU(I,J,1)

ENDDO

ENDIF

ENDDO

C**********************************************************

DO I=1,NX

DO J=1,NY

DO L=1,4

QQS=QU(I,J,L)

QU0(I,J,L)=QQS

ENDDO

ENDDO

ENDDO

SMN=1.E10

SMX=-1.E10

DO I=3,NX-2

DO J=3,NY-2

TTX=QU(I,J,1)

IF(TTX.GE.SMX) THEN

SMX=TTX

ENDIF

IF(TTX.LE.SMN) THEN

SMN=TTX

ENDIF

ENDDO

ENDDO

C Конец основного цикла

GOTO 777

ENDIF

C Запись результатов

IF(DT0.GT.DK) THEN

open(40,FILE='REZ')

WRITE(40,*) 'VARIABLES="X","Y","D","U","V"'

WRITE(40,*) 'ZONE T="XX",I=118,J=384,F=POINT'

DO I=3,NX-2

DO J=3,NY-2

TT=QU(I,J,1)

UU=QU(I,J,2)/TT

VV=QU(I,J,3)/TT

WRITE(40,*) X1(I),X2(J),TT,UU,VV

ENDDO

ENDDO

CLOSE(40)

CLOSE(33)

ENDIF

STOP

END

C Расчет средних по методу Рое

SUBROUTINE CSXK(NK)

PARAMETER (NX=388,NY=122)

COMMON/WRMX/QSR(NX,5)

COMMON/SWR/PRW(NX,5)

COMMON/GWGW/GAM1

NN=388

IF(NK.EQ.2) NN=NY

DO I=1,NN-1

RR0=PRW(I,1)

RR1=PRW(I+1,1)

RSL=SQRT(RR0)

RSR=SQRT(RR1)

ALF=RSL/(RSL+RSR)

ALG=1.-ALF

QSR(I,1)=RSL*RSR

QSR(I,2)=ALF*PRW(I,2)/RR0+ALG*PRW(I+1,2)/RR1

QSR(I,3)=ALF*PRW(I,3)/RR0+ALG*PRW(I+1,3)/RR1

QSR(I,4)=ALF*PRW(I,4)/RR0+ALG*PRW(I+1,4)/RR1

QSR(I,5)=ALF*PRW(I,5)+ALG*PRW(I+1,5)+GAM1/2.*

1 QSR(I,1)/(RSL+RSR)**2*

1 ((PRW(I+1,2)/RR0-PRW(I,2)/RR0)**2+

1 (PRW(I+1,3)/RR0-PRW(I,3)/RR0)**2)

ENDDO

RETURN

END

FUNCTION FF(X)

COMMON/DEP/EPS,MH

FF=ABS(X)

IF(MH.EQ.1.OR.MH.EQ.4) THEN

IF(FF.LT.2.*EPS) THEN

FF=(X**2/(4.*EPS))+EPS

ENDIF

ENDIF

RETURN

END

 

3. Приложение 2. Текст Норма-программы

MAIN PART GDST.

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

oijm:((i=1..Nx);(j=1..Ny);(m=1..4)).

ok1:(k=1..Nx). ok2:(k=1..Ny).

VARIABLE QU0,QU DEFINED ON oijm.

VARIABLE x1 DEFINED ON ok1.

VARIABLE x2 DEFINED ON ok2.

VARIABLE gamma,gam1,c0,dt0,an,fi,pi,dx REAL.

pi=3.1415926.

gamma=1.4.

gam1=gamma-1.0.

c0=0.5.

an=1.0/6.0.

fi=pi/3.0.

COMPUTE init(pi,gam1,fi,an

RESULT dx,x1 ON ok1, x2 ON ok2, QU0 ON oijm,dt0).

COMPUTE GD(pi,gamma,gam1,c0,fi,an,dx,dt0,x1 ON ok1, x2 ON ok2, QU0 ON oijm

RESULT QU ON oijm).

COMPUTE writer(x1 ON ok1, x2 ON ok2, QU ON oijm).

END PART.

PART GD.

! входные параметры для распределенной программы

pi,gamma,gam1,c0,fi,an,dx,dt0,x1,x2,QU0

! результаты вычисления раздела

RESULT QU

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

DISTRIBUTION INDEX i=1..6,j=1..5.

oi:(i =1..Nx). oi3nm2:(i=3..Nx-2). oi2nm2:(i=2..Nx-2).

oj:(j =1..Ny). oj3nm2:(j=3..Ny-2). oj2nm2: (j =2..Ny-2).

m4:(m =1..4). m5: (m =1..5).

! i=1,Nx j=1,Ny,m=1,4

oijm:(oi;oj;m4).

!i=3,Nx-2 j=3,Ny-2,m=1,4 i=3,Nx-2 j=3,Ny-2

oin:(oi3nm2;oj3nm2;m4). oijstep:(oi3nm2;oj3nm2).

!i=3,Nx-2 j=1,Ny,m=1,4

oin1:(oi3nm2;oj;m4).

!i=1,2 j=1,Ny,m=1,4 i=Nx-1,Nx j=1,Ny; m=1,4

og1:((i=1..2);oj;m4). og2:((i=Nx-1..Nx);oj;m4).

!i=3,Nx-2 j=1,2;m=1,4 i=3,Nx-2, j=Ny-1,Ny; m=1,4

og3:(oi3nm2;(j=1..2);m4). og4:(oi3nm2;(j=Ny-1..Ny);m4).

og33:(oi3nm2;(j=1..2)). og44:(oi3nm2;(j=Ny-1..Ny)).

!i=3,Nx-2 j=1,2;m=1,4 i=3,Nx-2, j=Ny-1,Ny; m=1,4

og5:(oi3nm2;(j=1..2);m4). og6:(oi3nm2;(j=Ny-1..Ny);m4).

!i=1,Nx j=1,Ny i=1,Nx j=1,Ny m=1,5

oprx:(oi;oj). oprw:(oprx;m5).

!i=2,Nx-2 j=3,Ny-2 m=1,4 i=3,Nx-2 j=2,Ny-2 m=1,4

oif:(oi2nm2;oj3nm2;m4). oig:(oi3nm2;oj2nm2;m4).

ok1:(k=1..Nx). ok2:(k=1..Ny).

VARIABLE QU0,QU,QU2 DEFINED ON oijm.

VARIABLE umin DEFINED ON oijstep.

VARIABLE x1 DEFINED ON ok1.

VARIABLE x2 DEFINED ON ok2.

VARIABLE f DEFINED ON oif.

VARIABLE g DEFINED ON oig.

VARIABLE prx DEFINED ON oprx.

VARIABLE prw DEFINED ON oprw.

VARIABLE pi,gam1,gamma,c0,dt,dt0,dt0it,dk,dx,dtr,an,fi REAL.

VARIABLE A1,A2,A9,B9,C9,D9 REAL.

A9=0.25. B9=0.75. C9= 0.33333333. D9=0.66666666. dk=0.15.

ITERATION QU,dt0it ON itime.

BOUNDARY

!i=1..2; j=1..Ny;m=1..4 & i=Nx-1..Nx; j=1..Ny;m=1..4

FOR og1;og2 ASSUME QU=QU0.

END BOUNDARY

INITIAL itime=0:

dt0it=dt0.

! i=3..Nx-2; j=1..Ny;m=1..4

FOR oin1 ASSUME QU=QU0.

END INITIAL

!i=3..Nx-2; j=3..Ny-2

FOR oijstep ASSUME

COMPUTE step1(gam1,gamma,c0,dx,QU[itime-1] ON oijm/(i=i,j=j) RESULT umin).

dt =MIN((oijstep)umin).

dt0it=dt0it[itime-1]+dt.

OUTPUT dt(FILE='tst').

dtr=dt/dx.

ITERATION QU2 ON irk.

BOUNDARY

! i=1..2; j=1..Ny;m=1..4 & i=Nx-1..Nx; j=1..Ny;m=1..4

! i=3..Nx-2; j=1..2;m=1..4 & i=3..Nx-2; j=Ny-1..Ny;m=1..4

FOR og1;og2;og3;og4 ASSUME QU2=QU[itime-1].

END BOUNDARY

INITIAL irk=0:

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME QU2=QU[itime-1].

END INITIAL

! i=1..Nx;j=1..Ny

FOR oprx ASSUME prx=gam1*(QU2[irk-1,m=4]

-( QU2[irk-1,m=2]**2+ QU2[irk-1,m=3]**2)/2.0/ QU2[irk-1,m=1]).

! i=1,Nx;j=1,Ny; m=1,3

FOR oprw/m=1..3 ASSUME prw=QU2[irk-1].

! i=1,Nx;j=1,Ny; m=4

FOR oprw/m=4 ASSUME prw=QU2[irk-1]+prx.

! i=1,Nx;j=1,Ny; m=5

FOR oprw/m=5 ASSUME prw=gamma*prx/QU2[irk-1,m=1].

! i=2..Nx-2;j=3..Ny-2;m=1..4

COMPUTE mhdx (gam1,prw ON oprw, prx ON oprx, QU2[irk-1] ON oijm

RESULT f ON oif).

! i=3..Nx-2;j=2..Ny-2;m=1..4

COMPUTE mhdy (gam1,prw ON oprw, prx ON oprx, QU2[irk-1] ON oijm

RESULT g ON oig).

COMPUTE coeff (irk, A9,B9,C9,D9 RESULT A1,A2).

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME QU2=A1*QU[itime-1]+A2*(QU2[irk-1]-dtr*(f-f[i-1]+g-g[j-1])).

EXIT WHEN (irk=3).

END ITERATION irk.

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME

COMPUTE absqu (QU2 ON oin/(i=i,j=j),i,j,an,x1 ON ok1

RESULT QU ON oin/(i=i,j=j)).

! i=3..Nx-2; j=1..2;m=1..4

FOR og33 ASSUME

COMPUTE downbound (QU ON og3/(i=i,j=3),i,x1 ON ok1,an,pi,gam1

RESULT QU ON og3/(i=i,j=j)).

! i=3..Nx-2; j= Ny-1,Ny;m=1..4

FOR og44 ASSUME

COMPUTE upbound (i,j,x1 ON ok1, x2 ON ok2,dt0it,fi,an,pi,gam1

RESULT QU ON og4/(i=i,j=j)).

EXIT WHEN (dt0it[itime]> dk).

END ITERATION itime.

END PART.

PART mhdx.

!входные параметры для распределенной программы

gam1,prw,prx,QU

!результаты вычисления раздела

RESULT f

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

DISTRIBUTED INDEX i=1..6,j=1..5.

oi: (i =1..Nx). oi1: (i=1..Nx-1). oi2nm2: (i =2..Nx-2). oi2: (i=1..Nx-2).

oj: (j =1..Ny). oj3nm2: (j =3..Ny-2).

m5: (m=1..5). m4: (m=1..4).

! i=1,Nx;j=1,Ny;m=1,4 i=2,Nx-2;j=3,Ny-2;m=1,4

oijm: (oi;oj;m4). of: (oi2nm2;oj3nm2;m4).

!i=1,Nx;j=1,Ny i=1,Nx;j=1,Ny;m=1,5

oprx:(oi;oj). oprw:(oprx;m5).

!i=1,Nx;j=3,Ny-2;m=1,4 i=1,Nx-1;j=3,Ny-2;m=1,5

oeq: (oi;oj3nm2;m4). oslp: (oi1;oj3nm2;m5).

oslp1: (oi1;oj3nm2).

oslp2: (oi1;oj3nm2;m4).

!i=1,Nx-2;j=3,Ny-2;m=1,4; i=1,Nx-1;j=3;Ny-2

oggf: (oi2;oj3nm2;m4). orsl: (oi1;oj3nm2).

oft,off:of/slp<>0.

EXTERNAL FUNCTION FF.

VARIABLE QU DEFINED ON oijm.

VARIABLE eq DEFINED ON oeq.

VARIABLE prw DEFINED ON oprw.

VARIABLE prx DEFINED ON oprx.

VARIABLE qsr,rsr,rsl,alf,alg DEFINED ON oslp.

VARIABLE slp DEFINED ON oslp2.

VARIABLE csr DEFINED ON oslp1.

VARIABLE rsr,rsl,alf,alg DEFINED ON oslp.

VARIABLE ggf DEFINED ON oggf.

VARIABLE as,gmx,xfi,ed,f DEFINED ON of.

VARIABLE gam1 REAL.

VARIABLE sl, stx DEFINED ON oggf.

!i=1,Nx;j=3,Ny-2;m=1,4

FOR oeq/m=1 ASSUME eq=QU[m=2].

FOR oeq/m=2 ASSUME eq=QU[m=2]**2/QU[m=1]+prx.

FOR oeq/m=3 ASSUME eq=QU[m=2]*QU/QU[m=1].

FOR oeq/m=4 ASSUME eq=(QU+prx)*QU[m=2]/QU[m=1].

!i=1,Nx-1;j=3,Ny-2;m=1,5

FOR oslp ASSUME rsl=SQRT(prw[m=1]);

rsr = SQRT(prw[i+1,m=1]);

alf=rsl/(rsl+rsr);

alg=1.0-alf.

FOR oslp/m=1 ASSUME qsr=rsl*rsr.

FOR oslp/m=2..4 ASSUME qsr=alf*prw/prw[m=1]+alg*prw[i+1]/prw[i+1,m=1].

FOR oslp/m=5 ASSUME qsr=alf*prw+alg*prw[i+1]

+gam1/2.0*qsr[m=1]/(rsl+rsr)**2*((prw[i+1,m=2]/prw[m=1]

-prw[m=2]/prw[m=1])**2+(prw[i+1,m=3]/prw[m=1]

-prw[m=3]/prw[m=1])**2).

FOR oslp1 ASSUME

COMPUTE cslx(gam1, QU ON oijm/(i=i+1,j=j),

QU ON oijm/(i=i,j=j), qsr ON oslp/(i=i,j=j)

RESULT slp ON oslp2/(i=i,j=j),csr).

!i=1,Nx-2;j=3,Ny-2;m=1,4

FOR oggf ASSUME sl=2.0*slp;

stx=SIGN(1.0,sl);

ggf=stx*AMAX1(0.0,AMIN1(ABS(sl),stx*2.0*slp[i+1],

stx*(slp+slp[i+1])/2.0)).

!i=2,Nx-2;j=3,Ny-2;m=1,4

FOR of/m=1 ASSUME as=qsr[m=2]-csr.

FOR of/m=2 ASSUME as=qsr.

FOR of/m=3 ASSUME as=qsr[m=2].

FOR of/m=4 ASSUME as=qsr[m=2]+csr.

FOR oft ASSUME gmx=FF(as,m)*(ggf-ggf[i-1])/slp/2.0.

FOR off ASSUME gmx=0.0.

FOR of ASSUME xfi=FF(as,m)*(ggf+ggf[i-1])/2.0-FF(as+gmx,m)*slp.

FOR of/m=1 ASSUME ed=xfi+xfi[m=2]+xfi[m=3]+xfi[m=4].

FOR of/m=2 ASSUME ed=(qsr-csr)*xfi[m=1]+qsr*xfi

+qsr*xfi[m=3]+(qsr+csr)*xfi[m=4].

FOR of/m=3 ASSUME ed=qsr*xfi[m=1]+(qsr+csr)*xfi[m=2]

+(qsr-csr)*xfi+qsr*xfi[m=4].

FOR of/m=4 ASSUME ed=(qsr-qsr[m=2]*csr)*xfi[m=1]

+((qsr[m=2]**2+qsr[m=3]**2)/2.0+qsr[m=3]*csr)*xfi[m=2]

+((qsr[m=2]**2+qsr[m=3]**2)/2.0-qsr[m=3]*csr)*xfi[m=3]

+(qsr+qsr[m=2]*csr)*xfi.

FOR of ASSUME f=(eq+eq[i+1]+ed)/2.0.

END PART.

PART mhdy.

!входные параметры для распределенной программы

gam1,prw,prx,QU

!результаты вычисления раздела

RESULT g

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

DISTRIBUTED INDEX i=1..6,j=1..5.

oj: (j =1..Ny). oj2: (j=1..Ny-2). oj2nm2: (j =2..Ny-2).

oj1: (j =1..Ny-1).

oi: (i =1..Nx). oi1: (i=1..Nx-1). oi3nm2:(i=3..Nx-2).

m5: (m=1..5). m4: (m=1..4).

! i=1,Nx;j=1,Ny;m=1,4 i=2,Nx-2 j=3,Ny-2 m=1,4

oijm: (oi;oj;m4). of: (oi3nm2;oj2nm2;m4).

!i=1,Nx;j=1,Ny; i=1,Nx;j=1,Ny;m=1,5

oprx:(oi;oj). oprw:(oprx;m5).

!i=3,Nx-2;j=1,Ny;m=1,4 i=3,Nx-2;j=1,Ny-1;m=1,5

oeq: (oi3nm2;oj;m4). oslp: (oi3nm2;oj1;m5).

oslp1: (oi3nm2;oj1).

oslp2: (oi3nm2;oj1;m4).

!i=3,Nx-2;j=1,Ny-2;m=1,4 i=3,Nx-2;j=1,Ny-1

oggf: (oi3nm2;oj2;m4). orsl: (oi3nm2;oj1).

oft,off:of/slp<>0.

EXTERNAL FUNCTION FF.

VARIABLE QU DEFINED ON oijm.

VARIABLE eq DEFINED ON oeq.

VARIABLE prw DEFINED ON oprw.

VARIABLE prx DEFINED ON oprx.

VARIABLE qsr,rsr,rsl,alf,alg DEFINED ON oslp.

VARIABLE csr DEFINED ON oslp1.

VARIABLE slp DEFINED ON oslp2.

VARIABLE ggf DEFINED ON oggf.

VARIABLE as,gmx,xfi,ed,g DEFINED ON of.

VARIABLE gam1 REAL.

VARIABLE sl, stx DEFINED ON oggf.

!i=1,Nx;j=3,Ny-2;m=1,4

FOR oeq/m=1 ASSUME eq=QU[m=3].

FOR oeq/m=2 ASSUME eq=QU[m=3]*QU[m=2]/QU[m=1].

FOR oeq/m=3 ASSUME eq=QU[m=3]**2/QU[m=1]+prx.

FOR oeq/m=4 ASSUME eq=(QU[m=4]+prx)*QU[m=3]/QU[m=1].

!i=1,Nx-1;j=3,Ny-2;m=1,5

FOR oslp ASSUME rsl=SQRT(prw[m=1]);

rsr = SQRT(prw[j+1,m=1]);

alf=rsl/(rsl+rsr);

alg=1.0-alf.

FOR oslp/m=1 ASSUME qsr=rsl*rsr.

FOR oslp/m=2..4 ASSUME qsr=alf*prw/prw[m=1]+alg*prw[j+1]/prw[j+1,m=1].

FOR oslp/m=5 ASSUME qsr=alf*prw+alg*prw[j+1]

+gam1/2.0*qsr[m=1]/(rsl+rsr)**2*((prw[j+1,m=2]/prw[m=1]

-prw[m=2]/prw[m=1])**2+(prw[j+1,m=3]/prw[m=1]

-prw[m=3]/prw[m=1])**2).

FOR oslp1 ASSUME

COMPUTE csly(gam1,QU ON oijm/(i=i,j=j+1),QU ON oijm/(i=i,j=j),

qsr ON oslp/(i=i,j=j)

RESULT slp ON oslp2/(i=i,j=j),csr).

!i=1,Nx-2;j=3,Ny-2;m=1,4

FOR oggf ASSUME sl=2.0*slp; stx=SIGN(1.0,sl); ggf=stx*AMAX1(0.0,AMIN1(ABS(sl),stx*2.0*slp[j+1],

stx*(slp+slp[j+1])/2.0)).

!i=2,Nx-2;j=3,Ny-2;m=1,4

FOR of/m=1 ASSUME as=qsr[m=3]-csr.

FOR of/m=2 ASSUME as=qsr[m=3].

FOR of/m=3 ASSUME as=qsr[m=3].

FOR of/m=4 ASSUME as=qsr[m=3]+csr.

FOR oft ASSUME gmx=FF(as,m)*(ggf-ggf[j-1])/slp/2.0.

FOR off ASSUME gmx=0.0.

FOR of ASSUME xfi=FF(as,m)*(ggf+ggf[j-1])/2.0-FF(as+gmx,m)*slp.

FOR of/m=1 ASSUME ed=xfi+xfi[m=2]+xfi[m=3]+xfi[m=4].

FOR of/m=2 ASSUME ed=qsr*xfi[m=1] +(qsr+csr)*xfi[m=2]

+(qsr-csr)*xfi[m=3]+qsr*xfi[m=4].

FOR of/m=3 ASSUME ed=(qsr-csr)*xfi[m=1]

+qsr*xfi[m=2]

+qsr*xfi[m=3]

+(qsr+csr)*xfi[m=4].

FOR of/m=4 ASSUME ed=(qsr-qsr[m=3]*csr)*xfi[m=1]

+((qsr[m=2]**2+qsr[m=3]**2)/2.0+qsr[m=2]*csr)*xfi[m=2]

+((qsr[m=2]**2+qsr[m=3]**2)/2.0-qsr[m=2]*csr)*xfi[m=3]

+(qsr+qsr[m=3]*csr)*xfi.

FOR of ASSUME g=(eq+eq[j+1]+ed)/2.0.

END PART.

 

4. Приложение 3. Текст Фортран-программ, используемых в Норма-программе