[Jury]

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:

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:

.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]

The files can be found here:



https://drive.google.com/drive/folde...6R?usp=sharing