Vectorization of an if condition with cirle function

1 view (last 30 days)
Steffen B.
Steffen B. on 18 Mar 2024
Edited: Steffen B. on 22 Mar 2024
Hello,
I have a question regarding vectorization of loops.
In the case the value hT is changing depending on the position of the nodes (x-z plaine) with respect to two concentric circles.
It's not difficult to script this with a few loops, but this is extremly slow.
The value of hT, which is changing due to the geometric position is afterwards used in function TRB. The function TRB can be easy vectorized, but the I have the problem that the value of hT is not variet of the index variable j and ij.
Is it possible to make the code leaner, and faster with vectorization of the if elseif else condtion?
Here is the Structure of the slow loop script:
for i=[1 ny]
if i==1
for ij=2:1:nz-1
for j=2:1:nx-1
if ((j*dx-x0)^2+(ij*dz-z0)^2)<(RLN2^2)
hT=h(m,1);
TinfT=TinfTLN2N2;
elseif (((j*dx-x0)^2+(ij*dz-z0)^2)>((RLN2)^2))&&(((j*dx-x0)^2+(ij*dz-z0)^2)<=((RN2)^2))
hT=hTN2;
TinfT=TinfTLN2N2;
else
hT=hTAir;
TinfT=TinfTAir;
end
%
TRB(i,j,ij)=(hT*dx*dy+ BlaBla*To(1,j-1,ij)+);
end
end
else
The loop is cheking if the node is in the inner circle, between the inner and outer circle or outside of the circles and for each case the value of hT is different. Afterwards the TRB is calculated depending on the hT value.

Sign in to answer this question.

Accepted Answer

