Решите дифференциальное уравнение Используя многосеточный предварительный формирователь на распределенной дискретизации

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

Этот пример продолжает темы, затронутые в Использовании Распределенные Массивы, чтобы Решить Системы Линейных уравнений с Итерационными методами. На основе [1], модели в качестве примера нагревают распределение в комнате при помощи уравнения Пуассона в форме, известной как гомогенное установившееся уравнение тепла. Устойчивое состояние означает, что тепло не меняется в зависимости от времени, и гомогенный означает, что нет никакого внешнего источника тепла.

$- Δ u = - (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + \frac{\partial^{2} u}{\partial z^{2}}) = 0$

В уравнении, $u$ представляет температуру в каждой точке $(x, y, z)$ из комнаты. Чтобы решить уравнение, вы сначала аппроксимируете его системой линейных уравнений с помощью метода дискретизации конечной разности. Затем вы используете предобусловленные методы сопряженных градиентов (pcg) метод, чтобы решить систему. Предварительное создание условий преобразовывает проблему улучшать производительность числового решателя. При помощи распределенных массивов можно усилить объединенную память о кластере машин и позволить более прекрасные дискретизации.

Узнать, как к:

Настройте дискретную 3-D сетку и граничные условия.
Задайте дискретизацию и многосеточный предварительный формирователь.
Примените предобусловленный числовой решатель, чтобы решить уравнение тепла через 3-D объем.

Дискретизируйте пространственные размерности

В этом примере, кубе стороны 1 моделирует комнату. Первый шаг должен дискретизировать его с помощью 3-D сетки.

Метод перед созданием условий в этом примере использует несколько сеток с разными уровнями гранулярности. Каждый уровень огрубляет сетку фактором 2 в каждой размерности. Задайте количество многосеточных уровней.

multigridLevels = 2;

Задайте число точек в каждой размерности, XY, и Z, из самой прекрасной сетки. Метод перед созданием условий требует, чтобы число точек в этой сетке было делимым 2^multigridLevels. В этом случае число точек должно быть делимым 4, потому что количеством многосеточных уровней является 2.

numPoints.X = 32;
numPoints.Y = 32;
numPoints.Z = 32;

Дискретизируйте пространственные размерности с 3-D сеткой при помощи meshgrid функция. Разделите каждую размерность однородно согласно числу точек при помощи linspace. Обратите внимание на то, что, чтобы включать контуры куба, необходимо добавить две дополнительных точки.

[X,Y,Z] = meshgrid(linspace(0,1,numPoints.X+2), ...
    linspace(0,1,numPoints.Y+2), ...
    linspace(0,1,numPoints.Z+2));

Задайте граничные условия

Предположим, что комната имеет окно и дверь. Стены и потолок имеют постоянную температуру 0 градусов, окно имеет постоянную температуру 16 градусов, и дверь имеет постоянную температуру 15 градусов. Пол в 0,5 градусах. Цель состоит в том, чтобы определить температурное распределение через внутреннюю часть комнаты.

Задайте координаты пола, окна и двери с помощью операторов отношения, и задайте температуру на этих граничных элементах. Контуры являются фасетами куба и, поэтому, один из XY, или Z должен быть 0 или 1. Установите остальную часть контура и внутренней части куба к 0.

floor  = (0.0 <= X & X <= 1.0) & (0.0 <= Y & Y <= 1)   & (Z == 0.0);
window = (X == 1)              & (0.2 <= Y & Y <= 0.8) & (0.4 <= Z & Z <= 0.6);
door   = (0.4 <= X & X <= 0.6) & (Y == 1.0)            & (0.0 <= Z & Z <= 0.6);

u = zeros(size(X));
u(floor)  = 0.5;
u(window) = 16;
u(door)   = 15;

Эти граничные условия задают постоянные значения, что решение должно взять с собой контур области. Этот тип граничного условия известен как граничное условие Дирихле.

Визуализируйте граничные условия с помощью slice функция. Используйте срезы, расположенные за пределами куба, которые показывают ненулевые граничные условия.

