Conditions aux limites

[turbulent]

Bonjour tout le monde. Merci de m'aider à résoudre mon problème de formulation des conditions aux limites. Pour votre information, le problème réside exclusivement dans la façon d'établir les conditions aux limites: Primo, parce que l'expression de F(p) marche impeccablement dans un autre code. Secondo, dans ce code j'ai affecté un scalaire à F(p) et malgré tout le problème persiste, c'est à dire le code n'affiche aucune erreur et tourne sans fin. Merci

clc;clear all;close all;
Nx=5;
Nz=5;
R=3;
h0=.4;
h1=.2;
deltax=R/Nx;
D=10;
B=.35;
L=100;
deltaz=L/Nz; 
deltax=R/Nx; 
a=.1;
dt=(a*deltaz^2)/D; 
b=B*(dt/(2*deltaz));
pz=[2:Nz];
px=[2:Nx];
Z=[0:deltaz:L]; 
x=[0:deltax:R]; 
colorstflex = 'bkgrmcy';
plotStyle = {'s','*','-','o','^','-+',':'};
%%%%%%Initial condition %%%%%%%%%%%%%%%%%%%
t=0; F(1,:)=h0*ones(1,Nz+1);              % 
Z=F(1,:);                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(x,Z,plotStyle{1},'color',colorstflex(1),'LineWidth',1.5,'MarkerSize',3),hold on
for i=2:3%7
  if i==2
    T=60/30;%60*60/12;
  end
  if i==3
    T=60/12;%60*60/6;
  end
  %  if i==4
  %    T=60*60/3;
  %  end
  %  if i==5
  %    T=60*60/2;
  %  end
  %  if i==6
  %    T=60*60;%1*60;
  %  end
  %  if i==7
  %    T=120*60;%2*60;3
  %    end
  
  Ndt=round(T/dt); % Nombre de pas de temps
  for n=1:Ndt % Incrément de temps
    t=t+dt;
    Fn=F;
    %%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%
    while (px = 1)                                   %
      F(1)=h0;                                         %
    endwhile                                         %
    F(1)=Fn(2);                                      %
    while (px = Nx)                                  % 
      AN(2:pz)=Fn(2:pz);                               %
    endwhile                                         %
    while (px = Nx+1)                                %
      BN(2:pz)=Fn(2:pz);                               %
      Fn(pz:round(h0-h1)/deltaz)=-pz*deltaz;           %
      Fn(round(h0-h1)/deltaz:h0/deltaz)=h1-pz*deltaz;  %
    endwhile                                         %
    AN(2:pz)=BN(2:pz);                               %
    F(Nz+1)=Fn(Nz);                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Expression de F(p) en fonction de Fn 
    Z=-F(p);
  end
  if  i > 1 && i <5%
    T=[00 5 10];%T=[00 5 10 20 30  60 120];
    plot(x,Z,plotStyle{i},'color',colorstflex(i),'LineWidth',1.5,'MarkerSize',3), hold on
    legendInfo{1} = ['T = ' num2str(00) "\ \mn"];
  end
end

legend(legendInfo,"location", "southeast")
xlabel('x Coordinates (m)')
ylabel('Z (m)')
xlim([0 R])
ylim([0 h0])
%  set(gca,"ydir","reverse")

-----

[umfred]

A quel moment tes px sont mis à jour dans tes boucles while( px= ...) ?? et un test d'égalité, il faut doubler le =, il me semble; donc le while (px=1) fait que tu affectes 1 à px et que ton "test" est toujours vrai.
et px est un array, non ? (px=[2:Nx]) ?

[turbulent]

J'ai apporté une légère modification au code, mais toujours le même problème qui persiste:[code]

