2D Finite Volume code in MATLAB

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 February 8, 2023, 13:01 2D Finite Volume code in MATLAB #1 New Member   Chris Join Date: Nov 2021 Posts: 12 Rep Power: 4 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

 February 8, 2023, 13:06 #2 New Member   Chris Join Date: Nov 2021 Posts: 12 Rep Power: 4 The files can be found here: https://drive.google.com/drive/folde...6R?usp=sharing

 Tags elliptic pde, finite volume method, matlab

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Ganesh FLUENT 15 November 18, 2020 06:09 mikmzo CFD Freelancers 3 September 25, 2018 14:30 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 AlmostSurelyRob OpenFOAM 3 June 24, 2011 13:06

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

 Contact Us - CFD Online - Privacy Statement - Top