Importing forward modeling results and performing ERT simulations

[ghost]

Dear pygimli developers,

I have gone through the documentation of pygimli, including the 2017 (Computers and Geosciences) paper and online help, but I could not find how to implement inverse modeling in pygimli using forward modeling results from a different code. I have forward modeling results for flow simulations from some other code that I would like to use to perform inverse ERT simulations in pygimli. Specifically, the process is similar to as described in section 4.3.3 of Wagner (2016) dissertation, which involves using changes in saturation to changes in the electrical bulk resistivity of the medium. It would be helpful if someone can explain how I can do this in pygimli. I would appreciate the help.

Best regards, L4student

[halbmy]

Maybe I am not getting it but do you basically mean how you can use pyGIMLi to do inversion using your own forward modelling code. There are inversion tutorials with very simple own-written forward classes. If it comes to coupling ERT with transport modelling, just have a look at the hydrogeophysical inversion example of the CG paper. It does not matter whether one of the objects or all of them are not pyGIMLi based, just the main frame (the class providing a response function) needs to be derived from the modelling base class. Does this answer your question? (Less an issue, more for the Q&A part)

[ghost]

Dear Thomas,

Thank you for trying to address my question. I mean to say that if I want to use forward modeling results from another code (let's say Dumux or some other fluid flow forward modeling simulator), then what is the appropriate way in pygimli (or bert) to import those results and use them in inverse ERT simulations? This is similar to the process described in section 4.3.3 of Wagner (2016) dissertation, except that I am not trying to optimize my parameter(s). If it's not clear yet then please let me know and I will be glad to clarify it further.

Thanks.

[ghost]

Dear @halbmy ,

This was the response I got earlier for a closed issue from 3 days ago. Although @florian-wagner suggested having this capability (and he implemented that in Chapter 4 of his dissertation), I can't find anything related to importing forward modeling results from another code that can be used to perform ERT inverse simulations on those imported results.

for a task as complex as CO2 migration, I would definitely recommend to use a benchmarked multi-phase reservoir simulator that offers the appropriate constitutive relations and physical processes that might play a role in the CO2 plume behavior (thermal, mechanical, geochemical). An open-source code that has been successfully applied for CO2 reservoir simulation is http://www.dumux.org/ for example. But there are others as well.

Once you have x, y, z, time, co2 saturation data, it is straightforward to use this for electrical resistivity simulations using pygimli. You might also want to have a look at Chapter 4 of my disseration (www.diss.fwagner.info), where CO2 saturations from the commercial ECLIPSE reservoir simulator were transferred to electrical resistivity using an experimentally established petrophysical relation for 3D-ERT simulations with pygimli.

Best wishes Florian

[Coastal0]

For what it's worth, I'm also in the midst of a similar problem, trying to convert FEFlow meshes to a style where BERT/Gimli can read and model the electrical responses.

[halbmy]

I see, combining pyGIMLi with other codes such as flow and transport modelling seems to be a great topic many people are working on. Technically, one has (in the response function) to generate input files, call the external modelling and read in the output, with Python it is easier. However, mostly processes work on very different scales so that for inversion some sort of interpolation (or pilot-point strategy) must be followed as done by Florian or Carsten. There is no ready solution but we are working on it. I assign it to @florian-wagner.

Importing FEflow meshes and results can be helpful, so we appreciate any input from you. Also other modelling codes (modflow, seawat, OGS, HGS) should be able to be worked with.

[ghost]

Thanks for deciding that adding this capability may be useful for many people.

There is one open-source e4d code that is integrated with pflotran (a transport and flow code) to perform ERT simulations, but it only allows forward modeling results for a tracer. I am interested in performing ERT on CO2 migration results (like Chapter 4 in @florian-wagner dissertation), so I thought pygimli/bert may be useful in case there is a capability to import the forward modeling results from some open-source codes like pflotran.

[florian-wagner]

Hi @L4student,

pygimli/BERT is definitely well suited for this task, no matter if you want to realize a fully-coupled inversion or synthetic ERT simultations based on previously simulated flow and transport scenarios.

However, it is a bit difficult to discuss importing, if you have not decided yet, which transport simulator you will be using for the CO2 simulations. In any case, you would have to parse the output file with Python and transform it to something like x y z saturation temperature, etc. Saturation (and potentially temperature or pressure) would need to be transformed to electrical properties using appropriate petrophysical relations.

If you have that, than you can easily interpolate x y z electrical_resistivity to a gimli mesh for forward simulations. A relevant tutorial is the mesh interpolation one (https://www.pygimli.org/_tutorials_auto/1_basics/plot_5-mesh_interpolation.html#sphx-glr-tutorials-auto-1-basics-plot-5-mesh-interpolation-py) for example.

If you have your CO2/resistivity distribution interpolated to a mesh, you can perform ERT modelling as presented in the paper examples: https://www.pygimli.org/_downloads/example-2_modelling.py

So two things: pygimli is well suited for your task and we are happy to help, yet we do not provide ready-made importers for all the transport simulators out there, so manual transformation of the specific output format is required.

Hope that answers your question.

Best wishes Florian

[ghost]

Hi @florian-wagner

I can use any of the available options for open-source flow/transport simulators that would make this task easy. I am using pflotran, but I didn't think it was important to mention, so I didn't specify that yet. It outputs results in 5 different formats, e.g. in vtk, hdf5, Tecplot etc. If you are planning on adding the import capability then I would recommend adding this transport simulator as an option. But, if you are not making changes in the pygimli code regarding import options then I would see if I can work around the problem using your suggested steps.

Thanks.

[carsten-forty2]

Hi L4student,

I don't think we plan to write a special importer for different flow/transport simulator until we need them for our research, ourself. Writing such an importer is a task that usually only have success if you need the results. ;)

However, we have a VTK importer that might not be complete but fits our needs until today. You can try it and report issues so we can try to improve it. I think a simple pg.load(VTKfileName) will give some insights.

Alternatively you can take a look at pygimli.meshtools.convertHDF5Mesh(...), maybe you can adopt this code for your needs.

If you find a stable solution it would by nice if you allow us to incorporate it into pygimli.

Cheers, Carsten

[ghost]

Thanks, @carsten-forty2. I appreciate your suggestions. I will try them and let you know if I have something to add or if I face any issues.

Cheers :)

[ghost]

@carsten-forty2: just a quick request. Can you point me to any help/example/doc that could allow me to look at what is returned by the functions you mentioned? I was looking at the code /io/load.py, but it's not apparent what the returned output is. I know I can keep looking through the repository to find that out, but it would be more convenient if I can take your experienced help on where to find appropriate help for the use of these specific functions. Thanks.

[carsten-forty2]

it should be a:

https://www.pygimli.org/gimliapi/classGIMLI_1_1Mesh.html

in the case its 2D you can directly take a look at with

mesh = pg.load(VTKFile)
print(mesh)
print(mesh.dataInfo())
pg.show(mesh)
# or
pg.show(mesh, mesh.data(DataName))

VTK supports Data arrays as well es the Mesh class does .. but I'm currently unsure if its already implemented by us.

[ghost]

@carsten-forty2 My mesh is 3D, and importing my mesh vtk file as you suggested didn't give anything:

print(mesh)
<class 'pygimli.core._pygimli_.RVector'> 0

print(mesh.dataInfo())
Traceback (most recent call last):

  File "<ipython-input-19-39eba9c878f3>", line 1, in <module>
    print(mesh.dataInfo())

AttributeError: 'RVector' object has no attribute 'dataInfo'

@florian-wagner: I do have x y z saturation (and other variables if necessary) parsed into python, but I am not sure how to adapt them to the pygimli's mesh object structure so they can be used in the ERT modeling as desired. I tried using pg.load(VTKFile) command for the mesh, but it didn't produce any result as shown above. Besides the mesh, I also would like to have the saturation (and other variables) read into pygimli data structure.

Thanks, L4student

[carsten-forty2]

seems the load auto-detection fails .. he selects regarding the file suffix. Your vtk file should have a .vtk suffix.

Btw: I'm not sure if 3D meshes actually work with the default ERT Manager based on pygimli ... this manager is mainly for 2D demonstration purposes. For advanced 3D Problems the pyBERT package is preferable. https://gitlab.com/resistivity-net/bert (also open source but still behind a registry wall)

[ghost]

@carsten-forty2: Yes, my vtk file has .vtk suffix. I created a gitlab account with @L4student username. Please allow me in to use pyBERT for 3D modeling as suggested. Thanks.

[carsten-forty2]

just to be sure .. you tried the following:

mesh = pg.load('yourfile.vtk')
print(mesh)

[ghost]

Yes.

[carsten-forty2]

hmm .. weird .. I guess a hidden error/problem then. I think need to see the file to find the problem .. If you can provide it .. I could take a look

[ghost]

I have created a small test file...please find it attached and the visualization in paraview of the same. mesh.tar.gz

[carsten-forty2]

ahh ... its a 'STRUCTURED_GRID' Dataset .. this is not yet implemented since our internal mesh container is a cell collection.

Currently we only support:

"UNSTRUCTURED_GRID" and "POLYDATA"

We can try to implement it .. should be no big task .. but its part of the c++core and it will need a little time to the next binary release.

The mesh does not contain any data so I only can implement the mesh structure and only hope the vtk file specification is modular enough for the data import.

edit: will open a new issue

[ghost]

@carsten-forty2: You suggested that the support to read the vtk file I sent earlier was added. However, I am still getting an error reading the same file I attached with my message earlier:

mesh = pg.load('mesh.vtk', verbose=True)

Traceback (most recent call last):
  File "/anaconda2/envs/pyg/lib/python3.6/site-packages/pygimli/io/load.py", line 113, in load
    return import_routines[suffix](fname)
RuntimeError: /home/wagner/miniconda3/conda-bld/pygimli_1506089658849/work/gimli/gimli/src/mesh_io.cpp: 1023            void GIMLI::Mesh::importVTK(const string&)  not yet implemented
 libgimli-0.9.7
Please send the messages above, the commandline and all necessary data to the author.
/home/wagner/miniconda3/conda-bld/pygimli_1506089658849/work/gimli/gimli/src/mesh_io.cpp: 1023          void GIMLI::Mesh::importVTK(const string&)  not yet implemented
 libgimli-0.9.7
Please send the messages above, the commandline and all necessary data to the author.
File extension .vtk seems to be not correct. Trying auto-detect.

@florian-wagner: The error suggests informing you.

[halbmy]

Carsten wrote

but its part of the c++core and it will need a little time to the next binary release.

so it is not automatically there, only after we provide the next conda package.

@florian-wagner can you please inform us all when it can be tested?

[Coastal0]

Hi Guys,

The example given by Florian earlier does not currently work.

I think this has been brought up before though...

If you have your CO2/resistivity distribution interpolated to a mesh, you can perform ERT modelling as presented in the paper examples: https://www.pygimli.org/_downloads/example-2_modelling.py

runfile('C:/Users/264401k/Downloads/example-2_modelling.py', wdir='C:/Users/264401k/Downloads')
Solve Darcy equation ... 
F:\WinPython-64bit-3.6.2.0Qt5\python-3.6.2.amd64\lib\site-packages\matplotlib\backend_bases.py:2453: MatplotlibDeprecationWarning: Using default event loop until function specific to this GUI is implemented
  warnings.warn(str, mplDeprecation)
Solve Advection-diffusion equation ...
Traceback (most recent call last):

  File "<ipython-input-43-2709a427e2a4>", line 1, in <module>
    runfile('C:/Users/264401k/Downloads/example-2_modelling.py', wdir='C:/Users/264401k/Downloads')

  File "F:\WinPython-64bit-3.6.2.0Qt5\python-3.6.2.amd64\lib\site-packages\spyder\utils\site\sitecustomize.py", line 710, in runfile
    execfile(filename, namespace)

  File "F:\WinPython-64bit-3.6.2.0Qt5\python-3.6.2.amd64\lib\site-packages\spyder\utils\site\sitecustomize.py", line 101, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/264401k/Downloads/example-2_modelling.py", line 62, in <module>
    scheme='PS', verbose=0)

  File "F:\WinPython-64bit-3.6.2.0Qt5\python-3.6.2.amd64\lib\site-packages\pygimli\solver\solverFiniteVolume.py", line 667, in solveFiniteVolume
    userData=kwargs.pop('userData', None))

  File "F:\WinPython-64bit-3.6.2.0Qt5\python-3.6.2.amd64\lib\site-packages\pygimli\solver\solverFiniteVolume.py", line 398, in diffusionConvectionKernel
    S = pg.RSparseMapMatrix(dof, dof, 0) + identity(dof, scale=b)

RuntimeError: C:/msys64/home/Guenther.T/py36/gimli/gimli/src/sparsematrix.h: 188        GIMLI::SparseMapMatrix<ValueType, IndexType>::SparseMapMatrix(const IndexArray&, const IndexArray&, const RVector&) [with ValueType = double; IndexType = long long unsigned int; GIMLI::IndexArray = GIMLI::Vector<long long unsigned int>; GIMLI::RVector = GIMLI::Vector<double>]  5637 != 0

[halbmy]

Thanks for reporting. This is obviously an error only occurring in the Windows version as the script runs nicely under Linux. I will create a new Windows version and check whether it persists. However, this is another issue so we leave this one closed.