О реализации конечно-элементного моделирования в задачах остеосинтеза на кластерных системах СГУ Текст научной статьи по специальности «Физика»

ИНФОРМАТИКА

УДК 519.63:004.272.26

0 РЕАЛИЗАЦИИ КОНЕЧНО-ЭЛЕМЕНТНОГО МОДЕЛИРОВАНИЯ В ЗАДАЧАХ ОСТЕОСИНТЕЗА НА КЛАСТЕРНЫХ СИСТЕМАХ СГУ

Д.К.Андрейченко1, П.В. Ирматов2, М.С. Ирматова2, М.Г.Щербаков2

1 Саратовский государственный университет,

кафедра математической кибернетики и компьютерных наук E-mail: [email protected]

2Саратовский государственный университет,

Поволжский региональный центр новых информационных технологий

E-mail: [email protected], [email protected]

В статье рассмотрена оптимизированная кластерная версия пакета конечно-элементного моделирования, использованная в проекте «Разработка вычислительно-информационных технологий компьютерного моделирования на параллельных вычислительных комплексах травматологических и операционных процессов для оперативной выработки диагностических и лечебных рекомендаций», выполняемом в рамках федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007-2012 годы» по Государственному контракту от 30 сентября 2009 года 02.514.11.4121.

Ключевые слова: математическое моделирование, метод конечных элементов, кластерные системы, параллельное программирование.

About Realisation of Finite-Element Modeling in Problems of an Osteosynthesis on Cluster Systems of SSU

D.K. Andreichenko1, P.V. Irmatov2, M.S. Irmatova2, M.G. Scherbakov2

1 Saratov State University,

Chair of Mathematical Cybernetics and Computer Sciences E-mail: [email protected]

2Saratov State University,

Volga Regional Center of New Informational Technologies E-mail: [email protected], [email protected]

In this article we describe the optimized cluster version of a package of finite-element modeling. This project titled «Development of Computational and information technologies using computer modeling on parallel computing complexes for traumatological and surgical evaluations to enable efficient diagnostic and medical recommendations». The project which commenced on September 30th 2009 was a «state commission» and part of the federal special-purpose programme titled «Research and development on priority directions of the Russian science & technology complex 2007-2012», Project 02.514.11.4121.

Key words: mathematical modeling, finite-element method, cluster systems, parallel programming.

ВВЕДЕНИЕ

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

НАУЧНЫЙ ОТДЕЛ

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

1. КОНЕЧНО-ЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЕ С ТОЧКИ ЗРЕНИЯ МЕХАНИКИ СПЛОШНОЙ СРЕДЫ

Движение сплошной среды, занимающей в пространстве в текущий момент времени t некоторый объем V, описывается т. н. уравнением импульсов [3]

pwk = pFk + V%aki, i, k = 1, 2, 3. (1)