xSlices = 1;
ySlices = 1;
zSlices = 0;
f = slice(X,Y,Z,u,xSlices,ySlices,zSlices,'nearest');
title('Constant nonzero boundary conditions'), xlabel('x'), ylabel('y'), zlabel('z');
colorbar, colormap cool;
shading interp;
set(f,'EdgeColor',[0 0 0]);

Дискретизируйте и решите дифференциальное уравнение

Этот пример дискретизирует дифференциальное уравнение в линейную систему с помощью метода приближения конечных разностей и использует многосеточный предварительный формирователь, чтобы улучшать производительность итеративного решателя. В данном примере используйте дискретизацию и предварительный формирователь в функциях поддержки discretizePoissonEquation и multigridPreconditioner. Для других проблем выберите дискретизацию и предварительный формирователь, которые подходят для вашего приложения.

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

Среди других методов этот предварительный формирователь применяет сглаживание, чтобы минимизировать ошибки с рядом приближений. Задайте количество сглаживания шагов. Используя большее количество шагов делает приближения более точными, но также и более в вычислительном отношении интенсивными. Затем дискретизируйте дифференциальное уравнение и настройте предварительный формирователь.

numberOfSmootherSteps = 1;
[A,b,multigridData] = discretizePoissonEquation(numPoints,multigridLevels,numberOfSmootherSteps,u);

Level 0: The problem is of dimension 32768 with 223232 nonzeros.
Level 1: The problem is of dimension 4096 with 27136 nonzeros.
Level 2: The problem is of dimension 512 with 3200 nonzeros.

preconditioner = setupPreconditioner(multigridData);

Решите линейную систему с помощью предобусловленных методов сопряженных градиентов.

tol = 1e-12;
maxit = numel(A);
pcg(A,b,tol,maxit,preconditioner);

pcg converged at iteration 45 to a solution with relative residual 5.4e-13.

Масштабируйте с распределенными массивами

Если вам нужно больше вычислительных ресурсов, таких как память, можно масштабировать, использование распределило массивы, не будучи должен изменить код. Распределенные массивы распределяют ваши данные на нескольких рабочих, и они могут усилить вычислительную производительность и память о кластере машин.

Запустите пул параллельных рабочих. По умолчанию, parpool использует ваш кластер по умолчанию. Проверяйте, что ваш кластерный профиль по умолчанию на вкладке MATLAB Home, в области Environment, параллельно> Выбирает Default Cluster.

parpool;

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 12).

Распределите температурную переменную u через память о рабочих в вашем кластере при помощи distributed функция.

distU = distributed(u);

Можно использовать тот же код как прежде; никакие изменения не требуются, потому что дискретизация и функции перед формирователем создают распределенные массивы, если вход является распределенным массивом. Много функций MATLAB улучшены для распределенных массивов, таким образом, можно работать с ними таким же образом, вы работаете с массивами в оперативной памяти.

Обратите внимание на то, что discretizePoissonEquation возвращает структуру, содержащую распределенные данные. Чтобы использовать распределенные данные в структуре распределенным способом, необходимо создать структуру в spmd блок. Необходимо также вызвать любую функцию, которая использует его в spmd блок.

spmd
    [A,b,multigridData] = discretizePoissonEquation(numPoints,multigridLevels,numberOfSmootherSteps,distU);
    preconditioner = setupPreconditioner(multigridData);
end

Lab  1: 
  Level 0: The problem is of dimension 32768 with 223232 nonzeros.
  Level 1: The problem is of dimension 4096 with 27136 nonzeros.
  Level 2: The problem is of dimension 512 with 3200 nonzeros.

Используйте pcg в spmd блокируйтесь, чтобы решить линейную систему распределенным способом.

spmd
    x = pcg(A,b,tol,maxit,preconditioner);
end

Lab  1: 
  pcg converged at iteration 45 to a solution with relative residual 5.4e-13.

Постройте результаты

Решение от решателя является вектором, который умещается в памяти. Отправьте данные от рабочих клиенту при помощи gather. Измените данные назад в трехмерный массив и переупорядочьте размерности, чтобы произвести конечное решение. Установите внутреннюю часть u к этому решению. Внешняя часть, контур, уже содержит значение граничных условий.