clc;clear all;close all;
Nx=5;
Nz=5;
R=3;
h0=.4;
h1=.2;
deltax=R/Nx;
D=10;
B=.35;
L=100;
deltaz=L/Nz; 
deltax=R/Nx; 
a=.1;
dt=(a*deltaz^2)/D; 
b=B*(dt/(2*deltaz));
pz=[2:Nz];
px=[2:Nx];
Z=[0:deltaz:L]; 
x=[0:deltax:R]; 
colorstflex = 'bkgrmcy';
plotStyle = {'s','*','-','o','^','-+',':'};
%%%%%%Initial condition %%%%%%%%%%%%%%%%%%%
t=0; F(1,:)=zeros(1,Nz+1);%h0*ones(1,Nz+1);              % 
Z=-F(1,:);                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(x,Z,plotStyle{1},'color',colorstflex(1),'LineWidth',1.5,'MarkerSize',3),hold on
for i=2:3%7
  if i==2
    T=60/30;%60*60/12;
  end
  if i==3
    T=60/12;%60*60/6;
  end
  %  if i==4
  %    T=60*60/3;
  %  end
  %  if i==5
  %    T=60*60/2;
  %  end
  %  if i==6
  %    T=60*60;%1*60;
  %  end
  %  if i==7
  %    T=120*60;%2*60;3
  %    end
  
  Ndt=round(T/dt); % Nombre de pas de temps
  for n=1:Ndt % Incrément de temps
    t=t+dt;
    Fn=F;
    %%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%
    while (px = 1)                                   %
      F(1)=h0;                                       %
    endwhile                                         %
    F(1)=Fn(2);                                      %
    while (px = Nx)                                  % 
      AN(2:pz)=F(2:pz);                              %
    endwhile                                         %
    while (px = Nx+1)                                %
      BN(2:pz)=F(2:pz);                              %
      F(pz:round(h0-h1)/deltaz)=-pz*deltaz;          %
      F(round(h0-h1)/deltaz:h0/deltaz)=h1-pz*deltaz; %
    endwhile                                         %
    AN(2:pz)=BN(2:pz);                               %
    F(Nz+1)=Fn(Nz);                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Expression de F(p) en fonction de Fn 
  end
  Z=-F(1,:);
  if  i > 1 && i <4%
    T=[00 5 10];%T=[00 5 10 20 30  60 120];
    plot(x,Z,plotStyle{i},'color',colorstflex(i),'LineWidth',1.5,'MarkerSize',3), hold on
    legendInfo{1} = ['T = ' num2str(00) "\ \mn"];
  end
end

legend(legendInfo,"location", "southeast")
xlabel('x (m)')
ylabel('Z (m)')
xlim([0 R])
ylim([0 h0])
set(gca,"ydir","reverse")

[turbulent]

[QUOTE]Après avoir doublé le =, j'ai reçu ce message:warning: colon arguments should be scalars
warning: called from
Cas3 at line 64 column 13
warning: colon arguments should be scalars
warning: called from
Cas3 at line 64 column 13
error: 'BN' undefined near line 64, column 64
error: called from
Cas3 at line 64 column 13
L'introduction de BN et AN, c'est juste une tentative de ma part pour établir l'égalité entre entre F(2z) à px=Nx et F(2z) à px=Nx+1. Apparemment, cette astuce ne marche pas, s'il y en a une autre merci de me la proposer.[QUOTE]

clc;clear all;close all;
Nx=5;
Nz=5;
R=3;
h0=.4;
h1=.2;
deltax=R/Nx;
D=10;
B=.35;
L=100;
deltaz=L/Nz; 
deltax=R/Nx; 
a=.1;
dt=(a*deltaz^2)/D; 
b=B*(dt/(2*deltaz));
pz=[2:Nz];
px=[2:Nx];
Z=[0:deltaz:L]; 
x=[0:deltax:R]; 
colorstflex = 'bkgrmcy';
plotStyle = {'s','*','-','o','^','-+',':'};
%%%%%%Initial condition %%%%%%%%%%%%%%%%%%%
t=0; F(1,:)=zeros(1,Nz+1);%h0*ones(1,Nz+1);              % 
Z=-F(1,:);                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(x,Z,plotStyle{1},'color',colorstflex(1),'LineWidth',1.5,'MarkerSize',3),hold on
for i=2:3%7
  if i==2
    T=60/30;%60*60/12;
  end
  if i==3
    T=60/12;%60*60/6;
  end
  %  if i==4
  %    T=60*60/3;
  %  end
  %  if i==5
  %    T=60*60/2;
  %  end
  %  if i==6
  %    T=60*60;%1*60;
  %  end
  %  if i==7
  %    T=120*60;%2*60;3
  %    end
  
  Ndt=round(T/dt); % Nombre de pas de temps
  for n=1:Ndt % Incrément de temps
    t=t+dt;
    Fn=F;
    %%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%
    while (px ==1)                                   %
      F(1)=h0;                                       %
    endwhile                                         %
    F(1)=Fn(2);                                      %
    while (px == Nx)                                  % 
      AN(2:pz)=F(2:pz);                              %
    endwhile                                         %
    while (px == Nx+1)                                %
      BN(2:pz)=F(2:pz);                              %
      F(pz:round(h0-h1)/deltaz)=-pz*deltaz;          %
      F(round(h0-h1)/deltaz:h0/deltaz)=h1-pz*deltaz; %
    endwhile                                         %
    AN(2:pz)=BN(2:pz);                               %
    F(Nz+1)=Fn(Nz);                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Expression de F(p) en fonction de Fn 
  end
  Z=-F(1,:);
  if  i > 1 && i <4%
    T=[00 5 10];%T=[00 5 10 20 30  60 120];
    plot(x,Z,plotStyle{i},'color',colorstflex(i),'LineWidth',1.5,'MarkerSize',3), hold on
    legendInfo{1} = ['T = ' num2str(00) "\ \mn"];
  end
