How can I ignor NaN value when I calculate in a loop?

Hello

I am calculating a loop with the Euler method. I put criteria for each cell when it reaches Tmax being vanished from my calculation and replace it with NaN. The reason is that I assume that the cell is not physically existing anymore.

The problem is that for the next time step, the whole matrix's Cells become NaN. This happens when it is expected that only the neighboring cells to the NaN gain more Temperature and meet Tmax, consequently disappearing.

Is there any way to avoid such a problem with using NaN?

I used here T_Mask to force the data that they includ NaN become zero with the aid of T_Mask(isnan(T_Mask))=0;

%The code below shows this section.

k_MLI = 250;                                        
cp = 890;                                       
rho = 2700;                                    
rhoLiq = 2375;                                 
den_MLI = 2700;                                      
Cp_MLI = 900;                                         
T0_MLI_L=6.9;                                      
%eps_MLI=0.04;                                       
T_deg= 700;                                    

k_glass = 0.8;
Cp_Glass = 1000;
rho_glass= 2500;
alpha = k_MLI/(cp*rho);                                  
alpha2  = k_glass/(Cp_Glass*rho_glass);
eps_rad = 0.9;
eps_MLI = 0.08;
eps_MLI_Last= 0.02;
eps_Spacer = 0.3;
sigma = 5.67e-8;

%% domain
Thick_MLI = 0.000009;     
Thick_Spacer=0.000691;   % (i-1)th and the (i)th radiative layer
layers=4;
Lz = 0.2;                %  width (m)
Ly = 0.2;                %  depth (m)
Lx = (layers*Thick_MLI)+((layers-1)*Thick_Spacer);             %  height (m)
Nz=50;                   
Ny=50;                   
Nx=layers;             
dz=Lz/Nz;               
dy=Ly/Ny;               

dx=1;

dxx= 1;
xc=(Nz+1)/2;            % center plane of the domain
yc=(Ny+1)/2;
zc=(Nx+1)/2;
dt = (1/(3*alpha*((1/dz^2)+(1/dy^2))));	% stable time step, sec for Euler ethod   % for the third direction

tmax=5000;  % max time, sec 
t=0;
T_0_1 = 0.3325*t+273.15;              % Initial temperature in all nodes in kelvin
T_0_2 = 800+273.15;

T_Heatexchanger = 273.15;
%% Initial conditions for Steady State
T = ones(Nz,Ny,Nx)*T_0_1;
Kg=0.04;
Ks=0.05

Ks = 0.0500

T(:,1,:)    =  T_0_1;   %T_left;
T(:,Ny,:)   =  T_0_1;   %T_right;
T(1,:,:)    =  T_0_1;   %T_front;
T(Nz,:,:)   =  T_0_1;   %T_back;
T(:,:,1)    = T_0_1 ;%T_bott;


BC2=1; % set to 0 to disable

%% Constants

iter = 1;     % counter of iterations
error  = 1; % initial error for iterations
k1 = (alpha*dt); 
k2 = dt/(cp*rho);
k3= (sigma*k2);
eps1 = (1/eps_rad)+(1/eps_MLI)-1;
eps2 = (1/eps_MLI)+(1/eps_MLI)-1;
eps3 = (1/eps_MLI_Last)+(1/eps_MLI_Last)-1;
iterTot= abs(round(tmax/dt));
%% plotting 
T_min=273;
T_max=1100;  %Kelvin
numisosurf=25; % number of isosurfaces
isovalues=linspace(T_min,T_max,numisosurf);
x=linspace(0,Lx,Nx);
y=linspace(0,Ly,Ny);
z= linspace(0,Lz,Nz);
[Z Y X]=meshgrid(z,y,x);
figure('units','pixels','position',[50 50 500 500])
Pixel = dz*dy;

%% Steady-State Solution
tic
i=2:Nz-1; %inner nodes along y
j=2:Ny-1; %inner nodes along x
k=2:Nx; %inner nodes along z
T_new=T;
T_Mask=T;

