Филатов Е.Ю., Ясинский Ф.Н. Математическое моделирование течений жидкостей и газов

Подождите немного. Документ загружается.

принадлежит целиком 2-му процессору. Размер сетки 42х22,

нумерация начинается с 0. 0-й и 41-й узлы по горизонтали

являются узлами границы, так же как и 0-й и 21-й узлы по

вертикали. Расчет полей ведется на внутренних узлах, а на

граничных узлах выполняются граничные условия. Узлы между

клиентами распределяются так: каждый клиент забирает 10

внутренних узлов (по горизонтали) плюс по узлу с каждой

стороны с соседних процессоров (на которых ставят граничные

условия первого рода) или с границ.

Библиотека Model, как уже было сказано, во многом

идентична предыдущей программе (п. 12.2, 12.3).

Точкой входа клиентского приложения Client служит

функция WinMain. Она инициализирует MPI, определяет

количество процессоров, номера собственного и соседних

процессоров, вызывает функцию Run и деинициализирует MPI.

Программа построена по технологии событийного управления.

Для взаимодействия процессоров используется топология

«труба». Для первого клиентского процессора предыдущим

процессором считается сервер (0-й процессор). В функции Run

каждый процессор (начиная с первого) ожидает команду от

предыдущего процессора, обрабатывает ее и передает

следующему (если он – не последний).

Обрабатываются следующие команды:

- mpicmdExit – выход из программы;

- mpicmdInit – начать работу (получить необходимые

данные напрямую от сервера и инициализировать объект

CModel);

- mpicmdDone – завершить работу (в действительности

эта команда не выполняет никаких действий);

- mpicmdFieldsToServer – отправить участки полей

данного процессора напрямую на сервер;

- mpicmdNext – выполнить шаг по времени и обменяться

смежными узлами с соседними процессорами;

- mpicmdError – вывести сообщение об ошибке.

Серверное приложение AirAttack построено по

архитектуре «Документ – Представление» технологии MFC,

которая предоставляет довольно обширный шаблонный код.

Рассмотрим лишь ту часть исходного кода данного проекта,

которая касается процедуры распараллеливания и

вычислительной стороны программного комплекса.

В класс CAirAttackDoc были добавлены следующие

поля: режим выполнения программы (однопроцессорный или

многопроцессорный вариант), количество клиентов, координаты

участков полей у клиентов Regions (каждый клиент занимает

по горизонтали узлы с Regions[i-1]+1 по Regions[i] плюс

2 узла с соседних процессоров или с границ), объект CModel,

размеры сетки, шаги по времени и по расстоянию, а также

переменная blWorking, определяющая режим работы

программы, и переменная blAutorefreshing, отвечающая за

автоматическую перерисовку экрана.

На окне программы имеются 3 кнопки: “play”, “stop” и

“next”. Кнопка “next” (нажатие на эту кнопку обрабатывает

функция OnCompNext) рассчитывает один шаг по времени,

кнопка “play” (нажатие на эту кнопку обрабатывает функция

OnCompPlay) запускает расчет шагов в фоновом режиме (при

этом переменной blWorking присваивается значение TRUE), а

кнопка “stop” (нажатие на эту кнопку обрабатывает функция

OnCompStop) останавливает расчет шагов в фоновом режиме

(при этом переменной blWorking присваивается значение

FALSE). Функция OnIdle реализует расчет шагов в фоновом

режиме (функция UpdateAllViews выполняет перерисовку).

Переменная blAutorefreshing определяет, будет ли

происходить перерисовка экрана каждый шаг по времени при

работе в фоновом режиме на нескольких процессорах (отказ от

перерисовки значительно ускоряет расчет). Функция

OnNewDocument выполняет инициализацию полей класса

CAirAttackDoc, посылает команду инициализации клиентам и

создает свой объект CModel.

Рис. 6

Кроме того, в класс CAirAttackDoc добавлены

следующие методы:

- SendCmd – отправить команду первому клиенту;

- SendKoefsToClients – отправить необходимые

