Simulate apparent resistivities on 2-D anisotropic resistivity distributions?

[geoale90]

Problem description

Hello pyGIMLi team,

I am trying to use pygimli.physics.ert.simulate(mesh = mesh, scheme = scheme, res = resis) to simulate apparent resistivities in 2-D when resis is an anisotropic array of resistivities. However, there is a problem with the cell markers and I get an error.

The code I am using is the following:

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert as ert

# Model dimensions and discretization
xmin = -100
xmax = 100
ymin = -20
ymax = 0

nx = 100
ny= 10

# Create grid
mesh = pg.createGrid(x=np.linspace(xmin, xmax,nx + 1 ),
                       y=np.linspace(ymin, ymax, ny + 1))

# Create anisotropic input array 
resis = [None]*mesh.cellCount()
for c in mesh.cells():
    resis[c.id()] = pg.solver.createAnisotropyMatrix(lon=1.0, trans=0.1,
                                               theta=90/180 * np.pi)

elecs_surface = np.linspace(start = xmin, stop = xmax, num = 48)
surface_scheme = ert.createData(elecs= elecs_surface,schemeName='dd')

rhoa = ert.simulate(mesh = mesh, scheme = surface_scheme, res = resis)

The error message is the following:

Traceback (most recent call last):

File ".........py", line 58, in <module>
    rhoa = ert.simulate(mesh = mesh, scheme =scheme, res = resis)
File "C:\Software\BERT\pygimli\physics\ert\ert.py", line 138, in simulate
    if any([mark not in mapMarkers for mark in meshMarkers]):

  File "C:\Software\BERT\pygimli\physics\ert\ert.py", line 138, in <listcomp>
    if any([mark not in mapMarkers for mark in meshMarkers]):

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

The compiled resis of shape (1000,2,2) (one (diagonal) resistivity tensor for each pixel) works fine when I use pg.solver.solveFiniteElements(mesh, a=resis, bc= bc):

neumann_boundary = {1: 0.0, 2: 0.0}
dirichlet_boundary = {3: 1.0, 4: 0.0}
bc = {'Dirichlet': dirichlet_boundary ,'Neumann': neumann_boundary}

u =  pg.solver.solveFiniteElements(mesh, a=resis, bc= bc) 

Expected behavior

I would like to compute apparent resistivities coming from a 2-D anisotropic distribution of resistivities.

Actual behavior

I get an error with the cell markers, which is not occurring when I use pg.solver.solveFiniteElements(mesh, a=resis, bc= bc)

Thank you in advance. Best, Alejandro

[halbmy]

Seems like ert.simulate is not yet working with resistivity tensors (and even less inversion does), and actually nobody did a through study using anisotropy so far. So I suggest going back on step using the reference ERT modelling implementation demonstrated in https://www.pygimli.org/_examples_auto/3_dc_and_ip/plot_04_ert_2_5d_potential.html I would really love to see some studies with anisotropic resistivity, in comparison with examples from the literature, so that we could finally include it into the ERT manager or any derived class.

[halbmy]

I agree that anisotropic resistivity is a very interesting subject of research, particularly for the typical targets in inductive electromagnetics (see also packages like empymod or custEM). However, for ERT it is still an edge case. The ERT module and manager is designed for routine analysis with isotropic resistivity and will not include anisotropy in the near future. Such a development could be made in other repositories like BERT.

Independent on the ERT module (using a C++ core forward operator) it should be possible to use the generic finite-element solver as you pointed out, but this means for 2D you will have to use the wavenumber solution pointed out in https://www.pygimli.org/_examples_auto/3_dc_and_ip/plot_04_ert_2_5d_potential.html

[geoale90]

Hello Thomas, Thanks for your two replies. I thought I had replied to the first message.

Yes, after I finish some work with the inversion of scalar conductivities, I am planning to build a forward operator of apparent resistivity data for tensorial distribution of conductivities. For this, I would simply iterate the wavenumber solution that you point out over the different electrode positions. Then, I can first test the forward operator on some scalar distributions of conductivities and compare the results against ert.simulate. Then, my plan is to use this forward operator within a stochastic inversion framework (MCMC).

If interested/useful for the repository, I will provide updates on this work.

Merry Christmas

[halbmy]

Would be great to hear about updates from your work!