Problem with mod command

Hi everyone,

I am running a 3d simulation of a concentration profile in a 1m^3 box with a sinkterm. I have added a line: if(mod(t,x) ==0) to update the plot for a certain time interval due to computations becomming to difficult to run for my pc. However, for x > 0.3 my 3d simulation will not display anything. Another problem I have encountered is that the frequency with which the simulation is updated changes. First it updates very fast, but as the simulation is running the update interval changes as displayed in the images. It seems like there a factor 4 with which the interval changes.

clear all
close all

%paramters
diff  = 4.058*10^-5;

%dimensions
Lx = 10;
Ly = 10;
Lz = 10;
Nx = 21; Nt  = 400000; %amount of iterations
Ny = 21;
Nz = 21;
dx = Lx/(Nx-1);
dy = Ly/(Ny-1);
dz = Lz/(Nz-1);

%CFL conditions, determines how big dt should be for the equation to
%converge
c = 1;
C=0.05; %C<1
dt = dx*C/c;

%field variables
cn=zeros(Nx,Ny,Nz); %concentration
x=linspace(0,Lx,Nx); %xdistance
y=linspace(0,Ly,Ny); %ydistance
z=linspace(0,Lz,Nz); %zdistance
[X,Y,Z]=meshgrid(x,y,z); %defining 3D grid

%Value of Diffusion
K = ones(Nx,Ny,Nz)+diff;
K([1 end],:,:) = 10^-15; %insulated problem which means that Diffusion is almost zero at the boundaries end does both sides
K(:,[1 end],:) = 10^-15;
K(:,:,[1 end]) = 10^-15;

%initial condition
cn(:,:,:) = 0.15*10^-9;
t=0;

for n=1:Nt
   cc = cn;
    t=t+dt; %new time
     
    %New temperature
for i=2:Nx-1
    for j=2:Ny-1
        for k=2:Nz-1
            cn(i,j,k) = cc(i,j,k)+dt*(K(i,j,k))*...
       ((cc(i+1,j,k) - 2*cc(i,j,k) + cc(i-1,j,k))/dx/dx +...
       (cc(i,j+1,k) - 2*cc(i,j,k) + cc(i,j-1,k))/dy/dy +...
       (cc(i,j,k+1) - 2*cc(i,j,k) + cc(i,j,k-1))/dz/dz);
        end
    end
end

%Insulation conditions
cn(1,:,:)=cn(2,:,:); %walls
cn(end,:,:) = cn(end-1,:,:); %walls
cn(:,1,:)=cn(:,2,:); %walls
cn(:,end,:) = cn(:,end-1,:); %walls 
cn(:,:,1)=cn(:,:,2); %floor
cn(:,:,end) = cn(:,:,end-1); %roof

%sink term Filter
cn(4,4,1) = 0; 

%Visualization 
    if(mod(t,0.06) == 0)
    slice(X,Y,Z,cn,5,5,0); colorbar; caxis([0 1.5*10^-10]);
    title(sprintf('time = %f minutes', t/60))
    pause(0.01);
    end
end

Answer 1

[Walter Roberson]

http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

[Yorick de]

Does this mean the error has to do with the sprintf command?

[Walter Roberson]

No, the error is in thinking that t+dt will reliably end up becoming exact integers when dt is not an integer or simple negative power of 2. Or, in your case, that t+dt will accumulate to become exact multiples of

>> fprintf('%.999g\n', 0.6)
0.59999999999999997779553950749686919152736663818359375

0.59999999999999997779553950749686919152736663818359375

Especially if Lx and Nx were ever altered to not have the exact ratio 1 : 2

C=0.05

C will not be exactly 1/20

>> fprintf('%.999g\n', 0.05)
0.05000000000000000277555756156289135105907917022705078125

If you want to update every 60 iterations of n, then instead of testing t, you should test your loop counter, n, like mod(n, 60)

[Yorick de]

Thank you so much! it worked