данные клиентам для инициализации объекта CModel;

- RecvFieldsFromClients – принять участки полей от

клиентов для формирования полей своего экземпляра класса

CModel.

Функция OnDraw вывода полученных полей находится в

классе CAirAttackView.

Результат работы этого программного комплекса

представлен на рис. 6.

Листинг 4. Проект Model, фрагменты файла Model.cpp

// Model.cpp : Defines the initialization routines for the DLL.

#include "stdafx.h"

#include <afxdllx.h>

#include <math.h>

#include "Model.h"

. . .

#define M_PI 3.14159265358979323846

static AFX_EXTENSION_MODULE ModelDLL = { NULL, NULL };

extern "C" int APIENTRY

DllMain(HINSTANCE hInstance, DWORD dwReason, LPVOID lpReserved)

{

. . .

}

//////////////////////////////////////////////////////////////////////

// Model Class

// Инициализация

//////////////////////////////////////////////////////////////////////

CModel::CModel()

{

for (int k=0; k<nFields; k++)

Fields[k]=NULL;

Cells=NULL;

} // CModel::CModel

CModel::~CModel()

{

FreeFields();

} // CModel::~CModel

// Выделить память для полей

void CModel::AllocFields()

{

int k, i;

for (k=0; k<nFields; k++)

if (!Fields[k]) {

Fields[k] = new double*[Nx+2];

for (i=0; i<=Nx+1; i++)

Fields[k][i] = new double[Ny+2]; }

if (!Cells) {

Cells = new char*[Nx+2];

for (i=0; i<=Nx+1; i++)

Cells[i] = new char[Ny+2]; }

} // CModel::AllocFields

// Освободить память для полей

void CModel::FreeFields()

{

int k, i;

for (k=0; k<nFields; k++)

if (Fields[k]) {

for (i=0; i<=Nx+1; i++)

delete []Fields[k][i];

delete []Fields[k];

Fields[k]=NULL; }

if (Cells) {

for (i=0; i<=Nx+1; i++)

delete []Cells[i];

delete []Cells;

Cells=NULL; }

} // CModel::FreeFields

// Инициализировать модель

void CModel::Init (int akProc, int akLeft, int akRight, int

akNuMethod,

int aNx, int aNy, double ah, double atau, RECT aBuild)

{

FreeFields();

kProc=akProc;

kLeft=akLeft; kRight=akRight;

kNuMethod=akNuMethod;

Nx=aNx; Ny=aNy;

h=ah; tau=atau;

AllocFields();

InitFields (aBuild);

} // CModel::Init

// Инициализировать коэффициенты и поля

void CModel::InitFields (RECT aBuild)

{

// Инициализируем константы

Time=0;

kMu=29*0.001/*кг/моль*/;

kR=8.31/*Дж/моль/К*/;

kMolNu=0.015;

kNu=1000*kMolNu; //в 1000 раз > молекулярной вязкости 0.015

kRo=kMu*kP/kR/kT/*1,2928/*кг/м3*/;

kP=100000/*Па*/;

kT=295/*К*/;

kUx0=20/*м/с*/; kUy0=0/*м/с*/;

kMinT=320; kMaxT=500; /*К*/

kA=1;

kQ_CRo=1;

kMaxPsi=kUx0*Ny*h;

kEpsPsi=0.05;

Build=aBuild;

F=FP=FFr=0;

// Инициализируем карту Cells и поля

int i, j;

for (i=0; i<=Nx+1; i++)

for (j=0; j<=Ny+1; j++) {

Cells[i][j]=kAir;

Ux[i][j]=0; Uy[i][j]=0;

Om[i][j]=0;

Psi[i][j] = -kMaxPsi/2+kMaxPsi*j/(Ny+1);

T[i][j]=kMinT;

Nu[i][j]=10*kMolNu; }

for (i=Build.left; i<=Build.right; i++)

for (j=Build.top; j<=Build.bottom; j++) {

Cells[i][j]=kBuildings;

Nu[i][j]=0; }

// Инициализируем поля на границах

for (j=0; j<=Ny+1; j++) {

if (kLeft==kBoundary) {

Ux[0][j]=kUx0; Uy[0][j]=kUy0; }

if (kRight==kBoundary) {

Ux[Nx+1][j]=kUx0; Uy[Nx+1][j]=kUy0; } }

for (i=0; i<=Nx+1; i++) {

Ux[i][0]=Ux[i][Ny+1]=kUx0;

Uy[i][0]=Uy[i][Ny+1]=kUy0; }

} // CModel::InitFields

