CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

2D Finite Volume code in MATLAB

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 8, 2023, 13:01
Default 2D Finite Volume code in MATLAB
  #1
New Member
 
Chris
Join Date: Nov 2021
Posts: 12
Rep Power: 4
Jury is on a distinguished road
Hey everyone,


I'm trying to implement a finite volume code in MATLAB for a ellitical dgl -u'' = f in 2d. The 1d case worked out for me:



Code:
function unum = solvFVM(mesh, prob)
% This function probves an elliptical problem -u''(x) = f(x)

% We need a matriz of size:
nloc = mesh.N-1;
A = sparse(nloc,nloc); % One cell, one line
hri = mesh.hri;
% Here we need a loop.
ii = 1;
A(ii,ii  ) =  hri(ii) + 1./(mesh.xr(1)-mesh.xg(1));
A(ii,ii+1) = -hri(ii);
for ii = 2:nloc-1
    A(ii,ii-1) =           -hri(ii-1);
    A(ii,ii  ) =  hri(ii) + hri(ii-1);
    A(ii,ii+1) = -hri(ii);
end
ii = nloc;
A(ii,ii-1) =                                  -hri(ii-1);
A(ii,ii  ) = 1./(mesh.xg(end)-mesh.xr(end))  + hri(ii-1);

% spy(A)
% RHS

fi = mesh.f(mesh.xr) .* mesh.hg; % Valor de la funcion tamaņo intervalo
fi(1  ) = fi(1  ) - 1./(mesh.xr(1)-mesh.xg(1)) * prob.u0;
fi(end) = fi(end) - 1./(mesh.xg(end)-mesh.xr(end)) * prob.u1;

%unum = sparse(A)\fi';

unum = [prob.u0; -sparse(A)\fi'; prob.u1];
For the 2d case my code is as follows:



Code:
.rtcContent { padding: 30px; }.lineNode {font-size: 10pt; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-style: normal; font-weight: normal; }function unum = solvFVM2D(mesh, prob)
% This function probves an elliptical problem -u''(x) = f(x)

% We need a matriz of size:
nxr = mesh.Nx-1;
nyr = mesh.Ny-1;

nloc = (nxr)*(nyr);

A = sparse(nloc,nloc); % One cell, one line
fi = zeros(1, nloc);

% Here we need a loop.
for jj = 1:nyr
    for ii = 1:nxr
        fi(ii + (jj-1)*nxr) = mesh.f(mesh.xr(ii), mesh.yr(jj)) * mesh.hxg(ii) * mesh.hyg(jj);
        if ii > 1
            W = mesh.hyg(jj)/mesh.hxr(ii-1);
        else
            W = mesh.hyg(jj)/(mesh.xr(ii)-mesh.xg(ii));
            fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hyg(jj)./(mesh.xr(ii)-mesh.xg(ii)) * mesh.usol(mesh.xp(ii), mesh.yp(jj));
        end

        if ii <= nxr-1
            E = mesh.hyg(jj)/mesh.hxr(ii);
        else
            if ii == 1
                E = mesh.hyg(jj)/(mesh.xg(ii)-mesh.xr(ii));
                fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hyg(jj)./(mesh.xg(ii)-mesh.xr(ii)) * mesh.usol(mesh.xp(ii), mesh.yp(jj));
            elseif ii < nxr
                E = mesh.hyg(jj)/(mesh.xg(ii)-mesh.xr(ii));
                fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hyg(jj)./(mesh.xg(ii)-mesh.xr(ii)) * mesh.usol(mesh.xr(ii), mesh.yr(jj));
            else
                E = mesh.hyg(jj)/(mesh.xg(ii+1)-mesh.xr(ii));
                fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hyg(jj)./(mesh.xg(ii+1)-mesh.xr(ii)) * mesh.usol(mesh.xg(ii+1), mesh.yg(jj+1));
            end
        end

        if jj > 1
            S = mesh.hxg(ii)/mesh.hyr(jj-1);
        else
            S = mesh.hxg(ii)/(mesh.yr(jj)-mesh.yg(jj));
            fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hxg(ii)./(mesh.yr(jj)-mesh.yg(jj)) * mesh.usol(mesh.xp(ii), mesh.yp(jj));
        end

        if jj <= nyr-1
            N = mesh.hxg(ii)/mesh.hyr(jj);
        else
            if jj == 1
                N = mesh.hxg(ii)/(mesh.yg(jj)-mesh.yr(jj));
                fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hxg(ii)./(mesh.yg(jj)-mesh.yr(jj)) * mesh.usol(mesh.xp(ii), mesh.yp(jj));
            elseif jj < nyr
                N = mesh.hxg(ii)/(mesh.yg(jj)-mesh.yr(jj));
                fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hxg(ii)./(mesh.yg(jj)-mesh.yr(jj)) * mesh.usol(mesh.xg(ii), mesh.yr(jj));
            else
                N = mesh.hxg(ii)/(mesh.yg(jj+1)-mesh.yr(jj));
                fi(ii + (jj-1)*nxr) = fi(ii + (jj-1)*nxr) - mesh.hxg(ii)./(mesh.yg(jj+1)-mesh.yr(jj)) * mesh.usol(mesh.xg(ii+1), mesh.yg(jj+1));
            end
        end
        P = N + S + E + W;

        indx = ii + (jj-1) * nxr;

        A(indx, indx) = P;

        if indx < nloc
            A(indx, indx + 1) = -E;
        end
        if indx > 1
            A(indx, indx - 1) = -W;
        end
        if indx <= nloc-nxr
            A(indx, indx + nxr) = -N;
        end
        if indx > nxr
            A(indx, indx - nxr) = -S;
        end
    end
end

%%
spy(A)

%unum = [prob.u0; sparse(A)\fi'; prob.u1];
sol = -sparse(A)\fi';
sol_mat = reshape(sol, [nxr, nyr]);

unum = [mesh.usol(mesh.xg(1), mesh.yp(2:end-1))+ zeros(1, nyr); sol_mat; mesh.usol(mesh.xg(end), mesh.yp(2:end-1)) + zeros(1, nyr)];
unum = [mesh.usol(mesh.xp(:), mesh.yp(1))+zeros(nxr+2, 1), unum, mesh.usol(mesh.xp(:), mesh.yp(end))+zeros(nxr+2, 1)];
The code only works for a quasi 1D case (only one cell in x or y direction), but fails when I'm trying to increase the cells (see picture below).




Can anyone give me a hint what the error could be?


Thanks in advance and regards!


Chris
Jury is offline   Reply With Quote

Old   February 8, 2023, 13:06
Default
  #2
New Member
 
Chris
Join Date: Nov 2021
Posts: 12
Rep Power: 4
Jury is on a distinguished road
The files can be found here:



https://drive.google.com/drive/folde...6R?usp=sharing
Jury is offline   Reply With Quote

Reply

Tags
elliptic pde, finite volume method, matlab

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
how to set periodic boundary conditions Ganesh FLUENT 15 November 18, 2020 06:09
Looking for freelancer to write finite difference Matlab code mikmzo CFD Freelancers 3 September 25, 2018 14:30
I need a finite volume matlab code olalekan1986 CFD Freelancers 7 April 12, 2017 01:51
[blockMesh] Errors during blockMesh meshing Madeleine P. Vincent OpenFOAM Meshing & Mesh Conversion 51 May 30, 2016 10:51
channelFoam for a 3D pipe AlmostSurelyRob OpenFOAM 3 June 24, 2011 13:06


All times are GMT -4. The time now is 08:37.