end

legend(legendInfo,"location", "southeast")
xlabel('x (m)')
ylabel('Z (m)')
xlim([0 R])
ylim([0 h0])
set(gca,"ydir","reverse")

[A voir en vidéo sur Futura]

[turbulent]

Après avoir substituer les boucles while par celles de if et en commentant momentanément les lignes BN et AN, le code commence à répondre, mais il donne pas le résultat escompté. Cela pourrait être à cause du manque de cette condition aux limites.

clc;clear all;close all;
Nx=5;
Nz=5;
R=3;
h0=.4;
h1=.2;
deltax=R/Nx;
D=10;
B=.35;
L=100;
deltaz=L/Nz; 
deltax=R/Nx; 
a=.1;
dt=(a*deltaz^2)/D; 
b=B*(dt/(2*deltaz));
pz=[2:Nz];
px=[2:Nx];
Z=[0:deltaz:L]; 
x=[0:deltax:R]; 
colorstflex = 'bkgrmcy';
plotStyle = {'s','*','-','o','^','-+',':'};
%%%%%%Initial condition %%%%%%%%%%%%%%%%%%%
t=0; F(1,:)=zeros(1,Nz+1);%h0*ones(1,Nz+1);              % 
Z=-F(1,:);                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(x,Z,plotStyle{1},'color',colorstflex(1),'LineWidth',1.5,'MarkerSize',3),hold on
for i=2:3%7
  if i==2
    T=60*60/12;
  end
  if i==3
    T=60*60/6;
  end
  %  if i==4
  %    T=60*60/3;
  %  end
  %  if i==5
  %    T=60*60/2;
  %  end
  %  if i==6
  %    T=60*60;%1*60;
  %  end
  %  if i==7
  %    T=120*60;%2*60;3
  %    end
  
  Ndt=round(T/dt); % Nombre de pas de temps
  for n=1:Ndt % Incrément de temps
    t=t+dt;
    Fn=F;
    %%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%
%    while (px ==1)                                   %
%      F(1)=h0;                                       %
%    endwhile                                         %
%    F(1)=Fn(2);                                      %
%    while (px == Nx)                                 % 
%      AN(2:pz)=F(2:pz);                              %
%    endwhile                                         %
%    while (px == Nx+1)                               %
%      BN(2:pz)=F(2:pz);                              %
%      F(pz:round(h0-h1)/deltaz)=-pz*deltaz;          %
%      F(round(h0-h1)/deltaz:h0/deltaz)=h1-pz*deltaz; %
%    endwhile                                         %
%    AN(2:pz)=BN(2:pz);                               %
%    F(Nz+1)=Fn(Nz);                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (px ==1)                                      %
      F(1)=h0;                                       %
    end                                              %
    F(1)=Fn(2);                                      %
%    if (px == Nx)                                   % 
%      AN(2:pz)=F(2:pz);                             %
%    end                                             %
    if (px == Nx+1)                                  %
    %  BN(2:pz)=F(2:pz);                             %
      F(pz:round(h0-h1)/deltaz)=-pz*deltaz;          %
      F(round(h0-h1)/deltaz:h0/deltaz)=h1-pz*deltaz; %
    end                                              %
    %AN(2:pz)=BN(2:pz);                              %
    F(Nz+1)=Fn(Nz);                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
  
    % Expression de F(p) en fonction de Fn 

    end
  Z=-F(1,:);
  if  i > 1 && i <4%
    T=[00 5 10];%T=[00 5 10 20 30  60 120];
    plot(x,Z,plotStyle{i},'color',colorstflex(i),'LineWidth',1.5,'MarkerSize',3), hold on
    legendInfo{1} = ['T = ' num2str(00) "\ \mn"];
  end
