Significant performance regression in transient advection diffusion from v2.8 to v3.0

[JiayanJI]

For the same transient advection diffusion problem, the solving speed of the model in version 3.0 is significantly slower than that of version 2.8.

The code for 2.8:

print('Initial boundary')
tad.set_IC(0)
print('Setup the transient algorithm settings')
tad.setup(t_scheme='cranknicolson', t_final=100, t_output=2, t_step=1, t_tolerance=1e-12)
tad.run()

The code for 3.0:

tspan = (0, 100)
saveat = 50
tad.run(x0=0, tspan=tspan, saveat=saveat)

Can you give me some advice?
Thank you very much!!!

[ma-sadeghi]

@mkaguer Could you please take a look?

[JiayanJI]

Hi, professer,
I further explored the code and found that the network was the main reason that affected the code operation.
I had tried different netwoeks and I found that when the spacing (paraments in the network) is small (1e-2), the code can not run. It shows the error.

When the spacing is big enough (1e-1), the code can run.
In the version 2.8, we can set the t_tolerance and t_step, I did not meet this problem.
I hope you can help me to solve this question.
Thank you very much!!!

[ma-sadeghi]

Hi @JiayanJI Please post a minimal working example (so I can just copy and run it on my machine) that produces the error.

[JiayanJI]

import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
print('Start Tuto_stepbystep')
print('Read libraries')
import numpy as np
import openpnm as op
import porespy as ps
import numpy as np
import openpnm as op
import porespy as ps
np.set_printoptions(precision=5)
pn = op.network.Cubic(shape=[10, 5, 3], spacing=1e-3)  # when we change the spacing to 1e-1, it will run.
print(pn)
geo =op.models.collections.geometry.spheres_and_cylinders
pn.add_model_collection(geo)
pn.regenerate_models()
print(pn)
print('Check health')
h = op.utils.check_network_health(network=pn)
print(h)
op.topotools.trim(pn, pores=h['disconnected_pores'])
h = op.utils.check_network_health(network=pn)
print(h)
print('Assign water phase and physics')
water = op.phase.Water(network=pn)
water['pore.diffusivity'] = 1e-9
water.add_model_collection(op.models.collections.phase.water)
water.add_model_collection(op.models.collections.physics.standard)
water.regenerate_models()
print('Performing Stokes flow')
sf = op.algorithms.StokesFlow(network=pn, phase=water)
sf.set_value_BC(pores=pn.pores('left'), values=50.0)
sf.set_value_BC(pores=pn.pores('right'), values=0)
sf.run()
water.update({'pore.pressure':sf['pore.pressure']})
fd = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water)
inlet  = pn.pores('left')
outlet = pn.pores('right')
fd.set_value_BC(pores=inlet, values=1.0)
f = op.models.physics.source_terms.power_law
water.add_model(propname='pore.reaction',
             model=f,
             X='pore.concentration',
             A1=1e-5,
             A2=1,
             A3=0.0,
             regen_mode='deferred')
fd.set_source(propname='pore.reaction', pores=outlet)
print(fd)
tspan = (0, 20)
saveat = 5
fd.run(x0=0, tspan=tspan, saveat=saveat)
print(fd)
print(water)
print('the value of outlet')
c = fd['pore.concentration']
print(c)

[ma-sadeghi]

