PDE convection-diffusion equation using the fully implicit method.

24 Apr 2022
0 Answers
Updated 23 Jan 2025
3 Views (30 days)

0 votes

I have wrote thd code below. While runing it I got the following errores that I need help to fix them:
Operator '+' is not supported for operands of type 'struct'.
Error in HW7>createMesh3D (line 394)
G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
Error in HW7>createMesh3D (line 408)
MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]);
Error in HW7 (line 350)
m = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the mesh
394 G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
The code:
%CP7
%Written by Viridiana Salazar
clear variables
close all
%% 1. Space discretization
Lx = 1.0; dx =0.125; N=Lx/dx+1; x=0:dx:Lx;
Ly = 1.0; dy =0.025; y=0:dy:Ly;
Lz = 1.0; dz =0.025; z=0:dz:Lz;
%% 2. Time discretization
tf =0.5; dt =0.001; M=tf/dt+1; t= 0:dt:tf;
%% Constants
Mu=0.5; r=Mu*dt/(dx)^2; q=-0.5*dt/dx;
%% Analytical Solution (We're gonna take from here the initial and
% boundary conditions)
UA = zeros(N,N,N,M);
for n=1:M %time
for i=1:N %x
for j=1:N %y
for k=1:N %z
UA(i,j,k,n)=(1+exp(x(i)/(2*Mu)+y(j)/(2*Mu)+z(k)/(2*Mu)-3*t(n)/(4*Mu)))^(-1);
end
end
end
end
%% Boundary Conditions
U = zeros(N,N,N,M);
% For x=1 and x=N
U(1,:,:,:)=UA(1,:,:,:); U(N,:,:,:)=UA(N,:,:,:);
%For y=1 and y=N
U(:,1,:,:)=UA(:,1,:,:); U(:,N,:,:)=UA(:,N,:,:);
%For z=1 and z=N
U(:,:,1,:)=UA(:,:,1,:); U(:,:,N,:)=UA(:,:,N,:);
%% Initial conditions
% For t=0
U(:,:,:,1)=UA(:,:,:,1);
%%
for n=1:M-1 %Time loop
% Initial guest at level n+1
for i=1:N %x
for j=1:N %y
for k=1:N %z
if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N %At the boundaries 1 and N, U is known
U(i,j,k,n+1)=U(i,j,k,n);
end
end
end
end
iteration=0;
EPS=1;
while EPS>0.0000000001 && iteration<20
%Jacobian Matrix
J=zeros(N^3,N^3);
% Vector f
F=zeros(N^3,1);
i=0;
for I=1:N^3
K=I;
% Fix the index i,j,k
if mod(I,N*N)==1
i=i+1;
j=0;
end
if mod(I,N)==1
k=1;
j=j+1;
else
k=k+1;
end
if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
F(K,1)=(q*(-3*U(i,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j,k-1,n+1)+U(i,j,k-1,n+1))+(-(6*r+1)*U(i,j,k,n+1)+r*(U(i+1,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j+1,k,n+1)+U(i,j-1,k,n+1)+U(i,j,k+1,n+1)+U(i,j,k-1,n+1)))+U(i,j,k,n));
end
%% derivatives for i=N and/or j=N and/or k=N
if i==N && j==N && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i==N && j~=N && k~=N
if i==N && j~=N && k~=N && j~=1 && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j==1 && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j~=1 && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j==N && k~=N
if i~=N && j==N && k~=N && i~=1 && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k~=N && i==1 && k~=1
J(K,I-1)=q*+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j==N && k~=N && i~=1 && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j~=N && k==N
if i~=N && j~=N && k==N && i~=1 && j~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j~=N && k==N && i==1 && j~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j~=N && k==N && i~=1 && j==1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j==N && k~=N
if i==N && j==N && k~=N && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j==N && k~=N && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j~=N && k==N
if i==N && j~=N && k==N && j~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q*+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k==N && j==1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i~=N && j==N && k==N
if i~=N && j==N && k==N && i~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k==N && i==1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%% derivatives for i=1 and/or j=1 and/or k=1
if i==1 && j==1 && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i==1 && j~=1 && k~=1 && j~=N && k~=N
if i==1 && j~=1 && k~=1 && j~=N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j==N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j~=N && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i~=1 && j==1 && k~=1
if i~=1 && j==1 && k~=1 && i~=N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i==N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i~=N && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=1 && j~=1 && k==1
if i~=1 && j~=1 && k==1 && i~=N && j~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i==N && j~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i~=N && j==N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i==1 && j==1 && k~=1
if i==1 && j==1 && k~=1 && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j==1 && k~=1 && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i==1 && j~=1 && k==1 && j~=N
if i==1 && j~=1 && k==1 && j~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k==1 && j==N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i~=1 && j==1 && k==
if i~=1 && j==1 && k==1 && i~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k==1 && i==N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
end
% Solving the system of equations J*dU=-F, we can use LU decomposition to save
% computational resources in robust system
dU=(J\(-F));
EPS=max(abs(dU(:,1)));
a=0;
for i=1:N %x
for j=1:N %y
for k=1:N %z
a=a+1;
if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
U(i,j,k,n+1)=U(i,j,k,n+1)+dU(a);
end
end
end
end
iteration=iteration+1;
end
end
%% Plots and subrotines for 3D plotting
m = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the mesh
Numerical=CellVariable(m, U(:,:,:,M));
Analytical=CellVariable(m, UA(:,:,:,M));
Error=CellVariable(m,abs(U(:,:,:,M)-UA(:,:,:,M)));
figure
visualizeCells3D(Numerical);
title('Numerical solution at t=0.25')
figure
visualizeCells3D(Analytical);
title('Analytical solution at t=0.25')
figure
visualizeCells3D(Error);
title('Absolute Error at t=0.25')
function visualizeCells3D(phi)
% Modified from 2012-2016 Ali Akbar Eftekhari
phi.value = phi.value(2:end-1,2:end-1,2:end-1);
[X,Y,Z]=meshgrid(phi.domain.cellcenters.y, phi.domain.cellcenters.x, ...
phi.domain.cellcenters.z);
phi.value(1)=phi.value(1)+eps; % to avoid an strange error for assigning color limits
Sx = [phi.domain.cellcenters.x(1) phi.domain.cellcenters.x(end)];
Sy = [phi.domain.cellcenters.y(1) phi.domain.cellcenters.y(end)];
Sz = [phi.domain.cellcenters.z(1) phi.domain.cellcenters.z(end)];
slice(X,Y,Z, phi.value, Sy, Sx, Sz);
xlabel('[y vlaues]'); % this is correct [matrix not rotated]
ylabel('[x vlaues]'); % this is correct [matrix not rotated]
zlabel('[z vlaues]');
axis equal tight
colorbar
end
function MS = createMesh3D(varargin)
% Modified from Ali Akbar Eftekhari
Nx=varargin{1};
Ny=varargin{2};
Nz=varargin{3};
Width=varargin{4};
Height=varargin{5};
Depth=varargin{6};
% cell size is dx
dx =0.125;
dy = 0.125;
dz = 0.125;
G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
CellSize.x= dx*ones(Nx+2,1);
CellSize.y= dy*ones(Ny+2,1);
CellSize.z= dz*ones(Nz+2,1);
CellLocation.x= [1:Nx]'*dx;
CellLocation.y= [1:Ny]'*dy;
CellLocation.z= [1:Nz]'*dz;
FaceLocation.x= [0:Nx]'*dx;
FaceLocation.y= [0:Ny]'*dy;
FaceLocation.z= [0:Nz]'*dz;
c=G([1,end], [1,end], [1, end]);
e1=G([1, end], [1, end], 2:Nz+1);
e2=G([1, end], 2:Ny+1, [1, end]);
e3=G(2:Nx+1, [1, end], [1, end]);
MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]);
end
% Used for 3D plotting
classdef HW7
% From 2012-2016 Ali Akbar Eftekhari
properties
dimension
dims
cellsize
cellcenters
facecenters
corners
edges
end
methods
function meshVar = Mesh(dimension, dims, cellsize, ...
cellcenters, facecenters, corners, edges)
if nargin>0
meshVar.dimension = dimension;
meshVar.dims = dims;
meshVar.cellsize = cellsize;
meshVar.cellcenters = cellcenters;
meshVar.facecenters = facecenters;
meshVar.corners= corners;
meshVar.edges= edges;
end
end
end
end

