MATLAB error when run

% Given parameters
Pr = 2.65;              % Average density (g/cm³)
Kr = 7.75e-3;           % Thermal conductivity (cal/cm s °C)
Cp = 0.197;             % Heat capacity (cal/g °C)
Tinit = 323;            % Initial temperature (K)
Tapplied = 423;         % Applied temperature (K)
tmax = 2;               % Simulation time (years)
reservoirExtent = 75;   % Reservoir extent (m)
% Convert units
Kr = Kr * 418.4;        % Convert thermal conductivity to (W/m K)
Cp = Cp * 4.184;        % Convert heat capacity to (J/kg K)
% Discretization parameters
nr = 100;               % Number of radial grid points
dr = reservoirExtent / (nr - 1); % Radial grid spacing
% Time discretization parameters
dt = 0.01;              % Time step size
nt = round(tmax * 365.25 * 24 * 60 * 60 / dt); % Number of time steps
% Initialize temperature matrix
T = zeros(nr, nt+1);

Error using zeros
Requested 100x6311520001 (4702.4GB) array exceeds maximum array size preference (30.9GB). This might cause MATLAB to become unresponsive.

T(:, 1) = Tinit;        % Set initial temperature
% Perform time-stepping
for i = 1:nt
    % Perform radial discretization
    for j = 2:nr-1
        % Calculate thermal diffusivity
        alpha = Kr / (Cp * Pr);
        
        % Calculate radial derivatives
        dT_dr = (T(j+1, i) - T(j-1, i)) / (2 * dr);
        d2T_dr2 = (T(j+1, i) - 2 * T(j, i) + T(j-1, i)) / (dr^2);
        
        % Update temperature using finite difference method
        T(j, i+1) = T(j, i) + alpha * dt * (d2T_dr2 + (1 / j) * dT_dr);
    end
    
    % Apply boundary conditions
    T(1, i+1) = T(2, i+1);                                % Symmetry boundary condition
    T(nr, i+1) = T(nr-1, i+1) + dr * (Tapplied - T(nr, i+1)); % Heat conduction at reservoir boundary
end
% (1) Reservoir temperature at a distance of 20 m after 2 years
distance = 20;
index = round(distance / dr) + 1;  % Add 1 to account for MATLAB indexing
temperature = T(index, end);
% (2) Energy needed to heat the reservoir at 20 m
mass = Pr * reservoirExtent * pi * dr^2;
energy = mass * Cp * (Tapplied - Tinit);
% (3) Viscosity change at 20 m (using a hypothetical correlation)
viscosityInitial = 10;   % Initial viscosity (cP)
correlationConstant = 0.5;
viscosityChange = correlationConstant * (temperature - Tinit);
% Plotting temperature profile
r = linspace(0, reservoirExtent, nr);
figure;
plot(r, T(:, end));
xlabel('Radial Distance (m)');
ylabel('Temperature (K)');
title('Temperature Profile in the Reservoir');
% Displaying the results
fprintf('Reservoir temperature at 20 m after 2 years: %.2f K\n', temperature);
fprintf('Energy needed to heat the reservoir at 20 m: %.2f J\n', energy);
fprintf('Viscosity change at 20 m: %.2f cP\n', viscosityChange);

Answer 1

[Walter Roberson]

You are asking to calculate and store data for 100 radial grid points for 100 samples per second for 2 complete years .

2 complete years is about 6.3 giga-samples, and you want 100 times that so about 630 giga-samples, each of which takes 8 bytes (double precision): it is not surprising that it wants 4702 gigabytes.

You have several choices:

1. find a computer system with at least 13 terabytes of RAM that you can use; or
2. provide your system with at least 13 terabytes of swap space, and turn off the preference that limits array size, and let your program run, which will take a long time; or
3. do all necessary work to convert to tall arrays; or
4. don't store the iterations for nearly as long -- do you truly need to save all of that to disk or truly need to draw a graph 6 billion pixels wide?? or
5. reduce the resolution by about a factor of 200