Здесь: и = {nkj}, nik = nki — симметричный тензор напряжений, заданный своими контравариантны-ми компонентами; p — плотность среды; F = {Fk}, k = 1, 2, 3, — вектор массовых сил, действующих на индивидуальные частицы сплошной среды, и заданный своими контравариантными компонентами; w = {wk}, k = 1, 2, 3, — ускорение индивидуальных точек сплошной среды, заданное своими контра-вариантными компонентами; V_, i = 1, 2, 3, — символ ковариантного дифференцирования в используемой системе координат (возможно, подвижной, криволинейной и деформируемой); по одноименному верхнему и нижнему индексу, как обычно в тензорном исчислении, выполняется суммирование [4]. С точки зрения Лагранжа, свойственной механике твердого деформируемого тела, индивидуализация точек сплошной среды выполняется посредством указания их координат r с Vo в начальный момент времени t = t0, когда среда заполняет объем V0. Если при этом вектор перемещений u отсчитывается от некоторой инерциальной системы координат, то векторы скорости v и ускорения w индивидуальных частиц сплошной среды суть v = du/dt, w = d2u/dt2. При необходимости в массовых силах F следует учесть переносные и кориолиссовы силы инерции [5]. Радиус-вектор r индивидуальной точки сплошной среды относительно используемой системы координат задается набором криволинейных координат {xk}, k = 1, 2, 3, а расстояние между двумя бесконечно близкими точками суть ds2 = dr2 = gkidxkdxi, k, i = 1, 2, 3; gki = gik, где gki — ковариантные компоненты метрического тензора используемой криволинейной системы координат. Любой вектор a может быть задан как своими ковариантными компонентами ak, k = 1, 2, 3, так и своими контравариантными компонентами ai, i = 1, 2, 3, причем ak = gkjai. Ковариантные, контравариантные и смешанные компоненты тензоров преобразуются по аналогичному закону [4]. Матрицы, составленные из ковариантных {gki} и контравариантных {gki} компонент метрического тензора являются взаимно-обратными, а смешанные компоненты метрического тензора образуют единичную матрицу. При преобразовании системы криволинейных координат xk = xk(С1, С2,С3), k = 1, 2, 3, компоненты векторов и тензоров в новой системе координат Ck (отмечены символом Л) преобразуются по закону

= dCL aj a. = dxxj a. Aik = dCL Aji A = dxj dx1

a = dxja, a_ = dCi aj, A = dxj dxi A , Aik = 6? d£k Ajl.

По правилам тензорного исчисления операции ковариантного дифференцирования над векторными и тензорными полями выполняются следующим образом [4]:

Viak = dak/dxi + aj j, Vi ak = dak/dxi - aj Yjki, Vi ak = gkj Viaj,

Vi Akj = dAkj /dx1 + Alj rj + Akl rji, Vi Aj = dAj/dxi - Aj Т^ - Akij,

V7 л _ „ „ v7 \lm -pk _ 1 Jks (%s . dgjs dgij\ k

ViAkj = gkl gjm ViA , Tkj = -g + 1т-Г -1— =Гji.

Ковариантные компоненты симметричного тензора деформаций е = {е^}, к,г = 1, 2, 3 связаны с ковариантными компонентами вектора перемещений и формулами Коши [3]:

еы = 2 (Ук Щ + УгПк + Ук Щ Угил) . (2)

Здесь вектор перемещений задан проекциями на оси криволинейной системы координат, соответствующей начальному состоянию деформированной среды, а ковариантное дифференцирование выполняется в метрике начального состояния. С точки зрения Лагранжа, метрические тензоры в начальном дк0) и актуальном состоянии суть = дк0 + 2е*г. Любой вектор а и тензор А могут быть заданы как своими компонентами ак, Агк в метрике начального состояния, так и своими компонентами аак, Агк в метрике актуального состояния, причем (6к — символ Кронекера)

ак = аг {5к + Уг ик), Акг = Ал (5к + у ик )($ + У иг).

С точки зрения Лагранжа, закон сохранения массы индивидуальной частицы среды принимает вид р0вУ0 = рвУ, т.е.

Ро = ^ . (3)

Здесь ро, р — плотность среды в начальном и актуальном состоянии, вУ0 и вУ — объем малой индивидуальной частицы среды в начальном и актуальном состоянии.

В механике деформируемого твердого тела связь тензоров напряжений (в метрике актуального состояния) и деформаций задается либо функционалом и = и[е], учитывающим историю деформирования индивидуальной частицы сплошной среды, либо некоторыми функциями компонент тензора деформаций и компонент тензора скоростей деформаций е = {еу}, j, I = 1, 2, 3, еу = Зеу /ЗЬ

и = ц(е,е). (4)

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

икг = и = и (е), (5)

ро Зекг

причем предполагается, что свободная энергия деформирования (упругий потенциал) и зависит лишь от компонент тензора деформаций е и, возможно, температуры Т. Если полагать, что относительные деформации достаточно малы (порядка 10-3 при упругом деформировании литых металлов), то [6]

