Pre-allocating cell containers

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');

[per isakson]

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

• Code Performance Profiler, improve performance and find potential problems in MATLAB® code
• profile, Profile execution time for function

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.

Answer 1

[per isakson]

"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.

[Chukwuemeka Enemanna]

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]

Thus, I was wrong about the main cause of the poor performance.

Have you tried running profile() ?