//////////////////////////////////////////////////////////////////////

// Расчет щага по времени

//////////////////////////////////////////////////////////////////////

// Сделать шаг по времени

void CModel::UpdateData()

{

int i, j, k;

UpdatePsi (TRUE); // Расчет Psi

for (i=1; i<=Nx; i++) // Расчет Ux, Uy

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) {

Ux[i][j] = (Psi[i][j+1]-Psi[i][j-1])/h/2;

Uy[i][j] = (Psi[i-1][j]-Psi[i+1][j])/h/2; }

UpdateNu (); // Расчет Nu

for (k=0; k<3; k++) // Расчет Omega

UpdateG (k, kOmModel, Om, kNu,

(kLeft==kBoundary)?2:1, 2, (kRight==kBoundary)?3:1, 2);

FindForce(); // Расчет силы F

Time+=tau;

} // CModel::UpdateData

// Сделать шаг по времени для поля Psi

void CModel::UpdatePsi (BOOL blOnlyAir)

{

int i, j;

double newPsi, e;

double ksi = pow ( (cos(M_PI/(Nx-1)) + cos(M_PI/(Ny-1))) / 2, 2);

double theta = 2 * (1-sqrt(1-ksi))/ksi;//1.85;

// Вычисляем Psi

do {

e=0;

for (i=1; i<=Nx; i++)

for (j=1; j<=Ny; j++)

if (!blOnlyAir || Cells[i][j]==kAir) {

//newPsi = (Psi[i+1][j]+Psi[i-1][j]+Psi[i][j+1]+Psi[i][j-1]

+ h*h*Om[i][j]) / 4;

newPsi = (1-theta)*Psi[i][j] + theta/4 * (Psi[i+1][j]

+Psi[i-1][j]+Psi[i][j+1]+Psi[i][j-1] + h*h*Om[i][j]);

if (e < fabs(Psi[i][j]-newPsi))

e = fabs(Psi[i][j]-newPsi);

Psi[i][j]=newPsi; }

} while (e>kEpsPsi);

// Вычисляем Psi на границах

for (i=0; i<=Nx+1; i++) {

Psi[i][0] = -kMaxPsi/2;

Psi[i][Ny+1] = kMaxPsi/2; }

for (j=1; j<=Ny; j++) {

if (kRight==kBoundary)

Psi[Nx+1][j] = Psi[Nx][j];

//if (kRight==kBoundary)

//Psi[Nx+1][j] = 2*Psi[Nx][j]-Psi[Nx-1][j];

if (kLeft==kBoundary)

Psi[0][j] = -kMaxPsi/2+kMaxPsi*j/(Ny+1); }

// Psi внутри здания

if (blOnlyAir)

for (i=2; i<=Nx-1; i++)

for (j=2; j<=Ny-1; j++)

if (Cells[i][j]==kBuildings) Psi[i][j]=0;

} // CModel::UpdatePsi

// Сделать шаг по времени для поля Nu

void CModel::UpdateNu()

{

if (kNuMethod==kLaminaryMeth) return;

} // CModel::UpdateNu

// Сделать шаг по времени (часть шага номер step) для модели kind

// (поле G, коэффициент alfa)

void CModel::UpdateG (int step, TModelKind kind, double **G,

double alfa, int G1, int G2, int G3, int G4)