@JiayanJI I spent some time with your code and I think the problem is not the spacing, but your boundary conditions: For Stokes flow, you set a pressure gradient, which induces a flow from left to right. However, for your transport algorithm, you only set the concentration at the inlet (plus a source term at the outlet). This implies that you have no-flux at the outlet, which to me seems to be unphysical (in OpenPNM, the default boundary condition is no-flux, so if you don't specify anything, it means nothing exists your boundaries).

In case it's not clear why this is probably unphysical, note that you have a net positive flow from left to right, but you don't allow any mass to exit the system (by not setting any boundary condition at the outlet for your transport algorithm). There are two options at the moment in OpenPNM which you can use: either use an "outflow" condition, which means no "gradient" (not necessarily no-flux), or just set a concentration boundary condition at the outlet.

One final note: whichever you choose, you can no longer add a source term to your outlet since it'll be considered a "boundary" face.

[JiayanJI]

Thank you very much for your help.
Follow your advice, I have changeed the boundary condition in outlet (outflow), and selected a pore as the source boundary.
But unfortunately, that problem still exists.
Here is my new code, thank you for your help:

import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
print('Start Tuto_stepbystep')
print('Read libraries')
import os
import pickle
import imageio
import scipy as sp
import numpy as np
import openpnm as op
import porespy as ps
import scipy.ndimage as spim
import matplotlib.pyplot as plt
from porespy.filters import find_peaks, trim_saddle_points, trim_nearby_peaks
from porespy.tools import randomize_colors
from skimage.segmentation import watershed
ps.visualization.set_mpl_style()
np.set_printoptions(precision=4)
np.random.seed(10)
np.set_printoptions(precision=5)
pn = op.network.Cubic(shape=[10, 5, 3], spacing=1e-1)  # when we change the spacing to 1e-1, it will run.
print(pn)
geo =op.models.collections.geometry.spheres_and_cylinders
pn.add_model_collection(geo)
pn.regenerate_models()
print(pn)
print('Check health')
h = op.utils.check_network_health(network=pn)
print(h)
op.topotools.trim(pn, pores=h['disconnected_pores'])
h = op.utils.check_network_health(network=pn)
print(h)

print('Assign water phase and physics')
water = op.phase.Water(network=pn)
water['pore.diffusivity'] = 1e-9
water.add_model_collection(op.models.collections.phase.water)
water.add_model_collection(op.models.collections.physics.standard)
water.regenerate_models()
print('Performing Stokes flow')
sf = op.algorithms.StokesFlow(network=pn, phase=water)
sf.set_value_BC(pores=pn.pores('left'), values=50.0)
sf.set_value_BC(pores=pn.pores('right'), values=0)
sf.run()
water.update({'pore.pressure':sf['pore.pressure']})
print('Performing Transient Advection Diffusion')
tad = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water)
inlet  = pn.pores('left')
outlet = pn.pores('right')
tad.set_value_BC(pores=inlet, values=1.0)
tad.set_outflow_BC(pores=outlet)
# add source term
pn.source_pores = np.arange(90)  # a pore not in the 'left' and 'right'  "new boundary condition"
f = op.models.physics.source_terms.power_law
water.add_model(propname='pore.reaction',
             model=f,
             X='pore.concentration',
             A1=1e-5,
             A2=1,
             A3=0.0,
             regen_mode='deferred')
tad.set_source(propname='pore.reaction', pores=pn.source_pores)
print(tad)
tspan = (0, 20)
saveat = 5
tad.run(x0=0, tspan=tspan, saveat=saveat)
print(tad)
print(water)
print('the value of outlet')
c = tad['pore.concentration']
print(c)

[ma-sadeghi]