end

legend(legendInfo,"location", "southeast")
xlabel('x (m)')
ylabel('Z (m)')
xlim([0 R])
ylim([0 h0])
set(gca,"ydir","reverse")

[umfred]

comme je l'avais dit, px est un array (un tableau) px=[2:Nx]; donc le message d'erreur qui apparait quand tu testes px==1 est logique: tu compares un tableau avec une valeur, tu veux surement tester le px courant de la boucle donc px(n) (donc px(n)==1)

[turbulent]

[QUOTE]En fait, je cherche à affecter des valeurs pour les conditions aux limites à px=1 et px=Nx+1 ainsi que à pz=1 et pz=Nz+1. Le px=1 et px=Nx+1 sont les indexes des coordonnées suivant l'axe des x à gauche et à droite du domaine du calcul. pz=1 et pz=Nz+1 sont les indexes des coordonnées suivant l'axe des z en haut et en bas du domaine du calcul. Donc, l'index ou le compteur px ne varie pas en fonction d'un autre compteur, celui du temps (n). Ce qui varie ici, c'est bien la fonction F. Merci pour l'intérêt que vous porter à ce code.

[umfred]

je re-re-redis que tu as défini px et pz en tant que tableaux
px=[2:Nx] est un tableau contenant 2,3,4,5=> px=[2,3,4,5] (idem pour pz)

tu as une boucle for n=1:Ndt, dans laquelle n prendra tour à tour les valeurs de 1 à Ndt; donc c'est peut-être ce n que tu dois tester au lieu de px .... Reposes le problème à plat sur papier et mathématiquement.

[turbulent]

Après avoir apporter des modifications profondes, le code ne signale aucune erreur mais il n'aboutit à aucun résultat (calcul sans fin). Pourriez-vous m'éclairer pourquoi.

clc;clear all;close all;
Nx=15;
Nz=15;
R=3;
h0=1.6;
h1=.2;
D=10;
B=.35;
L=1.00;
dz=h0/Nz; 
Nx=Nz;
dx=R/Nx; 
a=.1;
dt=(a*dz^2)/D; 
b=B*(dt/(2*dz));
pz=[2:Nz];
px=[2:Nx];
Z=[0:dz:h0];
x=[0:dx:R]; 
alpha=dt/(dz^2);
phi=.3;
ks=11e-5;
thetas=.3;
thetar=.01;
Ss=1;%5e-5;
alphaV=3.3;
nV=4.1;%1;%1.5;%.5%
mV=1-1/nV;

colorstflex = 'bkgrmcy';
plotStyle = {'s','-','+','o','^','-+',':'};

%%%%%%Initial condition %%%%%%%%%%%%%%%%%%%
t=0; H(Nx+1,Nz+1)=h0; H(1,:)=h0*ones(1,Nz+1); % 
F(1,:)=h0*ones(1,Nz+1); % 
h(Nx+1,Nz+1)=0; 
k(Nx+1,Nz+1)=100*ks; 
k(1,:)=ks*ones(1,Nz+1);
for K=1:Nx
  for j=1:Nz;
    h(1,:)=H(1,:)-j.*dz;                % 
  end
