How to stop time loop when steady state is reached?

I'm using unsteady case. so, it will reach a steady state at a certain time level. I fixed time at 'j'th column wise and it ran at maximum time level, but I want to stop the 'time loop' when steady state is reached. Even though I'm using a steady criteria i.e error >tolerance along dt>0 using 'while' loop, it didn't work and using dt>0, my code didn't stop.

clc; close all; clear all;
ymax=20; n=80; dy=ymax/n;        
xmax=1; m=20; dx=xmax/m;                       
tmax=100; nt=500; dt=tmax/nt; t=0:dt:tmax;      
%tol=1e-5; max_difference=1.0; 
UOLD=zeros(n,nt); VOLD=zeros(n,nt); 
TNEW=0; TOLD=ones(n,nt); 
A=zeros([1,m]); 
B=A;
C=A;
D=A;
T=TOLD; 
tic
%while dt>0 && max_difference>tol
for j=1:m
    for i=1:n
        if j>i
            C(i)=(dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
        elseif i>j
            A(i)=(-dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
        elseif i==j
            B(i)=1+(dt*UOLD(i,j)/2*dx)+((dt)/(dy^2));
        end
    end
end
for j=2:m
        for i=1:n
                D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)-TOLD(i,j)));
            
        end
    end
    T(:,j)=TriDiag(A,B,C,D);  %call tridiagonal 
    dt=0.2+dt;

Thanks for any advice.

[Mathieu NOE]

where does the TriDiag function comes from ?

[Mathieu NOE]

ok

I am not sure to understand all the details of your code , but something disturbs me

if we consider that the metric (error) that defines how your computation converge to the supposed steady state values

difference = (abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));

then if I plot it for the 80 iterations , we see the difference increasing and not decreasing .... so I wonder if you have 100% confidence in your code

ploted with y log scale to make it clear :

[Image Analyst]

Which loop, the one over i or the one over j, is the "time loop"?

I don't see any line like "error >tolerance". Where exactly is that line?

Answer 1

[Mathieu NOE]

hello

think I have understood

your difference computation must be indexed with j , so it must be inside the for loop

I used a break to exist the for loop when the difference is below tol

now, tol = 1e-5 will not be achieved with only 500 iteration

for the sake of the demo I put tol = 1e-2 and it will do only 103 iterations

plot of difference vs j with tol = 1e-5

plot of difference vs j with tol = 1e-2

code :

clc; close all; clear all;
ymax=20; n=80; dy=ymax/n;        
xmax=1; m=20; dx=xmax/m;                       
tmax=100; nt=500; dt=tmax/nt; t=0:dt:tmax;      
tol=1e-2; 
UOLD=zeros(n,nt); VOLD=zeros(n,nt); 
TNEW=0; TOLD=TNEW*ones(n,nt); TWALL=ones(1,length(t));
A=zeros([1,m]); 
B=A;
C=A;
D=A;
T=TOLD; 
difference(1) = 1;
tic
for j=1:nt
    for i=1:n
        if j>i
            C(i)=(dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
        elseif i>j
            A(i)=(-dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
        elseif i==j
            B(i)=1+(dt*UOLD(i,j)/2*dx)+((dt)/(dy^2));
        end
    end
end

for j=2:nt
    if j==1
        for i=1:n
            if i==1
                D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)-TOLD(i,j)))-((dt/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TNEW);
            elseif i==n
                D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TWALL(j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)))-(((-dt)/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TWALL);
            else
                D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)));
            end
        end
    else
        for i=1:n
            if i==1
                D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)+TOLD(i,j)))-((dt/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TNEW);
            elseif i==n
                D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TWALL(j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)))-(((-dt)/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TWALL(j));
            else
                D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)));
            end
        end
    end
    T(:,j)=TriDiag(A,B,C,D);  %call tridiagonal 
    dt=0.2+dt;
    TOLD=T;
    
    difference(j) = max(abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));
    
    if difference(j)<tol
        break
    end
end

[Mathieu NOE]

hmm

the data is constantly decreasing, so it does not follow your suggestion going first up then down

I cannot tell why is it so as this I think is pretty much related of your core computation (T)

maybe another metric insted of difference will show this inversed bell curve but I doubt as we can see on this plot that T has a uniform converging shape

obtained with this code line

contourf(log(T));colormap('jet');set(gca,'YDir','normal');colorbar('vert');

whatever slice you make in the x or y direction, I don't see how any metric based on T could show an inversed bell curve

also you can see the difference from one col to the previous one does not converge towards zero , there is still a slight trend

[Mathieu NOE]

Correction

maybe ploting the difference in loglog scale would show more clearly the initial bump in the convergence curve

it shows also that in order to reach tol = 1e-5 you probably need more than 50,000 iterations

clc; close all; clear all;
ymax=20; n=80; dy=ymax/n;        
xmax=1; m=20; dx=xmax/m;                       
tmax=100; nt=5000; dt=tmax/nt; t=0:dt:tmax;      
tol=1e-5; 
UOLD=zeros(n,nt); VOLD=zeros(n,nt); 
TNEW=0; TOLD=TNEW*ones(n,nt); TWALL=ones(1,length(t));
A=zeros([1,m]); 
B=A;
C=A;
D=A;
T=TOLD; 
difference(1) = NaN;
tic
for j=1:nt
    for i=1:n
        if j>i
            C(i)=(dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
        elseif i>j
            A(i)=(-dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
        elseif i==j
            B(i)=1+(dt*UOLD(i,j)/2*dx)+((dt)/(dy^2));
        end
    end
end

for j=2:nt
    if j==1
        for i=1:n
            if i==1
                D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)-TOLD(i,j)))-((dt/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TNEW);
            elseif i==n
                D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TWALL(j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)))-(((-dt)/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TWALL);
            else
                D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)));
            end
        end
    else
        for i=1:n
            if i==1
                D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)+TOLD(i,j)))-((dt/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TNEW);
            elseif i==n
                D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TWALL(j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)))-(((-dt)/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TWALL(j));
            else
                D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)));
            end
        end
    end
    T(:,j)=TriDiag(A,B,C,D);  %call tridiagonal 
    dt=0.2+dt;
    TOLD=T;
    
    difference(j) = max(abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));
    
    if difference(j)<tol
        break
    end
end

loglog(difference)

[Sam Chak]

@Gayathri, Can you mathematically define the desired steady-state tolerance in your case?

[Yanni]

@Sam Chak yes. steady state tolerance limit is between 1e-4 to 1e-6.