# How to write 2D finite volume solver in a functional approach?

I'm solving 2D poisson problem of a finite volume problem using TDMA.

The solution of the problem is given:

\[Phi][x_, y_] := x^2 y + x y^2 + 1


The definitions are the next:

(*Input parameters*)
nx = 12; x0 = -1.0; xl = 1.0;
ny = 10; y0 = -1.0; yl = 1.0;
\[CapitalGamma] = 1.;

(*Meshing parameters*)
dx = (xl - x0)/nx;
dy = (yl - y0)/ny;
dz = 1.0;
dv = dx dy dz;
x = Range[x0, xl, dx];
y = Range[y0, yl, dy];
xc = (x + RotateRight[x])[[2 ;;]]/2.;
yc = (y + RotateRight[y])[[2 ;;]]/2.;

(*Source terms*)
g0y = yc - yc^2 + 1;
gly = yc + yc^2 + 1;
g0x = xc - xc^2 + 1;
glx = xc + xc^2 + 1;

(*Surface elements*)
Se = Outer[dy*dz &, x[[2 ;;]], yc];
Sw = Outer[dy*dz &, x[[;; -2]], yc];
Sn = Outer[dx*dz &, xc, y[[2 ;;]]];
Ss = Outer[dx*dz &, xc, y[[;; -2]]];

(*Coeffincients *)
aE = \[CapitalGamma]*Se/dx;
aW = \[CapitalGamma]*Sw/dx;
aN = \[CapitalGamma]*Sn/dy;
aS = \[CapitalGamma]*Ss/dy;
aP = aE + aW + aN + aS;
sP = Outer[(2 #1 + 2 #2)*dv &, xc, yc];

(*Boundary values*)
aP[[All, -1]] = aP[[All, -1]] + aN[[All, -1]];
sP[[All, -1]] = sP[[All, -1]] + 2 aN[[All, -1]]*glx;
aN[[All, -1]] = 0;

aP[[All, 1]] = aP[[All, 1]] + aS[[All, 1]];
sP[[All, 1]] = sP[[All, 1]] + 2 aS[[All, 1]]*g0x;
aS[[All, 1]] = 0;

aP[[-1, All]] = aP[[-1, All]] + aE[[-1, All]];
sP[[-1, All]] = sP[[-1, All]] + 2 aE[[-1, All]]*gly;
aE[[-1, All]] = 0;

aP[[1, All]] = aP[[1, All]] + aW[[1, All]];
sP[[1, All]] = sP[[1, All]] + 2 aW[[1, All]]*g0y;
aW[[1, All]] = 0;

(*Gauss TDMA function*)
tDMA[a_, b_, c_, d_] := Module[{m},
m = SparseArray[{Band[{1, 1}] -> b, Band[{2, 1}] -> a[[2 ;;]],
Band[{1, 2}] -> c[[;; -2]]}, {Length[b], Length[b]}];
LinearSolve[m, d, "Method" -> "Krylov"]
]


As using 2D i have to map my solution. Finite volumes book give me the algoritms by using loops as follows:

(*Gauss TDMA by lines for the x line*)
For[k = 0, k < 100, k++,
sP2 = sP + aN*RotateLeft[Tb, {0, 1}] + aS*RotateRight[Tb, {0, 1}];
For[j = 1, j <= ny, j++,
Tb[[All, j]] =
tDMA[-aW[[All, j]], aP[[All, j]], -aE[[All, j]], sP2[[All, j]]]
]
]

(*Gauss TDMA by lines for the y line*)
For[k = 0, k < 100, k++,
sP2 = aE*RotateLeft[Tb, {1, 0}] + aW*RotateRight[Tb, {1, 0}] + sP;
For[i = 1, i <= nx, i++,
Tb[[i, All]] =
tDMA[-aS[[i, All]], aP[[i, All]], -aN[[i, All]], sP2[[i, All]]]
]
]


I think this is not efficient if I use a large mesh. So I would like to write the loops in a functional approach.

I tried somthin like this

Tb = Outer[1 &, xc, yc];
Tb = Nest[
tDMA, {-aS, aP, -aN,
aE*RotateLeft[#, {1, 0}] + aW*RotateRight[#, {1, 0}] + sP}] &,
Tb, 20];


But it only works for the y line, not for the x line

For the xlines i tried trasposing all the aP,aE,aW,aN,aS, and sP arrays but also I think is not the most efficient solution.

(*Line X*)

s2 = Transpose[
tDMA, {-aWT, aPT, -aET,
aNT*RotateLeft[#, {1, 0}] + aST*RotateRight[#, {1, 0}] +
sPT}] &, Transpose[Tb], 20]]; // AbsoluteTiming


Where aWT, aPT, aET, aNT, aST, sPT are transposed aP,aE,aW,aN,aS, and sP.

• I would just change the For to Do loop. No need to do anything else. Clarity of the algorithm is the most important thing, not how short of code you can get it to be. That is why all algorithms in textbooks are written in normal procedural fashion. Mar 17, 2023 at 1:09
• 1. “I think this is not efficient…” If the main concern is speed, since LinearSolve already involves in definition of tDMA, I don't think optimizing the loop will help much. (It's better to avoid For loop to make the code cleaner, of course. ) 2. Can you add a bit more background info? 3. Where's the definition of Tb? 4. To make your last approach work for x line, you just need Transpose. Mar 17, 2023 at 3:36
• I alredy add the initial definition of Tb. I dont thin transpose helps very well, thinking in making then the solver 3D. I alredy used some aproaches, I edit the question with more details
– kpaz
Mar 17, 2023 at 16:03
• Many, many, many moons back I wrote an introduction to FVM with code, you find the introduction here. And the code is in the same repo. Nowadays, I'd do thinks very differently, but the theory did not change. There is one specially about this, that's that we use shape functions here - kind of a cross between FEM and FVM. Just thought I pass this along, maybe you can find something useful. Mar 17, 2023 at 16:30
• Here is a talk I gave not quite as many moons back. This uses the FEM BUT, you'd only need to replace the discretization and then you'd have a pretty efficient solver. It also shows the usage of Compile for generating an efficient element. Should not be too hard to replace the very small FEM code with a very small FVM code. Mar 17, 2023 at 16:36