[mkhm]

Dear foamers,



I am struggling to find a way to have a better insight of convergence of simulations done with rhoCentralFoam as the latter does not print residuals. I tried three ways but although it can be seen in the paraview that the case is converged, the residuals obtained by these methods do to reflect this fact. The methods being used so far are:

1) Compare the initial values of rho, rhoU and rhoE with their final calculated values for all cell points and take the maximum value of it; for instance for rho we have: mag(max(rho_f-rho_0)).value() where rho_f and rho_0 refere, respectively, to the final calculated value and initial value for rho. The same is done for rhoU and rhoE.
This was not a very helpful. I can imagine that if the calculation is converged for all the points except one, this method will fail as the maximum of difference could remain the same. So, it is enough that you fix the outlet waveTransmissive boundary condition for pressure to something wrong, the whole domain converges except the final point where the converged value is different from the fixed value.

2) Calculate the difference of calculated variables and their initial values for all cell points and sum the absolute values over all cells; for rho we have: sum(mag(rho_f-rho_0)).value(). Again this method can fail for the same reason as the one mentioned above. If I am not wrong, this method is what is called as unscaled residual in the link below.

3) Calculate the difference of calculated variables and their initial values for all cell points and sum the absolute values over all cells and divide the sum of absolute values per the sum of absolute calculated values; for rho we have: (sum(mag(rho_f-rho_0))/sum(mag(rho_f))).value(). Again this method can fail for the same reason as the one mentioned above. If I am not wrong, this method is what is called as unscaled residual in the link below.


link: https://www.afs.enea.it/project/nept...ug/node812.htm


Could you help me to figure it out how to solve this important issue ?

Thanks in advance,
Mary

[dlahaye]

1/ Do you mean mag(rho_f).value() / mag(rho_0).value() ?

(alternative order of brackets)

2/ Can you see residual norms on the screen? Capture them in log-file and plot with appropriate tools?

3/ Do you want to see residual fields in paraview?

Cheers, D.

[mkhm]

Originally Posted by dlahaye

1/ Do you mean mag(rho_f).value() / mag(rho_0).value() ?

(alternative order of brackets)

2/ Can you see residual norms on the screen? Capture them in log-file and plot with appropriate tools?

3/ Do you want to see residual fields in paraview?

Cheers, D.

Dear dlahaye,
1. I think whether you do mag(something_1).value() / mag(something_2).value() or (mag(something_1/something_2)).value(), they are both equivalent. However, I don't understand why this division should reflect the convergence of the solution.
2. Seeing residuals is not a problem. I am using gnuplot.
3. I don't know how to see the residuals in paraview. But, don't forget that here, there is no calculation of residual in the openFoam native solver of rhoCentralFoam. So first, we need to find an appropriate way to calculate residuals and after I'll learn how to do it with paraview.

[dlahaye]

Possibly we are lost in translation.

However, mag(vec1/vec2) <> mag(vec1) / mag(vec2).

Left-hand side involves a division of vectors. This operation is typically not defined.

Right-hand side involves a division of numbers.

D.

[usv001]

Hi Mary,

I am not sure why you're taking the difference between the initial and final values. It may be something particular to your application. My understanding of residuals is that you need to take the difference between subsequent time steps to determine how much difference there is in one time step. When your solution has converged, this should be very small.

So, you may try to replace this code in rhoCentralFoam:

solve(fvm::ddt(rho) + fvc::div(phi));

with this:

const volScalarField rho0(rho); // @ previous time step

solve(fvm::ddt(rho) + fvc::div(phi));

const scalar massError
(
    sum(mag(rho - rho0)*mesh.V()).value()
);

Info<< "Mass error: " << massError << endl;

to compute the mass error. That should give you an indication of the convergence of your results.

Hope this helps.

Cheers,
USV

[mkhm]

Originally Posted by usv001

Hi Mary,

I am not sure why you're taking the difference between the initial and final values. It may be something particular to your application. My understanding of residuals is that you need to take the difference between subsequent time steps to determine how much difference there is in one time step. When your solution has converged, this should be very small.

So, you may try to replace this code in rhoCentralFoam:

solve(fvm::ddt(rho) + fvc::div(phi));

with this:

const volScalarField rho0(rho); // @ previous time step

solve(fvm::ddt(rho) + fvc::div(phi));

const scalar massError
(
    sum(mag(rho - rho0)*mesh.V()).value()
);

Info<< "Mass error: " << massError << endl;

to compute the mass error. That should give you an indication of the convergence of your results.

Hope this helps.

Cheers,
USV

Dear USV,
I do consider what you call as the subsequent time steps as follows :

volScalarField rho_0("rho_0", rho);
solve(fvm::ddt(rho) + fvc::div(phi));
volScalarField rho_f("rho_f", rho);
scalar rho_Res = sum(mag(rho_f-rho_0)).value();

Is it correct what I do ? Why in your code, there is "*mesh.V()"? What does that mean ?
If it is correct what I am doing, I don't understand why it does not reflect the convergence of the simulation. If it is not correct, could you please indicate how I can correct what is wrong ?


Thanks in advance

[usv001]

Hi Mary,

"mesh.V()" is the cell volume. So, when multiplied by density, the residual has units of mass. I think you can try to divide it by "sum(mesh.V()).value()" to get back density.

By the way, what you're doing should also be valid as it measures the change between two time steps. Why do you say that it does not reflect convergence? It should if your solution has converged.

USV

[mkhm]

Originally Posted by usv001

Hi Mary,

"mesh.V()" is the cell volume. So, when multiplied by density, the residual has units of mass. I think you can try to divide it by "sum(mesh.V()).value()" to get back density.

By the way, what you're doing should also be valid as it measures the change between two time steps. Why do you say that it does not reflect convergence? It should if your solution has converged.

USV

Thanks USV for your reply. I attach the plot of residuals. My supervisor had told me that the residuals should usually decrease for 3 order of magnitude. You can see in the attached file that is not the case. However, I am sure that it is not within 10000 iterations that the case is converged, it should be after. And I am sure that the case is converged because if I plot results on three different axis of the geometry, nothing changes after a while. So I am pretty sure that the case is converged but the plot of residuals and more precisely my residual calculation does not reflect that.