@JiayanJI I think I have figured the problem. Your source term is very huge, causing the system of equations to become very stiff. Typically, reaction rate is per unit volume (in case it's homogeneous reaction), or per unit area (in case in heterogeneous, like catalytic reactions), but in either cases, your rate constant needs to be multiplied by the pore volume or area, which shrinks it down quite a bit, making the system of equations relatively stable. The reason you don't have this problem when using a spacing of 1e-1, is that (1e-1)**2 and (1e-1)**3 (for area and volume, respectively) are not drastically differen than 1e-1 itself, whereas for relatively small spacings, i.e., 1e-3, the area and volume are multiple orders of magnitudes smaller, hence the rest of the story!

Here's a minimal script that works:

import matplotlib.pyplot as plt
import numpy as np

import openpnm as op

np.random.seed(0)
op.visualization.set_mpl_style()

shape = [10, 5, 3]
pn = op.network.Cubic(shape=[10, 5, 3], spacing=1e-3)
geo = op.models.collections.geometry.spheres_and_cylinders
pn.add_model_collection(geo)
pn.regenerate_models()

h = op.utils.check_network_health(network=pn)
op.topotools.trim(pn, pores=h['disconnected_pores'])
h = op.utils.check_network_health(network=pn)

water = op.phase.Water(network=pn)
water['pore.diffusivity'] = 1e-9
water.add_model_collection(op.models.collections.phase.water)
water.add_model_collection(op.models.collections.physics.standard)
water.regenerate_models()

sf = op.algorithms.StokesFlow(network=pn, phase=water)
sf.set_value_BC(pores=pn.pores('left'), values=50.0)
sf.set_value_BC(pores=pn.pores('right'), values=0)
sf.run()

water.update({'pore.pressure': sf['pore.pressure']})
fd = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water)
inlet = pn.pores('left')
outlet = pn.pores('right')
fd.set_value_BC(pores=inlet, values=1.0)
fd.set_outflow_BC(pores=outlet)
f = op.models.physics.source_terms.power_law
k0 = 2.0
water.add_model(propname='pore.reaction',
                model=f,
                X='pore.concentration',
                A1=pn["pore.volume"]*k0,
                A2=1,
                A3=0.0,
                regen_mode='deferred')

non_boundary = pn.pores(["left", "right"], mode="not")
Ps_rxn = np.random.choice(non_boundary, size=90, replace=False)
fd.set_source(propname='pore.reaction', pores=Ps_rxn)

tspan = (0, 1)
saveat = None
# integrator = op.integrators.ScipyLSODA()
integrator = op.integrators.ScipyRK45()
# integrator = op.integrators.ScipyRadau()
fd.run(x0=0, tspan=tspan, saveat=saveat)

c = fd.soln['pore.concentration']

fig, ax = plt.subplots()

for t in np.linspace(*tspan, 5):
    c_x = c(t).reshape(shape).mean(axis=(1, 2))
    ax.plot(c_x, label=f't={t:.2f}')

ax.legend()

PS. Since you're using a "source" term, not a sink term, you see concentrations greater than 1.

[JiayanJI]

Thank you very much for your help.
This problems have been solved by your help.

[ma-sadeghi]

Thank you very much for your help. This problems have been solved by your help. Now, I meet a new problem. For the code you provided, if I replace the network with a new network (35,000 pores, 63026 throats), the calculation speed is very slow, even if I do not calculate the chemical reaction source term, only the transient advection diffusion. This process (transient advection diffusion) allows for fast calculations on version 2.8.

Could you let me know how slow are we talking about? (ex. how long does it take to solve transient diffusion @ `tspan = (0,1)` with and without source term?)

[JiayanJI]

Thanks for your attention.
I calculate my new network that with the boundary condition (stokesflow: left: rate, right: pressure). After 12 hours calculating, it still can not get the results. That without the source term, just the 'stokesflow' and 'Transient advection diffusion'.
In the version 2.8 (with the same network, phase, and physics), It will take me five minutes to get the result.
If I change the network to a simple cubic neitwork (Other Settings remain unchanged), It can get the results in one minutes.
Or, I only change the boundary condition of stokesflow, it also can get the results in one minutes.
Here is my code:

from pint import quantity
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
print('Start Tuto_stepbystep')
print('Read libraries')
import os
import pickle
import imageio
import scipy as sp
import numpy as np
import openpnm as op
import porespy as ps
import scipy.ndimage as spim
import matplotlib.pyplot as plt
from porespy.filters import find_peaks, trim_saddle_points, trim_nearby_peaks
from porespy.tools import randomize_colors
from skimage.segmentation import watershed
ps.visualization.set_mpl_style()
np.set_printoptions(precision=4)
np.random.seed(10)
op.visualization.set_mpl_style()

