Flux calculation in Simple diffusion does not work!

Numerical methods and mathematical models of Elmer
jfaraudo
Posts: 17
Joined: 11 May 2018, 11:50
Antispam: Yes

Flux calculation in Simple diffusion does not work!

Post by jfaraudo »

Hi,

We have tried a very simple example of a very simple diffusion system (modified from an official ELMER test example) to check that ELMER correctly calculates diffusion fluxes, but we get very annoying results.

The example is simply a 2D domain with an adsorbing surface at the botom and constant concentration at the top and lateral periodic boundaries .
Looking at the development of concentration profile with time using paraview it looks as expected (a linear concentration profile develops).

But we asked ELMER to calculate and save to a data file magnitudes such as diffusion flux, total mass inside the system and convective flow.
All of them are completelly off. There is an unphysical convective flow in one of the surface (but velocity was imposed to be zero, convection was not calculated and NS equations were not solved!) which is even larger than the diffusive flux. The values of the diffusion flux are not what theory predicts. In the steady state the flux adsorbing on the bottom is different from the flux leaving the top surface. The total mass in the system (obtained by integration of the concentration) is two orders of magnitude off...

We have posted our mesh file, .sif file, some results and some discussion here:

https://github.com/jfaraudo/ELMER_EXAMP ... _diffusion

We are completelly upset by these results, since diffusion is present in many elmer applications.

Is there something wrong with the .sif file (which is based on an elmer example)? Any ideas ?
Any help will be much appreciated.

jordi
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Flux calculation in Simple diffusion does not work!

Post by mzenker »

Hi,

I am by no means a fluidics/duffusion expert and have never used that solver myself.
So just two very trivial things:

1. Did you check that the original Elmer test runs OK and gives valid results?
2. If I understand your sif file correctly, you have a time step of 20 seconds. That seems large to me. But indeed diffusion processes are not that quick. Nevertheless, does it change anything to use smaller timestep?

I hand over to the people with more expertise in that field... ;)

Matthias
jfaraudo
Posts: 17
Joined: 11 May 2018, 11:50
Antispam: Yes

Re: Flux calculation in Simple diffusion does not work!

Post by jfaraudo »

Yes, of course we did these tests.

1.- We reproduce all reported example data (but warning: the flux data produced in the example but not shown in the elmer wiki is also unphysical!).
2.- We tried many different steps, even reducing three orders of magnitude the time step we obtain exactly the same result (but with lots of useless data and large computer time). From my experience in diffusion and the time scales in this particular problem , this time scale has to be perfectly OK, this is what I used in my own home-made codes. In any case, the obtained steady-state concentration profile looks fine but the steady state calculated fluxes are wrong, so this looks completelly unrelated to time steps. Even more surprisingly the unphysical convection flux at the adsorbing surface (convection is not calculated, velocity was fixed to zero and NS equations not solved) is independent of the time step.

If you can forward these comments to anyone expert in the field it would be great.
We are also trying to reproduce published results with Elmer to check whether the problems are there or not.

Jordi
mzenker wrote:Hi,

I am by no means a fluidics/duffusion expert and have never used that solver myself.
So just two very trivial things:

1. Did you check that the original Elmer test runs OK and gives valid results?
2. If I understand your sif file correctly, you have a time step of 20 seconds. That seems large to me. But indeed diffusion processes are not that quick. Nevertheless, does it change anything to use smaller timestep?

I hand over to the people with more expertise in that field... ;)

Matthias
raback
Site Admin
Posts: 4828
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Flux calculation in Simple diffusion does not work!

Post by raback »

Hi Jordi,

I don't think there is anything wrong except for confusing documentation and funny defaults.

Convective flux actually uses the given "variable" as the convection and the given "coefficient" as the multiplier. Unfortunately when the variable is of wrong dimension (as here 1 instead of 2) it does not complain by error. The result is integration over given variable i.e. concentration over the BC which of course does not make sense. The intended use for this has been to check mass conservation when velocity is computed by N-S equation.

