How can I speed up the multiplication of matrices?

3 ビュー (過去 30 日間)
Jose Salazar
Jose Salazar 2018 年 4 月 24 日
コメント済み: Walter Roberson 2018 年 4 月 30 日
Hello all,
I am having problems trying to improve the speed of my code. I wrote a code to solve the pressure and water saturation of a reservoir. I have to iterate up to 500 days (with increments of 1 day), solving the pressure implicitly and the water saturation explicitly. See the code please:
clc; clear;
tic
warning('off', 'MATLAB:singularMatrix');
%%P R O P E R T I E S & C O M P U T A T I O N O F M A T R I C E S
% Calls the function inputvalueSN for the input files
[reservoir, fluid, numerical, well, pressure, rperm, Sw] = inputProject2_2;
% Computes the matrices T, B, Q, G
Sw_hyst = zeros(size(reservoir.permx));
[T, Tw, d12, d11, D, J, Q, Qwprod, G, gamma_w] = TDJQG2(reservoir, fluid, numerical, well, pressure, rperm, Sw, Sw_hyst);
%%N U M E R I C A L S O L U T I O N S
% BEFORE 500 DAYS
P_old = sparse(pressure.Po);
P_new = sparse(zeros(size(pressure.Po)));
Po_total = zeros(length(pressure.Po), numerical.final_t/numerical.dt + 1);
Po_total(:,1) = pressure.Po;
Sw_total = zeros(length(Sw), numerical.final_t / numerical.dt + 1);
Sw_total(:,1) = Sw;
t = 1;
for i = 0:numerical.dt:15 %numerical.switch - 1
P_new = (T + J + D) \ (D * P_old + Q + G); % It computes the oil pressure
Sw_new = (Sw_total(:,t) + d12 \(-d11 * (P_new - P_old) + Qwprod...
- Tw * (P_new - gamma_w .* reservoir.depth - pressure.Pc)));
indexzero = (reservoir.depth == 0);
Sw_new(indexzero) = 0; P_new(indexzero) = 0;
[pressure, Sw_hyst2] = Pr_hyst(Sw_total(:,t), Sw_new, rperm, P_new, pressure, Sw_hyst); % Update capillary pressure and Pw
Sw_hyst = Sw_hyst2;
% Updating the matrices and vector dependent of pressure and saturation
[T, Tw, d12, d11, D, J, Q, Qwprod, G, gamma_w] = TDJQG2(reservoir, fluid, numerical, well, pressure, rperm, Sw_new, Sw_hyst);
P_old = P_new;
Po_total(:, t+1) = P_new;
Sw_total(:, t+1) = Sw_new;
t = t + 1;
end
However, the matrices T, J, D, Q and G are pressure dependent, i.e., I have to update them every time step within the loop. The function that computes them is: TDJQG2. The script is the following:
function [T, Tw, d12, d11, D, J, Q, Qwprod, G, gamma_w] = TDJQG2(reservoir, fluid, numerical, well, pressure, rperm, Swater, Sw_hyst)
% This function computes the transmisibility matrix T, accumulation matrix
% B, source vector Q, G and J
%%FVF of oil and water
fluid.Bo = fluid.Bo_sat * (1 - fluid.co .* (pressure.Po - fluid.Pb));
fluid.Bw = fluid.Bw_ref * (1 - fluid.cw .* (pressure.Pw - pressure.ref_Bw));
% fluid.Bw = fluid.Bw_ref * (1 - fluid.cw .* (pressure.Pw - pressure.ref - pressure.Pe));
%%Update densities of oleic phase
oleic_dens = (fluid.density_o + fluid.Rs * (0.0458171 /5.615)) ./ fluid.Bo;
%%TDJQG
Tw = zeros(6000);
To = zeros(6000);
D = zeros(6000);
d11 = zeros(6000);
d12 = zeros(6000);
d21 = zeros(6000);
d22 = zeros(6000);
Jw = zeros(6000);
Jo = zeros(6000);
Qw = zeros(6000, 1);
Qo = zeros(6000, 1);
Qwprod = zeros(6000, 1);
for L = 1:6000
if rem(L, numerical.Nx) ~= 1 % If the block is NOT on left edge
[tw, to] = inter4(L, L-1, reservoir, fluid, numerical, pressure, rperm, Swater, oleic_dens); % Computation of transmissibilities
% Water Transmissibilities on the left
Tw(L, L-1) = -tw; Tw(L, L) = Tw(L, L) - Tw(L, L-1);
% Oil Transmissibilities on the left
To(L, L-1) = -to; To(L, L) = To(L, L) - To(L, L-1);
end
%%%%%%%%MORE CODE THAT WORKS %%%%%%%%%%%%%%
end
d12 = sparse(d12); d11 = sparse(d11); d22 = sparse(d22); D = sparse(D);
C = isnan(D);
D(C) = 0;
Tw = sparse(Tw * 6.33e-3); To = sparse(To * 6.33e-3);
T = ((-d22 * d12^-1) * Tw + To) * 6.33e-3;
Qw = sparse(Qw); Qo = sparse(Qo);
Q = (-d22 * d12^-1) * Qw + Qo;
Qwprod = sparse(Qwprod);
Jw = sparse(Jw); Jo = sparse(Jo);
J = ((-d22 * d12^-1) * Jw + Jo) * 6.33e-3;
gamma_o = ((fluid.density_o + (fluid.Rs * fluid.density_g * (1/5.615))) ./ fluid.Bo)./144;
gamma_o = sparse(gamma_o);
gamma_w = (fluid.density_w ./ fluid.Bw) ./144;
gamma_w = sparse(gamma_w);
prePc = sparse(pressure.Pc);
depths = sparse(reservoir.depth);
G = ((-d22 * d12^-1) * (Tw * prePc) + diag(-d22 / d12) .*...
(Tw * (gamma_w .* reservoir.depth))+ (To * (gamma_o .* depths)));
G = sparse(G);
end
My main issue is that the following lines consume 57% of all the computation time. (I clic "Run and Time" for 15 days).
G = ((-d22 * d12^-1) * (Tw * prePc) + diag(-d22 / d12) .*... (Tw * (gamma_w .* reservoir.depth))+ (To * (gamma_o .* depths)));
d12 = sparse(d12); d11 = sparse(d11); d22 = sparse(d22); D = sparse(D);
T = ((-d22 * d12^-1) * Tw + To) * 6.33e-3;
Q = (-d22 * d12^-1) * Qw + Qo;
J = ((-d22 * d12^-1) * Jw + Jo) * 6.33e-3;
Is there any way to improve the computation of those matrices mentioned above? d22, d12, d11 are diagonal matrices, Q only has 10 entries, and J only 10 entries. G is a column vector. I already tried the backslash \ operator to compute those value but I get the message that the matrices are baldy scalded. The solutions are correct, but the is poor in performance.
Thank you in advance.
PD: these are the speed performance results:
  1 件のコメント
Walter Roberson
Walter Roberson 2018 年 4 月 30 日
Please do not close questions that have an answer.

サインインしてコメントする。

回答 (1 件)

James Tursa
James Tursa 2018 年 4 月 24 日
Q1: Why do you pre-allocate a bunch of 6000x6000 full matrices for d11, d12, and d22 only to downstream use sparse( ) on them and tell us they are diagonal? Why not just create the diagonals to begin with and work with that?
Q2: Are you sure things stay diagonal throughout your calculations, or do small off-diagonal elements creep in? This might make a big difference in how the linear algebra calculations are done for the sparse matrices.
Q3: Why do you do the (-d22 * d12^-1) calculation repeatedly? Why not do it only once?
Q4: Are any of your large calculations involving a mix of sparse and full matrices? This could cause calculations to be done in full and slow things down (because the 0's are explicitly multiplied).
  1 件のコメント
Jose Salazar
Jose Salazar 2018 年 4 月 25 日
Hello James,
Q1: I am now working with the diagonals only. I am using spdiags.
Q2: There are some values where the main diagonal is just full of zeros. Like 20%.
Q3: I have not thought about! It is indeed helpful.
Q4: Should I multiply everything in sparse?
Thank you!

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeAxis Labels についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by