print('Read image')
path = 'D:/Marie file/Modelisation_copy/Modelisation/Test_PNM/Images/N03/'
file_format = '.tif'
file_name = 'N03_bin_AV_TRAV'
file = file_name + file_format
fetch_file = os.path.join(path, file)
im = imageio.mimread(fetch_file, memtest='2000MiB')
im = ~np.array(im, dtype=bool)

print('Calculate porosity')
print(ps.metrics.porosity(im))

print('Import Net')
# import dict
with open("N03_bin_AV_TRAV.net.pkl", "rb") as tf:
    net = pickle.load(tf)

print('Update network')
pn = op.network.Network()
pn.update(net)
prj = pn.project

print('Resolution_label')
resolution_label=0.01
op.topotools.label_faces(pn,resolution_label)
geo =op.models.collections.geometry.spheres_and_cylinders
pn.add_model_collection(geo)
pn.regenerate_models()
print(pn)

#####
#If I change the network to a cubic network, it is ok.
# shape = [10, 5, 3]
# pn = op.network.Cubic(shape=[10, 5, 3], spacing=1e-3)
# geo = op.models.collections.geometry.spheres_and_cylinders
# pn.add_model_collection(geo)
# pn.regenerate_models()
#
###

h = op.utils.check_network_health(network=pn)
op.topotools.trim(pn, pores=h['disconnected_pores'])
h = op.utils.check_network_health(network=pn)

water = op.phase.Water(network=pn)
water['pore.diffusivity'] = 1e-9
water.add_model_collection(op.models.collections.phase.water)
water.add_model_collection(op.models.collections.physics.standard)
water.regenerate_models()

sf = op.algorithms.StokesFlow(network=pn, phase=water)
sf.set_rate_BC(pores=pn.pores('left'), rates=7e-9) #rate boundary  （m^3/s）
sf.set_value_BC(pores=pn.pores('right'), values=101325)  #pressure boundary (pa)
####
# If I set this boundary condition, it also can get the results quickly.
#sf.set_value_BC(pores=pn.pores('left'), values=50.0)
#sf.set_value_BC(pores=pn.pores('right'), values=0)
#
#####
sf.run()
water.update({'pore.pressure': sf['pore.pressure']})
fd = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water)
inlet = pn.pores('left')
outlet = pn.pores('right')
fd.set_value_BC(pores=inlet, values=10.0)
fd.set_outflow_BC(pores=outlet)
tspan = (0, 1)
saveat = None
# integrator = op.integrators.ScipyLSODA()
integrator = op.integrators.ScipyRK45()
# integrator = op.integrators.ScipyRadau()
fd.run(x0=0, tspan=tspan, saveat=saveat)

c = fd.soln['pore.concentration']
print(c)

[JiayanJI]

files.zip
Here are my files, Thank you very much.

[ma-sadeghi]

Thanks for sharing the script. I'll look into it, but it has to wait as I'm busy with other stuff. I'll reply on this thread when I have something.

In the meantime, just a general comment: the old approach was not very robust, simply because the time step was specified by the user. The error of a transient solver is tightly linked with the time steps you choose, so if you set a fixed time step, you absolutely have no control over the error. In the new approach, the time step is automatically set by the solver itself to maintain a specified atol and rtol.

Long story short, just because the old solver "runs" doesn't mean that the results you're getting are accurate. So, I advise against using the old solver as you never know when your results are accurate and when they are not.

[JiayanJI]

OK, Thank you very much.
Look forward to your reply.

[JiayanJI]

Hi, professer,
I measured the calculation time of the code, and found that the problem of slow computing was not the occurrence of the source item, but the flow transpor equation.
Here is my test code（and the results are in the comments）：
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
print('Start Tuto_stepbystep')
print('Read libraries')
from timeit import default_timer as timer
tic = timer()
import os
import pickle
import imageio
import scipy as sp
import phreeqpython
import numpy as np
import openpnm as op
import porespy as ps
from scipy.integrate import odeint
import scipy.ndimage as spim
import matplotlib.pyplot as plt
from porespy.filters import find_peaks, trim_saddle_points, trim_nearby_peaks
from porespy.tools import randomize_colors
from skimage.segmentation import watershed

