Курс, 524 группа Адаптация комплекса программ M2dgd для работы на мвс с использованием среды параллельного программирования ost павлухин Павел

Вид материалаДокументы

Содержание


4.Реализация алгоритма в OST
Подобный материал:
1   2   3

4.Реализация алгоритма в OST


Текущая версия OST [1] написана на языке Python, поэтому прикладной объект также пишется на этом языке, но т. к. Python — язык интерпретируемый, то скорость работы программ на нем будет меньше скорости аналогичных программ, написанных на компилируемых языках. Поэтому, чтобы свести эту разницу к минимуму, вычислительное ядро алгоритма и обходы ячеек пишутся на C/fortran, компилируются в виде динамической библиотеки, и она подключается к объекту на Python, в самом объекте реализуются удаленные вызовы к соседям, вызовы setXYZT, задается топология задачи.


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


#define NVAL 33


typedef struct {double x, y;} point;

typedef struct {double f[NVAL];} cell;

typedef struct {int N,M; point *p; cell *c;} grid;

typedef struct {int typeb,L;point *p; cell *c;} border;


point — структура узла сеточного разбиения с координатами x и y; cell — содержит осредненное значение вектора-решения на ячейке (33 компоненты); grid описывает параметры блока: N, M — число ячеек в двух сторонах прямоугольного блока, p и c — указатели на массивы соответствующих структур; border описывает структуру границы для блока: typeb определяет, что соответствующая граница — физическая (если typeb>1) с соответствующим индексом или определяет статус ячеек, полученных от соседа — необсчитаны (typeb=0) или обсчитаны (typeb=1), L задает число ячеек в граничном массиве.


В соответствии с пользовательским описанием определяется топология для назначения связей между объектами (на Python):


class OST_Plane_Point:

def __init__(self, x, y):

self.coord = None

self.neighbors = []


self.coord = [x,y]

self.time = 0


self.neighbors.append( { "coord":[x, y-1], "obj":None } )

self.neighbors.append( { "coord":[x+1, y], "obj":None } )

self.neighbors.append( { "coord":[x, y+1], "obj":None } )

self.neighbors.append( { "coord":[x-1, y], "obj":None } )


Все объекты будут определяться по двум своим координатам, а ссылки на соседей будут отличаться на единицу в одной из координат (монитор OST установит эти ссылки на фактических соседей, распределенным по процессорам в МВС в процессе инициализации задачи).


Сам прикладной объект наследуется от определенного в OST и в нем определяются все необходимые поля и методы для счета задачи:


class appObject(OST_Object_Abstract):

def __init__(self):

OST_Object_Abstract.__init__(self)

self.relate={} #словарь с записями номер соседа:номер гр стороны в соседе

self.threads={} #хранилище потоков ввода-вывода

self.s={} #хранилище значений f для гр ячеек для передачи соседу

self.f={} #хранилище значений f для гр ячеек для принятия из соседа

self.g={} #хранилище геометрии для ячеек,получаемых от соседа

self.tmpgeom={} #вспомогательное хранилище для геометрии

self.forw={} #вспом для forward

self.flag={} #флаги для проверки получ ячеек от соседей


Каждому направлению в блоке приписывается индекс: «снизу» - 0, «справа» - 1 «сверху» - 2, «слева» - 3. В self.relate будут храниться соответствия направления соседа и стороны в соседе, для которой отсылаются ячейки (например в сосед «справа» отсылаются ячейки для той стороны, которая является «левой» в соседе, следовательно, это соответствие 1 — 3).


В init определяются координаты объекта и его принадлежность к классу (типу) black и white:


def init(self,x,y,kind):

self.kind=kind #тип объекта 0-black,1-white

self.point = OST_Plane_Point(x,y)


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


class point(Structure):#структура узла сетки

_fields_ = [("x", c_double),("y", c_double)]


class cell(Structure):#структура вектора-решения на ячейке

_fields_ = [("f",c_double*NVAL)]

class grid(Structure):#структура,содерж инф о блоке(коорд узлов и f)

_fields_ = [("N",c_int),("M",c_int),("p",POINTER(point)),("c",POINTER(cell))]


class border(Structure):#структура границы-у блока их будет 4

#L-длина границы,type - физ гр(>=2) или гр из соседа(=0,1)

_fields_= [("typeb",c_int),("L",c_int),("p",POINTER(point)),("c",POINTER(cell))]


