Help indexing: How to index at the interfaces of cells and not the cells

2 views (last 30 days)
Jamie Al on 7 May 2021
Edited: per isakson on 9 May 2021
I am struggling to understand how to setup some indexing problem in MATLAB. If my entire domain is defined from 1 to N with a padding of two cells. So say my internal domain (domain without padded cells) from 2 to N-2, then if I want to calculate some values at the "interfaces" of the cells instead at the cells, how can I setup this on MATLAB? For example I have the following code for the internal domain that I would like to change my indicies to the interfaces?
My set up is like:
% calculate Flux at j+1/2
%
% j+1/2 Cell's grid:
% | wL| |
% | /|wR | 1 2 3 4 N-2 N-1 N
% | / |\ | {x=0} |-o-|-o-|-o-|-o-| ... |-o-|-o-|-o-| {x=L}
% |/ | \ | 1 2 3 4 5 N-1 N N+1
% | | \|
% | | | NC: Here cells 1 and N are ghost cells%
% j j+1
I keep getting the error:
Index in position 2 exceeds array bounds
The code for example:
function [ U_L, U_R] = MUSCL(U,N,NEQ )
size_U = size(U);
U_L = zeros(size_U);
U_R = zeros(size_U);
delta = 1e-6; % epsilon
for i= 1:NEQ
for j = 3:N-1
denom_L = (U(i,j) - U(i,j-1));
denom_R = (U(i,j+2) - U(i,j+1)); % error here since index exceeds the bounds of the array N (1001 here)
denom_L = sign(denom_L+1e-64)*max(delta, abs(denom_L));
denom_R = sign(denom_R+1e-64)*max(delta, abs(denom_R));
r_L = (U(i,j) - U(i,j-1))/denom_L;
r_R = (U(i,j+2) - U(i,j+1))/denom_R;
theta_L = Limiter(r_L);
theta_R = Limiter(r_R);
end
end
for j = 2:N-2
U_L(:,j) = U(:,j)+0.5*theta_R.*(U(:,j)-U(:, j-1));
U_R(:,j) = U(:,j+1)-0.5*theta_L.*(U(:,j+2)-U(:,j+1));
end
end
How can I find U_L, U_R, r_L, r_R at j+1/2 ??

per isakson on 8 May 2021
Edited: per isakson on 8 May 2021
Index in position 2 exceeds array bounds
U(:,j+2) is most likely the problem.
I assume that the size of U is NEQxN. The maximum value of j is N-1 (for j = 3:N-1). Thus, the maximum value of j+2 is N+1, which exceeds the array bound, N
per isakson on 9 May 2021
Your code is about "flux", "cells" and "interfaces". The name, MUSCL, does that stand for "muscle"? However, I neither understand the physics nor the code. Futhermore, I cannot run the code, since I neither have the function, Limiter(), nor a realistic sample of U.
It looks strange to me that none of the left-hand-side variables in the double loop is indexed.
%% A sample matrix
M = 40 .* membrane( 1, 25 );
%%
to_the_interfaces = @(row) interp1( (1:numel(row)), row, (1.5:1:numel(row)-0.5) );
row = (1:10);
to_the_interfaces(row)
ans = 1×9
1.5000 2.5000 3.5000 4.5000 5.5000 6.5000 7.5000 8.5000 9.5000
%%
M_at_the_interfaces = nan(size(M)-[0,1]);
for jj = 1 : size(M,1)
M_at_the_interfaces(jj,:) = to_the_interfaces(M(jj,:));
end
imagesc( M_at_the_interfaces ) DGM on 8 May 2021
Edited: DGM on 8 May 2021
Just interpolate.
x = 0:10;
y = x.^2
xhalf = 0.5:9.5;
yhalf = interp1(x,y,xhalf,'linear')
plot(x,y,'rd',xhalf,yhalf,'bo') Jamie Al on 8 May 2021
I have no clue how to apply this in my example. Can you show something for my for loop?