Diffusivity calculation in Version 2.3 produces unexpected/unreliable result, not consistent with Version 2.0.3

[taylrcawte]

When running the Gostick _Et _al.__ (2007) paper recreation on OpenPNM v. 2.3 the single phase diffusivity of air is greater than the bulk diffusivity which doesn't make physical sense. The same example was run on OpenPNM v. 2.0.3 and a different answer was produced. The difference in diffusivity between each version (using the exact same code) was 3 orders of magnitude (i.e. v. 2.0.3 D_eff = 10^-6, v. 2.3 D_eff = 10^-3).

[jgostick]

That does seem bad. We've obviously changed the way something is calculated and the example didn't notice. I'm guessing we have a #NBVAL_IGNORE_OUTPUT in the cell?

One bug we recently fixed is that re-using an algorithm within a for-loop does not work. It runs, but is erroneous since not all coefficients in the A matrix were being purged. Can you check the rate of material entering one face (ie. 'left') and leaving the other (i.e. 'right')? If this is the bug you're experiencing then these values won't match.

Also, I assume you're using the version pip? What do you get if you downloaded the 'dev' version from github?

[taylrcawte]

The same issue persists using the 'dev' version and when checking the rates the bc matrix shows all NaNs.

[jgostick]

That must mean the conductance values are not being calculated correctly. Do you see any error / warning messages? If not, have you changed the log level?

FYI, we are working on a massive revamp of the examples, and part of this will be to remove the paper recreations (at least for now) since they are huge and very difficult to maintain. Not exactly a solution, but it is an admission that I'm not surprised that they don't work as they should since we haven't been able to maintain them properly. We have enough trouble just keeping up with maintaining the actual code!

[taylrcawte]

There are no errors or warnings and the log levels have not been changed.

I am using the paper recreation as an example and have run into the same problem in my own work using the package. We first noticed the issue when playing around with boundary conditions, in v. 2.0.3 changing the boundary conditions did not have an impact on the diffusivity of the network while in version 2.3 when the boundary conditions were changed the impact on the diffusivity was noticeable.

[jgostick]

Let's simplify this problem down to bare basics to isolate the issue. When I run the following code, I get the same effective diffusivity no matter what boundary conditions I apply, which tells me that the solver and code are fine. The problem must be coming from your conductance values.

import openpnm as op

pn = op.network.Cubic(shape=[10, 10, 10], spacing=1e-4)
geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
air = op.phases.Air(network=pn, name='air')
phys_air = op.physics.Standard(network=pn, phase=air, geometry=geo)

fd = op.algorithms.FickianDiffusion(network=pn)
fd.setup(phase=air)
fd.set_value_BC(pores=pn.pores('left'), values=100)
fd.set_value_BC(pores=pn.pores('right'), values=0)
fd.run()
fd.calc_effective_diffusivity()

So as you add complexity, at which point does the problem start? Do the problems only occur when doing multiphase diffusion? If so how are you calculating the conductance of water filled pores? Also, generally, have you checked for negative throat lengths? Are there nans in the 'throat.diffusive_conductance'?

[raymond-guan1]

One bug we recently fixed is that re-using an algorithm within a for-loop does not work. It runs, but is erroneous since not all coefficients in the A matrix were being purged. Can you check the rate of material entering one face (ie. 'left') and leaving the other (i.e. 'right')? If this is the bug you're experiencing then these values won't match.

It appears that the reason for the issue that Taylr is encountering is similar to the one you mentioned. Running multiple algorithms results in erroneous values.

import openpnm as op

pn = op.network.Cubic(shape=[10, 10, 10], spacing=1e-4)
geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
air = op.phases.Air(network=pn, name='air')
phys_air = op.physics.Standard(network=pn, phase=air, geometry=geo)

fd = op.algorithms.FickianDiffusion(network=pn)
fd.setup(phase=air)
fd.set_value_BC(pores=pn.pores('left'), values=100)
fd.set_value_BC(pores=pn.pores('right'), values=0)
fd.run()
diff = fd.calc_effective_diffusivity()

import openpnm as op

pn = op.network.Cubic(shape=[10, 10, 10], spacing=1e-4)
proj = pn.project
geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
air = op.phases.Air(network=pn, name='air')
phys_air = op.physics.Standard(network=pn, phase=air, geometry=geo)

sf = op.algorithms.StokesFlow(network=pn)
sf.setup(phase=air)
sf.set_value_BC(pores=pn.pores('left'), values=100)
sf.set_value_BC(pores=pn.pores('right'), values=0)
sf.run()
sf.calc_effective_permeability()
proj.purge_object(sf)

fd = op.algorithms.FickianDiffusion(network=pn)
fd.setup(phase=air)
fd.set_value_BC(pores=pn.pores('left'), values=100)
fd.set_value_BC(pores=pn.pores('right'), values=0)
fd.run()
diff = fd.calc_effective_diffusivity()

In the first code-block, the simulation predicts an effective diffusivity that has a reasonable magnitude (9.58E-7), However, we run into issues in the second code block, where we run StokesFlow prior to FickianDiffusion. Despite purging the StokesFlow object, we obtain a completely different effective diffusivity (1.74E-4). Like you said, the A matrix is probably not being purged correctly in OpenPNM version 2.3.1.

[jgostick]

Wow, this is a very strange error. Will investigate immediately. Thanks for posting this code snippet. It's very helpful when an error can be reproduced in as few lines as possible.

[jgostick]

OK, this is not a bug, it's a 'feature' :-) In version 2.3 we added some automatic property regeneration, to support our reactive transport and non-linear solvers. We are working on non-newtonian flow, concentration dependent diffusion, and other fancy applications, so this feature was essential.

For example, if you have a viscosity model that is a function of temperature, and your algorithm somehow changes temperature, then the flow algorithm will automatically iterate (by updating the viscosity) until it converges on a stable solution. HOWEVER, in order to do this, the algorithm must write the new temperature to the PHASE object, so that all the other phase properties (and dependent physics properties) can "see" the new temperature and use it in the update.

So in your case, the StokesFlow algorithm was updating the air pressure, and this was causing all properties that depend on pressure to get updated, including the diffusion coefficient of air, which then impacted the diffusive conductance values.

This behavior caught me by surprize. It is one of the principles of object-oriented programming that data should be encapsulated, to prevent this exact surprize. The fact that the algorithm 'leaks' data back to the phase is an issue, and we have an open discussion about this (#1354). We are not exactly sure how to solve it, but we'll put it on the agenda for tomorrows dev telecon.

In the meantime, the solution is: don't reuse a phase unless you WANT the results of previous algorithms to impact simulations down the line. This is not the official solution (at least not yet), but will solve your current problems.