end
k(1,:)=ks; 
F(1,:)=h0*ones(1,Nz+1);
plot(x,F(1,:),plotStyle{1},'color',colorstflex(1),'LineWidth',1.5,'MarkerSize',3),hold on
legendInfo{1} = ['T = ' num2str(00) "\ \mn"];
for i=2:3%7
  if i==2
    T=60/6;% T=60*60/12;
  end
  if i==3
    T=300;%/3;%T=60*60/6;
  end
  %  if i==4
  %    T=60*60/3;
  %  end
  %  if i==5
  %    T=60*60/2;
  %  end
  %  if i==6
  %    T=60*60;%1*60;
  %  end
  %  if i==7
  %    T=120*60;%2*60;3
  %    end
  Ndt=round(T/dt); % Nombre de pas de temps
  for n=1:Ndt % Incrément de temps
    t=t+dt;
    for K=2:Nx
      for j=2:Nz;
        
        Hn=H;          
        h_new=h;
        H_new=H;
        F_new=F;
        %%%%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
        H(1,j) = h0;     %AB
        h(1,j) = h0-j*dz;     %AB
        H(K,Nz+1) = H_new(K,Nz); % BC
        h(K,Nz+1)= h_new(K,Nz)-dz; % BC
        H(1, K) =H(2, K);% AF
        h(1, K) =h_new(2, K)+dz;% AF
        H(Nx+1, 2:round(h1/dz)+1)=h1*ones(1,round(h1/dz));% FE
        h(Nx+1, 2:round(h1/dz)+1) =(h1-j*dz)*ones(1,round(h1/dz));% FE
        for j=round(h1/dz)+1:round(h1/dz)+3
          H(Nx, j) =dz*j;%.*ones(1,j);%j*dz;% ED
          h(Nx+1, j)=0;%zeros(1,j);%0;% ED
        end
        for j=round(h1/dz)+3:round(h0/dz)
          H(Nx,j)=H(Nx+1,j);% DC
          h(Nx+1,j)=h(Nx,j);% DC
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        c(K,j)=alphaV*(thetas-thetar)*(nV-1)*(alphaV*abs(min(h(K,j),0))).^nV.^(nV-1)./...
        ((1+(alphaV*abs(min(h(K,j),0))).^nV).^((2*nV-1)/nV));
        k(K,j)=ks*(1-(abs(min(h(K,j),0)).^(nV-1)).*(1+abs(min(h(K,j),0)).^nV).^(-mV)).^2./...
        ((1+abs(min(h(K,j),0)).^nV).^(mV/2));
        
        h(K,j)=h_new(K,j)+alpha.*((k(K,j)+k(K+1,j)).*(h_new(K+1,j)+h(K,j))-(k(K,j)+k(K-1,j)).*(h_new(K-1,j)+...
        h_new(K,j))+(k(K,j)+k(K,j+1)).*(h_new(K,j+1)+h(K,j))-(k(K,j)+k(K,j-1)).*(h_new(K,j-1)+h(K,j))+(k(K,j+1)-...
        k(K,j-1)))./(2.*Ss.*phi+c(K,j));
        if h(K,j)==0;
          F_new(K,j)=j*dz; 
          if  i > 1 && i <4%
            T=[00 5 10];%T=[00 5 10 20 30  60 120];
            
            plot(x,F(1,:),plotStyle{i},'color',colorstflex(i),'LineWidth',1.5,'MarkerSize',3),,hold on
            legendInfo{i} = ['T = ' num2str(T(i )) "\ \mn"];
          end
        end
      end      
      
    end
    
  end
end
legend(legendInfo,"location", "southeast")
xlabel('x (m)')
ylabel('Z (m)')
xlim([0 R])
ylim([0 h0])
%set(gca,"ydir","reverse")

[umfred]

peut-être en parti à cause du fait que tu réutilises j dans les boucles de ta partie des conditions initiales, ce qui fait changer le j de la boucle au dessus (for j=2:Nz) et du fait que tu fais 2'636'719*15 + 263'672*15 boucles, ce qui prend du temps ....

[MissJenny]

k(1,:)=ks*ones(1,Nz+1); 
for K=1:Nx   
  for j=1:Nz;    
     h(1,:)=H(1,:)-j.*dz;                %    
  end 
end
k(1,:)=ks;

cette partie est bizarre : tu commences par mettre quelque-chose dans la variable k, puis tu fais un calcul où k n'intervient pas, puis tu remets autre-chose dans k.

[umfred]

Je me demande si turbulent comprend le code qu'il a écrit
et effectivement, dans le code que tu relèves, on modifie Nx*Nz fois la même partie de la matrice soit h(1, : ) pour y mettre au final H(1, : )-Nz.*dz

[turbulent]

Non, K majuscule est différent de k minuscule. K est un indice et k est une variable

[turbulent]

Effectivement, vous avez raison. Je devais supprimer la ligne : k(1,=ks;

[turbulent]

[QUOTE]Je reconnais que je ne maîtrise pas la programmation, mais je comprends parfaitement la problématique sur laquelle je travaille. C'est pour cette raison que je sollicite ce forum. Merci pour les remarques que vous m'adressez.

[umfred]

As-tu décomposé sur papier les différentes étapes simples que le programme doit effectuer?
Par exemple, qu'est censé faire le bout de code relevé par MissJenny ? parce que là actuellement, la boucle ne sert à rien car au final, h(1, : ) = H(1, : )-Nz.*dz
par K n'est jamais utilisé et on est toujours sur les mêmes éléments pour H et h

[turbulent]

[QUOTE]J'ai supprimé la boucle de K pour la cdt initiale et j'ai décomposé la double boucle des cdts aux limites en celle de K là où il n'y a variation que de K et en boucle de j là où il n'y a variation que de j. Je crois que cela va prendre énormément de temps. Pour être plus clair, il s'agit de résoudre une EDP 2D instationnaire: 2D ( les 2 variables d'espaces) et instationnaire (c.a.d variable dans le temps).

