Pre-allocating cell containers

3 visualizaciones (últimos 30 días)
Chukwuemeka Enemanna
Chukwuemeka Enemanna el 7 de Abr. de 2021
Comentada: per isakson el 14 de Abr. de 2021
I am setting up a dual porosty model with a well defined control volume. I will be using influx and recovery from this control volume to compare performance of different transfer functions. In calculating influx/recovery (oil_influx/recovery_Top, Left, Bottom, Right) from this control volume, I concatenated flux from all the control volume boundary faces. Implying influx/recovery continually increase after each iteration loop.
My challenge
  1. It takes too much time to run - over 24 hours and
  2. I get this warning message - The variable appears to change in every time loop iteration. Consider pre-allocating for speed
  3. I have pre-allocated memory for the these influx/recovery (oil_influx/recovery_Top, Left, Bottom, Right) variables but still having same challenges as in 1 and 2
Please help. Any tips will be appreciated.
%% Two-phase water injection in an oil-saturated dual-porosity reservoir.
% Water flows quickly in the fracture system, while transfer to the matrix
% happens via imbibition.
mrstModule add ad-blackoil ad-core ad-props dual-porosity
%% Set up grid
G = cartGrid([100, 1, 100], [10, 1, 10]*meter);
G = computeGeometry(G);
%% Set up rock properties
rock_matrix = makeRock(G, 10*milli*darcy, 0.2);
rock_fracture = makeRock(G, 100*milli*darcy, 0.1);
%% Define Model Boundary faces
eps = 0.01;
Ztop_faces2 = find(G.faces.centroids(:, 3)> (0.0 - eps) &...
G.faces.centroids(:, 3)< (0.0 + eps) &...
G.faces.centroids(:, 1)> (0) &...
G.faces.centroids(:, 1)< 10);
Zbtm_faces2 = find(G.faces.centroids(:, 3)> (10 - eps) & ...
G.faces.centroids(:, 3)< (10 + eps)&...
G.faces.centroids(:, 1)> (0) &...
G.faces.centroids(:, 1)< 10);
XLHS_faces2 = find(G.faces.centroids(:, 1)> (0 - eps) &...
G.faces.centroids(:, 1)< (0 + eps) &...
G.faces.centroids(:, 3)> (0) &...
G.faces.centroids(:, 3)< 10);
XRHS_faces2 = find(G.faces.centroids(:, 1)> (10 - eps) &...
G.faces.centroids(:, 1)< (10 + eps)&...
G.faces.centroids(:, 3)> (0) &...
G.faces.centroids(:, 3)< 10);
%% C. Define Control Volume and boudaries faces for flux calculation in the control volume
CntrV = find(G.cells.centroids(:, 1)> 4.00-eps &...
G.cells.centroids(:, 1)< 6.00+eps &...
G.cells.centroids(:, 3)> 4.00-eps &...
G.cells.centroids(:, 3)< 6.00+eps);
CntrV_Zface_top = find(G.faces.centroids(:, 1)> (4.00-eps) &...
G.faces.centroids(:, 1)< (6.00+eps) &...
G.faces.centroids(:, 3)> (4.00-eps) &...
G.faces.centroids(:, 3)< (4.00+eps));
CntrV_Zface_btm = find(G.faces.centroids(:, 1)> (4.00-eps) &...
G.faces.centroids(:, 1)< (6.00+eps) &...
G.faces.centroids(:, 3)> (6.00-eps) &...
G.faces.centroids(:, 3)< (6.00+eps));
CntrV_Xface_LHS= find(G.faces.centroids(:, 3)> (4.00-eps) &...
G.faces.centroids(:, 3)< (6.00+eps) &...
G.faces.centroids(:, 1)> (4.00-eps) &...
G.faces.centroids(:, 1)< (4.00+eps));
CntrV_Xface_RHS= find(G.faces.centroids(:, 3)> (4.00-eps) &...
G.faces.centroids(:, 3)< (6.00+eps) &...
G.faces.centroids(:, 1)> (6.00-eps) &...
G.faces.centroids(:, 1)< (6.00+eps));
CtrV_faces = [CntrV_Zface_top CntrV_Zface_btm CntrV_Xface_LHS CntrV_Xface_RHS];
%% Plot grid
figure(3);
plotGrid(G, 'FaceColor', 'none', 'FaceAlpha', 0);
view([0.143826069396107 0.170936586318325]); hold on
plotGrid(G, CntrV, 'FaceColor', 'b'); hold on
xlabel('X');
ylabel('Y');
zlabel('Z');
%% Set up fluid structure and Mobility fields
fluid_matrix = initSimpleADIFluid('phases', 'WO', ...
'rho', [1000, 800]*kilogram/meter^3,...
'c', [1e-20, 1e-20]/barsa);
fluid_fracture = initSimpleADIFluid('phases', 'WO', ...
'rho', [1000, 800]*kilogram/meter^3,...
'c', [1e-20, 1e-20]/barsa);
Pe = 0.15*barsa;
fluid_matrix.pcOWm = @(swm)Pe*swm;
fluid_matrix.krWm_h = @(swm_av)swm_av;
fluid_matrix.krOm_h = @(swm_av) (1-swm_av);
%% Set the DP model.
% Here, a two-phase model is used. Rock and fluid are sent in the constructor as a cell array {fracture,matrix}
model = TwoPhaseOilWaterModelDP(G, {rock_fracture,rock_matrix},{fluid_fracture,fluid_matrix});
%% Initializing state
state = initResSol(G, 2.5*psia, [eps, (1-eps)]);
state.wellSol = initWellSolAD([], model, state);
state.pressure_matrix = state.pressure;
state.sm = state.s;
%% Setting transfer function.
% model.transfer_model_object = KZ_TF('KZ_SF', [5,5,5]);
% model.transfer_model_object = GK_TF('GK_SF',[5,5,5]);
% model.transfer_model_object = QS_TF('QS_SF',[5,5,5]);
% model.transfer_model_object = G_TF('G_SF',[5,5,5]);
% model.transfer_model_object = G_TF_NoVermeulen('G_SF',[5,5,5]);
model.transfer_model_object = G_TF_TrblShoot('G_SF',[5,5,5]);
%% 5. Pore Volume calculation & Boundary conditions
pv_matrix = sum(model.operators.pv_matrix);
pv_fracture = sum(model.operators.pv);
tpv = pv_matrix + pv_fracture;
tmax =1000*day;
injRate =100* tpv/tmax;
bc = fluxside([], G, 'ZMin', injRate, 'sat', [1,0]);
bc = pside(bc, G, 'ZMax', 5*psia, 'sat', [0,1]);
%% Solver & Model Validation
solver = NonLinearSolver();
model = model.validateModel();
n = 500; % Total number of timesteps
dt = (tmax/n); % Unit timestep
states = cell(n+1, 1); % Pre-allocate cell containers for states
states{1} = state; % define states{1} as initail state
solver.verbose = true;
%% 8. time loop
for i = 1:n
state = solver.solveTimestep(states{i}, dt, model, 'bc', bc);
states{i+1} = state;
end
%% 9B. Filter out control volume boundary face fluxes
Ztopflux = states{i}.flux(CntrV_Zface_top,:);
Zbtmflux = states{i}.flux(CntrV_Zface_btm,:);
XLHSflux = states{i}.flux(CntrV_Xface_LHS,:);
XRHSflux = states{i}.flux(CntrV_Xface_RHS,:);
%% 9C. Preallocate memory and Initialize phase fluxes - influx and recovery
% CntrV Top
oil_influx_T = cell((n+1), 1);
oil_influx_T{1}= 0;
oil_recovery_T = cell((n+1),1);
oil_recovery_T{1}= 0;
water_influx_T = cell((n+1),1);
water_influx_T{1}= 0;
water_recovery_T =cell((n+1),1);
water_recovery_T{1}= 0;
%% 9D. Boundary Fluxes and flux sign interpretation
% i. Top Flux
for i=2:n+1
CntrV_Ztopflux = states{i}.flux(CntrV_Zface_top,:);
oil_influx_T_i = 0;
oil_recovery_T_i = 0;
water_influx_T_i = 0;
water_recovery_T_i = 0;
for j=1:length(CntrV_Zface_top)
% Oil Flux
if CntrV_Ztopflux(j,1) > 0
oil_influx_T_i = oil_influx_T_i + CntrV_Ztopflux(j,1)*dt;
else
oil_recovery_T_i = oil_recovery_T_i + abs(CntrV_Ztopflux(j,1))*dt;
end
% Gas Flux
if CntrV_Ztopflux(j,2) > 0
water_influx_T_i = water_influx_T_i + CntrV_Ztopflux(j,2)*dt;
else
water_recovery_T_i = water_recovery_T_i + abs(CntrV_Ztopflux(j,2))*dt;
end
end
oil_influx_T = [oil_influx_T oil_influx_T(end) + oil_influx_T_i];
oil_recovery_T = [oil_recovery_T oil_recovery_T(end) + oil_recovery_T_i];
water_influx_T = [water_influx_T water_influx_T(end) + water_influx_T_i];
water_recovery_T = [water_recovery_T water_recovery_T(end) + water_recovery_T_i];
end
%% LHS Boundary Flux
% CntrV LHS - Preallocate memory and Initialize phase fluxes - influx and recovery
oil_influx_L = cell((L+1),1);
oil_influx_L{1} = 0;
oil_recovery_L = cell((L+1),1);
oil_recovery_L{1} = 0;
water_influx_L = cell((L+1),1);
water_influx_L{1} = 0;
water_recovery_L = cell((L+1),1);
water_recovery_L{1} = 0;
for i=2:n+1 % iterate for every time step 1- n where n = 500
CntrV_XLHSflux = states{i}.flux(CntrV_Xface_LHS,:); % All CntrV XLHS fluxes for time step i
oil_influx_L_i = 0;
oil_recovery_L_i = 0;
water_influx_L_i = 0;
water_recovery_L_i = 0;
for j=1:length(CntrV_Xface_LHS)
% Oil Flux
if CntrV_XLHSflux(j,1) > 0
oil_influx_L_i = oil_influx_L_i + CntrV_XLHSflux(j,1)*dt;
else
oil_recovery_L_i = oil_recovery_L_i + abs(CntrV_XLHSflux(j,1))*dt;
end
% Gas Flux
if CntrV_XLHSflux(j,2) > 0
water_influx_L_i = water_influx_L_i + CntrV_XLHSflux(j,2)*dt;
else
water_recovery_L_i = water_recovery_L_i + abs(CntrV_XLHSflux(j,2))*dt;
end
end
oil_influx_L = [oil_influx_L oil_influx_L(end) + oil_influx_L_i];
oil_recovery_L = [oil_recovery_L oil_recovery_L(end) + oil_recovery_L_i];
water_influx_L = [water_influx_L water_influx_L(end) + water_influx_L_i];
water_recovery_L = [water_recovery_L water_recovery_L(end) + water_recovery_L_i];
end
%% Btm Boundary Flux
% CntrV Btm- Preallocate memory and Initialize phase fluxes - influx and recovery
oil_influx_B = cell((L+1),1);
oil_influx_B{1} = 0;
oil_recovery_B = cell((L+1),1);
oil_recovery_B{1} = 0;
water_influx_B = cell((L+1),1);
water_influx_B{1} = 0;
water_recovery_B = cell((L+1),1);
water_recovery_B{1} = 0;
for i=2:n+1 % iterate for every time step 1- n where n = 1000
CntrV_Zbtmflux = states{i}.flux(CntrV_Zface_btm,:); % All CntrV Ztop fluxes for time step i
oil_influx_B_i = 0;
oil_recovery_B_i = 0;
water_influx_B_i = 0;
water_recovery_B_i = 0;
for j=1:length(CntrV_Zface_btm)
% Oil Flux
if CntrV_Zbtmflux(j,1) > 0
oil_recovery_B_i = oil_recovery_B_i + CntrV_Zbtmflux(j,1)*dt;
else
oil_influx_B_i = oil_influx_B_i + abs(CntrV_Zbtmflux(j,1))*dt;
end
% Gas Flux
if CntrV_Zbtmflux(j,2) > 0
water_recovery_B_i = water_recovery_B_i + CntrV_Zbtmflux(j,2)*dt;
else
water_influx_B_i = water_influx_B_i + abs(CntrV_Zbtmflux(j,2))*dt;
end
end
oil_influx_B = [oil_influx_B oil_influx_B(end) + oil_influx_B_i];
oil_recovery_B = [oil_recovery_B oil_recovery_B(end) + oil_recovery_B_i];
water_influx_B = [water_influx_B water_influx_B(end) + water_influx_B_i];
water_recovery_B = [water_recovery_B water_recovery_B(end) + water_recovery_B_i];
end
%% RHS Boundary Flux
% CntrV RHS- Preallocate memory and Initialize phase fluxes - influx and recovery
oil_influx_R = cell((L+1),1);
oil_influx_R{1} = 0;
oil_recovery_R = cell((L+1),1);
oil_recovery_R{1} = 0;
water_influx_R = cell((L+1),1);
water_influx_R{1} = 0;
water_recovery_R = cell((L+1),1);
water_recovery_R{1} = 0;
for i=2:n+1 % iterate for every time step 1- n where n = 1000
CntrV_XRHSflux = states{i}.flux(CntrV_Xface_RHS,:); % All CntrV XLHS fluxes for time step i
oil_influx_R_i = 0;
oil_recovery_R_i = 0;
water_influx_R_i = 0;
water_recovery_R_i = 0;
for j=1:length(CntrV_Xface_RHS)
% Oil Flux
if CntrV_XRHSflux(j,1) > 0
oil_recovery_R_i = oil_recovery_R_i + CntrV_XRHSflux(j,1)*dt;
else
oil_influx_R_i = oil_influx_R_i + abs(CntrV_XRHSflux(j,1))*dt;
end
% Gas Flux
if CntrV_XRHSflux(j,2) > 0
water_recovery_R_i = water_recovery_R_i + CntrV_XRHSflux(j,2)*dt;
else
water_influx_R_i = water_influx_R_i + abs(CntrV_XRHSflux(j,2))*dt;
end
end
% Control Volume Total Fluxes
oil_influx_R = [oil_influx_R oil_influx_R(end) + oil_influx_R_i];
oil_recovery_R = [oil_recovery_R oil_recovery_R(end) + oil_recovery_R_i];
water_influx_R = [water_influx_R water_influx_R(end) + water_influx_R_i];
water_recovery_R = [water_recovery_R water_recovery_R(end) + water_recovery_R_i];
end
%% Net Control Volume Influx and Recovery
CntrV_W_influx_case1 = water_influx_T + water_influx_R + water_influx_B + water_influx_L;
CntrV_O_influx_case1 = oil_influx_T + oil_influx_R + oil_influx_B + oil_influx_L;
CntrV_cum_influx_case1 = CntrV_W_influx_case1 + CntrV_O_influx_case1;
CntrV_W_Recovery_case1 = water_recovery_T + water_recovery_R + water_recovery_B + water_recovery_L;
CntrV_O_Recovery_case1 = oil_recovery_T + oil_recovery_R + oil_recovery_B + oil_recovery_L;
CntrV_cum_Recovery_case1 = CntrV_W_Recovery_case1 + CntrV_O_Recovery_case1;
CntrV_cum_Netflux_case1 = CntrV_cum_influx_case1 - CntrV_cum_Recovery_case1;
%% Transfer Function
Net_W_flux_case1 = CntrV_W_influx_case1 - CntrV_W_Recovery_case1;
Net_O_flux_case1 = CntrV_O_influx_case1- CntrV_O_Recovery_case1;
%% Plot outputs
% Initialize plot variables
time = cell((n+1),1);
time{1} = 0;
inj_W_Model = cell((n+1),1);
inj_W_Model{1} = 0;
inj_W_CntrV = cell((n+1),1);
inj_W_CntrV{1} = 0;
prod_W_Model = cell((n+1),1);
prod_W_Model{1} = 0;
prod_O_Model = cell((n+1),1);
prod_O_Model{1} = 0;
prod_cum_Model = cell((n+1),1);
prod_cum_Model{1} = 0;
CntrV_tpv = (2*2*1*0.2); % Compute CntrV pv = CntrV Volume(Lx*Ly*Lz) * Matrix porosity (0.3)
%% Cummulative fluid injection and production from model
for i = 2:n+1
time = [time time(end)+ dt];
inj_W_Model = [inj_W_Model inj_water_GD200_Pe015_Kf100md(end) + sum(states{i}.flux(Ztop_faces,2))*dt];
inj_W_CntrV = [inj_W_CntrV inj_W_CntrV(end) + sum(states{i}.flux(CntrV_Zface_top,2))*dt];
prod_W_Model = [prod_W_Model prod_W_Model(end) + sum(states{i}.flux(Zbtm_faces,2))*dt];
prod_O_Model = [prod_O_Model prod_O_Model(end) + sum(states{i}.flux(Zbtm_faces,1))*dt];
prod_cum_Model = [prod_cum_Model prod_cum_Model(end) + sum(states{i}.flux(Zbtm_faces,1))*dt...
+ sum(states{i}.flux(Zbtm_faces,2))*dt];
end
%% Model Recoveries
figure(8);
subplot(2,1,1);
grid on;
hold on;
plot(time/day, prod_O_Model /tpv, 'DisplayName', 'RF_O(Δρ=200)', 'LineWidth',2);
hold on;
plot(time/day, prod_W_Model/tpv, 'DisplayName', 'RF_G(Δρ=200)', 'LineWidth',2);
xlabel('Time [days]');
ylabel('Recovery Factor');
subplot(2,1,2);
grid on;
hold on;
plot(time/day, CntrV_O_Recovery_case1/CntrV_tpv,'DisplayName','F-M_O Trnsf(Δρ=200)','LineWidth',2);
hold on
plot(time/day, CntrV_W_Recovery_case1/CntrV_tpv, 'DisplayName','F-M_W Trnsf(Δρ=200)', 'LineWidth',2);
xlabel('Time [days]');
ylabel('F-M Fluid Transfer-Outflux');
%%
figure(9);
subplot(2,1,1);
grid on;
hold on
plot(time/day, inj_W_Model/tpv, 'DisplayName','inj PV(Δρ=200)', 'LineWidth',2);
hold on;
plot(time/day, prod_cum_Model/tpv, 'DisplayName','Cum RF(Δρ=200)', 'LineWidth',2);
xlabel('Time [days]');
ylabel('PV Injected/Recovered - Model');
%
subplot(2,1,2);
grid on;
hold on;
plot(time/day, prod_O_Model /tpv, 'DisplayName', 'RF_O(Δρ=200)', 'LineWidth',2);
hold on;
plot(time/day, prod_W_Model/tpv, 'DisplayName', 'RF_G(Δρ=200)', 'LineWidth',2);
xlabel('Time [days]');
ylabel('Recovery Factor');
% subplot(2,2,3:4);
% grid on;
% hold on;
% plot(time/day, Net_O_flux_case1 /CntrV_tpv,'DisplayName','F-M NetOil Trnsf-CntrV 1D','LineWidth',2);
% hold on
% plot(time/day, Net_water_flux_GD200_Pe015_Kf100md/CntrV_tpv, 'DisplayName','F-M NetGas Trnsf-CntrV 1D', 'LineWidth',2);
% xlabel('Time [days]');
% ylabel('Normalised Net Gas Transfer');
%% CntrV Water Saturation Distribution Plot
% After Day 1
figure(10);
subplot(2,5,1);
plotCellData(G, states{1}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('Gas Sat in CntrV After Day 1 ');
% After Day 100
% figure(100);
subplot(2,5,2);
plotCellData(G, states{28}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 100 Days');
% After Day 250
% figure(11);
subplot(2,5,3);
plotCellData(G, states{69}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 250 Days');
% After 500 Days
% figure(12);
subplot(2,5,4);
plotCellData(G, states{137}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 500 Days');
% After Day 1000
% figure(13);
subplot(2,5,5);
plotCellData(G, states{274}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 1000 Days');
% After 1500 Days
% figure(14);
subplot(2,5,6);
plotCellData(G, states{411}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 1500 Days');
% After Day 2000
% figure(15);
subplot(2,5,7);
plotCellData(G, states{548}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 2000 Days');
% After 2500 Days
% figure(16);
subplot(2,5,8);
plotCellData(G, states{685}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 2500 Days');
% After 3000 Days
% figure(16);
subplot(2,5,9);
plotCellData(G, states{822}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 3000 Days');
% After 3650 Days
% figure(16);
% subplot(2,5,10);
plotCellData(G, states{1001}.s(:,1), CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
view([0.143826069396107 0.170936586318325]);
title('After 3650 Days');
%% CntrV Pressure Distribution Plot
% After Day 1
figure(11);
subplot(2,5,1);
plotCellData(G, states{1}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('Pressure Distr in CntrV After Day 1');
% After Day 1
% figure(170);
subplot(2,5,2);
plotCellData(G, states{28}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 100 Days');
% After 250 Days
subplot(2,5,3);
plotCellData(G, states{69}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 250 Days');
% After 500 Days
subplot(2,5,4);
plotCellData(G, states{137}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 500 Days');
% After 1000 Days
subplot(2,5,5);
plotCellData(G, states{274}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 1000 Days');
% After 1500 Days
subplot(2,5,6);
plotCellData(G, states{411}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 1500 Days');
% After 2000 Days
subplot(2,5,7);
plotCellData(G, states{548}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 2000 Days');
% After 2500 Days
subplot(2,5,8);
plotCellData(G, states{685}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 2500 Days');
% After 3000 Days
subplot(2,5,9);
plotCellData(G, states{822}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 3000 Days');
% After 3650 Days
subplot(2,5,10);
plotCellData(G, states{1001}.pressure(:,1)/mega, CntrV, 'EdgeColor', 'none', 'LineStyle', 'none');
% caxis([9.956, 10.0]);
view([0.143826069396107 0.170936586318325]);
title('After 3650 Days');
  1 comentario
per isakson
per isakson el 14 de Abr. de 2021
Its difficult to spot the cause of the poor performance by inspecting the code of your question.
Matlab provides good tools to analyze performance. They require running the code. See
For us to help you we need a code that we can run. Now the code of the question only produces "Undefined function or variable" errors.

Iniciar sesión para comentar.

Respuestas (1)

per isakson
per isakson el 8 de Abr. de 2021
Editada: per isakson el 8 de Abr. de 2021
"I have pre-allocated memory". After pre-allocating memory for a variable you must not change its size. If you do, the effect of pre-allocationg is lost.
>> checkcode('cssm.m') % cssm.m contains your code
L 170 (C 5-16): The variable 'oil_influx_T' appears to change size on every loop iteration. Consider preallocating for speed.
L 171 (C 5-18): The variable 'oil_recovery_T' appears to change size on every loop iteration. Consider preallocating for speed.
L 172 (C 5-18): The variable 'water_influx_T' appears to change size on every loop iteration. Consider preallocating for speed.
L 173 (C 5-20): The variable 'water_recovery_T' appears to change size on every loop iteration. Consider preallocating for speed.
... % altogether 22 warnings of the same type
All warnings refer to statements of this type, horizontal concatenations.
oil_influx_T = [oil_influx_T oil_influx_T(end) + oil_influx_T_i];
oil_recovery_T = [oil_recovery_T oil_recovery_T(end) + oil_recovery_T_i];
water_influx_T = [water_influx_T water_influx_T(end) + water_influx_T_i];
water_recovery_T = [water_recovery_T water_recovery_T(end) + water_recovery_T_i];
Search "oil_influx_T" with the editor, Notepad++
Line 138: oil_influx_T = cell((n+1), 1); % pre-allocate cell array
Line 139: oil_influx_T{1}= 0;
Line 151: oil_influx_T_i = 0; % assign a double scalar
Line 159: oil_influx_T_i = oil_influx_T_i + CntrV_Ztopflux(j,1)*dt;
Line 170: oil_influx_T = [oil_influx_T oil_influx_T(end) + oil_influx_T_i];
...
I don'f fully understand your code. (Maybe because I use R2018b.) Why doesn't the statement below throw an error? Concatenating cell array and double shouldn't work.
oil_influx_T = [oil_influx_T oil_influx_T(end) + oil_influx_T_i];
Anyhow, the statement changes the size of oil_influx_T, thus the warning and the poor performance.
  2 comentarios
Chukwuemeka Enemanna
Chukwuemeka Enemanna el 13 de Abr. de 2021
Hey Isakson,
Thanks for your prompt and detailed response. I have effected the corrections -
  • Correctly pre-allocated memory for oil_influx_T variables and
oil_influx_T = zeros(1, (n+1));
oil_recovery_T = zeros(1, (n+1));
water_influx_T = zeros(1, (n+1));
water_recovery_T = zeros(1, (n+1));
  • And corrected the concatenated output
oil_influx_T_increament = oil_influx_T(end) + oil_influx_T_i;
oil_influx_T(1, i) = oil_influx_T_increament;
clear oil_influx_T_increament;
oil_recovery_T_increament = oil_recovery_T(end) + oil_recovery_T_i;
oil_recovery_T(1, i) = oil_recovery_T_increament;
clear oil_recovery_T_increament
water_influx_T_increament = water_influx_T(end) + water_influx_T_i;
water_influx_T(1, i) = water_influx_T_increament;
clear water_influx_T_increament
water_recovery_T_increament = water_recovery_T(end) + water_recovery_T_i;
water_recovery_T(1, i) = water_recovery_T_increament;
clear water_recovery_T_increament
Currently running a case since about 11pm last night and still running till now.
per isakson
per isakson el 14 de Abr. de 2021
Editada: per isakson el 14 de Abr. de 2021
Thus, I was wrong about the main cause of the poor performance.
Have you tried running profile() ?

Iniciar sesión para comentar.

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by