кг

= Зи/Зекг, и = и (е). (6)

В частности, если свободная энергия деформирования и является квадратичной (положительно определенной) формой компонент тензора деформаций

и (е) = 2 СкгЛ екг ел, СкгЛ = Сл'кг = Сгкл = Скг'л, (7)

то связь напряжений и деформаций становится линейной

икг = с кгЛ ел

кг = с кгЛ ел. (8)

Для уравнений импульса (1) требуется поставить начальные и граничные условия. Начальные условия требуются лишь в том случае, когда в левой части уравнения импульса (1) учитываются компоненты (относительного) ускорения и состоят в задании в начальный момент времени Ь0 координат и скоростей индивидуальных точек сплошной среды, занимающей объем У0. При задании

граничных условий требуется учитывать, что на некоторой части $1 (рис. 1) граничной поверхности $ = дУ = $1 и могут задаваться смещения и(0) точек деформируемого твердого тела, а на другой

части $2 граничной поверхности могут задаваться действующие на деформируемое твердое тело внешние поверхностные силы р(0), т. е.

пк = 40) (£), г е $1, к = 1, 2, 3.

(9)

Пусть п — единичный вектор нормали к элементарной площадке, проходящей через данную точку сплошной среды. Он задан своими ковариантными компонентами п = {пг}, г = 1, 2, 3. В частности, это может быть и единичный вектор внешней нормали к части $2 граничной поверхности. Известно, что контравариантные компоненты рк развивающегося на данной элементарной площадке вектора р поверхностных сил (для элементарных площадок, проходящих через внутренние точки объема V, р есть вектор полного внутреннего напряжения) связаны с компонентами тензора напряжений и т.н. формулами Коши [3, 6] рк = икгпг. Следовательно,

Рис. 1. К постановке граничных условий

ЦкгПг = рко),

г е $2,

(10)

где р(0) суть контравариантные компоненты заданного на $2 вектора внешних поверхностных сил р(0). Возможна ситуация, когда на некоторых фрагментах $* поверхности $ для некоторых значений индекса к = 1, 2, 3 задаются краевые условия (9), а для остальных значений индекса к задаются краевые условия (10).

Вместе с тем, возможна ситуация, когда требуется отдельно рассмотреть деформирование каждой из частей V/ и V// объема V = V/UVп (например, в силу существенной неоднородности конструкции). Полагаем: дVI = и $3, дVII = Бц и $3, $ = и Бц. В каждом из объемов V/ и V// должно выполняться уравнение импульсов (1), для которого на соответствующих фрагментах поверхности $1 П и $1 П потребуется задать краевые условия типа (9), а на фрагментах поверхности $2 П и $2 П — краевые условия типа (10). На разделяющей объемы V/ и V// лежащей внутри V поверхности $з должны выполняться условия

и(/) = и(//), и (!) пг = а(11 )пг?

г е $з, г е $з.

(11) (12)

Далее, на поверхности контакта $** твердых деформируемых тел I и II (рис. 2) должны совпадать нормальные составляющие векторов перемещений и поверхностных сил

к к И(/) Пк = И(Л) Пк,

г е $ *

г е $ *

(13)

кг _ кг с**

Рис. 2. Контактная поверхность и(/)ПгПк = и(И)ПгПк' г е $ . (14)

При отсутствии трения на контактной поверхности на стороне I и на стороне II отсутствуют касательные составляющие векторов касательных сил, т. е. достаточно потребовать выполнения любых двух равенств из трех равенств (15) и любых двух равенств из трех равенств (16)

_кг

и(/у

кг

)Пг - (иПгП)пк =0, к = 1, 2, 3, г е $** И {II) пг - (и^ )пг Щ )пк = 0, к = 1, 2, 3, г е $ *

(15)

(16)

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

