Different values ​​between set_rate_BC input and calculated rate value in the advection-diffusion algorithm

[avlmachado]

When running the transient advection-diffusion algorithm the defined rate in the boundary is not assured. I don't know if I'm making a mistake or it is a bug. Shouldn't the defined and calculated value be the same?

The test code is:

import openpnm as op 
import numpy as np 
np.random.seed(0) 

pn = op.network.Cubic(shape=[5,5,5], spacing=1e-3) 
geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts) 
water = op.phases.Water(network=pn) 
phys = op.physics.Standard(network=pn, phase=water, geometry=geo) 

sf = op.algorithms.StokesFlow(network=pn, phase=water) 
inlet=pn.pores('left') 
outlet=pn.pores('right') 
sf.set_value_BC(pores=inlet, values=1) 
sf.set_value_BC(pores=outlet, values=0.0) 
sf.run() 
water.update(sf.results()) 

mod = op.models.physics.ad_dif_conductance.ad_dif 
phys.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw') 
ad = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water) 
ad.set_rate_BC(pores=inlet, total_rate=1e-13) 
ad.set_outflow_BC(pores=outlet) 
ad.set_IC(0) 

ad.setup(t_scheme='implicit', t_final=50, t_output=1, t_step=10, t_tolerance=1e-12) 
ad.run() 

print(ad.rate(pores=inlet))

[ma-sadeghi]

Thanks for reporting this. This is almost certainly a bug. We'll investigate and report back here.

[mkaguer]

Hi @avlmachado. @ma-sadeghi asked me to look into this. I think the problem is that your system is not at steady state. Only at steady state will the net rate calculated be equal to the total_rate specified. It looks like your problem has a very long time scale. First try setting t_scheme to 'steady'. You should notice that the calculated rate now equals the rate specified. After that, try using really long t_final times. I tried t_final = 1000000s and got just 8.99155046e-14 so try something even larger! Or try increasing your pressure drop to increase flow rate and hence reduce the time scale of your problem. Does that help? If you don't understand why your problem needs to be at steady state I can clarify that as well.

[avlmachado]

Hi @mkaguer and @ma-sadeghi , thank you for the quick answer.

Indeed, testing higher injection velocities the rate values start to get close.

I agree that the concentration rate in the inlet and the outlet should be the same only in the steady-state. However, I think the Neumann boundary condition used in the inlet and the calculated value in the inlet should be the same at all times.

Checking the sum of ad['pore.bc_rate'][inlet] values, it is the same as defined, i.e., 1e-13. Therefore, I think the algorithm solution is correct, the problem is only when using the function "rate" to return the calculated value.

[mkaguer]

Hi @avlmachado. It's important to note that set_rate_BC is not a neumann boundary condition. I understand that can be confusing but the rate BC is better thought of as a source term rather than a Neumann boundary condition. The rate BC adds a constant rate of either consumption or generation to the specified pores. Does that help?

[ma-sadeghi]

@avlmachado One issue that I see with your simulation is that you're pumping mass into your domain from "left" pores, and forcing zero gradient at "right" pores, but don't have any source terms. Not sure if it's relevant though.

[ma-sadeghi]

@avlmachado I also checked with COMSOL (a commercial FEM solver), and indeed, when you impose a flux on a face, the flux changes over time in a transient simulation. This is because the flux right before x = 0 (call it J-) and the one right after x = 0 (call it J+) only match at steady state. In a transient simulation, they always differ by dc/dt, and so, when you ask for the rate via alg.rate(), it's always going to give you a smaller number, since part of the influx at x = 0 is contributing to dc/dt. Once steady state is reached, then alg.rate() will indeed match the flux that was originally set at x = 0.

[ma-sadeghi]

@avlmachado I'm closing this issue, feel free to reopen in case you still have questions.