{

double *L, *M;

L = new double [max(Nx,Ny)+2];

M = new double [max(Nx,Ny)+2];

double tau_hh = tau/h/h, tau_h = tau/h;

double aij, bij, cij, fij, rij;

int i, j, i1, j1;

#define NU kNu

switch (step) {

case 0: // Первая часть шага

if (G1==1) i1=1; else i1=0;

for (j=1; j<=Ny; j++) { // По всем строкам

// Вычисляем прогоночные коэффициенты

// в начале

switch (G1) {

case 1: L[1]=0; M[1]=G[0][j]; break;

case 2: L[1]=1; M[1]=0; break; }

// везде

for (i=1; i<=Nx; i++) if (Cells[i][j]==kAir) {

// рядом со зданием

if (Cells[i-1][j]==kBuildings) {

L[i]=0; M[i]=(Psi[i][j]-Psi[i+1][j])/h/h; }

// вычисляем aij, bij, cij, fij

if (Ux[i][j]>0) aij = -tau_h*Ux[i][j]-NU*tau_hh;

else aij = -NU*tau_hh;

bij = 1 + tau_h*fabs(Ux[i][j]) + 2*NU*tau_hh;

if (Ux[i][j]>0) cij = -NU*tau_hh;

else cij = tau_h*Ux[i][j]-NU*tau_hh;

fij = G[i][j];

// результат

L[i+1] = -cij/(bij+aij*L[i]);

M[i+1] = (fij-aij*M[i])/(bij+aij*L[i]); }

// Вычисляем новое значение поля G

// в конце

switch (G3) {

case 1: NewG[Nx+1][j] = G[Nx+1][j]; break;

case 2: NewG[Nx+1][j] = -M[Nx+1]/(L[Nx+1]-1); break;

case 3: NewG[Nx+1][j] = -((L[Nx]-2)*M[Nx+1]+M[Nx]) / ((L[Nx]-

2)*L[Nx+1]+1); break; }

// везде

for (i=Nx; i>=i1; i--) {

if (Cells[i][j]==kAir)

NewG[i][j] = L[i+1]*NewG[i+1][j]+M[i+1];

// рядом со зданием

if ((Cells[i][j]==kBuildings) && (Cells[i-1][j]==kAir))

NewG[i][j] = (Psi[i-1][j]-Psi[i-2][j])/h/h; } }

// Сохраняем новое значение поля G

for (i=i1; i<=Nx+1; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) G[i][j]=NewG[i][j];

break;

case 1: // Вторая часть шага

if (G1==1) j1=1; else j1=0;

for (i=1; i<=Nx; i++) { // По всем столбцам

// Вычисляем прогоночные коэффициенты

// в начале

switch (G2) {

case 1: L[1]=0; M[1]=G[i][0]; break;

case 2: L[1]=1; M[1]=0; break; }

// везде

for (j=1; j<=Ny; j++) if (Cells[i][j]==kAir) {

// рядом со зданием

if (Cells[i][j-1]==kBuildings) {

L[j]=0; M[j]=(Psi[i][j]-Psi[i][j+1])/h/h; }

// вычисляем aij, bij, cij, fij

if (Uy[i][j]>0) aij = -tau_h*Uy[i][j]-NU*tau_hh;

else aij = -NU*tau_hh;

bij = 1 + tau_h*fabs(Uy[i][j]) + 2*NU*tau_hh;

if (Uy[i][j]>0) cij = -NU*tau_hh;

else cij = tau_h*Uy[i][j]-NU*tau_hh;

fij = G[i][j];

// результат

L[j+1] = -cij/(bij+aij*L[j]);

M[j+1] = (fij-aij*M[j])/(bij+aij*L[j]); }

// Вычисляем новое значение поля G

// в конце

switch (G4) {

case 1: NewG[i][Ny+1] = G[i][Ny+1]; break;

case 2: NewG[i][Ny+1] = -M[Ny+1]/(L[Ny+1]-1); break;

case 3: NewG[i][Ny+1] = -((L[Ny]-2)*M[Ny+1]+M[Ny]) / ((L[Ny]-

2)*L[Ny+1]+1); break; }

// везде

for (j=Ny; j>=j1; j--) {

if (Cells[i][j]==kAir)

NewG[i][j] = L[j+1]*NewG[i][j+1]+M[j+1];

// рядом со зданием

if ((Cells[i][j]==kBuildings) && (Cells[i][j-1]==kAir))

NewG[i][j] = (Psi[i][j-1]-Psi[i][j-2])/h/h; } }

// Сохраняем новое значение поля G

for (i=1; i<=Nx; i++)

for (j=j1; j<=Ny+1; j++)

if (Cells[i][j]==kAir) G[i][j]=NewG[i][j];

if ((G1!=1) && (G2!=1)) G[0][0]=(G[0][1]+G[1][0])/2;

if ((G1!=1) && (G4!=1)) G[0][Ny+1]=(G[0][Ny]+G[1][Ny+1])/2;

if ((G3!=1) && (G2!=1)) G[Nx+1][0]=(G[Nx+1][1]+G[Nx][0])/2;

if ((G3!=1) && (G4!=1)) G[Nx+1][Ny+1]=(G[Nx+1][Ny]+G[Nx][Ny+1])/2;

break;

case 2: // Третья часть шага (Метод расщепления)

for (i=1; i<=Nx; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) {

rij=0;

G[i][j] += tau*rij; }

break;

}

delete []L; delete []M;

} // CModel::UpdateG