pp = phreeqpython()

pp = phreeqpython.PhreeqPython(database='phreeqc.dat')
ps.visualization.set_mpl_style()
np.set_printoptions(precision=5)
np.random.seed(10)
op.visualization.set_mpl_style()

np.random.seed(0)
op.visualization.set_mpl_style()

shape = [100, 50, 30]
pn = op.network.Cubic(shape=[100, 50, 6], spacing=1e-6)
geo = op.models.collections.geometry.spheres_and_cylinders
pn.add_model_collection(geo)
pn.regenerate_models()

h = op.utils.check_network_health(network=pn)
op.topotools.trim(pn, pores=h['disconnected_pores'])
h = op.utils.check_network_health(network=pn)

print('Assign water phase and physics')
water = op.phase.Water(network=pn)
water['pore.diffusivity'] = 1e-9
water.add_model_collection(op.models.collections.phase.water)
water.add_model_collection(op.models.collections.physics.standard)
water.regenerate_models()

print('Performing Stokes flow')
sf = op.algorithms.StokesFlow(network=pn, phase=water)
sf.set_value_BC(pores=pn.pores('left'), values=10) #Solution time: (10pa,1385.74s; 50pa,1508.75s; 100pa,1345.62s; 500pa,1495.70s;
#1000pa,1569.67s; 5000pa,1220.78s; 10000pa,1298.79s; 50000pa,1345.95s;100000pa,1362.59s;500000pa,2347.01)

sf.set_rate_BC(pores=pn.pores('left'), rates=7e-9) # rate boundary(m^3/s）:A day of calculations yielded no results

sf.set_value_BC(pores=pn.pores('right'), values=0)
sf.run()
water.update({'pore.pressure':sf['pore.pressure']})
print('Performing Transient Advection Diffusion')
mod = op.models.physics.ad_dif_conductance.ad_dif # calculate methods
water.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw')

ad = op.algorithms.AdvectionDiffusion(network=pn, phase=water) # 稳态

tad = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water) #瞬态
print('Adding boundary conditions')
inlet = pn.pores('left')
outlet = pn.pores('right')
tad.set_value_BC(pores=inlet, values=0.1249)
tad.set_outflow_BC(pores=outlet)

when the injection pressure is 10 pa, the solution time with source_terms is 1221.40s (without source_term =1385.74s)

we can know that the problem is not in the source_term

# add source term

g=op.models.physics.source_terms.general_symbolic # fast

water['pore.a'] = 0.00273

water['pore.b'] = -0.47923

water['pore.c'] = -0.000738433

water['pore.d'] = 0.15527

water['pore.x'] = np.random.random(water.Np)

y = 'ax**b + cx**d'

arg_map = {'a':'pore.a', 'b':'pore.b', 'c':'pore.c','d':'pore.d'}

water.add_model(propname='pore.general',

model=g,

eqn=y, x='pore.x', **arg_map)

non_boundary = pn.pores(["left", "right"], mode="not")

tad.set_source(propname='pore.reaction', pores=non_boundary)

print(tad)

tspan = (0,1)
saveat = 1
integrator = op.integrators.ScipyRK45()
tad.run(x0=0, tspan=tspan, saveat=saveat)

print(tad)
print('the value of outlet')
c = tad['pore.concentration']
print(c)
toc = timer()
messuer_time=toc-tic
print('messure time :', messuer_time)

op.io.project_to_xdmf(prj,'final_code')

[ma-sadeghi]

@JiayanJI Just wanted to update you. I was hoping to further look into this when I get the chance, but it seems that I'm only getting busier. Just wanted to let you know that I probably won't have the time. Sorry about that.