% method='Euler';
while t<=tmax  %error >= tol_ss % by tmax or by tolerance (for Staeady State)
    
    %% von Neuman B.C. (fluxes at walls)
    % grad T as central difference
    % for q=1:layers
    % if BC2==1
    %      
    
    
    for k= 2:Nx

        if k<Nx
            for kk= k:k
                Kg(:,:,kk)=Kg; %This will be replace by a funtion
                Ks(:,:,kk)=Ks;
            end
        end
        if k>=3 && k<Nx
            if anynan(T(i,j,k-1))
                for iiii= i(1):i(end)
                    for jjjj= j(1):j(end)
                        if anynan(T(iiii,jjjj,k-1))
                            T_new(iiii,jjjj,k) = T(iiii,jjjj,k) + k1 * (((T_Mask(iiii-1,jjjj,k)-(2*T_Mask(iiii,jjjj,k))+T_Mask(iiii+1,jjjj,k))/(dz*dz)) + ((T_Mask(iiii,jjjj-1,k)-(2*T_Mask(iiii,jjjj,k))+T_Mask(iiii,jjjj+1,k))/(dy*dy))) ...
                                + (k3/eps1)*(T(iiii,jjjj,1)^4-T(iiii,jjjj,k)^4)...
                                + (k3/eps2)*(T(iiii,jjjj,k+1)^4-T(iiii,jjjj,k)^4) ... 
                                + Kg(iiii,jjjj,k) *k2 *(T(iiii,jjjj,k+1)-T(iiii,jjjj,k))/(dx*dxx) ...
                                + Ks(iiii,jjjj,k) *k2 *(T(iiii,jjjj,k+1)-T(iiii,jjjj,k))/(dx*dxx);
                                
                        else
                            T_new(iiii,jjjj,k) = T(iiii,jjjj,k) + k1 * (((T_Mask(iiii-1,jjjj,k)-(2*T_Mask(iiii,jjjj,k))+T_Mask(iiii+1,jjjj,k))/(dz*dz))+((T_Mask(iiii,jjjj-1,k)-(2*T_Mask(iiii,jjjj,k))+T_Mask(iiii,jjjj+1,k))/(dy*dy))) ...
                                +(k3/eps2)*(T(iiii,jjjj,k-1)^4-T(iiii,jjjj,k)^4) ...
                                +(k3/eps2)*(T(iiii,jjjj,k+1)^4-T(iiii,jjjj,k)^4) ...
                                + Kg(iiii,jjjj,k-1)*k2*(T(iiii,jjjj,k-1)-T(iiii,jjjj,k))/(dx*dxx) + Ks(iiii,jjjj,k-1)*k2*(T(iiii,jjjj,k-1)-T(iiii,jjjj,k))/(dx*dxx) ...
                                + Kg(iiii,jjjj,k)*k2*(T(iiii,jjjj,k+1)-T(iiii,jjjj,k))/(dx*dxx) + Ks(iiii,jjjj,k)*k2*(T(iiii,jjjj,k+1)-T(iiii,jjjj,k))/(dx*dxx) ;
                            
                        end
                    end
                end
                
            else
                T_new(i,j,k) = T(i,j,k) + k1 * (((T_Mask(i-1,j,k)-(2*T_Mask(i,j,k))+T_Mask(i+1,j,k))/(dz*dz))+((T_Mask(i,j-1,k)-(2*T_Mask(i,j,k))+T_Mask(i,j+1,k))/(dy*dy))) ...
                    +(k3/eps2)*(T(i,j,k-1)^4-T(i,j,k)^4) ...
                    +(k3/eps2)*(T(i,j,k+1)^4-T(i,j,k)^4) ...
                    + Kg(i,j,k-1)*k2*(T(i,j,k-1)-T(i,j,k))/(dx*dxx) + Ks(i,j,k-1)*k2*(T(i,j,k-1)-T(i,j,k))/(dx*dxx) ...
                    + Kg(i,j,k)  *k2*(T(i,j,k+1)-T(i,j,k))/(dx*dxx) + Ks(i,j,k)  *k2*(T(i,j,k+1)-T(i,j,k))/(dx*dxx) ;
                end
        end
        
        if k==2
            T_new(i,j,k) = T(i,j,k)+  k1*(((T_Mask(i-1,j,k)-(2*T_Mask(i,j,k))+T_Mask(i+1,j,k))/(dz*dz)) + ((T_Mask(i,j-1,k)-(2*T_Mask(i,j,k))+T_Mask(i,j+1,k))/(dy*dy))) ...
                + (k3/eps1)*(T(i,j,1)^4-T(i,j,k)^4) ...
                + (k3/eps2)*(T(i,j,k+1)^4-T(i,j,k)^4) ...
                + Kg(i,j,k)*k2*(T(i,j,k+1)-T(i,j,k))/(dx*dxx) ...
                + Ks(i,j,k)*k2*(T(i,j,k+1)-T(i,j,k))/(dx*dxx);
        end
        if k == Nx
            T_new(i,j,k) = T(i,j,k)+ k1 *(((T_Mask(i-1,j,k)-(2*T_Mask(i,j,k))+T_Mask(i+1,j,k))/(dz*dz)) + ((T_Mask(i,j-1,k)-(2*T_Mask(i,j,k))+T_Mask(i,j+1,k))/(dy*dy))) ...
                + (k3/eps2)*(T(i,j,k-1)^4-T(i,j,k)^4) ...
                + (k3/eps3)*(T_Heatexchanger^4-T(i,j,k)^4) ...
                + Kg(i,j,k-1)*k2*(T(i,j,k-1)-T(i,j,k))/(dx*dxx) + Ks(i,j,k-1)*k2*(T(i,j,k-1)-T(i,j,k))/(dx*dxx);  Kg(i,j,k-1)*k2*(273.15-T(i,j,k))/(5*dx*dxx) ...
                
        end
        
        if max(max(T_new(i,j,k)>T_deg))
            for iii= i(1):i(end)
                for jjj= j(1):j(end)
                    
                    if T_new(iii,jjj,k)>T_deg
                        T_new(iii,jjj,k)=NaN;
                    end
                end 
            end
            
        end
        T_Mask(i,j,k)=T_new(i,j,k);  
        T_Mask(isnan(T_Mask))=0;
        
        T(i,j,k)=T_new(i,j,k);
        if mod(iter,5)==0
            
            Area(k,iter)= sum (sum (isnan (T(i,j,k))))*Pixel;
        end 
    end
    
    iter =iter+1;
    t    =iter*dt; % pointless for Steady State
    %     error=max(max(max(abs(T_new-T))));    
    % end
    if t<2408
        T(1:Nz,1:Ny,1) = 0.3325*t+273.15;
    else
        T(1:Nz,1:Ny,1) = T_0_2;
    end
    if mod(iter,100)==0  % e.g. plot every 5 iter
        plot3D(Z,Y,X,Lz,Ly,Lx,dz,dy,Thick_Spacer,T,T_min,T_max,isovalues,iter,t); 
        
    end
    
    
end

Unrecognized function or variable 'K_g_mod'.

%% 
toc 


disp('COMPLETE!');

Answer 1

[Naman Kaushik]

Hi Ali,

As per my understanding, you want to omit or ignore the NaN values in the loop.

One workaround for this, is to use the “isnan” function in MATLAB, which returns a logical array containing 1 (true) where the elements of A are NaN, and 0 (false) where they are not.

Consider the code shown below:

if(isnan(variable_name)) 
    %continue or perform the operation that you want 
else  
    %continue or perform the operation that you want 
end 

For more information on the “isnan” function, you can refer to the following documentation:

https://in.mathworks.com/help/matlab/ref/isnan.html

[Ali Hariri]

Yes, I have already solve the problem. Your method also works

Thank you