Case Studies [plot_hydrogeophysical_modeling.py] Report an error

[bravegr]

Problem description

Thank you very much for reading my question, I am a student studying convective dispersion and solute transport motion of contaminants in soil and groundwater. I want to use PyGIMLI to do finite element calculations of convective dispersion equations and solute transport equations in soil and groundwater module, only forward, and calculate the concentration of contaminants in soil and groundwater without inversion. So I found the case: plot_hydrogeophysical_modeling.py. I downloaded it and ran it with pycharm, and it prompted SparseEfficiencyWarning: splu requires CSC matrix format warn('splu requires CSC matrix format', SparseEfficiencyWarning) I would like to ask you two questions: plot_hydrogeophysical_modeling (1).zip

1. Is it appropriate for me to use PyGIMLI and do 2D and 3D orthorectified modeling in the model plot_hydrogeophysical_modeling.py?
2. What is the effect of this warning and how to remove it?

Translated with www.DeepL.com/Translator (free version)

Your environment

Operating system: Windows, Python version: 3.8 pyGIMLi version:1.2.4 Way of installation: Conda package

Steps to reproduce

Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error: C:\Users\zly_w\AppData\Roaming\Python\Python38\site-packages\scipy\sparse\linalg\dsolve\linsolve.py:318: SparseEfficiencyWarning: splu requires CSC matrix format warn('splu requires CSC matrix format', SparseEfficiencyWarning)

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.ert as ert
import pygimli.physics.petro as petro
from pygimli.physics import ERTManager
world = mt.createWorld(start=[-20, 0], end=[20, -16], layers=[-2, -8],
                       worldMarker=False)

# Create a heterogeneous block
block = mt.createRectangle(start=[-6, -3.5], end=[6, -6.0], marker=4,
                           boundaryMarker=10, area=0.1)

# Merge geometrical entities
geom = world + block
pg.show(geom, boundaryMarker=True)
# Create a mesh from the geometry definition
mesh = mt.createMesh(geom, quality=32, area=0.2, smooth=[1, 10])
pg.show(mesh)
kMap = [[1, 1e-8], [2, 5e-3], [3, 1e-4], [4, 8e-4]]
# Map conductivity value per region to each cell in the given mesh
K = pg.solver.parseMapToCellArray(kMap, mesh)
pg.show(mesh, data=K, label='Hydraulic conductivity $K$ in m$/$s', cMin=1e-5,
        cMax=1e-2, logScale=True, grid=True)
# Dirichlet conditions for hydraulic potential
left = 0.75
right = 0.0
pBound = {"Dirichlet": {1: left, 2: left, 3: left, 5: right, 6: right, 7: right}}
# Solve for hydraulic potential
p = pg.solver.solveFiniteElements(mesh, a=K, bc=pBound)
# Solve velocity as gradient of hydraulic potential
vel = -pg.solver.grad(mesh, p) * np.asarray([K, K, K]).T
ax, _ = pg.show(mesh, data=pg.abs(vel) * 1000, cMin=0.05, cMax=0.15,
                label='Velocity $v$ in mm$/$s', cMap='YlGnBu', hold=True)
ax, _ = pg.show(mesh, data=vel, ax=ax, color='k', linewidth=0.8, dropTol=1e-5,
                hold=True)
S = pg.Vector(mesh.cellCount(), 0.0)
# Fill injection source vector for a fixed injection position
sourceCell = mesh.findCell([-19.1, -4.6])
S[sourceCell.id()] = 1.0 / sourceCell.size()  # g/(l s)
# Choose 800 time steps for 6 days in seconds
t = pg.utils.grange(0, 6 * 24 * 3600, n=800)
# Create dispersitivity, depending on the absolute velocity
dispersion = pg.abs(vel) * 1e-2
# Solve for injection time, but we need velocities on cell nodes
veln = mt.cellDataToNodeData(mesh, vel)
c1 = pg.solver.solveFiniteVolume(mesh, a=dispersion, f=S, vel=veln, times=t,
                                 scheme='PS', verbose=0)
