PDE convection-diffusion equation using the fully implicit method.
8 views (last 30 days)
Show older comments
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
7 Comments
kmaonme
on 12 Sep 2022
thanks bro! More details for all
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
Answers (0)
See Also
Categories
Find more on Boundary Conditions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!