8 Comments
Show 6 older comments Hide 6 older comments

Torsten
Torsten on 26 Apr 2022
Seems there is still some error in your code.
You'd better focus on the numerical solution method rather than the graphical representation of the results.
%CP7
%Written by Viridiana Salazar
%clear variables
%close all
%% 1. Space discretization
Lx = 1.0; dx =0.125; N=Lx/dx+1; x=0:dx:Lx;
Ly = 1.0; dy =0.025; y=0:dy:Ly;
Lz = 1.0; dz =0.025; z=0:dz:Lz;
%% 2. Time discretization
tf =0.5; dt =0.001; M=tf/dt+1; t= 0:dt:tf;
%% Constants
Mu=0.5; r=Mu*dt/(dx)^2; q=-0.5*dt/dx;
%% Analytical Solution (We're gonna take from here the initial and
% boundary conditions)
UA = zeros(N,N,N,M);
for n=1:M %time
for i=1:N %x
for j=1:N %y
for k=1:N %z
UA(i,j,k,n)=(1+exp(x(i)/(2*Mu)+y(j)/(2*Mu)+z(k)/(2*Mu)-3*t(n)/(4*Mu)))^(-1);
end
end
end
end
%% Boundary Conditions
U = zeros(N,N,N,M);
% For x=1 and x=N
U(1,:,:,:)=UA(1,:,:,:); U(N,:,:,:)=UA(N,:,:,:);
%For y=1 and y=N
U(:,1,:,:)=UA(:,1,:,:); U(:,N,:,:)=UA(:,N,:,:);
%For z=1 and z=N
U(:,:,1,:)=UA(:,:,1,:); U(:,:,N,:)=UA(:,:,N,:);
%% Initial conditions
% For t=0
U(:,:,:,1)=UA(:,:,:,1);
%%
for n=1:M-1 %Time loop
% Initial guest at level n+1
UAn = reshape(UA(:,:,:,n),[N,N,N]);
Un = reshape(U(:,:,:,n),[N,N,N]);
disp('Error between analytical and numerical solution:')
norm(UAn(:)-Un(:))
for i=1:N %x
for j=1:N %y
for k=1:N %z
if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N %At the boundaries 1 and N, U is known
U(i,j,k,n+1)=U(i,j,k,n);
end
end
end
end
iteration=0;
EPS=1;
while EPS>0.0000000001 && iteration<20
%Jacobian Matrix
J=zeros(N^3,N^3);
% Vector f
F=zeros(N^3,1);
i=0;
for I=1:N^3
K=I;
% Fix the index i,j,k
if mod(I,N*N)==1
i=i+1;
j=0;
end
if mod(I,N)==1
k=1;
j=j+1;
else
k=k+1;
end
if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
F(K,1)=(q*(-3*U(i,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j,k-1,n+1)+U(i,j,k-1,n+1))+(-(6*r+1)*U(i,j,k,n+1)+r*(U(i+1,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j+1,k,n+1)+U(i,j-1,k,n+1)+U(i,j,k+1,n+1)+U(i,j,k-1,n+1)))+U(i,j,k,n));
end
%% derivatives for i=N and/or j=N and/or k=N
if i==N && j==N && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1);
%DF/U(i,j,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i==N && j~=N && k~=N
if i==N && j~=N && k~=N && j~=1 && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j==1 && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j~=1 && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j==N && k~=N
if i~=N && j==N && k~=N && i~=1 && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1);
%DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k~=N && i==1 && k~=1
J(K,I-1)=q*+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j==N && k~=N && i~=1 && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j~=N && k==N
if i~=N && j~=N && k==N && i~=1 && j~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1);
%DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j~=N && k==N && i==1 && j~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j~=N && k==N && i~=1 && j==1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j==N && k~=N
if i==N && j==N && k~=N && k~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1);
%DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j==N && k~=N && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j~=N && k==N
if i==N && j~=N && k==N && j~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q*+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k==N && j==1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i~=N && j==N && k==N
if i~=N && j==N && k==N && i~=1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k==N && i==1
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%% derivatives for i=1 and/or j=1 and/or k=1
if i==1 && j==1 && k==1
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i==1 && j~=1 && k~=1 && j~=N && k~=N
if i==1 && j~=1 && k~=1 && j~=N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j==N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j~=N && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i~=1 && j==1 && k~=1
if i~=1 && j==1 && k~=1 && i~=N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i==N && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i~=N && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=1 && j~=1 && k==1
if i~=1 && j~=1 && k==1 && i~=N && j~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i==N && j~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i~=N && j==N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i==1 && j==1 && k~=1
if i==1 && j==1 && k~=1 && k~=N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j==1 && k~=1 && k==N
J(K,I-1)=q+r; %DF/U(i,j,k-1)
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i==1 && j~=1 && k==1 && j~=N
if i==1 && j~=1 && k==1 && j~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k==1 && j==N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I-N)=q+r; %DF/U(i,j-1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i~=1 && j==1 && k==
if i~=1 && j==1 && k==1 && i~=N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I+N*N)=r; %DF/U(i+1,j,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k==1 && i==N
J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
J(K,I+1)=r; %DF/U(i,j,k+1)
J(K,I+N)=r; %DF/U(i,j+1,k)
J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
end
% Solving the system of equations J*dU=-F, we can use LU decomposition to save
% computational resources in robust system
dU=(J\(-F));
EPS=max(abs(dU(:,1)));
a=0;
for i=1:N %x
for j=1:N %y
for k=1:N %z
a=a+1;
if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
U(i,j,k,n+1)=U(i,j,k,n+1)+dU(a);
end
end
end
end
iteration=iteration+1;
end
end
%% Plots and subrotines for 3D plotting
%m = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the mesh
%Numerical=CellVariable(m, U(:,:,:,M));
%Analytical=CellVariable(m, UA(:,:,:,M));
%Error=CellVariable(m,abs(U(:,:,:,M)-UA(:,:,:,M)));
%figure
%visualizeCells3D(Numerical);
%title('Numerical solution at t=0.25')
%figure
%visualizeCells3D(Analytical);
%title('Analytical solution at t=0.25')
%figure
%visualizeCells3D(Error);
%title('Absolute Error at t=0.25')
end
function visualizeCells3D(phi)
% Modified from 2012-2016 Ali Akbar Eftekhari
phi.value = phi.value(2:end-1,2:end-1,2:end-1);
[X,Y,Z]=meshgrid(phi.domain.cellcenters.y, phi.domain.cellcenters.x, ...
phi.domain.cellcenters.z);
phi.value(1)=phi.value(1)+eps; % to avoid an strange error for assigning color limits
Sx = [phi.domain.cellcenters.x(1) phi.domain.cellcenters.x(end)];
Sy = [phi.domain.cellcenters.y(1) phi.domain.cellcenters.y(end)];
Sz = [phi.domain.cellcenters.z(1) phi.domain.cellcenters.z(end)];
slice(X,Y,Z, phi.value, Sy, Sx, Sz);
xlabel('[y vlaues]'); % this is correct [matrix not rotated]
ylabel('[x vlaues]'); % this is correct [matrix not rotated]
zlabel('[z vlaues]');
axis equal tight
colorbar
end
function MS = createMesh3D(varargin)
% Modified from Ali Akbar Eftekhari
Nx=varargin{1};
Ny=varargin{2};
Nz=varargin{3};
Width=varargin{4};
Height=varargin{5};
Depth=varargin{6};
% cell size is dx
dx =0.125;
dy = 0.125;
dz = 0.125;
G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
CellSize.x= dx*ones(Nx+2,1);
CellSize.y= dy*ones(Ny+2,1);
CellSize.z= dz*ones(Nz+2,1);
CellLocation.x= [1:Nx]'*dx;
CellLocation.y= [1:Ny]'*dy;
CellLocation.z= [1:Nz]'*dz;
FaceLocation.x= [0:Nx]'*dx;
FaceLocation.y= [0:Ny]'*dy;
FaceLocation.z= [0:Nz]'*dz;
c=G([1,end], [1,end], [1, end]);
e1=G([1, end], [1, end], 2:Nz+1);
e2=G([1, end], 2:Ny+1, [1, end]);
e3=G(2:Nx+1, [1, end], [1, end]);
MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]);
end
% Used for 3D plotting
%classdef HW7
% % From 2012-2016 Ali Akbar Eftekhari
% properties
% dimension
% dims
% cellsize
% cellcenters
% facecenters
% corners
% edges
% end
% methods
% function meshVar = Mesh(dimension, dims, cellsize, ...
% cellcenters, facecenters, corners, edges)
% if nargin>0
% meshVar.dimension = dimension;
% meshVar.dims = dims;
% meshVar.cellsize = cellsize;
% meshVar.cellcenters = cellcenters;
% meshVar.facecenters = facecenters;
% meshVar.corners= corners;
% meshVar.edges= edges;
% end
% end
% end
%end
Hana Bachi
Hana Bachi on 27 Apr 2022
I listed above the errors I got after running the code. what are the other errors so I can fix them?
Torsten
Torsten on 27 Apr 2022
Edited: Torsten on 27 Apr 2022
what are the other errors so I can fix them?
I don't know. I just ran your code without graphics and found that the error between analytical solution and numerical approximation grows from timestep to timestep.
For your assignment, the graphical routines of the code are not necessary. According to your description, you should only compute the 2-norm of the difference between numerical and analytical solution at t = 0.5.
Hana Bachi
Hana Bachi on 1 May 2022
While trying to run the code to get the 3D plots, I'm getting this error message:
Unrecognized function or variable 'CellVariable'.
Error in Cp72nd (line 344)
Numerical=CellVariable(m, U(:,:,:,M));
How can this be fixed?
the code I'm using as a separate .m file is the following
classdef MeshStructure
%CELLVARIABLE Summary of this class goes here
% Detailed explanation goes here
% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file
properties
dimension
dims
cellsize
cellcenters
facecenters
corners
edges
end
methods
function meshVar = MeshStructure(dimension, dims, cellsize, ...
cellcenters, facecenters, corners, edges)
if nargin>0
meshVar.dimension = dimension;
meshVar.dims = dims;
meshVar.cellsize = cellsize;
meshVar.cellcenters = cellcenters;
meshVar.facecenters = facecenters;
meshVar.corners= corners;
meshVar.edges= edges;
end
end
end
end
Torsten
Torsten on 1 May 2022
Edited: Torsten on 1 May 2022
You copied an incomplete code.
The necessary functions are here:
https://de.mathworks.com/matlabcentral/fileexchange/46637-a-simple-finite-volume-solver-for-matlab
But all you need for visualization can be done using MATLAB's "slice" command:
https://de.mathworks.com/help/matlab/ref/slice.html
Mary Barber
Mary Barber on 6 Sep 2023
By discretizing the equation in both time and space, we can construct a system of equations that can be solved iteratively to obtain the solution at each time step. This method is particularly useful for cases where stability is a concern, but it comes with a higher computational cost https://de.mathworks.com/matlabcentral/fileexchange/46637-a-simple-finite-volume-solver-for-matlab/uno online
Alex
Alex on 23 Jan 2025
Discretization is a powerful tool for solving complex problems in many fields such as physics, engineering, and computer science. However, it also comes with a high computational cost. The choice of grid size and method of solving the system of equations will directly affect the accuracy and efficiency of the results.https://de.mathworks.com/matlabcentral/fileexchange/46637-a-simple-finite-volume-solver-for-matlab/ heardle

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on Electromagnetics in Help Center and File Exchange

Tags

Asked:

Hana Bachi
Hana Bachi
on 24 Apr 2022

Commented:

Alex
Alex
on 23 Jan 2025

Community Treasure Hunt

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

Start Hunting!