В аналитической механике уравнения движения (равновесия) механических систем обычно получаются при помощи т. н. вариационного принципа возможных перемещений Эйлера - Лагранжа [5], согласно которому, работа всех действующих на индивидуальные точки механической системы внешних и внутренних физических сил, а также сил инерции, на любом виртуальном перемещении тождественно обнуляется. Под виртуальным (т. е. возможным) перемещением понимается любое бесконечно малое перемещение, совместное со связями, наложенными на механическую систему. При этом реакции связей по определению не совершают никакой работы на виртуальных перемещениях. Пусть ¿и — векторное поле виртуальных перемещений точек твердого деформируемого тела. Выполнение краевых условий (9) приводит к следующему ограничению, накладываемому на поле виртуальных перемещений:

¿и = 0, г е Бь (17)

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

¿£ы = 2 (V*,5щ + Vi5пк). (18)

Принцип Эйлера - Лагранжа применительно к механике деформируемого твердого тела принимает

вид

У [р(^к - )5пк - 5еы] (V + £р\0)5пк (Б = 0. (19)

В силу инвариантности операций тензорного анализа подынтегральные выражения в (19) могут вычисляться в любой подходящей криволинейной системе координат, однако элементы объема и площади соответствуют актуальному состоянию. Из (17) следует, что интеграл по фрагменту поверхности Б1 в левой части (19) тождественно обнуляется. Если на некоторых фрагментах Б * поверхности Б для некоторых значений индекса к = 1, 2, 3 ставятся краевые условия (9), а для остальных значений индекса к задаются краевые условия (10), то в поверхностном интеграле по Б * в левой части (19) выполняется суммирование по соответствующему подмножеству значений индекса к. Можно показать, что из принципа Эйлера - Лагранжа (19) и произвольности поля виртуальных перемещений ¿и следуют уравнение импульса (1), а также краевые условия (9) и (10). Далее, если выполняются кинематические условия (11), из принципа Эйлера - Лагранжа (19) следует краевое условие (12). Если предположить, что на контактной поверхности Б** отсутствуют силы трения, то исходя из (19), можно показать, что выполнение условия (13) влечет выполнение краевых условий (14)-(16). Если деформируемая среда является упругой либо гиперупругой, то из (3), (5) следует ак¿е^ (V = ¿и(е)(УЬ, и уравнение (19) принимает вид

J к - тк ^ (V - ¿у и (е) (V + ^ рк{0^Пк (Б = 0, к = 1, 2, 3. (20)

V Ус 52

Здесь V0 — объем, занимаемый средой в недеформированном состоянии. При численном решении (1) сплошную среду обычно наделяют конечным числом степеней свободы, например, аппроксимируя поля упругих смещений отрезками фурье-разложений по той или иной полной системе базисных функций Фк3-(г) = Фк3-(ж1,ж2,ж3), удовлетворяющих краевым условиям (9),

N

Пк(г, А) Як3 (*)Фк3. (г), к = 1, 2, 3, г е V). (21)

¿=1

При использовании метода конечных элементов [7-9] индекс 3 в (21) можно сопоставить с номером текущего узла конечно-элементной сетки, а базисные функции Ф^. (г) отличны от нуля лишь в пределах тех конечных элементов, в которые входит текущий узел 3. Особенностью метода конечных элементов является тот факт, что выполнение краевых условий типа (9) сводится к заданию значений обобщенных координат для узлов на граничной поверхности Б1. При этом в (21) функции (А), 3 = N +1,..., N, к = 1, 2, 3, становятся известными функциями времени и тем самым исключаются

из числа искомых функций, а дальнейшему определению подлежат функции (к; (£), з = 1, 2,..., N, к = 1, 2, 3. В качестве базисных функций метода конечных элементов [7-9] обычно используют базисные функции т. н. элементов сирендипова семейства, либо базисные функции т. н. элементов лагранжева семейства (где базисные функции в том или ином смысле представляют собой произведения интерполяционных полиномов Лагранжа локальных координат конечного элемента). Поскольку варьируются лишь искомые обобщенные координаты, из (21) следует

N

5пк (г,*) = 5]Фк. (г )8дк., к = 1, 2, 3, г е Ус.

¿=1

Уравнение Эйлера - Лагранжа (19) должно выполняться при любых значениях вариаций обобщенных координат б(к;. После группировки слагаемых, приравнивания нулю сомножителей при б(к; и преобразования соответствующих уравнений к метрике начального состояния находим

| (бк + V,- ик )уг Фк; (г)^0+

+ / [Р(С) - (бк + VикП] - ^12)/(^11)^22) - (^ )2)Фк; (г)й$(С) +

о(0) Д2

+ 1 Рс (^к - ик )Ф(; (г)^Ус = 0, к = 1, 2, 3, з = 1, . (22)

Здесь величины и операции в метрике актуального состояния отмечены символом Из (21), (2), (4)-(7) следует, что (22) представляет собой записанную в неявной форме систему обыкновенных дифференциальных уравнений второго порядка (весьма большой размерности) относительно набора обобщенных координат , к = 1, 2, 3, з = 1, 2, . Фактически (22) представляет собой некоторый вариант проекционного метода Галеркина [10]. Сходимость приближенных решений (21) к точному решению исходной краевой задачи (в пренебрежении относительными ускорениями частиц сплошной среды) частично исследована в [11].

Элементы матрицы при обобщенных ускорениях (к, к = 1, 2, 3, з = 1, 2, , называемой

обычно матрицей инерции, задаются выражениями рс(г)дкг(г)Ф(; (г)Ф^ (г)йУс, представляющими собой скалярные произведения по начальному объему Ус базисных функций с весовыми коэффициентами рс (г)дь(г). Поскольку весовые коэффициенты образуют положительно определенную симметричную матрицу, матрица инерции будет симметричной положительно определенной. Так как базисные функции отличны от нуля лишь в пределах конечных элементов, содержащих текущий узел, матрица инерции будет разреженной.

Предположим, что сплошная среда является упругой либо гиперупругой, и рассмотрим (20) в пренебрежении (относительным) ускорением частиц сплошной среды. Если материал не является резиноподобным, то в пределах упругих деформаций компоненты тензора деформаций можно считать малыми. Далее, если упруго деформируемое тело не является тонкостенным, то величины упругих смещений также можно считать малыми. Следовательно,

£ы = 2 Vи + ViUk) (23)

и с точностью до малых высшего порядка можно полагать, что положение упруго деформируемого тела в актуальном состоянии практически совпадает с положением в начальном состоянии, т. е. V ~ УС, $ « $(с). Если полагать, что компоненты вектора массовых Е = к} и поверхностных р(С) = (ркС)} сил не зависят от поля смещений и, то (20) принимает вид

бФ = 0, Ф = Ф[и] = У р^кик ¿V - У и(е) ¿V + у рС) ик (24)

V V ^2

т. е. решение уравнений равновесия сплошной среды доставляет экстремальное (а фактически — минимальное) значение функционалу (24), что и составляет содержание известного вариационного

метода Ритца [б]. Подстановка приближенного решения (21) в (23) и далее в (24) позволяет явно выразить величину функционала Ф через значения обобщенных координат q = {qk }, k = 1, 2, 3, j = 1,2,...,N

Ф = Ф(Я), <№ = ¿ f дф^. (25)

k=1 j = 1 dqkj

Из (24), (25) и произвольности вариаций обобщенных координат Sqkj следует, что приближенному решению (21) уравнений равновесия сплошной среды соответствует выполнение требований

дфМ=о, k = 1, 2, 3, j = 1,2, ...,Nk. (2б)

dqkj

Если теперь полагать, что деформируемая среда является линейно-упругой, т. е. справедливы формулы (7), то в этом случае функционал Ритца (24) будет квадратичной функцией обобщенных координат:

N 2 1 NN 3 3

Ф = Bvk qkv - 2 y^A vmkj qkv qjm ,

v=1 k=1 v=1 m=1 k=1 j = 1 (27)

Avmü = J Ckjl VkФг^ ■ Vj Ф^ dV = Amvlt, Bv* = J pFlФг^ dV + J p(o) Фг^ dS.

V V S2

В выражениях для Avm¿l и суммирование по i и l отсутствует. Матрица коэффициентов Avmkj квадратичной формы в правой части (27), называемая обычно матрицей жесткости, очевидно, симметрична. Поскольку квадратичная форма в правой части (27) получается посредством интегрирования по объему V некоторой квадратичной формы, положительно определенной в каждой точке объема V, матрица жесткости также будет положительно определенной. Поскольку базисные функции, типичные для метода конечных элементов, отличны от нуля лишь в пределах конечных элементов, содержащих текущий узел, матрица Avmkj будет разреженной. Решение исходной краевой задачи сводится к решению системы линейных уравнений

3 Ni 3 N

ЕЕ Avmkj j = Bvk £ Avmk, j

l = 1 m=1 l = 1 m=Ni +1

с разреженной симметричной положительно определенной матрицей.

Заметим далее, что методы, рассмотренные в п. 1, автоматически реализуются в большинстве коммерческих [12] и свободно распространяемых [13] пакетов конечно-элементного моделирования.

2. ОПТИМИЗИРОВАННАЯ ВЕРСИЯ КОНЕЧНО-ЭЛЕМЕНТНОГО РЕШАТЕЛЯ

За основу оптимизированной кластерной редакции конечно-элементного решателя был взят исходный код параллельной редакции конечно-элементного решателя (версия 5.5.5), входящего в состав свободно распространяемого пакета конечно-элементного моделирования Elmer [13] (www.csc.fi/english/pages/elmer). Технология MPI [14] (www.mpi-forum.org), основанная на обмене сообщениями между параллельно работающими процессами, является стандартной для параллельной версии пакета Elmer. При сборке оптимизированной версии конечно-элементного решателя были выбраны библиотеки стандартных программ MPICH2 (версия 2-1.2 свободно доступна в сети Интернет, http://www.mcs.anl.gov/research/projects/mpich2/index.php). Сборка оптимизированной кластерной редакции конечно-элементного решателя (Solverá) требует использования обновленного комплекта компиляторов GNU для Linux х8б_б4 (в частности, версии 4.4.3 и выше компиляторов gcc, g++ и gfortran). В связи с этим были исправлены некоторые ошибки в скриптах компиляции и сборки исходного кода пакета Elmer, не позволявшие автоматически отследить изменение используемой версии комплекта компиляторов GNU.

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

фективности выполнения базовых операций по факторизации разреженных матриц. При расчете глобальных конечно-элементных матриц (т. е. при сборке конечно-элементной модели) любой конечно-элементный решатель выполняет большое число базовых операций линейной алгебры (перемножение матриц, умножение матриц на векторы, значительно реже — решение систем линейных уравнений) над плотно заполненными матрицами разного размера (от достаточно малого до достаточно большого). По указанной причине при сборке оптимизированной редакции решателя входящие в состав базовой версии пакета Elmer неоптимизированные версии пакетов LAPACK и т. н. Reference BLAS были заменены на оптимизированную редакцию указанных пакетов в составе библиотеки GotoBLAS2 (http://www.tacc.utexas.edu/resources/software), позволяющие использовать специфические машинные команды процессора Xeon, ресурсы кэш-памяти и т.д.. Для того, чтобы избежать конфликтов технологий MPI и OpenMP в пределах узлов кластера, содержащих по два 4-ядерных процессора Xeon, генерировались последовательные (не многопоточные) библиотеки GotoBLAS2.

Значительные размеры разреженных матриц, возникающих при математическом моделировании на основе метода конечных элементов (характерный размер от нескольких десятков тысяч до нескольких десятков миллионов при среднем количестве ненулевых элементов в строке/столбце от нескольких сотен до нескольких тысяч), требуют использования специальных алгоритмов при решении систем линейных уравнений с разреженными матрицами. В частности, это могут быть итерационные методы решения систем линейных уравнений. Однако относительно медленная сходимость итерационных методов (а в некоторых случаях, например, при наличии сгущения конечно-элементной сетки — отсутствие гарантированной сходимости) вынуждает по возможности использовать т. н. прямые методы решения линейных уравнений, связанные с выполнением факторизации разреженных матриц. Большинство свободно распространяемых, а также доступных для некоммерческого использования пакетов прямых методов численного решения систем линейных уравнений с разреженными матрицами (UMFPACK [15],CHOLMOD [16], MA57 [17], PARDISO [18]) не используют технологию MPI и, следовательно, будут непригодны для вычислительных систем с кластерной архитектурой. По указанной причине кластерные редакции большинства пакетов конечно-элементного моделирования с открытым исходным кодом (например, Elmer) ограничиваются использованием итерационных методов решения систем линейных уравнений с разреженными матрицами, ограничивая использование прямых методов лишь локальными версиями пакетов. К сожалению, это приводит к ощутимому снижению производительности при решении задач конечно-элементного моделирования на кластерных системах. Если говорить применительно к кластерным архитектурам, то одним из немногих свободно распространяемых пакетов, реализующих прямые методы решения линейных уравнений с разреженными матрицами, является пакет MUMPS [19] (MUltifrontal Massively Parallel Solver — «Мультифронтальный масштабируемый параллельный решатель»), исходный код которого после регистрации доступен с сайта производителя (http://mumps.enseeiht.fr либо http://graal.ens-lyon.fr/MUMPS). Пакет выполняет факторизацию достаточно широкого класса разреженных матриц (симметричных и эрмитовых положительно определенных, симметричных неопределенных, несимметричных действительных и комплексных). Исходный код пакета представлен модулями на языках Fortran77/9O и C. Пакет позволяет использовать в пределах отдельных узлов кластеров оптимизированные версии библиотеки BLAS (в частности, упомянутый выше GotoBLAS2). Для сборки пакета MUMPS дополнительно необходимы пакеты BLACS и ScaLAPACK (www.netlib.org), представляющие собой расширения пакетов BLAS и LAPACK на кластерные архитектуры, а также пакет PARMETIS (http://glaros.dtc.umn.edu/gkhome/metis/metis/download), предназначенный для анализа графов и переупорядочения элементов разреженных матриц на кластерных системах, что позволяет достичь требуемой производительности.

Оптимизированная версия конечно-элементного решателя функционирует на типовых кластерах ПРЦНИТ СГУ (содержащих примерно по 10 узлов, с двумя 4-ядерными процессорами Intel Xeon и 16 ГБ оперативной памяти в каждом узле). При ее создании к базовой поставке исходного кода пакета Elmer были добавлены пакеты MPICH2, GotoBLAS2, BLACS, SCALAPACK, PARMETIS и MUMPS, а также было выполнено некоторое исправление скриптов компиляции и сборки, позволяющее использовать обновленные версии комплекта компиляторов GNU (gcc/g++/gfortran). В результате удалось сократить характерное время выполняемого на кластере СГУ в ходе проведения виртуальной опе-

рации остеосинтеза типового этапа конечно-элементных расчетов с нескольких десятков минут до нескольких десятков секунд, что вполне приемлемо для клинической практики. Характерный пример моделирования (поле эквивалентных напряжений, МПа) в системе «бедренная кость - фиксатор» после отображения пгш¥ченных пянных и ппотгтпттеооппе ппикепен ня пио. 3

Vonmises 4.54е+006

Рис. 3. Поле эквивалентных напряжений

Библиографический список

1. Мюллер, М.Е. Руководство по внутреннему остео-синтезу / М.Е. Мюллер, М. Алльговер, Р. Шнайдер, Х. Виллингер. - М.: Springer-Verlag, 1996. - 750 с.

2. Соловьёв, В.М. Технология построения твердотельных моделей бедренных костей на основе данных компьютерной томографии / В.М. Соловьёв, П.В. Ирматов, М.С. Ирматова, М.Г. Щербаков // Изв. Сарат. ун-та. Нов. сер. - 2010. - Т. 10. - Сер. Математика. Механика. Информатика, вып. 2. - С. 81-87.

3. Седов, Л.И. Механика сплошной среды: в 2 т. / Л.И. Седов. - М.: Наука, 1976. - Т. 1. - 536 с.

4. Новиков, С.П. Элементы дифференциальной геометрии и топологии / С.П. Новиков, А.Т. Фоменко. - М.: Наука, 1987. - 432 с.

5. Четаев, Н.Г. Теоретическая механика / Н.Г. Четаев. - М.: Наука, 1987. - 368 с.

6. .Лурье, А.И. Теория упругости / А.И. Лурье. - М.: Наука, 1970. - 939 с.

7. Галлагер, Р. Метод конечных элементов. Основы / Р. Галлагер. - М.: Мир,1984. - 428 с.

8. Зенкевич, О. Метод конечных элементов в технике / О. Зенкевич. - М.: Мир, 1975. - 541 с.

9. Zienkiewich, O.C. The Finite Element Method for Solid and Structural Mechanics / O.C. Zienkiewich, R.L. Taylor. - Elsevier, 2005. - 631 p.

10. Флетчер, К. Численные методы на основе метода Галеркина / К. Флетчер. - М.: Мир, 1988. - 352 с.

11. Оден, Дж. Конечные элементы в нелинейной механике сплошных сред / Дж. Оден. - М.: Мир, 1976. -464 с.

12. ANSYS [Электронный ресурс]. URL: http://www. ansys.com (дата обращения: 20.04.2010).

13. Elmer Solver Manual /CSC - IT Center for Science.-Finland, 2010. — 119 p. [Электронный ресурс] URL: http://www.nic.funet.fi/pub/sci/physics/elmer/doc/Elmer SolverManual.pdf. (дата обращения: 20.04.2010).

14. Антонов, А.С. Параллельное программирование с использованием технологии MPI (www.parallel.ru) / А.С. Антонов. - М.: Изд-во Моск. ун-та, 2004. - 72 с.

15. Davis, T.A. UMFPACK v.5.4.0 User Guide. Dept. of Computer and InformationScience and Engineering Univ. of Florida / T.A. Davis. - Gainesville, Florida, 2009. — 133 p. [Электронный ресурс] URL: http:// www.cise.ufl.edu/research/sparse/umfpack/current/UMF PACK/Doc/UserGuide.pdf (дата обращения: 20.04.2010).

16. Davis, T.A. User Guide for CHOLMOD: a sparse Cholesky factorization and modification package. Ver. 1.7. Dept. of Computer and Information Science and Engineering Univ. of Florida / T.A. Davis. - Gainesville, Florida, 2008. - 140 p. [Электронный ресурс] URL: http:// www.cise.ufl.edu/research/sparse/cholmod/CHOLMOD/ Doc/UserGuide.pdf (дата обращения: 20.04.2010).

17. The HSL Mathematical Software Library [Электронный ресурс]. URL: http://hsl.rl.ac.uk (дата обращения: 20.04. 2010).

18. Intel® Math Kernel Library. Reference Manual. Document Number: 630813-031US. - 2009. - 3680 p. [Электронный ресурс]. URL: http://software.intel.com/ sites/products/documentation/hpc/mkl/mklman.pdf (дата обращения: 20.04.2010).

19. MUltifrontal Massively Parallel Solver (MUMPS 4.9.2) User Guide. - 2009. - 46 p. [Электронный ресурс]. URL: http://graal.ens-lyon.fr/MUMPS/doc/ userguide_4.9.2.pdf. (дата обращения: 20.04.2010).