CC BY
16
УДК 621.321
Разработка итерационного решателя библиотеки конечно-элементного моделирования магнитного поля и способов ускорения его работы
А.С. Шмелев, А.И. Тихонов, Л.Н. Булатов ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
г. Иваново, Российская Федерация Е-mail: [email protected]
Авторское резюме
Состояние вопроса: Одним из главных направлений развития современных САПР является разработка и совершенствование систем для решения инженерных задач моделирования и анализа физических процессов в технических устройствах. Существующие пакеты для решения подобных задач отличаются закрытостью, сложностью в использовании и высокой стоимостью. Альтернативный подход состоит в использовании программных компонентов и библиотек с открытым кодом для решения полевых задач. В связи с этим актуальной является проблема разработки эффективных решателей с открытым кодом для этих библиотек.
Материалы и методы: Для расчета магнитного поля и ускорения работы итерационного решателя использованы библиотека конечно-элементного моделирования EMLib и технология параллельных вычислений Cuda на графических процессорах nVidia.
Результаты: Разработан математический аппарат итерационного решателя библиотеки конечно-элементного моделирования магнитного поля EMLib. Предложены способы ускорения работы итерационного решателя с использованием экспоненциальной экстраполяции для поиска ожидаемого решения и путем распараллеливания алгоритма в технологии Cuda.
Выводы: Разработанный итерационный решатель библиотеки конечно-элементного моделирования целесообразно использовать при решении задач расчета физических полей на мелких сетках, а также в трехмерной постановке.
Ключевые слова: конечно-элементная модель, магнитное поле, прямые и итерационные решатели систем моделирования физических полей, технология параллельных вычислений Cuda.
Developing an iterative problem solver of the finite-element magnetic field modelling library and speeding up its operation
A.S. Shmelev, A.I. Tikhonov, L.N. Bulatov Ivanovo State Power Engineering University, Ivanovo, Russian Federation
E-mail: [email protected]
Abstract
Background: One of the main directions in modern CAD development is creating and improving systems of engineering simulation and analysis of physical processes in technical devices. The existing packages for such tasks are rather closed, difficult to use and expensive. An alternative approach is to use open-source software components and libraries to solve the field problems. Therefore, it is urgent to develop effective open source problem solvers for these libraries. Materials and methods: The finite-element modeling library EMLibIn and Cuda technology of parallel computing on graphics processors nVidia are used to calculate the magnetic field and to speed up the iterative solver. Results: We have developed mathematical tools of the iterative solver of the magnetic field finite element modeling library EMLib and suggested the ways of speeding up the iterative solver operation by using exponential extrapolation to find expected solutions and by parallelizing algorithms in Cuda technology.
Conclusions: The developed iterative solver of the finite element modeling library should be used in solving problems of physical field calculation in fine grids, as well as in a three-dimensional setting.
Key words: finite-element model, magnetic field, direct and iterative solvers for physical field modeling systems, parallel computing technology, Cuda.
Одним из главных направлений развития современных САПР является разработка и совершенствование CAE-систем (Computer-aided engineering), предназначенных для решения различных инженерных задач, в частности, задач моделирования и анализа физических процессов в технических устройствах. Так, среди CAE-систем, используемых при проектировании электрических машин и аппаратов, наибольшую популярность в России получили пакеты ANSYS, Maxwell, Comsol, ElCut. Однако высокая конкуренция на рынке вынудила их разработчи-
ков отойти от базовых принципов создания САПР - открытости, совместимости и развития. К недостаткам этих систем можно также отнести их высокую стоимость.
Определенную альтернативу данным системам могут составить программные компоненты с открытым кодом, которые могут быть внедрены в разрабатываемое приложение, обеспечивая его функционалом для решения полевых и цепных задач. Так, в настоящее время в открытом пользовании имеется достаточное количество библиотек, обеспечивающих разра-
ботчика программного обеспечения средствами численной математики, с помощью которых можно создавать собственные наукоемкие программные продукты, в частности, модели физических процессов в технических устройствах. Анализ тенденции роста доступности библиотек численной математики показывает, что приходит время, когда доступными для открытого использования при разработке собственных программных средств станут средства для решения задач математической физики.
Например, в ИГЭУ на протяжении многих лет развивается технология численного моделирования электромеханических устройств, основанная на комбинировании компонентов для расчета физических полей и цепей на базе открытых математических пакетов или разрабатываемых пользовательских приложений [1, 2]. Базовым компонентом для реализации данной технологии является библиотека конечно-элементного моделирования магнитного поля БМИЬ, разработанная и развиваемая в ИГЭУ [3, 4, 5]. С помощью данной библиотеки были созданы программные средства для численного моделирования физических процессов в торцевых синхронных генераторах [6], коллекторных машинах [7], асинхронных двигателях [8, 9], реакторах из ленты [10], силовых трансформаторах [11], магнитожидкостных герметизаторах [12] и т.п.
Библиотека БМИЬ поставляет в базовое приложение функции, с помощью которых может быть сгенерирована и исследована в соответствии с заданным программным кодом двухмерная конечно-элементная модель квазистационарного магнитного поля в плоскопараллельной или осесимметричной постановке с учетом нелинейности магнитных характеристик ферромагнитных сред. Задача минимизации функционала в расчетной области сводится к решению системы нелинейных алгебраических уравнений, решаемых методом Ньютона, на каждом шаге итераций которого решается методом Холесского линейная система уравнений. Развитие данной библиотеки ведется как в сторону оптимизации алгоритма, так и в сторону увеличения количества узлов конечно-элементной сетки (количества уравнений).
Из теории известно, что при использовании прямых методов решения систем алгебраических уравнений с ростом количества уравнений время решения и требуемый для хранения матрицы коэффициентов объем памяти возрастают в квадратичной зависимости. Это накладывает ограничение на количество узлов конечно-элементной сетки. Поэтому при количестве узлов в миллион и более целесообразно переходить на итерационные решатели.
Наиболее просто итерационная формула в методе конечных элементов реализуется при решении задачи Лапласа, т.е. при отсутствии источников поля.
В методе конечных элементов расчетную область разбивают на треугольные элементы таким образом, чтобы границы раздела сред совпадали со сторонами треугольников. При этом для каждого /'-го узла конечно-элементной сетки справедлив закон полного тока - циркуляция напряженности магнитного поля Н по замкнутому контуру Е равна потоку вектора плотности тока у через поверхность 5:
§ шЕ = \ рё. (1)
Е ё
При отсутствии источников тока выражение (1) принимает вид
^шЕ = о,
(2)
к=1
где Мд/ - количество треугольных элементов вокруг /-го узла конечно-элементной сетки.
Если в качестве контура интегрирования Е выбрать окружность с произвольным радиусом г (рис. 1), то уравнение (2) примет вид
М
чв,., г
в
ви ,к
(3)
—,—к а ,к =г Ъ—к а ,к = 0 к=1 к=1 ,к
где а/к - угол к-го элемента при /-м узле; В±,к -
составляющая магнитной индукции в к-м элемента при /-м узле, перпендикулярная вектору г/к, соединяющему /-й узел с центром дуги.
Рис. 1. Фрагмент конечно-элементной сетки
Принимая допущение о линейном распределении векторного магнитного потенциала в области треугольного элемента и одинаковом значении индукции внутри элемента, составляющую магнитной индукции, перпендикулярную вектору г,к, можно найти как плотность магнитного потока, пронизывающего контур, образованный данным вектором вдоль оси г, перпендикулярной плоскости расчетной области:
в,1 ,к = 1 (А /,к- А) = 1 (А 2 А ',к,3 - А,к 1
где А,К1 = А;
(4)
А '¡,к,2 = А/,к,1 + ~ (А,к,2 - Л,к,1);
Ч,к,2
А '/,к,3 = А
,к,1
' ¡,к,3
"(А,к,3 - А,к,1).
(5)
(6)
Здесь гк2 и 7к3 - длины сторон к-го элемента при /-м узле.
После подстановки (5), (6) в (4) получаем
В
Т/,к
А '/ ,к,2 г/ ,к,2
А'
/,к,3
7 ,к,3
- А ,к,'
1
1
'/,к,2 '/ ,к,3
(7)
Подставляя (7) в (4), получаем
I
к=1
Чк
Н-/,к
А '/,к,2 г/,к,2
А'/,к ,3 О ,к,3
- А/ ,к,'
1
1 /,к,2
'/,к,3
= 0.
(8)
Уравнение (8) можно преобразовать к виду
( А I ^
А/ = А/,к,1 = мд/ (
II
к=1
V ,!
"7, у-1 Н-/, у-1
а
/, У , У
*/' ,к
Н-/ ,к
1
1
= 0,
(9)
'/ ,к,2 '/,к,3
где М - количество соседних узлов вокруг /-го узла.
В числителе (9) стоит сумма по узлам, окружающим /-й узел, а в знаменателе - сумма по элементам вокруг /-го узла.
При наличии источников тока задача усложняется. Опуская промежуточные выкладки, уравнение (1) можно преобразовать к виду
М/
А; = Р,-Х(АЛ,у), (10)
У=1
где А/ - векторный магнитный потенциал /-го узла; А к - векторный магнитный потенциал к-го узла, соседствующего с /-м узлом; М - количество узлов вокруг 1-го узла;
(, У-11 ,а К УI
'г ,Р
к = У -А/, у-1
У II (К,к)
Н-/, у5д/, У
, з»: у
р =__1 к=1
,кБд/,к
Р =
к=1 ^/к^д/к г +1 при г +1 < 3,
3 I (5/,к )г,г к=1 ^/к^д/к
г - 2 при г +1 > 3,
Я =
(11)
[г + 2 при г + 2 < 3,
г -1 при г + 2 > 3,
(12)
где Мд/ - количество элементов вокруг /-го узла;
ц/к - магнитная проницаемость в к-м элементе
вокруг /-го узла; у1к - плотность тока в к-м эле-
менте вокруг -го узла; г - номер, под которым /-й узел встречается в описании к-го элемента
вокруг /-го узла (нумерация узлов элемента осуществляется против часовой стрелки); в нумерации коэффициентов 5 первый индекс соответствует номеру /-го узла, второй индекс -номеру элемента вокруг этого /-го узла (данный номер совпадает с номером У узла, образующего с -м узлом ребро, от которого данный элемент выстраивается против часовой стрелки,
(/—1 )-й номер соответствует элементу, построенному от данного ребра по часовой стрелке).
Коэффициенты (8/к)гр и ^¡к)^ являются элементами матрицы
( <2 2 , , . . Л
С,
" /,к,2 С * ,с,
(13)
которая вычисляется как 5 /,к = В[кВ/,к, (14)
где матрица производных для к-го элемента вокруг /-го узла
В
Ь/,к,1 Ь/,к,2
/,к
,к,3
"/',к,1
С/,к,2
" /,к ,3
(15)
находится из условия
1
2Б
'д /,к
"В/',к А/,к
дА
дХ
у
/,к
(16)
дА
/,к дУ,
где матрица значений векторного магнитного потенциала имеет вид
А ,к,3 )Г. (17)
А ,к,2
А/,к = (А ,к,1
Частные производные
дА дХ
и
/,к
_дА
дУ
в к-м
/,к
элементе вокруг /-го узла в (16) ввиду линейной конечно-элементной аппроксимации области имеют одинаковое значение во всех точках элемента.
Элементы матрицы В в (15) вычисляются по формулам:
а/,к,г = Х),к,РУ/,к,Я - Х1кЯУ1к,Р,
Ь ¡,к,г = у ,к,Р - у,к,Я, (18)
С/кг = Х/,к,я - Х/,к,р,
где Хщг и У/кг - координаты г-го узла к-го треугольного элемента вокруг -го узла; коэффициенты р и я вычисляются по (12).
Площадь к-го треугольного элемента вокруг /-го узла Бд/к в (11) и (16) рассчитывается по формуле
Бд/,к =
1 Х/,к,1 У',к,1
1 Х/,к,2 У',к,2 , (19)
1 Х/,к,3 У',к,3
Алгоритм решения систем уравнений (9)-(10) основан на итерационной формуле
( м/ Л
А/,Г+1 = А/,г + а/,Г
-1(А,уК,У) - А,Г
У=1
(20)
где t - номер итерации; а, - коэффициенты релаксации, выбор которых зависит от принятой схемы итерационного процесса.
Задача решается методом итераций по узлам. Данный алгоритм реализован в программном коде и протестирован на нескольких задачах, в том числе с учетом нелинейности кривой намагничивания стали.
Б
к
Главным недостатком итерационных методов решения систем алгебраических уравнений является медленная сходимость расчетного процесса. Для ускорения сходимости был использован прием, заключающийся в аппроксимации значений векторного магнитного потенциала в каждом /'-м узле экспонентой вида
А ,, = а,, •е
где
^ 2„ N
к/, (Г-Го)
а = А - Ь.
(21)
(22)
Рис. 2. Аппроксимация значений векторного магнитного потенциала в /-м узле экспонентой
к
1 А • ' А, ! ,С 1 'ч—; / : 1
Рис. 3. Реализация скачков итерационного процесса
1
Рис. 4. Ход итерационного процесса при ускорении с использованием экспонент
Так как экспонента для каждого /'-го узла в процессе итераций асимптотически стремится к
значению Ь (рис. 2), то после определения значений коэффициентов (22) можно реализовать скачок
А,г+1 = Ь, (23)
что существенно ускоряет процесс. Одного скачка недостаточно, поэтому итерационный процесс в итоге складывается из нескольких интервалов, на которых итерации осуществляются естественным путем с последующим вычислением коэффициентов (22), после чего осуществляется скачок (рис. 3). В результате приближение значения векторного магнитного потенциала для -го узла к точному значению принимает вид, показанный на рис. 4.
Для устранения «шумов» вычисление производных первого и второго порядков в (22) осуществляется по усредненным значениям. Для этого интервал естественных итераций делится на пять подинтервалов, на которых усредняются значения А,. По усредненным значениям рассчитываются коэффициенты (22), после чего осуществляется скачок (23).
Тестирование представленного алгоритма итерационного решателя осуществлялось в сравнении с прямым методом решения системы линейных уравнений методом Холесского, реализованного в библиотеке БМПЬ.
Как следует из теории и как показывает исследование, итерационные алгоритмы решения систем уравнений работают существенно медленнее, чем прямые, но только при малом количестве уравнений N. С ростом N время решения прямыми методами растет квадратично, в то время как время решения прямыми методами растет практически линейно. Поэтому, например, на тестовой задаче расчета магнитного поля проводника с током время решения обоими решателями становится одинаковым при N = 140000 (рис. 5).
N
О 50000 100000 150000 200000 25000
Рис. 5. Зависимости времени решения тестовой задачи прямым (1) и итерационным (2) методами от количества уравнений в системе
Однако многое при этом определяется постановкой задачи. Так, в тестовом примере, когда проводник с током окружен ферромагнитным сердечником, итерационный алгоритм существенно уступил по скорости прямому, и
время решения обоими решателями сравнялось только при N = 1800000. Так как библиотека EMLib не позволяет организовать решение прямым методом систем уравнений при таких N, кривая T = f(N) выстраивалась путем экстраполяции экспериментально полученной зависимости.
Все это говорит о том, что итерационные решатели становятся конкурентоспособными, по сравнению с прямыми, только при условии ускорения расчетов с использованием современных технологий параллельных вычислений. В настоящее время для этого широко используется технология Cuda, позволяющая задействовать для математических операций графические ускорители на многопроцессорных видеокартах nVidia. Распараллеливание прямых методов решения систем алгебраических уравнений в данной технологии не приводит к существенному выигрышу во времени. Выигрыш получается именно при использовании итерационных алгоритмов.
Главная проблема, возникающая при распараллеливании алгоритмов в технологии Cuda, состоит в опасности одновременного обращения к одним и тем же переменным из разных нитей (процессов). Для устранения данной проблемы помимо векторов значений векторного магнитного потенциала A в узлах и значений коэффициентов (11) в оперативной памяти на графической карте формируются дополнительные матрицы, в частности, матрица значений Ay в соседних узлах, что позволяет рассчитать A, по (9)—(10) в рамках адресного пространства одной нити. Матрица Ay заполняется после каждой итерации новыми значениями, получаемыми из вектора A, (рис. 6). При этом, для того чтобы не допустить одновременного обращения к одной и той же ячейке вектора A из разных нитей, на этапе подготовки данных для итерационного решателя формируется дополнительная матрица номеров соседних узлов. Особенность этой матрицы в том, что в каждой ее строке присутствуют уникальные номера узлов конечно-элементной сетки. По этим номерам обращение к элементам вектора A может осуществляться одновременно из разных нитей. Описанная технология распараллеливания итерационного алгоритма позволяет достичь ускорения в kNp раз, где k = 0,1-0,5 - коэффициент, зависящий от оптимальности программного кода; Np - количество графических процессоров на видеокарте (на современных видеокартах Np = 100-1000).
iV C*J щ
ж 5. S|SE'SS5cc
А, Аз А4 Ai А6 Ат Ад А? Аю А„
Afj.2 An A*,
V"" -""^^вднвпот^нцияпоя
с V У 0 Am i Atir 0 An.I - A, 0
А, Аз Ag A„ A,f Am Aid Aim A. A„ 0 At A,«: Акця
- .
f i
0 0 0 О A, 0 0 0 A, 0 0 0 0 0
Массив номеров сосодей
12 1 6 21 33 4 147 Э -1 N-1 N «
2 3 23 12 34 10 14B a 12 -1 1 1021 102G
■1 -1 -1 -t 9 -1 ■1 -1 7 -1 ■1 -1 -1
Рис. 6. Распределение адресного пространства на графической карте в целях ускорения работы итерационного решателя
Заключение
Разработанный итерационный решатель библиотеки конечно-элементного моделирования может быть использован при решении задач расчета магнитного поля при количестве узлов конечно-элементной сетки от миллиона и выше. Разработка данного решателя позволит распространить функционал библиотеки на решение задач расчета трехмерных полей. В настоящее время осуществляется оптимизация программного кода в целях ускорения вычислений в технологии С^а.
Список литературы
1. Тихонов А.И. Разработка моделей и методов анализа и синтеза решений в автоматизированном проектировании электромеханических устройств: дис. ... д-ра техн. наук. - Иваново, 2007.
2. Технология численного исследования электрических машин с использованием библиотеки конечно-элементного моделирования магнитного поля / А.И. Тихонов, С.Ю. Кучеров, И.М. Лашманов, Д.В. Рубцов // Вестник ИГЭУ. - 2006. - Вып. 3. - С. 5-8.
3.Свидетельство о государственной регистрации программы для ЭВМ № 2006610323. Программа для ЭВМ: Библиотека полевого моделирования / А.И. Тихонов. - М.: Федеральная служба по интеллектуальной собственности, патентам и товарным знакам, 2006.
4. Свидетельство о государственной регистрации программы для ЭВМ № 2011614852. Платформонезависи-мая библиотека конечно-элементного моделирования магнитного поля / А.И. Тихонов, Л.Н. Булатов. - М.: Федеральная служба по интеллектуальной собственности, патентам и товарным знакам, 2011.
5. Тихонов А.И., Казаков Ю.Б. Программирование численного эксперимента с использованием конечно-элементной модели магнитного поля в объектах электромеханики / Иван. гос. энерг. ун-т. - Иваново, 2008. - 80 с.
6. Тихонов А.И., Шишкин В.П., Рубцов Д.В. Поиск оптимальной конструкции торцевого ветрогенератора с использованием динамической полевой модели // Вестник ИГЭУ. - 2005. - Вып. 3. - С. 43-47.
7. Тихонов А.И., Казаков Ю.Б., Лашманов И.М. Численная модель коммутации коллекторных электрических машин в динамических режимах работы // Вестник ИГЭУ. - 2010. - Вып. 4. - С. 40-43.
8. Булатов Л.Н., Тихонов А.И. Разработка динамической модели асинхронной машины с использованием результатов конечно-элементного расчета // Вестник ИГЭУ. -
2012. - Вып. 4. - С. 32-34.
9. Булатов Л.Н., Тихонов А.И. Динамическая модель асинхронной машины с прямым обращением к конечно-элементной модели магнитного поля // Вестник ИГЭУ. -
2013. - Вып. 1. - С. 40-44.
10. Тихонов А.И., Иванов А.В. Использование библиотеки моделирования магнитного поля при расчете реактора из алюминиевой ленты // Вестник ИГЭУ. - 2009. -Вып. 3. - С. 25-28.
11. Тихонов А.И., Попов Г.В., Еремин И.В. Особенности методики расчета холостого хода трансформатора с сердечником из аморфной стали // Вестник ИГЭУ. -2013. - Вып. 4. - С. 32-36.
12. Казаков Ю.Б., Тихонов А.И. Расчет поля скоростей течения псевдопластической магнитной жидкости в зазоре магнитожидкостного герметизатора с неоднородным магнитным полем, мощности потерь и момента трения при максимальном перепаде давления // XIII Междунар. конф. по нанодисперсным магнитным жидкостям (сентябрь, 2008, Плес, Россия): сб. науч. тр. - Иваново, 2008. - С. 394-400.
References
1. Tikhonov, A.I. Razrabotka modeley i metodov analiza i sinteza resheniy v avtomatizirovannom proektirovanii elek-tromekhanicheskikh ustroystv. Diss. dokt. tekhn. nauk [Developing models and techniques of solution analysis and synthesis in computer-aided design of electromechanical devices. Dr. tech. sci. diss.]. Ivanovo, 2007.
2. Tikhonov, A.I., Kucherov, S.Yu., Lashmanov, I.M., Rubtsov, D.V. Vestnik IGEU, 2006, issue 3, pp. 5-8.
3. Tikhonov, A.I. Svidetel'stvo o gosudarstvennoy regis-tratsii programmy dlya EVM «Biblioteka polevogo modeliro-vaniya», no. 2006610323 [Certificate of state registration of computer program «Field modelling library», no. 2006610323]. Moscow, Federal'naya sluzhba po intellektual'noy sobstven-nosti, patentam i tovarnym znakam, 2006.
4. Tikhonov, A.I., Bulatov, L.N. Svidetel'stvo o gosu-darstvennoy registratsii programmy dlya EVM «Platformonezav-isimaya biblioteka konechno-elementnogo modelirovaniya mag-nitnogo polya», no. 2011614852 [Certificate of state registration of computer program «Platform-independent library of finite-element modeling of magnetic field», no. 2011614852]. Moscow, Federal'naya sluzhba po intellektual'noy sobstven-nosti, patentam i tovarnym znakam, 2011.
5. Tikhonov, A.I., Kazakov, Yu.B. Programmirovanie chislennogo eksperimenta s ispol'zovaniem konechno-elementnoy modeli magnitnogo polya v ob"ektakh elektromek-haniki [Numerical experiment programming using finite-element models of magnetic field in electromechanical objects]. Ivanovo, 2008. 80 p.
6. Tikhonov, A.I., Shishkin, V.P., Rubtsov, D.V. Vestnik IGEU, 2005, issue 3, pp. 43-47.
7. Tikhonov, A.I., Kazakov, Yu.B., Lashmanov, I.M. Vestnik IGEU, 2010, issue 4, pp. 40-43.
8. Bulatov, L.N., Tikhonov, A.I. Vestnik IGEU, 2012, issue 4, pp. 32-34.
9. Bulatov, L.N., Tikhonov, A.I. Vestnik IGEU, 2013, issue 1, pp. 40-44.
10. Tikhonov, A.I., Ivanov, A.V. Vestnik IGEU, 2009, issue 3, pp. 25-28.
11. Tikhonov, A.I., Popov, G.V., Eremin, I.V. Vestnik IGEU, 2013, issue 4, pp. 32-36.
Kazakov, Yu.B., Tikhonov, A.I. Raschet polya skorostey techeniya psevdoplasticheskoy magnitnoy zhidkosti v zazore magnitozhidkostnogo germetizatora s neodnorodnym magnit-nym polem, moshchnosti poter' i momenta treniya pri maksi-mal'nom perepade davleniya [Calculating pseudoplastic magnetic fluid velocity field in the gap of a magnetic fluid sealer with inhomogeneous magnetic field, loss power and friction torque under the maximum change of pressure]. Sbornik nauchnykh tpudov XIII Mezhdunarodnoy konferentsii po nanodispersnym magnitnym zhidkostyam [Scientific works of the XIII-th international conference on nanodisperse magnetic fluids]. Ivanovo, 2008, pp. 394-400.
Шмелев Александр Сергеевич,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», аспирант,
e-mail: [email protected] Тихонов Андрей Ильич,
ФГБОУВПО «Ивановский государственный энергетический университет имени В. И. Ленина», доктор технических наук, доцент, заведующий кафедрой физики, e-mail: [email protected]
Булатов Леонид Николаевич,
ФГБОУВПО «Ивановский государственный энергетический университет имени В. И. Ленина», аспирант,
e-mail: [email protected]