[umfred]

        %%%%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
        H(1,j) = h0;     %AB
        h(1,j) = h0-j*dz;     %AB
[..]
        H(1, K) =H(2, K);% AF
        h(1, K) =h_new(2, K)+dz;% AF
[..]
        H(Nx+1, 2:round(h1/dz)+1)=h1*ones(1,round(h1/dz));% FE
        h(Nx+1, 2:round(h1/dz)+1) =(h1-j*dz)*ones(1,round(h1/dz));% FE
        for j=round(h1/dz)+1:round(h1/dz)+3
          H(Nx, j) =dz*j;%.*ones(1,j);%j*dz;% ED
          h(Nx+1, j)=0;%zeros(1,j);%0;% ED
        end
        for j=round(h1/dz)+3:round(h0/dz)
          H(Nx,j)=H(Nx+1,j);% DC
          h(Nx+1,j)=h(Nx,j);% DC
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ici si j = K, les 2ères lignes ne servent à rien ?
les 2 for à la fin, font varier j , le même j de la grande boucle dans laquelle ce code est inclus, et ne s'exécute pas le nombre de fois demandé (bien ou mal, aucune idée, mais problème de logique ou de nom de variable)

[turbulent]

Normalement, j n'est pas forcément égal à K. Ce sont les indices des 2 variables indépendants (d'espace). Voici, ma dernière modification:

clc;clear all;close all;
Nx=15;
Nz=15;
R=3;
h0=1.6;
h1=.2;
%deltax=R/Nx;
D=10;
B=.35;
L=1.00;
dz=h0/Nz;%/round((h0-h1)/deltaz);%h0/Nz; 
Nx=Nz;
dx=R/Nx; 
a=.1;
dt=(a*dz^2)/D; 
b=B*(dt/(2*dz));
pz=[2:Nz];
px=[2:Nx];
%Z=[0:deltaz:h0]; 
Z=[0:dz:h0];
x=[0:dx:R]; 
alpha=dt/(dz^2);
phi=.3;
ks=11e-5;
thetas=.3;
thetar=.01;
Ss=1;%5e-5;
alphaV=3.3;
nV=4.1;%1;%1.5;%.5%
mV=1-1/nV;

colorstflex = 'bkgrmcy';
plotStyle = {'s','-','+','o','^','-+',':'};
%%%%%%Initial condition %%%%%%%%%%%%%%%%%%%
t=0; H(Nx+1,Nz+1)=h0; H(1,:)=h0*ones(1,Nz+1); % 
F(1,:)=h0*ones(1,Nz+1); % 
%H_new(Nx+1,Nz+1)=h0; H_new(1,:)=h0;
h(Nx+1,Nz+1)=0; 
k(Nx+1,Nz+1)=100*ks; 
k(1,:)=ks*ones(1,Nz+1);
%for K=1:Nx
for j=1:Nz;
  h(1,:)=H(1,:)-j.*dz;                % 