Voss
Voss on 18 Mar 2024
Ran in:
Use logical indexing. Something like as follows:
% make up missing variable values:
nz = 20;
nx = 24;
dz = 0.1;
dx = 0.1;
z0 = 1;
x0 = 1.2;
RLN2 = 0.4;
RN2 = 0.8;
h = 333;
m = 1;
hTN2 = 22;
hTAir = 1;
ij = 2:1:nz-1;
j = 2:1:nx-1;
xx2 = (j*dx-x0).^2;
zz2 = (ij.'*dz-z0).^2;
is_inner = xx2 + zz2 < RLN2^2
is_inner = 18×22 logical array
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
is_outer = xx2 + zz2 <= RN2^2
is_outer = 18×22 logical array
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
hT = hTAir*ones(size(is_inner));
hT(is_inner) = h(m,1);
hT(~is_inner & is_outer) = hTN2;
disp(hT)
1 1 1 1 1 1 1 1 1 1 22 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 22 22 22 22 22 22 22 1 1 1 1 1 1 1 1 1 1 1 1 1 22 22 22 22 22 22 22 22 22 22 22 1 1 1 1 1 1 1 1 1 1 22 22 22 22 22 22 22 22 22 22 22 22 22 1 1 1 1 1 1 1 1 1 22 22 22 22 22 22 333 22 22 22 22 22 22 1 1 1 1 1 1 1 1 22 22 22 22 22 333 333 333 333 333 22 22 22 22 22 1 1 1 1 1 1 1 22 22 22 22 333 333 333 333 333 333 333 22 22 22 22 1 1 1 1 1 1 1 22 22 22 22 333 333 333 333 333 333 333 22 22 22 22 1 1 1 1 1 1 22 22 22 22 333 333 333 333 333 333 333 333 22 22 22 22 22 1 1 1 1 1 1 22 22 22 22 333 333 333 333 333 333 333 22 22 22 22 1 1 1 1 1 1 1 22 22 22 22 333 333 333 333 333 333 333 22 22 22 22 1 1 1 1 1 1 1 22 22 22 22 22 333 333 333 333 333 22 22 22 22 22 1 1 1 1 1 1 1 1 22 22 22 22 22 22 22 22 22 22 22 22 22 1 1 1 1 1 1 1 1 1 22 22 22 22 22 22 22 22 22 22 22 22 22 1 1 1 1 1 1 1 1 1 1 22 22 22 22 22 22 22 22 22 22 22 1 1 1 1 1 1 1 1 1 1 1 1 1 22 22 22 22 22 22 22 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 22 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
That makes hT a matrix of size nz-2 by nx-2. If it should be nx-2 by nz-2, swap which of j and ij is transposed in the expressions for xx2 and zz2. If it should be nx by nz or nz by nx, other adjustments need to be made.
Also note the different behavior between this and your original when a node is exactly on the RLN2 circle (xx2+zz2==RLN2^2); in the original, hT would have the value hTAir in that case.
  2 Comments
Steffen B.
Steffen B. on 19 Mar 2024
Edited: Steffen B. on 19 Mar 2024
Thanks for the answer,
it's helping me a lot. Unfortunatley the use of the matrix for hT and Tinf in the vectorized Function TRB is causing some trouble.
% hT Matrix
ij=2:1:nz-1;
j=2:1:nx-1;
xx2=(j*dx-x0).^2;
zz2=(ij.'*dz-z0).^2;
is_inner=xx2+zz2<RLN2^2;
is_outer=xx2+zz2<=RN2^2;
hT=hTAir*ones(size(is_inner));
hT(is_inner)=h(m,1);
hT(~is_inner & is_outer)=hTN2;
% Tinf Matrix
ij=2:1:nz-1;
j=2:1:nx-1;
xx2=(j*dx-x0).^2;
zz2=(ij.'*dz-z0).^2;
is_inner=xx2+zz2<RLN2^2;
is_outer=xx2+zz2<=RN2^2;
TinfT=TinfTAir*ones(size(is_inner));
TinfT(is_inner)=TinfTLN2N2;
TinfT(~is_inner & is_outer)=TinfTLN2N2;
%
TRB(1,2:nx-1,2:nz-1)=(hT*dx*dz*TinfT+lambda*dy*dz/(2*dx)*To(1,1:nx-2,2:nz-1)+...
lambda*dy*dz/(2*dx)*To(1,3:nx,2:nz-1)+lambda*dy*dx/(2*dz)*To(1,2:nx-1,1:nz-2)+...
lambda*dy*dx/(2*dz)*To(1,2:nx-1,3:nz)+lambda*dx*dz/(dy)*To(2,2:nx-1,2:nz-1))/...
(hT*dx*dz+lambda*dy*dz/(2*dx)+lambda*dy*dz/(2*dx)+lambda*dy*dx/(2*dz)+lambda*dy*dx/(2*dz)+lambda*dx*dz/(dy));
%
TRB is the Temperature value Matrix of the top surface at the Simulation domaine with the size T(1,2:nx-1,2:nz-1),
but this is not working
Error using /
Arguments must be 2-D, or at least one argument must be scalar. Use RDIVIDE (./) for elementwise right division.
Error in Test_Vect_Red_Algo_Vanila_FDMcnBCAlphaFinal (line 261)
lambda*dy*dx/(2*dz)*To(1,2:nx-1,3:nz)+lambda*dx*dz/(dy)*To(2,2:nx-1,2:nz-1))/...
When I try to fix it wit ./
TRB(1,2:nx-1,2:nz-1)=(hT*dx*dz*TinfT+lambda*dy*dz/(2*dx)*To(1,1:nx-2,2:nz-1)+...
lambda*dy*dz/(2*dx)*To(1,3:nx,2:nz-1)+lambda*dy*dx/(2*dz)*To(1,2:nx-1,1:nz-2)+...
lambda*dy*dx/(2*dz)*To(1,2:nx-1,3:nz)+lambda*dx*dz/(dy)*To(2,2:nx-1,2:nz-1))./...
(hT*dx*dz+lambda*dy*dz/(2*dx)+lambda*dy*dz/(2*dx)+lambda*dy*dx/(2*dz)+lambda*dy*dx/(2*dz)+lambda*dx*dz/(dy));
The next error message arise
Unable to perform assignment because the size of the left side is 1-by-19-by-19 and the size of the right side is 19-by-19-by-19.
Error in Test_Vect_Red_Algo_Vanila_FDMcnBCAlphaFinal (line 259)
TRB(1,2:nx-1,2:nz-1)=(hT*dx*dz*TinfT+lambda*dy*dz/(2*dx)*To(1,1:nx-2,2:nz-1)+...
Voss
Voss on 19 Mar 2024
Ran in:
In addition to using ./ you need to use .* for all the multiplications involving matrices, to get element-wise multiplications (unless they really should be matrix multiplications).
Also, since it appears that To is of size n by nx by nz (for some n>1), it's convenient for hT and TinfT to be nx-2 by nz-2 rather than nz-2 by nx-2. (I mentioned the possibility of this ambiguity in my answer.) So I've swapped the transposes in the definitions of ij and j.
Now, the error is because To has an extra dimension in front compared to hT and TinfT. To account for that you can make a permuted copy of each of hT and TinfT that are of size compatible with To for element-wise multiplication.
Here's an example of the basic idea:
M = ones(1,3,4); % size: 1x3x4
P = ones(3,4); % size: 3x4
try
result = M.*P; % this can't be done
catch e
disp(e.message) % so you get this error:
end
Arrays have incompatible sizes for this operation.
% permute P to add another dimension in front to match the size of M, and
% now the mulitplication can be done:
P_permuted = permute(P,[3 1 2]); % size: 1x3x4 (same as M)
result = M.*P_permuted; % no problem
And here it is applied to your TRB calculation:
% made up stuff:
nz = 20;
nx = 24;
dz = 0.1;
dx = 0.1;
z0 = 1;
x0 = 1.2;
RLN2 = 0.4;
RN2 = 0.8;
h = 333;
m = 1;
hTN2 = 22;
hTAir = 1;
TinfTAir = 99;
TinfTLN2N2 = 55;
To = rand(2,nx,nz);
lambda = 0.345;
dy = 0.1;
% common variables used in both hT and Tinf calculation:
ij=2:1:nz-1;
j=2:1:nx-1;
xx2=(j.'*dx-x0).^2;
zz2=(ij*dz-z0).^2;
is_inner=xx2+zz2<RLN2^2;
is_outer=xx2+zz2<=RN2^2;
% hT Matrix
hT=hTAir*ones(size(is_inner));
hT(is_inner)=h(m,1);
hT(~is_inner & is_outer)=hTN2;
% Tinf Matrix
TinfT=TinfTAir*ones(size(is_inner));
TinfT(is_inner)=TinfTLN2N2;
TinfT(~is_inner & is_outer)=TinfTLN2N2;
% TRB
hT_p = permute(hT,[3 1 2]);
TinfT_p = permute(TinfT,[3 1 2]);
TRB_2(1,2:nx-1,2:nz-1) = ...
(hT_p*dx*dz.*TinfT_p+lambda*dy*dz./(2*dx).*To(1,1:nx-2,2:nz-1)+...
lambda*dy*dz/(2*dx)*To(1,3:nx,2:nz-1)+lambda*dy*dx/(2*dz).*To(1,2:nx-1,1:nz-2)+...
lambda*dy*dx/(2*dz)*To(1,2:nx-1,3:nz)+lambda*dx*dz/(dy).*To(2,2:nx-1,2:nz-1))./...
(hT_p*dx*dz+lambda*dy*dz/(2*dx)+lambda*dy*dz/(2*dx)+lambda*dy*dx/(2*dz)+lambda*dy*dx/(2*dz)+lambda*dx*dz/(dy));

Sign in to comment.

More Answers (1)

Steffen B.
Steffen B. on 22 Mar 2024
Edited: Steffen B. on 22 Mar 2024
Thank you for the answer,
I already solved the problem by myself.
I embedded the Tinf and hT matrixes in a suitable one filled with ones.
%
Ln=ones(nx-2,ny-2);
LL(1,:,:)=Ln;
hTneu(1,:,:)=hT;
TinfTneu(1,:,:)=TinfT;
%
TRB(1,2:nx-1,2:nz-1)=(hTneu.*TinfTneu*dx*dz+lambda*dy*dz/(2*dx)*To(1,1:nx-2,2:nz-1)+...
lambda*dy*dz/(2*dx)*To(1,3:nx,2:nz-1)+lambda*dy*dx/(2*dz)*To(1,2:nx-1,1:nz-2)+...
lambda*dy*dx/(2*dz)*To(1,2:nx-1,3:nz)+lambda*dx*dz/(dy)*To(2,2:nx-1,2:nz-1))./...
(hTneu*dx*dz+lambda*dy*dz/(2*dx)+lambda*dy*dz/(2*dx)+lambda*dy*dx/(2*dz)+lambda*dy*dx/(2*dz)+lambda*dx*dz/(dy));
Now it’s working quite well, and a lot faster than before.
Thanks again for your help.
Have a nice day.

Categories

Find more on Logical in Help Center and File Exchange

Tags

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!