Question about cfdAssembleIntoGlobalMatrixFaceFluxes

[zongy17]

I am a graduate student studying uFVM with the book FVM in CFD, and I find something puzzle.

In function cfdAssembleIntoGlobalMatrixFaceFluxes (in the file src/assemble/cfdAssembleIntoGlobalMatrixFaceFluxes.m), lines include that

...
bc(own)                       = bc(own)                       - theFluxes.FluxTf(iFace);
...
bc(nei)                       = bc(nei)                       + theFluxes.FluxTf(iFace);
....
bc(own) = bc(own) - theFluxes.FluxTf(iBFace);
....

I suppose that bc(...) here represents the right hand side (RHS) of the algebraic equations (just like vector b of system Ax=b). But as the book writes in equation (8.80) in Chapter 8,

bc = Qc*Vc + non-orthogonal terms

where non-orthogonal terms should be FluxVf calculated in function cfdAssembleDiffusionTerm.

So I am wondering why this is not adding or subtracting the term theFluxes.FluxVf(...) but rather FluxTf(...) ? I notice that FluxTf(...) calculated in function cfdAssembleDiffusionTerm represents the total flux through the face, and could not understand why it should be included in bc(...).

[mhamadmahdialloush]

The assembly in uFVM is based on the residual form. Instead of solving the system Ax=b, we solve the system Ax'=b-Ax, so, in uFVM, b is in fact "b-Ax", which is the residual of the equation. Now, note that we use face fluxes (fluxCf, fluxFf, fluxVf, and fluxTf) to assemble the global system (Ax'=b-Ax*) for terms such as advection and diffusion, so we store fluxCf and fluxFf the coefficients multiplied by the variables (this appears right after discretisation and/or linearisation) while we store in fluxVf the remaing part. fluxTf is the total flux across a face. This is explained here next. Imagine this equation involving phi:

f(phi) = 0,

which can be linearised to be

a1phi_C+a2phi_F+a3 = 0

The fluxes are evaluated as follows: fluxCf = a1 fluxFf = a2 fluxVf = a3 fluxTf = a1phi_C+a2phi_F+a3

Given that fluxTf is the total flux across the face, we can make use of it to assemble the residual form system Ax'=b-Ax*. Thus, the global assembly will have the following: b -= fluxTf

[zongy17]

The assembly in uFVM is based on the residual form. Instead of solving the system Ax=b, we solve the system Ax'=b-Ax, so, in uFVM, b is in fact "b-Ax", which is the residual of the equation. Now, note that we use face fluxes (fluxCf, fluxFf, fluxVf, and fluxTf) to assemble the global system (Ax'=b-Ax*) for terms such as advection and diffusion, so we store fluxCf and fluxFf the coefficients multiplied by the variables (this appears right after discretisation and/or linearisation) while we store in fluxVf the remaing part. fluxTf is the total flux across a face. This is explained here next. Imagine this equation involving phi:

f(phi) = 0,

which can be linearised to be

a1_phi_C+a2_phi_F+a3 = 0

The fluxes are evaluated as follows: fluxCf = a1 fluxFf = a2 fluxVf = a3 fluxTf = a1_phi_C+a2_phi_F+a3

Given that fluxTf is the total flux across the face, we can make use of it to assemble the residual form system Ax'=b-Ax*. Thus, the global assembly will have the following: b -= fluxTf

Wow, I really appreciate your quick reply! It helps me a lot in understanding the codes! Thank you so much!

[shanjierenyidp]

Hi Alloush, if the bC is actually the residual of the equation. When you calculate the equation residual of the current iteration, shouldn't the residual be: res = bC - A*dphi However, in the code (uFVM/src/solve/cfdComputeResidualsArray) it says:

for iElement=1:numberOfElements    
    residual(iElement) = bc(iElement) - ac(iElement)*phi(iElement);    
    for nNeighbour=1:length(cconn{iElement})
        iNeighbour = cconn{iElement}(nNeighbour);
        residual(iElement) = residual(iElement) - anb{iElement}(nNeighbour)*phi(iNeighbour);
    end    
end

Which basically indicates: res = bC - A*phi

Here, the definition of bC goes back to RHS in the book. Could you help me understand this?

[mhamadmahdialloush]

Hi,

bC is already the residual of the equation in fact. In ufvm, we assemble equations in residual form, meaning that thr right-hand-side of the system (bC) is the residual of the system.

In the linear solver, which contains the function you mentioned (cfdComputeResidualsArray), the residual is the residual of the smoothed system (solved system by the specified solver i.e. DILU), we use this residual to evaluate the performance of the linear solver as well as for other objectives.

[shanjierenyidp]

Hi,

bC is already the residual of the equation in fact. In ufvm, we assemble equations in residual form, meaning that thr right-hand-side of the system (bC) is the residual of the system.