# Solve without injection starting with last result
c2 = pg.solver.solveFiniteVolume(mesh, a=dispersion, f=0, vel=veln, u0=c1[-1],
                                 times=t, scheme='PS', verbose=0)
# Stack results together
c = np.vstack((c1, c2))
# We can now visualize the result:
for ci in c[1:][::200]:
    pg.show(mesh, data=ci * 0.001, cMin=0, cMax=3, cMap="magma_r",
            label="Concentration c in $g/l$")
ertScheme = ert.createERTData(pg.utils.grange(-20, 20, dx=1.0), schemeName='dd')

meshERT = mt.createParaMesh(ertScheme, quality=33, paraMaxCellSize=0.2,
                            boundaryMaxCellSize=50, smooth=[1, 2])
pg.show(meshERT)
# Select 10 time frame to simulate ERT data
timesERT = pg.IVector(np.floor(np.linspace(0, len(c) - 1, 10)))
# Create conductivity of fluid for salt concentration $c$
sigmaFluid = c[timesERT] * 0.1 + 0.01
# Calculate bulk resistivity based on Archie's Law
resBulk = petro.resistivityArchie(rFluid=1. / sigmaFluid, porosity=0.3, m=1.3,
                                  mesh=mesh, meshI=meshERT, fill=1)
# apply background resistivity model
rho0 = np.zeros(meshERT.cellCount()) + 1000.
for c in meshERT.cells():
    if c.center()[1] < -8:
        rho0[c.id()] = 150.
    elif c.center()[1] < -2:
        rho0[c.id()] = 500.
resis = pg.Matrix(resBulk)
for i, rbI in enumerate(resBulk):
    resis[i] = 1. / ((1. / rbI) + 1. / rho0)
# Initialize and call the ERT manager for electrical simulation:
ERT = ERTManager(verbose=False)
# Run  simulation for  the apparent resistivities
rhoa = ERT.simulate(meshERT, res=resis, scheme=ertScheme, 
                    returnArray=True, verbose=False)
for i in range(4):
    ERT.showData(ertScheme, vals=rhoa[i]/rhoa[0], cMin=1e-5, cMax=1)
### Expected behavior

Tell us what should happen or what you want to achieve.

### Actual behavior

Tell us what happens instead and/or provide the output of your script.

Paste output here.

> If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

#274 

[carsten-forty2]

Hello,

to 1.: Dispersion/Diffusion and Advection/Convection is solved by a very basic Finite Volume formulation tested for 2D meshes, 3D might work or not, but is untested. Note, there is a high numerical dispersion which can affect your results if you want to calculate with small diffusion coefficents. Higher Order Finite Element Solution for these problem classes are on TODO but we can not predict when they come.

to 2.: The Warning just says that the scipy sparse solver will be a little faster or use less memory if the internal matrix format would be in CSC instead of the used CSR. The results are okay so the warning can be ignored. A fix for this warning is already on TODO.

[bravegr]

Hello, I found after a few days of research about above case that the permeability coefficient kMap takes a value of 1e-4 even smaller,the area's flow line disappears，the final concentration diffusion is also not obvious.I would like to ask for advice， is this what you said in 1：there is a high numerical dispersion which can affect my results if I want to calculate with small diffusion coefficents?

In reality, the permeability coefficient of soil and groundwater is generally [1e-4,1e-7] m/s, so that they cannot be calculated using the finite element + finite volume method described above case? Do you have any good advice? I would be grateful if you could reply to me.

[halbmy]

We should rather use the term hydraulic conductivity (in m/s) as permeability is in m^2.

What do you mean with

the area's flow line disappears

Is it just a visualization problem? As far as I remember, it should work well (if an according discretization is used) for much smaller values.

At any rate, it would be interesting to find the limits of the solver and how numerical dispersion depends on the mesh.

[bravegr]

OK.Thank you very much, I will continue to work hard.

[carsten-forty2]

When the "flow lines disappears", then probably the "dropTol" keyword argument act as "Don't draw stream lines with velocity lower than drop tolerance."

[bravegr]

Okay, thank you very much for your answer, I probably understand.

[florian-wagner]

I am closing this, since it seems all questions have been answered.