Далее продолжается инициализация объекта:


self.lib=CDLL('./flib.so') #подключается библиотека с с- и fortran функциями

self.localtime=c_int(-1) # параметр,определяющий обход ячеек в блоке

self.grid=grid() #структура блока

self.border = (border*4)() #4 полосы-границы,инициализируются после

self.dt=c_double(0.0001) #квант времени для процедур fortran

self.dtpseudo=c_double(10000.) #псевдовремя для процедур fortran

self.sleeptime=0.0005 #интервал проверки флагов

self.forslope=c_int(1) #параметр для slope

self.initGrid() #инициализация блока ячеек

self.initBorder() #инициализация границ

self.initValues() #передаются ссылки на блок ячеек и границы,тип объекта, localtime в библиотечные функции (в код на c/fortran)

self.initRelate() #задается соответствие self.relate для всех сторон, по которым у объекта есть соседи

self.geomPrepare() #подготовка геометрии гр ячеек для отсылки соседям


self.point.time+=1

self.setXYZT() #для подтв того,что все соседи подготовили свои гр ячейки


self.getGeoms() # собирает гр ячейки от соседей

self.lib.input_()

self.lib.setinitvector_()


В initGrid() вызывается библиотечная функция для генерации блока ячеек, в качестве параметров ей передаются путь к файлу, содержащему геометрию всех блоков и граничные условия, ссылка на блок ячеек (self.grid) и координаты объекта. В initBorder() также генерируются массивы с граничными ячейками, дополнительно передается ссылка на self.border. Конструкция

self.point.time+=1

self.setXYZT()

увеличивает время объекта и отправляет заявку в монитор на продвижение времени для того, чтобы ячейки, за которыми обращается последующий метод getGeoms(), были подготовлены во всех соседях. Чтобы во время счета не передавать между объектами-соседями каждый раз геометрию ячеек, которая остается неизменной, она передается один раз перед началом счета в getGeoms(). Во всех последующих обменах ячейками будут передаваться только вектора f.


GetGeoms() вызывает getGeom() для каждого фактического объекта-соседа:


def getGeom(self,i): #забир у соседа соотв геометрию гр ячеек и уст-ет указатель на них в библиотеке

self.g[i]=self.point.neighbors[i]["obj"].takeGeom(self.relate[i])

self.lib.setgeompointer(self.g[i],i)


На этом примере показано, как обращаться к соседу: через простой вызов его методов, как если бы он (сосед) был локальным объектом с именем self.point.neighbors[i]["obj"], хотя фактически он может располагаться на другом вычислительном узле (распределение счетных объектов — одна из функций OST); ссылка на полученную геометрию ячеек(g[i]) передается в код на с/fortran. Вызываемый метод takeGeom() в удаленном объекте необходимо определить, а так как все счетные объекты (в том числе вызывающий и вызываемый) принадлежат одному классу appObject, то этот метод нужно определять в том же классе:


def takeGeom(self,i): #заполняет геометрию гр ячеек и отдает их вызвавшему соседу

return self.tmpgeom[i]


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


Self.lib.input_() и self.lib.setinitvector_() - процедуры fortran, которые подготавливают физические параметры задачи и начальные значения на расчетной области.


Для передачи ячеек между соседями (векторов f) используется self.setBorders(), которая для каждого соседа вызывает setBorder:


def setBorder(self,i): #заполняет гр ячейки и передает их вызываемому соседу

self.lib.fillborder(self.s[i],i)

self.point.neighbors[i]["obj"].acceptBorder(self.s[i],self.relate[i])


В удаленном объекте вызывается acceptBorder(), передаются заполненные fillborder ячейки — s[i] для соответствующей (self.relate[i]) стороны соседа. AcceptBorder() обновляет ссылки на присланные ячейки и выставляет флаг о получении ячеек:


def acceptBorder(self,income,i): #принимает ячейки от соседа

self.f[i]=income

self.lib.setpointer(self.f[i],i)

self.flag[i]=1 #значение флага - ячейки получены


Fillborder() для нижней стороны блока выглядит так:


void fillborder(cell* c,int i){

int n,j,index;

switch(i){

case(0):{//из нижней стороны блока

for (n=0;nc[n].f[j];

return;

}






Для других сторон изменяются только индексы ячеек gridArray->c[index].


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


AcceptBorder и acceptBorderForw являются односторонними функциями, т.е. вызывающий их объект не дожидается от соседа результата выполнения этих функций, а выполняет последующие операции, не тратя время на ожидание результата указанных функций.


Для получения данных о ячейках в вычислительном ядре используется вызов getnewcell, который через параметры-ссылки передает флаг R (отвечает за обсчитанность/необсчитанность ячейки или определяет индекс физической границы), g (геометрия ячейки), f (осредненный вектор ячейки). Для выбора обхода используется параметр localtime (который может изменяться и из кода Python):


void getnewcell_(int *R, double *g, double *f){

switch(*localtime){

case (1):{ simplesweep(R,g,f); return;} //omega

case (2):{ innersweep(R,g,f); return;} //slope

case (3):{ ringsweep(R,g,f); return;} //slope

case (4):{ innersweep(R,g,f); return;} //predicval

case (5):{ ringsweep(R,g,f); return;} //predicval

case (6):{if (*kind==black){//forward

smallinnersweep(R,g,f); return;

}

if (*kind==white){//forward

upinnersweep(R,g,f); return;

}

}






Все обходы из алгоритма реализуются одноименными функциями на C, например:


void simplesweep(int *R, double *g, double *f){//обход по ряду слева направо,переход на след верхний ряд

*R=1;

if ((icell==-1)&&(jcell==-1)) {icell=0;jcell=0;getgf(g,f,0,0);return;}//случай первой ячейки

if ((icell==N-1)&&(jcell==M-1)) {icell=0; jcell=0; *R=-1; return;}//последн ячейка в блоке

if (icell==N-1){//последн ячейка в строке

jcell++;

}

icell=(icell+1)%N;

getgf(g,f,icell,jcell);

return;

}


После определения индексов текущей ячейки на обходе вызывается getgf, которая непосредственно заполняет данные в g и f:


void getgf(double *g, double *f,int i,int j){//определение g и f по jcell и icell

int n,index=(N+1)*j;


g[0]=gridArray->p[index+i].x;//нижн лев

g[1]=gridArray->p[index+i].y;

g[2]=gridArray->p[index+i+1].x;//нижн прав

g[3]=gridArray->p[index+i+1].y;

g[4]=gridArray->p[index+N+i+2].x;//верхн прав

g[5]=gridArray->p[index+N+i+2].y;

g[6]=gridArray->p[index+N+i+1].x;//верхн лев

g[7]=gridArray->p[index+N+i+1].y;

index=N*j+i;

for (n=0; nc[index].f[n];

}


Для получения данных от ячеек-соседей в вычислительном ядре вызывается getnewneig(), которая, в отличие от getnewcell(), передает еще и данные о векторе общей стороны для текущей и соседней ячейки. Обход всех соседних ячеек осуществляется в порядке 0, 1, 2, 3 (индексы сторон, определенные выше), за индекс текущей стороны отвечает параметр neigh, при этом рассматриваются три случая: ячейка-сосед лежит внутри блока, либо получена от объекта-соседа, либо является физической границей для блока. Определение нижнего соседа:


void getnewneig_(int *R,double *gn,double *e,double *fn){

int n;

switch(neigh){

case(0):{//сосед снизу

if (jcell>0) {//текущ ячейка не в нижней полосе

getgf(gn,fn,icell,jcell-1);//значения из нижней ячейки

*R=treat[icell][jcell-1];//посчитан или нет сосед

e[0]=gn[6];//вектор общей грани(против часовой)

e[1]=gn[7];

e[2]=gn[4];

e[3]=gn[5];

break;

}

if (jcell==0){//нижняя полоса в блоке

*R=borders[0].typeb;//тип ячейки,смежной к граничной

for (n=0; n
if ((*R==0)||(*R==1)) {//ячейки - не физ граница

gn[0]=borders[0].p[icell].x;//нижн левая

gn[1]=borders[0].p[icell].y;

gn[2]=borders[0].p[icell+1].x;//нижн правая

gn[3]=borders[0].p[icell+1].y;

gn[4]=borders[0].p[icell+N+2].x;//верхняя правая

gn[5]=borders[0].p[icell+N+2].y;

gn[6]=borders[0].p[icell+N+1].x;//верхняя левая

gn[7]=borders[0].p[icell+N+1].y;


}

else{//ячейки - физ граница

gn[0]=borders[0].p[icell].x;//"нижн" левая

gn[1]=borders[0].p[icell].y;

gn[2]=borders[0].p[icell+1].x;//"нижн" правая

gn[3]=borders[0].p[icell+1].y;

gn[4]=gn[2];//верхняя правая

gn[5]=gn[3];

gn[6]=gn[0];//верхняя левая

gn[7]=gn[1];


}

e[0]=gn[6];//вектор общей грани(против часовой)

e[1]=gn[7];

e[2]=gn[4];

e[3]=gn[5];

}

break;

}







Для аписи обновленных значений текущей ячейки из вычислительного ядра вызывается setcurval():


void setcurval_(double* f){

int n,index=N*jcell+icell;

for (n=0; nc[index].f[n]=f[n];

}


Для записи значений в ячейку-сосед вызывается setneival(), рассматриваются случаи соседа в самом блоке и полученного из объекта-соседа:(для соседа снизу)


void setneival_(double* f){

int n;

int ng=neigh-1;//номер текущего соседа

switch(ng){

case(0):{//сосед снизу

if (jcell>0) {//текущ ячейка не в нижней полосе

for (n=0; nc[N*(jcell-1)+icell].f[n]=f[n];

}

if (jcell==0){//текущ ячейка в нижней полосе

if (borders[0].typeb!=0) printf("error in setneival:typeb!=0!!!\n");

for (n=0; n
}

return;

}

case(1):{//сосед справа

if (icell
for (n=0; nc[N*(jcell)+icell+1].f[n]=f[n];

}

if (icell==N-1){//текущ ячейка в левой полосе

if (borders[1].typeb!=0) printf("error in setneival:typeb!=0!!!\n");

for (n=0; n
}

return;

}






Наконец, сама функция счета:


def solve(self): #функция счета

self.localtime.value=1 #выбирается обход simplesweep

for timestep in range (0,100): #задается число шагов по времени(100)

for stepiniter in range (0,2): #число итераций на один шаг(2)

for i in self.relate:

self.border[i].typeb=0 #для всех нефизических границ ячейки объявляются необсчитанными

self.flag[i]=0 #сброс флагов (отвечающих за получение ячеек от объектов-соседей)


#Цикл вида for i in self.relate осуществляется по тем индексам, для которых имеются фактические объекты-соседи.


self.lib.omega_(byref(self.dt)) #запуск omega из библиотеки

self.localtime.value=2 # выбирается обход для следующей процедуры - innersweep

#заявка напродвижение времени:

self.point.time+=1

self.setXYZT()


#------------------------SLOPE-----------------------------------------

self.setBorders() #отсылка ячеек всем объектам-соседям

self.lib.slope_(byref(self.forslope)) #slope на innersweep

for i in self.relate:

while self.flag[i] != 1: # если ячейки от соседа еще не готовы для счета

time.sleep(self.sleeptime)

self.localtime.value=3 #выбирается обход ringsweep

self.lib.slope_(byref(self.forslope)) #slope на ringsweep

for i in self.relate:

self.flag[i] = 0 # сброс флагов получения гр ячеек от соседей

self.localtime.value=4 #выбор innersweep

self.point.time+=1

self.setXYZT()


#-----------------------------PREDICVAL---------------------------------------

#действия аналогичны slope

self.setBorders()

self.lib.predicval_(byref(self.dt)) #счет внутренней части

for i in self.relate:

while self.flag[i] != 1: # если ячейки от соседа еще не готовы для счета

time.sleep(self.sleeptime)


self.localtime.value=5 #выбирается обход ringsweep

self.lib.predicval_(byref(self.dt)) #счет граничной части

for i in self.relate:

self.flag[i] = 0 # сброс флагов получения гр ячеек от соседей

self.localtime.value=6

self.point.time+=1

self.setXYZT()


#------------------------------FORWARD_1---------------------------------------

self.lib.resettreat() #сброс статуса всех ячеек в состояние «необсчитаны»

# случай объекта black

if (self.kind==1):#black

self.setBorders() #отправка ячеек соседям

self.lib.forward_(byref(self.dt),byref(self.dtpseudo)) #forward на smallinnersweep

self.localtime.value=7 #выбор bigringsweep


#случай объекта white

if (self.kind==0):#white

self.lib.forward_(byref(self.dt),byref(self.dtpseudo)) #forward на upinnersweep

for i in self.relate:

while self.flag[i] != 1: # если ячейки от соседа еще не готовы для счета

time.sleep(self.sleeptime)

self.localtime.value=7 #ringsweep

self.logger.startExp("ring forward_1 in white:")

self.lib.forward_(byref(self.dt),byref(self.dtpseudo)) # forward на ringsweep

for i in self.relate:

self.flag[i] = 0 # сброс флагов получения гр ячеек от соседей

self.setBordersForw() #отсылка граничных ячеек из блока и обновленных копий, полученных от соседей, обратно соседям

self.localtime.value=8 #downinnersweep


#----------------------------------FORWARD_2----------------------------------------

#случай объекта black

if (self.kind==1):#black

for i in self.relate:

self.border[i].typeb=1 #присылаемые ячейки имеют статус «обсчитаны»

while self.flag[i] != 2: # если ячейки от соседа еще не готовы для счета

time.sleep(self.sleeptime)

for i in self.relate:

self.flag[i] = 0 # сброс флагов получения гр ячеек от соседей

#обновление из полученных копий граничных ячеек (невязки)

for tm in self.forw:

self.lib.resforward(self.forw[tm],tm)

self.lib.forward_(byref(self.dt),byref(self.dtpseudo)) #forward по bigringsweep

self.localtime.value=8 #выбор bigringsweepback

#продвижение по времени:

self.point.time+=1

self.setXYZT()


#случай объекта white

if (self.kind==0):#white

self.lib.forward_(byref(self.dt),byref(self.dtpseudo)) #forward по downinnersweep

self.localtime.value=9 #выбор downinnersweepback

#продвижение по времени

self.point.time+=1

self.setXYZT()


#----------------------BACKWARD_1--------------------------------------------

self.lib.resettreat() #сброс всех ячеек блока в статус «необсчитаны»

for i in self.relate:

self.border[i].typeb=0 #сброс ячеек границ в статус необсчитаны


#случай объекта black

if (self.kind==1):#black

self.lib.backward_(byref(self.dt),byref(self.dtpseudo)) #backward на bigringsweepback

self.setBorders() #отсылка ячеек white-соседям

self.localtime.value=9 # выбор smallinnersweepback


#случай объекта white

if (self.kind==0): #white

self.lib.backward_(byref(self.dt),byref(self.dtpseudo)) #backward по downinnersweepback

self.localtime.value=10 #выбор ringsweepback и upinnersweepback


#------------------------BACKWARD_2----------------------------------------

#случай объекта black

if (self.kind==1):#black

self.lib.backward_(byref(self.dt),byref(self.dtpseudo)) #backward по smallinnersweepback

self.lib.updateiter_() #обновление компонентов f с обходом simplesweep

self.localtime.value=1 #выбор simplesweep для omega.


#случай объекта white

if (self.kind==0):#white

for i in self.relate:

self.border[i].typeb=1 #принимаемые ячейки от соседей уже посчитаны

while self.flag[i] != 1: # если ячейки от соседа еще не готовы для счета

time.sleep(self.sleeptime)

for i in self.relate:

self.flag[i] = 0 # сброс флагов получения гр ячеек от соседей

self.lib.backward_(byref(self.dt),byref(self.dtpseudo))

self.lib.updateiter_() #обновление компонентов f с обходом simplesweep

self.localtime.value=1 #выбор simplesweep для omega.


#обновление компонент f после всех итераций на данном шаге по времени

self.lib.updatetime_()


# обработка f после завершения счета

self.lib.setfinvector_()


Для «раскрутки» задачи прикладной программист также пишет класс OST_Main_Appobj, в котором определяется число параллельно считаемых объектов и их разбиение на black и white:


class OST_Main_Appobj( OST_Main_Class_Abstract ):

'''Класс раскрутки задачи M2DGD'''


# Задание, количества объектов

NN = 4 #число объектов в горизонтальном направлении

MM = 1 #число объектов в вертикальном направлении

def __init__(self):

'''Функция-конструктор объекта раскрутки'''

OST_Main_Class_Abstract.__init__(self)


def init(self):

'''Функция инициализации счетной модели'''


for n in range( 0, self.NN):

for m in range( 0,self.MM):

obj = self.createObject(appObject) #создается прикладной объект класса appObject

obj.init(n,m,(n+m)%2) #так задаваемые координаты и тип объекта(black-0,white-​1) образуют на пространстве блоков необходимое «шахматное» распределение по типам