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

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

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

multigridLevels = 2;

Определите число точек в каждой размерности, X, Y, и Z, из лучших сетки. Метод предварительного обусловления требует, чтобы число точек в этой сетке было делено на 2^multigridLevels. В этом случае число точек должно быть делено на 4, потому что количество уровней multigrid 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 степени. Цель состоит в том, чтобы определить распределение температуры по внутренней части помещения.

Определите координаты пола, окна и двери с помощью реляционных операторов и определите температуру на этих краевых элементах. Контуры являются гранями куба и, следовательно, одной из X, Y, или 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, в области окружение, в Parallel > Select a Default Cluster.

parpool;

Starting parallel pool (parpool) using the 'MyCluster' 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

Analyzing and transferring files to the workers ...done.
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] Донгарра, Дж., М. А. Херу и П. Лющек. HPCG Benchmark: A New Metric for Ranking High Performance Computing Systems (неопр.) (недоступная ссылка). Ноксвилл, TN: Университет Теннесси, 2015.

[2] Elman, H. C., D. J. Silvester, and A. J. Wathen. Конечные элементы и быстрые итеративные решатели: с приложениями в несжимаемой динамике жидкости. Оксфорд, Великобритания: Oxford University Press, 2005, Section 2.5.