x3D = reshape(gather(x),numPoints.X,numPoints.Y,numPoints.Z);
u(2:end-1,2:end-1,2:end-1) = permute(x3D, [2, 1, 3]);

Визуализируйте решение с помощью slice функция. Добавьте дополнительные срезы, чтобы построить температуру в кубе. Можно использовать Rotate инструмент, или варьируются положение срезов, чтобы исследовать решение.

xSlices = [.5,1];
ySlices = [.5,1];
zSlices = [0,.5];
f = slice(X,Y,Z,u,xSlices,ySlices,zSlices,'nearest');
title('Heat distribution'), xlabel('x'), ylabel('y'), zlabel('z');
colorbar, colormap cool;
shading interp;

Можно попробовать различные значения numPoints в этом примере, чтобы протестировать разные уровни дискретизации. Используя большие повышения стоимости разрешение, но требует большей памяти. Кроме того, больший multigridLevels больше памяти, эффективной, предварительный формирователь. Однако больший multigridLevels подразумевает менее точный предварительный формирователь, поскольку огрубление уменьшает точность на каждом уровне. В результате решателю может быть нужно больше итераций, чтобы достигнуть того же уровня точности.

Задайте предварительный формирователь

Задайте многосеточный предварительный формирователь для использования с предобусловленным методом сопряженных градиентов. Этот тип предварительного формирователя использует несколько сеток дискретизации с разными уровнями гранулярности, чтобы аппроксимировать решение системы линейных уравнений более эффективно. Метод перед созданием условий в этом примере основан [2] и следует за этими основными этапами:

Предсглаженное использование метода приближения Зейделя Гаусса.
Вычислите остаточное решение на более грубом уровне.
Рекурсивно предусловие на более грубом уровне, или решают непосредственно если на самом грубом уровне.
Обновите решение с более грубым решением для сетки.
Постсглаженное использование метода приближения Зейделя Гаусса.

function x = multigridPreconditioner(mgData,r,level)

if(level < mgData(level).MaxLevel)
    x = zeros(size(r),'like',r);
    
    % Presmooth using Gauss-Seidel
    for i=1:mgData(level).NumberOfSmootherSteps
        x = mgData(level).Matrices.LD \ (-mgData(level).Matrices.U*x + r);
        x = mgData(level).Matrices.DU \ (-mgData(level).Matrices.L*x + r);
    end
    
    % Compute residual on a coarser level
    Axf = mgData(level).Matrices.A*x;
    rc = r(mgData(level).Fine2CoarseOperator)- Axf(mgData(level).Fine2CoarseOperator);
    
    % Recursive call until coarsest level is reached
    xc = multigridPreconditioner(mgData, rc, level+1);
    
    % Update solution with prolonged coarse grid solution
    x(mgData(level).Fine2CoarseOperator) = x(mgData(level).Fine2CoarseOperator)+xc;
    
    % Postsmooth using Gauss-Seidel
    for i = 1:mgData(level).NumberOfSmootherSteps
        x = mgData(level).Matrices.LD \ (-mgData(level).Matrices.U*x + r);
        x = mgData(level).Matrices.DU \ (-mgData(level).Matrices.L*x + r);
    end
else
    % Obtain exact solution on the coarsest level
    x = mgData(level).Matrices.A \ r;
end

end

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

function preconditioner = setupPreconditioner(multigridData)

if ~isempty(multigridData)
    preconditioner = @(x,varargin) multigridPreconditioner(multigridData,x,1);
else
    preconditioner = [];
end

end

Ссылки

[1] Dongarra, J., М. А. Херукс и П. Ласзкзек. "Сравнительный тест HPCG: новая метрика для рейтинга высокоэффективных вычислительных систем". Ноксвилл, TN: университет Теннесси, 2015.

[2] Элмен, H. C. Д. Дж. Сильвестер и А. Дж. Уотэн. Конечные элементы и быстро итеративные решатели: с приложениями в несжимаемой гидрогазодинамике. Оксфорд, Великобритания: Издательство Оксфордского университета, 2005, разделяет 2.5.

Документация