Два варианта организации вычислений - подпрограмма 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').

5.1.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 ofASSUME as= ;

xfi=;

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

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

5.1.2.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.

5.1.2.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 должны быть заданы индексы распределения в распределенных разделах, вызываемых из распределенных разделов - эта информация используется компилятором для необходимой организации параллельных вычислений.

5.1.2.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

5.1.2.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.1.2.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.

5.1.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=...,...

5.1.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.1.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.

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

В итоге, на основании исходной последовательной программы на Фортране (Приложение 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].

7 Литература

А.Н. Андрианов, А.Б. Бугеря, К.Н. Ефимкин, И.Б. Задыхайло. Норма. Описание языка. Рабочий стандарт. Препринт ИПМ им. М.В.Келдыша РАН, N 120, 1995, 50 с.

И.Б. Задыхайло, К.Н. Ефимкин. Содержательные обозначения и языки нового поколения. Информационные технологии и вычислительные системы, 1996, N 2, с.46-58.

http: //www.keldysh.ru/norma

Информационно-аналитический центр PARALLEL.RU. http://parallel.ru

А.Н. Андрианов, Г.Н. Гусева, И.Б. Задыхайло. Применение языка Норма для расчета дозвукового течения вязкого газа. Математическое моделирование, т.11, N 9, 1999, с. 45-53.

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.

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.

А.Н.Андрианов. Использование языка Норма для решения вычислительных задач на нерегулярных сетках. Тез. Всеросс. научн. конф. "Фундаментальные и прикладные аспекты разработки больших распределенных программных комплексов", Новороссийск, 1998, с.120-123.

А.Н. Андрианов, А.В. Жохова, Б.Н. Четверушкин. Использование параллельных алгоритмов для расчетов газодинамических течений на нерегулярных сетках. В сб. Прикладная математика и информатика, М., Диалог-МГУ, 2000, с. 68-76.

H.C. Yee. Construction of explicit and implicit symmetric TVD schemes and their applications, J. Comput. Physics, v.68, 1987, p.151.

P.L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Physics, v. 43, 1981, p. 357.

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.

S.-W. Shu. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II, J. Comput. Physics, v. 83, 1989, p. 32.

Бугеря А.Б. Компиляция и запуск программ на языке Норма. Препринт ИПМ им. М.В.Келдыша РАН, N 34, 2001, 12 с.

8 Приложение 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

9 Приложение 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-2j=3,Ny-2m=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.

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

C $DIR CONFIG INFORMATION

C $DIR PART init

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

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 ENDC $DIR CONFIG INFORMATIONC $DIR PART writerC $DIR proc:C $DIR send:C $DIR END OF CONFIG INFORMATION

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

C $DIR CONFIG INFORMATION

C $DIR PART step1

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

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

C $DIR CONFIG INFORMATION

C $DIR PART cslx

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

SUBROUTINE cslx(gam1,QU1,QU2,qsr,slp,csr)

real QU1(4),QU2(4),qsr(5),slp(4),csr

QR1=QU1(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

C $DIR CONFIG INFORMATION

C $DIR PART csly

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

SUBROUTINE csly(gam1,QU1,QU2,qsr,slp,csr)

real QU1(4),QU2(4),qsr(5),slp(4),csr

QR1=QU1(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)=((TMY+TM2)*QR1-TMC*QR2-(CS2+TMD)*QR3+GMC*QR4)/2.

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

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

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

RETURN

END

C $DIR CONFIG INFORMATION

C $DIR PART FF

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

FUNCTION FF(X,MH)

EPS=0.05

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

C $DIR CONFIG INFORMATION

C $DIR PART coeff

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

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

C $DIR CONFIG INFORMATION

C $DIR PART upbound

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

SUBROUTINE upbound(i,j,x1,x2,dt0,fi,an,pi,gam1,QU4)

INCLUDE 'TASK.PAR'

real QU4(4)

real x1(Nx),x2(Ny)

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

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

QU4(1)=8.

QU4(2)=8.25*SIN(PI/4.)*QU4(1)

QU4(3)=-8.25*COS(PI/4.)*QU4(1)

QU4(4)=116.5/GAM1+(QU4(2)**2+QU4(3)**2)/2./QU4(1)

ELSE

QU4(1)=1.4

QU4(2)=0.

QU4(3)=0.

QU4(4)=1./GAM1

ENDIF

RETURN

END

C $DIR CONFIG INFORMATION

C $DIR PART downbound

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

SUBROUTINE downbound(QU3,i,x1,an,pi,gam1,QU4)

INCLUDE 'TASK.PAR'

real QU4(4),QU3(4)

real x1(Nx)

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

QU4(1)=QU3(1)

QU4(2)=QU3(2)

QU4(3)=0.

QU4(4)=QU3(4)

ELSE

QU4(1)=8.

QU4(2)=8.25*SIN(PI/4.)*QU4(1)

QU4(3)=-8.25*COS(PI/4.)*QU4(1)

QU4(4)=116.5/GAM1+(QU4(2)**2+QU4(3)**2)/2./QU4(1)

ENDIF

RETURN

END

C $DIR CONFIG INFORMATION

C $DIR PART absqu

C $DIR proc:

C $DIR send:

C $DIR END OF CONFIG INFORMATION

SUBROUTINE absqu(QU2,i,j,an,x1,QU02)

INCLUDE 'TASK.PAR'

real QU2(4),QU02(4)

real x1(Nx)

DO m=1,4

QU02(m)=QU2(m)

ENDDO

IF(X1(I).GT.AN.AND.J.EQ.3) THEN

QU02(3)=abs(QU2(3))

ENDIF

RETURN

END

Файл TASK.PAR

PARAMETER(NX=388,NY=122,ipoints=65,jpoints=25)

Cipoints=[Nx/iproc], jpoints=[NY/jproc]

PARAMETER(iproc=6,jproc=5)

INTEGER ITASKJ(jproc*kproc)

COMMON /NORMAMESSPASS/ITASKJ, IPRWWW, JPRWWW, IVVVVJ,IWWWWJ