Error using / Matrix dimensions must agree.

Hi there, I'm trying to piece together this 2-D implicit heat code. However I am having trouble creating kappa which can be used in further loops when implicitly solving the equation. Firstly this does not work, if I create a matrix (Km) and run the code it says it cannot run as the sizes do not match on either side.

I have attached the code below:

% Solution of 2D temperature equation on a staggered grid 
% for a non-moving medium with constant conductivity;
% conversion to implicit formulation.


% Clean all variables
clear all;
% Clear all figures
clf;

% Numerical model parameters
% Boundary conditions = constant temperature
% setup correspond to a squre temperature wave
% Model size, m
xsize=100000; % Horizontal
ysize=100000; % Vertical
% Numbers of nodes
Nx=51; % Horizontal
Ny=31; % Vertical
% Grid step
dx=xsize/(Nx-1); % Horizontal
dy=ysize/(Ny-1); % Vertical
% Defining coordinates of staggered points
% P-Points
xpr=-dx/2:dx:xsize+dx/2; % Horizontal
ypr=-dy/2:dy:ysize+dy/2; % Vertical

% Defining the number of timesteps in that the loop iterates through
ntimesteps=20; % number of timesteps

% Material properties
K=3; % Thermal conductivity, W/m/K
% Defining density in xy-points
RHOCP=zeros(Ny+1,Nx+1); % Volumetric heat capacity, J/K/m^3

% Making vectors for nodal points positions (basic nodes)
x=0:dx:xsize; % Horizontal
y=0:dy:ysize; % Vertical

% Timestep
dt=min(dx,dy)^2/(4*K/min(min(RHOCP))); % Timestep, s

% Creating temperature arrays (profiles) for iteration of temperature changes through time
TDT=zeros(Ny+1,Nx+1); % New temperature, K
T0=zeros(Ny+1,Nx+1); % Old temperature, K

% Populating T0 using assigned temperatures and RHOCP at eulerian nodes
for j=1:1:Nx+1
    for i=1:1:Ny+1
        d=((xpr(j)-xsize/2)^2+(ypr(i)-ysize/2)^2)^0.5; % distance to the centre from the given nodal point
        if(d<20000)
            % Plume
            RHOCP(i,j)=3.2e+6;
            T0(i,j)=1873;
        else
            % Medium
            RHOCP(i,j)=3.3e+6;
            T0(i,j)=1573;
        end
        % Sticky air layer 
        if(ypr(i)<ysize*0.2)
            RHOCP(i,j)=3.5e+6;
            T0(i,j)=273;
        end
    end
end

% Formulating kappa (κ) to simplify the implicit equation:
% Equation: κ = k/RHO*Cp
% Thereby formulated as:
kappa(i,j)=K/RHOCP;

Error using /
Matrix dimensions must agree.

% Beginning Time-Cycle for implicit solution
timesum=0; % Elapsed time
for t=1:1:ntimesteps

    % Matrix of coefficients initialization for implicit solving
    L=sparse(Nx*Ny,Nx*Ny);
    % Vector of right part initialization for implicit solving
    R=zeros(Nx*Ny,1);
    
    % Implicitly solving for the 2-D temperature equation using 5-point cross:
    % RHOCP*dT/dt=K(d2T/dx^2+d2T/dy^2)
    %       |               ∆x
    %  _____|________________|___________________
    %       |                -
    %       |               To2
    %       |                -
    %       |                -
    %  ∆y---|   --To1-------To3-------To5--
    %       |                -
    %       |                -
    %       |               To4
    %       |                -
    %       |
    % 
    % Composing matrix of coefficients L()
    % and vector (column) of right parts R()
    % Process all grid points
    for i=1:1:Ny
        for j=1:1:Nx
            % Global index
            k=(j-1)*Ny+i;
            % Boundary nodes
            if(i==1 || i==Ny || j==1 || j==Nx)
                % Upper boundary
                if(i==1)
                    L(k,k)=1;
                    R(k,1)=273;                
                end
                % Lower boundary
                if(i==Ny)
                    L(k,k)=1;
                    R(k,1)=1573;               
                end
                % Left boundary
                if(j==1 && i>1 && i<Ny)
                    L(k,k)=1;
                    R(k,1)=0;                 
                end
                % Right boundary
                if(j==Nx && i>1 && i<Ny)
                    L(k,k)=1;
                    R(k,1)=0;                   
                end
            % Internal nodes
            else
                % dT/dt=kappa*(d2T/dx2+d2T/dy2)
                % TDT(i)/dt-kappa*(TDT(i,j-1)-2*TDT(i,j)+TDT(i,j+1))/dx^2-kappa*(TDT(i-1,j)-2*TDT(i,j)+TDT(i+1,j))/dy^2=T0(i,j)/dt
                L(k,k-Ny)=-kappa/dx^2; % coefficient for T(i,j-1)
                L(k,k-1)=-kappa/dy^2; % coefficient for T(i-1,j)
                L(k,k)=1/dt+2*kappa/dx^2+2*kappa/dx^2; % coefficient for T(i,j)
                L(k,k+1)=-kappa/dy^2; % coefficient for T(i+1,j)
                L(k,k+Ny)=-kappa/dx^2; % coefficient for T(i,j+1)
                R(k,1)=T0(i,j)/dt;
            end
        end
    end

    % Obtaining solution
    S=L\R;
    % Reloading solution to the new temperature array
    for i=1:1:Ny
        for j=1:1:Nx
            % Global index
            k=(j-1)*Ny+i;
            % Reload solution
            TDT(i,j)=S(k);
        end
    end
    
    % Reloading TDT to T0 for the interation of the next timestep
    T0=TDT;
   
    % Openning figure
    figure(1);
    
    % Visualization of geometrical arrays OMEGA(), PSI(), vx(), vy()
    figure(1);clf;colormap('Jet')
    % Drawing RHOCP()
    subplot(1,2,1)
    pcolor(xpr,ypr,RHOCP)
    shading flat
    colorbar
    axis ij image;
    title('RHOCP, J/m^3/K')
    % Plotting implicit solution (TDT) for every timestep
    subplot(1,2,2);
    pcolor(xpr,ypr,TDT);
    shading interp;
    colorbar;
    axis ij image;
    title(['Implicit: step=',num2str(t),' time,Myr=',num2str(timesum/(1e+6*365.25*24*3600))]);
    
    % Stop for .1 second
    pause(0.2);
    
    % Add time counter
    timesum=timesum+dt;
    
    % Reassigning temperature profile (TDT) for the next step as base
    % temperature profile (T0)
    T0=TDT;
end

The error code is as follows:

Error using  / 
Matrix dimensions must agree.

Any help or a prod in a direction will be greatly appreciated.

Regards,

Aleks

Answer 1

[Torsten]

kappa=K./RHOCP;

instead of

kappa(i,j)=K/RHOCP;

[Alex]

@Torsten

That seems to have worked however this now causes the following error:

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in Explicit_Implicit_2D (line 130)

L(k,k-Ny)=-kappa./dx^2; % coefficient for T(i,j-1)

[Torsten]

kappa is a matrix. You cannot assign a matrix (-kappa./dx^2) to a scalar (L(k,k-Ny)).

Most probably, all the "kappa" have to be replaced by "kappa(i,j)" in this section of your code.

[Alex]

super thank you very much