end
%end
%h_new(Nx+1,Nz+1)=0; h_new(1,:)=0;
%k_new(Nx,Nz)=ks; k_new(1,:)=ks;                                %
F(1,:)=h0*ones(1,Nz+1);
plot(x,F(1,:),plotStyle{1},'color',colorstflex(1),'LineWidth',1.5,'MarkerSize',3)
%[X,V]=meshgrid(x,Z);
%figure(1)
%    contourf(x,Z,h)
%    clabel(contourf(x,Z,h))
%    colorbar
%    colormap(jet)
,hold on
legendInfo{1} = ['T = ' num2str(00) "\ \mn"];
for i=2:3%7
  if i==2
    T=60/6;% T=60*60/12;
  end
  if i==3
    T=300;%/3;%T=60*60/6;
  end
  %  if i==4
  %    T=60*60/3;
  %  end
  %  if i==5
  %    T=60*60/2;
  %  end
  %  if i==6
  %    T=60*60;%1*60;
  %  end
  %  if i==7
  %    T=120*60;%2*60;3
  %    end
  Ndt=round(T/dt); % Nombre de pas de temps
  for n=1:Ndt % Incrément de temps
    t=t+dt;
    for j=round(h1/dz)+1:round(h1/dz)+3
      H(Nx, j) =dz*j;%.*ones(1,j);%j*dz;% ED
      h(Nx+1, j)=0;%zeros(1,j);%0;% ED
    end
    for j=round(h1/dz)+3:round(h0/dz)
      
    end
    %Hn=H;          
    h_new=h;
    H_new=H;
    %F_new=F;
    %%%%%%Boundary conditionns %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
    for j=2:Nz;
      H(1,j) = h0;     %AB
      h(1,j) = h0-j*dz;     %AB
      h(Nx+1, 2:round(h1/dz)+1) =(h1-j*dz)*ones(1,round(h1/dz));% FE
    end      
    
    for K=2:Nx
      H(K,Nz+1) = H_new(K,Nz); % BC
      h(K,Nz+1)= h_new(K,Nz)-dz; % BC
      H(1, K) =H(2, K);% AF
      h(1, K) =h_new(2, K)+dz;% AF
    end
    H(Nx,j)=H(Nx+1,j);% DC
    h(Nx+1,j)=h(Nx,j);% DC
    H(Nx+1, 2:round(h1/dz)+1)=h1*ones(1,round(h1/dz));% FE
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for K=2:Nx
      for j=2:Nz;
        c(K,j)=alphaV*(thetas-thetar)*(nV-1)*(alphaV*abs(min(h(K,j),0))).^nV.^(nV-1)./...
        ((1+(alphaV*abs(min(h(K,j),0))).^nV).^((2*nV-1)/nV));
        %Ss(i,j)=5e-5;%(1+(alphaV*abs(min(h(i,j),0))).^nV).^(-mV);
        k(K,j)=ks*(1-(abs(min(h(K,j),0)).^(nV-1)).*(1+abs(min(h(K,j),0)).^nV).^(-mV)).^2./...
        ((1+abs(min(h(K,j),0)).^nV).^(mV/2));
        
        h(K,j)=h_new(K,j)+alpha.*((k(K,j)+k(K+1,j)).*(h_new(K+1,j)+h(K,j))-(k(K,j)+k(K-1,j)).*(h_new(K-1,j)+...
        h_new(K,j))+(k(K,j)+k(K,j+1)).*(h_new(K,j+1)+h(K,j))-(k(K,j)+k(K,j-1)).*(h_new(K,j-1)+h(K,j))+(k(K,j+1)-...
        k(K,j-1)))./(2.*Ss.*phi+c(K,j));
        if abs(h(K,j)) ==0;%<.1%,
          %if h(K,j)==0;
          F=j*dz; % 
          if  i > 1 && i <4%
            T=[00 5 10];%T=[00 5 10 20 30  60 120];
            
            plot(x,F(1,:),plotStyle{i},'color',colorstflex(i),'LineWidth',1.5,'MarkerSize',3),,hold on
            legendInfo{i} = ['T = ' num2str(T(i )) "\ \mn"];
          end
        end
      end
    end
    %    figure(1)
    %    contourf(x,Z,h)
    %    clabel(contourf(x,Z,h))
    %    colorbar
    %    colormap(jet)
    
  end
end%;hold on
legend(legendInfo,"location", "southeast")
xlabel('x (m)')
ylabel('Z (m)')
xlim([0 R])
ylim([0 h0])
%set(gca,"ydir","reverse")

[umfred]

je réitère que ce code

for j=1:Nz;
  h(1,:)=H(1,:)-j.*dz;                % 
end

revient à

h(1,:)=H(1,:)-Nz.*dz

puisque j ne sert que pour la multiplication d'une matrice sans qu'il y ai une réutilisation d'autre chose.