The diffusive flux is instantaneous flux. Hence if you want to compare it to the total concentration in the domain you need to integrate over it. In this case the integral is roughly "sum(20*result(:,4))" and compare this to "result(:,2)". I think you will find the results to be at least reasonably close.

The typical challenge for flux evaluation comes from poor mesh quality. Even here there can be compromised accuracy at the start when the is a sharp gradient. Increasing mesh density of setting element to p:2, for example (quadratic p-element) might help.

The most accurate way of computing fluxes is combination of "calculate loads" and "boundary sum".

It might be useful to have the convective flux to complain when called in a wrong way, and improve on its documentation. Things to do.

-Peter
jfaraudo
Posts: 17
Joined: 11 May 2018, 11:50
Antispam: Yes

Re: Flux calculation in Simple diffusion does not work!

Post by jfaraudo »

Sorry Peter, but many of your comments indicate that you didn't take the time to really look in detail into our comments and files.
We can simply forget and move to other software, but we do really want to understand what happens and to help if possible developers and other users.

The example that we are reporting is extremelly simple and a trivial steady state was reached in a very short time. Everything we argue is calculated in the steady state.

In the steady state the flux leaving the top surface has to be equal to the flux entering the bottom surface, but Elmer report different values in the steady state (so mass conservation is NOT verified).

The linear concentration profile in the steady steate is calculated correctly by Elmer, but the value of the steady state flux reported by Elmer is incorrect (this is constant over time so don't blame the time dependency!!). The diffusive flux in this steady state is homogeneous over all the space , so don't blame the discretization or integration problems.

The total mass is trivial to integrate with this very simple linear concentration profile reported by Elmer (so don't blame the grid) but the result of the integral given by Elmer is orders of magnitude wrong.

Also it is amazing to see that the program reports a large convective flow (4 orders of magnitude larger than the diffusive flow) when the fluid velocity reported by Elmer is zero.

Our preliminary study of published data shows that the problem is also in at least a published example in which Elmer was employed. In my opinion (as a scientist and also as member of editorial board of 2 journals) this sounds as a problem to be SERIOUSLY analyzed.

We give our .sif and grid files. If our files have errors, please indicate which ones. We were based on Elmer examples, so everything we did is what is supposed to be done. I think I know what I am doing, I wrote my first program to solve advection-diffusion equation in 1995 . Since that time, I have spend all my academic career doing simulations, using my own codes or codes written by other people in the scientific community. I know by experience that developing code is hard, and that it may involve bugs...

We were thinking as using Elmer as an alternative to calculations that we currently do with our own programs or with COMSOL, but given the problems that we are facing, I am not sure which way we will take.

We are really disappointed to see that nobody is taking this issue seriously, and instead of looking in detail into the calculations the only responses are "canned" responses such as blaming the time dependency, the spatial resolution, generic mentions of the user manual.

I challenge anybody to produce the simplest possible example of diffusion with Elmer , to demonstrate that the program is able to report fluxes and adsorbing mass in this very simple problem: a boundary with a fixed concentration and another boundary with zero concentration (absorbing boundary). This can be calculated analytically with very simple equations, so we can compare and test. Until that, for me it is pretty clear that Elmer does not work as intended. Developing code is hard and difficult, and as code grows it could be difficult to check for bugs. For this reason I think it is a great idea to test simple situations with known analyticall solutions, to help identify possible bugs.

Jordi
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Flux calculation in Simple diffusion does not work!

Post by mzenker »

Hi,

I dare to do another try to perhaps find a possible cause of the problem. I apologize in advance if this is perceived as yet another "canned", incompetent or otherwise disappointing or inappropriate response...
The problem seems to be related to quantities like mass and flux which only make sense in 3D. Since you do your calculation in 2D, you need an assumption for the 3rd dimension. From the 2D simulations I did, I remember that Elmer assumes this 3rd dimension to be 1 m. So this has to be taken into account when interpreting the results.

HTH,

Matthias
jfaraudo
Posts: 17
Joined: 11 May 2018, 11:50
Antispam: Yes

Re: Flux calculation in Simple diffusion does not work!

Post by jfaraudo »

Thanks Matthias for taking your time to help us. Your answer does not solve completelly the problem, but it helps a lot.
Now, the situation is as follows:

1) Your point explains the value of the total mass in the stationary state. Assuming this "extra dimension" (not described in the manual, as far as I see) the value calculated integrating the concentration is now correct. I have updated the example at github:
https://github.com/jfaraudo/ELMER_EXAMP ... _diffusion

2) The values of diffusive fluxes in the stationary state are close to the predictions of Fick's law but not identical (see the Github readme), the deviation seems significant for such a very simple geometry and maybe more importanly, in the stationary state, the -time constant- flux of mass leaving the surface at constant concentration is different from the -time constant- flux of mass entering the absorbing surface. The difference is not large, but nevertheless indicates problems with the mass balance.