// Вычислить F

void CModel::FindForce()

{

if (Build.right==-1) return;

F=FP=FFr=0;

} // CModel::FindForce

Листинг 5. Проект Model, файл Model.h

// Model.h : main header file for the MODEL library

// "Процессорность": однопроцессорный выриант, распараллеливание

enum TProcKind { kOneProc, kServerProc, kClientProc };

// Модель: скорости, вихря, тока, температуры, вязкости

enum TModelKind { kUModel, kOmModel, kPsiModel, kTModel, kNuModel };

// Тип ячейки: воздух (жидкость), здание (тело)

enum TCellKind { kAir, kBuildings };

// Модель вычисления вязкости

enum TNuMethodKind { kLaminaryMeth, kPKMethod, kSekundovMeth };

// Стороны: граница области, соседний проц

enum TSideKind { kBoundary=1, kNextProc=0 };

#define nFields 7

class AFX_EXT_CLASS CModel //: public CObject // Класс "Модель"

{

public:

int kProc, kNuMethod; // "Процессорность", метод выч-я вязкости

int kLeft, kRight; // Стороны (типа TSideKind)

int Nx, Ny; // Размеры сетки

double h, tau, Time; // Шаги по расстоянию и по времени, время

double // Коэффициенты и экстр.значения:

kMu, kMolNu, kNu, kR, // молярная масса, молекулярная вязкость,

вязкость и унив. газовая пост. воздуха

kRo, kP, kT, // плотность, давление, температура воздуха

kUx0, kUy0, // начальная скорость ветра

kMinT, kMaxT, // мин. и макс. температуры

kA, kQ_CRo, // коэффициенты A и (q)/c/ro

kMaxPsi, kEpsPsi; // макс. Psi и ошибка Psi

union { // Поля

(double **)Fields[nFields];

struct { double **Ux, **Uy, **Om, **Psi, **T, **Nu, **NewG; };

};

char **Cells; // Карта типов ячеек

double F, FP, FFr; // Силы: общая, трения и давления

// POINT Area = { (Nx+1)*h, (Ny+1)*h } // Координаты рассматриваемой

области в метрах

RECT Build; // <every_coord>*h // Координаты здания (тела) в метрах

// Инициализация

CModel ();

virtual ~CModel();

void AllocFields(); // Выделить память для полей

void FreeFields(); // Освободить память для полей

void Init (int akProc, int akLeft, int akRight, int akNuMethod,

int aNx, int aNy, double ah, double atau, RECT aBuild); //

Инициализировать модель

void InitFields (RECT aBuild); // Инициализировать коэф-ты и поля

// Расчет шага

void UpdateData(); // Сделать шаг по времени