0
$\begingroup$

I want help with my Octave code regarding the unsteady heat equation.

My geometry and mesh are generated with FreeFEM++, so there is no problem with that (I tried it with the steady problem with no for loop and it worked just fine).

But when I introduce the for loop to iterate with the time step and try to visualize the results at multiple steps, the results remain the same for the entire loop with no change in figures or in Temperature data.

I tried varying the time step, the maximum number of iterations, the source term, and alpha (implicit/explicit Euler and Crank-Nicolson)

the formula I used

Results Here is the script :

% Visual aspects
set(0,'defaultAxesFontSize',20);
set(0,'defaultAxesLineWidth',2);

% Physical parameters
kappa=0.005;
amplitude_s=0.0; % 0.0
Thole=5;
Tborder=2;
Tcolormap=[2 5];

% Read the informations from the geometry
x=readarray('x.dat');
y=readarray('y.dat');
mesh_info=load('meshinfo.dat');

% Read matrices from the Finite Element Formulation
% Read mass matrix
M=readmat('M.dat');
[n,m]=size(M);
%pcolor(M);
%
K=readmat('K.dat');
%
% Get the position of the Dirichlet masks
BC_Dir_hole=readarray('BCdir5.dat');
ii_Dir_hole=find(abs(BC_Dir_hole)>0);
BC_Dir_borders=readarray('BCdir1234.dat');
ii_Dir_borders=find(abs(BC_Dir_borders)>0);
BC_Dir=readarray('BCdir.dat');
ii_int=find(BC_Dir==0);
%
% Build corresponding matrices
T=zeros(n,1);
alpha=0;
%alpha=1
%alpha=1/2
N=100;
delta_t=0.001;
A=1/delta_t*M+alpha*kappa*K;
Ared=A(ii_int,ii_int);
s=amplitude_s*ones(n,1);
s_new=s;
B=(1/delta_t*M-(1-alpha)*kappa*K)*T+M*(alpha*s_new+(1-alpha)*s);
% SOLVE
Td=zeros(n,1);
Td(ii_Dir_hole)=Thole;
Td(ii_Dir_borders)=Tborder;
%
for nt=1:N
C=B-A*Td;
Cred=C(ii_int);
Tred=Ared\Cred;
%
T=Td;
T(ii_int)=Tred;
%
if mod(nt,10)==0;
    FIG=visu_2D(x,y,mesh_info,T,'none','off');
    %FIG=visu_2D_visibility(x,y,mesh_info,T,'none','on');
     caxis(Tcolormap);
     %
     savename=['PIC/PIC_' sprintf('%05d',not)];
     print(FIG,'-r200','-djpg',savename);
     close all;
  end
 endfor

Thanks in advance.

$\endgroup$
1
  • 1
    $\begingroup$ Can you write down the differential equations, boundary and initial conditions for your problem? $\endgroup$
    – nicoguaro
    Commented Nov 2, 2022 at 0:04

1 Answer 1

3
$\begingroup$

Reading quickly through your code, it seems that you never update B. So the term $T^n$ of the previous step remains unchanged, hence you are effectively computing the same step over and over again...

There might be other issues but I cannot test the code at the moment.

$\endgroup$
3
  • $\begingroup$ Indeed, how did I miss that! a simple copy paste of the formula inside the loop did the trick. $\endgroup$
    – mibo27
    Commented Nov 1, 2022 at 15:30
  • $\begingroup$ Thanks you so much $\endgroup$
    – mibo27
    Commented Nov 1, 2022 at 15:30
  • $\begingroup$ @mibo27 You're welcome. Can you accept the answer so that this question is marked as "resolved" ? $\endgroup$
    – Laurent90
    Commented Nov 1, 2022 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.