Je réitère que à un moment j et K vont valoir la même valeur (d'autant plus avec ces boucles) et donc que h(1,j) et H(1,j) vont se faire écraser par les valeurs h(1,K) et H(1,K) de la boucle suivante.

for j=2:Nz;
      H(1,j) = h0;     %AB
      h(1,j) = h0-j*dz;     %AB
      h(Nx+1, 2:round(h1/dz)+1) =(h1-j*dz)*ones(1,round(h1/dz));% FE
    end      
    
    for K=2:Nx
      H(K,Nz+1) = H_new(K,Nz); % BC
      h(K,Nz+1)= h_new(K,Nz)-dz; % BC
      H(1, K) =H(2, K);% AF
      h(1, K) =h_new(2, K)+dz;% AF
    end

[turbulent]

En fait, H(1,: )= Nz.*dz=h0
du coup, on peut écrire

for j=1:Nz;
  %h(1,:)=H(1,: )-j.*dz;
  h(1,:)=(Nz-j).*dz;  % 
end

qui est complètement différent de

h(1,:)=H(1,:)-Nz.*dz

car h(1,: ) dépend de j.*dz.
Concernant la 2° remarque, là, il s'agit des conditions aux limites différentes: une concerne des conditions aux limites sur un axe vertical (j) et l'autre sur un un axe horizontal (K). Maintenant, numériquement, je ne sais pas s'il faut les écrire d'une autre façon ou pas.

[umfred]

h(1,: )=(Nz-j).*dz; c'est un produit scalaire (je crois que c'est le terme): tu multiplies le vecteur dz par un nombre et tu mets le résultat dans la 1ère rangée de h à chaque tour, donc au final, tu n'auras que le résultat du dernier tour soit (Nz-Nz).*dz = zeros(size(dz)) et donc ce qui est bien égal à H(1,: )-Nz.*dz puisque selon toi H(1, =Nz.*dz.
(et encore j et dz étant des scalaires, des nombres, j.*dz==j*dz)

la 2nde remarque c'est pareil, tu accèdes aux mêmes éléments h(1,j) == h(1,K) pour j=K (vu que tu boucles pour j de 2 à Nx(=15) et pour K de 2 à Nz(=15) tu boucles sur les mêmes éléments, donc tu modifies les mêmes éléments)

Si tu retires les points virgules, tu verras en direct les résultats des opérations et ainsi voir si c'est conforme à ce que tu attends

[turbulent]

[QUOTE]Chaque point (K, j) dans le domaine du calcul à sa propre h qui dépend surtout de sa position verticale (j). Par définition, H(K,j)=h(K,j)+j*dz, dz c'est le pas et non pas un vecteur. Mon objectif dans ce code est de tracer les courbes F(t) constituées par les points qui ont h=0.
Pour la 1° remarque:
Initialement (t=0), la 1° courbe F correspondante à h nulle est une droite horizontale qui coupe l'axe des z à h0. Maintenant, la boucle affecte des 0 à h, c'est juste, mais au même temps je dirai qu'elle est inutile, du moment que j'ai posé F(1,=h0*ones(1,Nz+1). Ainsi, je supprime carrément du code la boucle:

%for K=1:Nx
%for j=1:Nz;
  %h(1,:)=H(1,:)-j.*dz;
  h(1,:)=Nz.*dz-j.*dz;  % 
end
%end

ainsi que la ligne

H(1,:)=h0*ones(1,Nz+1);

Pour la 2° remarque, laisse moi un peu de temps, il faut que j'y réfléchisse. Merci pour ces remarques précieuses.

[JPL]

Turbulent, cela fait au moins la quatrième fois que je vois un [QUOTE] isolé qui ne correspond à rien au début de tes messages. Il y a quelque chose que tu ne semble pas bien maîtriser.

[umfred]

Met le "end "du for aussi en commentaire et il faut donner une valeur à j. Et comme tu veux lui donner la valeur Nz, cela affecte 0, donc pourquoi "se casser les bonbons" à faire un calcul.
(.* permet de multiplier une valeur à chaque élément d'une matrice, or ici Nz, dz, j sont des valeurs simples donc un simple multiplication suffirait)
(j'avais évoqué dz comme vecteur à cause de l'opérateur .*, mais en effet, ensuite j'ai remarqué que c'était un un simple nombre)

[MissJenny]

au fait, c'est quel langage?

[umfred]

Je pense que c'est du matlab (vu que % est utilisé pour les commentaires)

[turbulent]

GNU Octave

[umfred]

GNU Octave est très très proche de Matlab (GNU Octave est open source et gratuit, Matlab ne l'est pas et payant)

[turbulent]

Du moment que la boucle mène à l'écrasement des valeurs par la dernière valeur de la boucle, j'ai inséré ces 2 boucles (celles des conditions aux limites) dans la boucle de calcul et du plot. Apparemment, même cette procédure est insuffisante pour aboutir à un résultat.
une question: Si size(F(t=0)) est différente de size(F(t>0)) le code ne va pas marcher pas, non?