3) According to what Peter said in the forum, the convective flow is incorrectly computed by Elmer in this particular case by some reason that I do not really understand (somehting like it is incorrectly computed because the calculation is incorrect, or something?).

Any further ideas or help???

In the meantime, we keep trying with different options (changing mesh resolution, etc).


Jordi
Last edited by jfaraudo on 12 Jun 2018, 15:05, edited 1 time in total.
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Flux calculation in Simple diffusion does not work!

Post by mzenker »

Hi Jordi,

so we are one step further... :)
The difference in diffusive flux between the two boundaries is much less than 0.1%, so I would not worry about that too much. I would expect it to decrease with a finer mesh, and maybe also with a smaller timestep or if you do a steady state simulation.

I also don't understand the factor of 2 in the diffusive flux. Just taking Fick's law as you find it in a textbook, I get the same result as you do, so either we both miss something or Elmer does something different in the diffusive flux calculation.
For the convective flux, I don't understand either what is happening.

Peter, could you shed some light on this?

Matthias
jfaraudo
Posts: 17
Joined: 11 May 2018, 11:50
Antispam: Yes

Re: Flux calculation in Simple diffusion does not work!

Post by jfaraudo »

Hi Matthias,

Thanks a lot!!!

I corrected the github data, there was a mess with the geometry of the system.
The corrected fluxes (see the github) are similar to those predicted by Fick's law, but still we have this 1% difference.
What really worries me is that the value of the flux in the steady state is different on top or on the bottom surface for such a very simple system. This means that if we use the outputs for further calculations, we have an incorrect mass balance.

Also, this problem with the convective flux is rather worrying. We need to be able to rely on the values reported by the convective fluxes. This is for me the biggest problem that remains so far. If Peter can help us it will be absolutely great.

We keep trying with changing the mesh or time steps, to see whether this value reduces or not. For us is problematic because we are interested in balances of fluxes, and a missbalance of flux between boundaries is a big issue for us.

We will report further soon, I hope. In the meantime, any ideas will be appreciated.

I think we are now very close!!

Jordi
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Flux calculation in Simple diffusion does not work!

Post by mzenker »

Indeed, I overlooked that the distance between the boundaries is 1 cm and not 5 mm... :oops:

The difference between the diffusive fluxes is 0.058%. FEM in general does not enforce mass conservation, so this small difference is not that surprising. Do you really need such a high precision in your results? And do other FEM tools with the same type of mesh give much better precision?

There is still the convective flux which remains to be understood.
Until Peter finds the time to reply on that matter, you might take a look into the code of the SaveData solver yourself to maybe understand what is going on. It is FORTAN and not that difficult to read for someone with coding and simulation knowledge...

HTH,

Matthias
Post Reply