In the linear solver, which contains the function you mentioned (cfdComputeResidualsArray), the residual is the residual of the smoothed system (solved system by the specified solver i.e. DILU), we use this residual to evaluate the performance of the linear solver as well as for other objectives.

Thanks for you quick reply! I really appreciate that. Yes, you are right. But what i meant was that since the smooth solver is solving residual form denoted as Ax=b_C. The unknown vector (x) should be dphi. As I tested in the Matlab code, in every iteration, the initial guess of dphi=0. The smoother solver takes dphi=0 and solve the equation and updates dphi with some values (say dhpi'). So if you want to calculate the residual of the smoothed system, shouldn't it be: b_C-A dhpi' ? But the code didn't use dphi, it directly used phi. That's what I am confused about.

[mhamadmahdialloush]

dphi is the correction brought by the iterative solver. It is in fact the correction of the correction to phi. If you want to calculate the smoothed residual, you will need to add dphi first to x, which in other words:

Res = bC - A(x+dphi) = bC - Ax - Adphi

And again this is the smoothed system residual.

[shanjierenyidp]

dphi is the correction brought by the iterative solver. It is in fact the correction of the correction to phi. If you want to calculate the smoothed residual, you will need to add dphi first to x, which in other words:

Res = bC - A(x+dphi) = bC - Ax - Adphi

And again this is the smoothed system residual.

And here bC is Qc*Vc + non-orthogonal terms like zongy17 mentioned right?
Also I know why I was wrong in the code! I thought the varible phi we used in the function is the actual scalar field phi. But it is actually dphi itself. We fetched the updated dphi and named it phi.

function residual = cfdComputeResidualsArray(theCoefficients)
%--------------------------------------------------------------------------
%
%  Written by the CFD Group @ AUB, Fall 2018
%  Contact us at: cfd@aub.edu.lb
%==========================================================================
% Routine Description:
%   Compute residuals array: res = b_C - sum(a_F*phi_F)
%--------------------------------------------------------------------------

ac = theCoefficients.ac;
anb = theCoefficients.anb;
bc = theCoefficients.bc;
cconn = theCoefficients.cconn;
phi = theCoefficients.dphi;   <-----------------Here,  we fetched the updated dphi and named it phi. 

numberOfElements = length(ac);
residual = zeros(size(ac));
for iElement=1:numberOfElements    
    residual(iElement) = bc(iElement) - ac(iElement)*phi(iElement);    
    for nNeighbour=1:length(cconn{iElement})
        iNeighbour = cconn{iElement}(nNeighbour);
        residual(iElement) = residual(iElement) - anb{iElement}(nNeighbour)*phi(iNeighbour);
    end    
end

So this cfdComputeResidualsArray is calcuating the smoothed residual res = b-Ax-Adphi like you mentioned, where b stands for b =Qc*Vc + non-orthogonal terms. BTW, do you have any reference book or papers that show the residual form solving process? I didn't find it in the FVM book (The Finite Volume Method in Computational Fluid Dynamics by F. Moukalled). Also I want to thank you for writing up the matlab code. It saves so much time compared to reading the C++ code. Appreciate the work a lot. Do you mind add me on linkedin? I work in machine learning for fluid dynamics. Hopefully we can exchange ideas more in the future.

[MarcusAW]

This code has been great for helping me understand the inner workings of CFD solvers more thoroughly. Adding on from @shanjierenyidp , I am curious what the reasoning was for solving the residual form. The book indicates, "In this form, numerical errors during the solution of the equation are slightly less than those associated with the standard form for cases when small variations are expected for large values of phi." Unfortunately, I could not find more information in the literature and would also love to read more about it.

Thanks, Marcus

[mhamadmahdialloush]

with a residual form, you can ensure smaller steps to the solution (smaller values), and from a practical point of view, it also provides more flexibility at the discretisation level.

On Wed, Sep 4, 2024 at 7:28 PM Marcus Ang @.***> wrote:

This code has been great for helping me understand the inner workings of CFD solvers more thoroughly. Adding on from @shanjierenyidp https://github.com/shanjierenyidp , I am curious what the reasoning was for solving the residual form. The book indicates, "In this form, numerical errors during the solution of the equation are slightly less than those associated with the standard form for cases when small variations are expected for large values of phi." Unfortunately, I could not find more information in the literature and would also love to read more about it.

Thanks, Marcus

— Reply to this email directly, view it on GitHub https://github.com/mhamadmahdialloush/uFVM/issues/6#issuecomment-2329508763, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFLISFJDFPRKSIXH6VVLYJLZU4YJXAVCNFSM6AAAAABNUWE5XOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMRZGUYDQNZWGM . You are receiving this because you commented.Message ID: @.***>