Time-lapse 3D survey
To invert time-lapse survey there are two modes: reg_mode == 1 or
reg_mode == 2 (see R3t documentation). - reg_mode == 1 is
refered to as ‘background constrained’ inversion where the subsequente
survey are penalized from departure from the reference survey. -
reg_mode == 2 (default in ResIPy) is refered to as ‘difference
inversion’ but, for 3D, this requires to modify the protocol.dat to
compute \(d - d_0 + f(m_0)\) where \(d_0\) is the observed
transfer resitance of the reference survey, \(d\) is the observed
transfer resistances of the survey to be inverted, \(f(m_0)\)
represent the forward response (in terms of transfer resistances) of the
inverted reference model. While R2 and cR2 do d-d0+f(m0) internally,
R3t does not. But ResIPy has this build-in, hence no special initiative
is required from the user if you want to use reg_mode == 2. You will
see the --- Computing d-d0+f(m0) ---- when calling
Project.invert().
import warnings
warnings.filterwarnings('ignore')
import os
import sys
sys.path.append((os.path.relpath('../src'))) # add here the relative path of the API folder
testdir = '../src/examples/'
from resipy import Project
import pyvista as pv # for 3D plots
import numpy as np # to specify np.infty
API path = /media/jkl/data/phd/tmp/resipy/src/resipy
ResIPy version = 3.2.0
cR2.exe found and up to date.
R3t.exe found and up to date.
cR3t.exe found and up to date.
# difference inversion - see d-d0+f(m0) in the log
k = Project(typ='R3t')
k.createTimeLapseSurvey([testdir + 'dc-3d-timelapse-protocol/data/protocol3D-1.dat',
testdir + 'dc-3d-timelapse-protocol/data/protocol3D-2.dat'],
ftype='ProtocolDC')
k.importElec(testdir + 'dc-3d-timelapse-protocol/elec/electrodes3D-1.csv')
k.createMesh()
k.invert()
# by default `reg_mode == 2` (difference inversion)
print(k.param['reg_mode'])
should_run_async will not call transform_cell automatically in the future. Please pass the result to transformed_cell argument and any exception that happen during thetransform in preprocessing_exc_tuple in IPython 7.17 and above.
Working directory is: /media/jkl/data/phd/tmp/resipy/src/resipy
clearing dirname
2/2 imported
Creating tetrahedral mesh...fmd in gmshWrap.py: 16.343534
writing .geo to file completed, save location:
/media/jkl/data/phd/tmp/resipy/src/resipy/invdir
Reading mesh3d.msh
Gmsh version == 3.x
reading node coordinates...
Determining element type...Tetrahedra
Reading connection matrix...
ignoring 7581 elements in the mesh file, as they are not required for R2/R3t
Finished reading .msh file
interpolating topography onto mesh using triangulate interpolation...done
Done
ResIPy Estimated RAM usage = 0.164347 Gb
done
Writing .in file and protocol.dat... Matching quadrupoles between surveys for difference inversion...644 in common...done in 0.014881s
done!
------------ INVERTING REFERENCE SURVEY ---------------
>> R 3 t E R T M o d e l v 2.30 <<
>> Date: 19-01-2021
>> My beautiful 3D survey
>> I n v e r s e S o l u t i o n S e l e c t e d <<
>> A d v a n c e d M e s h I n p u t <<
>> T e t r a h e d r a l E l e m e n t M e s h <<
>> Reading mesh file
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading resistivity model from res0.dat
>> L o g - D a t a I n v e r s i o n <<
>> N o r m a l R e g u l a r i s a t i o n <<
>> Memory estimates:
For 1000 measurements the memory needed is: 0.250 Gb
For 2000 measurements the memory needed is: 0.491 Gb
For 5000 measurements the memory needed is: 1.213 Gb
For 10000 measurements the memory needed is: 2.417 Gb
>> Forming roughness matrix
>> Number of measurements read: 644
>> Total Memory required is: 0.164 Gb
Processing frame 1 - output to file f001.dat
Iteration 1
Initial RMS Misfit: 2.41 Number of data ignored: 0
Alpha: 67.701 RMS Misfit: 0.54 Roughness: 3.208
Step length set to 1.000
Final RMS Misfit: 0.54
Final RMS Misfit: 1.01
Solution converged - Outputing results to file
Calculating sensitivity map
End of data: Terminating
>> Program ended normally
----------------- Computing d-d0+f(m0) ---------------
Writing .in file and mesh.dat... done!
Writing protocol.dat... done!
Running forward model...
>> R 3 t E R T M o d e l v 2.30 <<
>> Date: 19-01-2021
>> My beautiful 3D survey
>> F o r w a r d S o l u t i o n S e l e c t e d <<
>> A d v a n c e d M e s h I n p u t <<
>> T e t r a h e d r a l E l e m e n t M e s h <<
>> Reading mesh file
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading resistivity model from resistivity.dat
>> Memory estimates:
For 1000 measurements the memory needed is: 0.006 Gb
For 2000 measurements the memory needed is: 0.006 Gb
For 5000 measurements the memory needed is: 0.006 Gb
For 10000 measurements the memory needed is: 0.006 Gb
>> Number of measurements read: 644
>> Total Memory required is: 0.006 Gb
>> Program ended normally
0/644 reciprocal measurements found.
0/644 reciprocal measurements found.
Forward modelling done.Matching quadrupoles between surveys for difference inversion...644 in common...done in 0.020808s
--------------------- MAIN INVERSION ------------------
>> R 3 t E R T M o d e l v 2.30 <<
>> Date: 19-01-2021
>> My beautiful 3D survey
>> I n v e r s e S o l u t i o n S e l e c t e d <<
>> A d v a n c e d M e s h I n p u t <<
>> T e t r a h e d r a l E l e m e n t M e s h <<
>> Reading mesh file
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading resistivity model from Start_res.dat
>> L o g - D a t a I n v e r s i o n <<
>> D i f f e r e n c e R e g u l a r i s a t i o n <<
>> Memory estimates:
For 1000 measurements the memory needed is: 0.250 Gb
For 2000 measurements the memory needed is: 0.491 Gb
For 5000 measurements the memory needed is: 1.213 Gb
For 10000 measurements the memory needed is: 2.417 Gb
>> Forming roughness matrix
>> Number of measurements read: 644
>> Total Memory required is: 0.164 Gb
Processing frame 1 - output to file f001.dat
Iteration 1
Initial RMS Misfit: 1.40 Number of data ignored: 0
Alpha: 69.072 RMS Misfit: 0.16 Roughness: 0.871
Step length set to 1.000
Final RMS Misfit: 0.16
Final RMS Misfit: 1.00
Solution converged - Outputing results to file
Calculating sensitivity map
End of data: Terminating
>> Program ended normally
2/2 results parsed (2 ok; 0 failed)
2
k.showResults(index=1, pvslices=[[10,20,30,40],[0],[]], attr='difference(percent)',
color_map='bwr', vmin=-10, vmax=10)
should_run_async will not call transform_cell automatically in the future. Please pass the result to transformed_cell argument and any exception that happen during thetransform in preprocessing_exc_tuple in IPython 7.17 and above. Converting np.character to a dtype is deprecated. The current result is np.dtype(np.str_) which is not strictly correct. Note that np.character is generally deprecated and 'S1' should be used. Converting np.character to a dtype is deprecated. The current result is np.dtype(np.str_) which is not strictly correct. Note that np.character is generally deprecated and 'S1' should be used.
In background constrained inversion (reg_mode == 1), the
difference(percent) attribute is not computed by R3t. In order to
have it computed by ResIPy, it is needed to tell R3t to output the full
mesh (so that ResIPy can match the mesh elements between reference and
other surveys). To do this set Project.param['num_xy_poly'] = 0 and
set zmin and zmax to \(-\infty\) and \(+\infty\)
respectively, just before Project.invert(). If you don’t change
that, R3t will crop the mesh by default and ResIPy will output a message
saying it couldn’t compute ‘difference(percent)’.
Note that setting num_xy_poly to 0 and changing zmin/zmax
does not change the vizualisation. It just changes the .vtk and .dat
files produced by R3t.
# background constrained - see d-d0+f(m0) in the log
k = Project(typ='R3t')
k.createTimeLapseSurvey([testdir + 'dc-3d-timelapse-protocol/data/protocol3D-1.dat',
testdir + 'dc-3d-timelapse-protocol/data/protocol3D-2.dat'],
ftype='ProtocolDC')
k.importElec(testdir + 'dc-3d-timelapse-protocol/elec/electrodes3D-1.csv')
k.createMesh()
k.param['reg_mode'] = 1
k.param['num_xy_poly'] = 0 # tells R3t to no crop the mesh
k.param['zmin'] = -np.inf
k.param['zmax'] = np.inf
k.invert()
Working directory is: /media/jkl/data/phd/tmp/resipy/src/resipy
clearing dirname
2/2 imported
Creating tetrahedral mesh...fmd in gmshWrap.py: 16.343534
writing .geo to file completed, save location:
/media/jkl/data/phd/tmp/resipy/src/resipy/invdir
Reading mesh3d.msh
Gmsh version == 3.x
reading node coordinates...
Determining element type...Tetrahedra
Reading connection matrix...
ignoring 7581 elements in the mesh file, as they are not required for R2/R3t
Finished reading .msh file
interpolating topography onto mesh using triangulate interpolation...done
Done
ResIPy Estimated RAM usage = 0.164347 Gb
done
Writing .in file and protocol.dat... done!
------------ INVERTING REFERENCE SURVEY ---------------
>> R 3 t E R T M o d e l v 2.30 <<
>> Date: 19-01-2021
>> My beautiful 3D survey
>> I n v e r s e S o l u t i o n S e l e c t e d <<
>> A d v a n c e d M e s h I n p u t <<
>> T e t r a h e d r a l E l e m e n t M e s h <<
>> Reading mesh file
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading resistivity model from res0.dat
>> L o g - D a t a I n v e r s i o n <<
>> N o r m a l R e g u l a r i s a t i o n <<
>> Memory estimates:
For 1000 measurements the memory needed is: 0.250 Gb
For 2000 measurements the memory needed is: 0.491 Gb
For 5000 measurements the memory needed is: 1.213 Gb
For 10000 measurements the memory needed is: 2.417 Gb
>> Forming roughness matrix
>> Number of measurements read: 644
>> Total Memory required is: 0.164 Gb
Processing frame 1 - output to file f001.dat
Iteration 1
Initial RMS Misfit: 2.41 Number of data ignored: 0
Alpha: 67.701 RMS Misfit: 0.54 Roughness: 3.208
Step length set to 1.000
Final RMS Misfit: 0.54
Final RMS Misfit: 1.01
Solution converged - Outputing results to file
Calculating sensitivity map
End of data: Terminating
>> Program ended normally
----------------- Computing d-d0+f(m0) ---------------
Writing .in file and mesh.dat... done!
Writing protocol.dat... done!
Running forward model...
>> R 3 t E R T M o d e l v 2.30 <<
>> Date: 19-01-2021
>> My beautiful 3D survey
>> F o r w a r d S o l u t i o n S e l e c t e d <<
>> A d v a n c e d M e s h I n p u t <<
>> T e t r a h e d r a l E l e m e n t M e s h <<
>> Reading mesh file
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading resistivity model from resistivity.dat
>> Memory estimates:
For 1000 measurements the memory needed is: 0.006 Gb
For 2000 measurements the memory needed is: 0.006 Gb
For 5000 measurements the memory needed is: 0.006 Gb
For 10000 measurements the memory needed is: 0.006 Gb
>> Number of measurements read: 644
>> Total Memory required is: 0.006 Gb
>> Program ended normally
0/644 reciprocal measurements found.
0/644 reciprocal measurements found.
Forward modelling done.
--------------------- MAIN INVERSION ------------------
>> R 3 t E R T M o d e l v 2.30 <<
>> Date: 19-01-2021
>> My beautiful 3D survey
>> I n v e r s e S o l u t i o n S e l e c t e d <<
>> A d v a n c e d M e s h I n p u t <<
>> T e t r a h e d r a l E l e m e n t M e s h <<
>> Reading mesh file
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading resistivity model from Start_res.dat
>> L o g - D a t a I n v e r s i o n <<
>> B a c k g r o u n d R e g u l a r i s a t i o n <<
>> Memory estimates:
For 1000 measurements the memory needed is: 0.250 Gb
For 2000 measurements the memory needed is: 0.491 Gb
For 5000 measurements the memory needed is: 1.213 Gb
For 10000 measurements the memory needed is: 2.417 Gb
>> Forming roughness matrix
>> Number of measurements read: 644
>> Total Memory required is: 0.164 Gb
Processing frame 1 - output to file f001.dat
Iteration 1
Initial RMS Misfit: 1.40 Number of data ignored: 0
Alpha: 69.072 RMS Misfit: 0.86 Roughness: 3.886
Step length set to 1.000
Final RMS Misfit: 0.86
Final RMS Misfit: 1.00
Solution converged - Outputing results to file
Calculating sensitivity map
End of data: Terminating
>> Program ended normally
2/2 results parsed (2 ok; 0 failed)
k.showResults(index=1, pvslices=[[10,20,30,40],[0],[]], attr='difference(percent)',
color_map='bwr', vmin=-10, vmax=10)
should_run_async will not call transform_cell automatically in the future. Please pass the result to transformed_cell argument and any exception that happen during thetransform in preprocessing_exc_tuple in IPython 7.17 and above. Converting np.character to a dtype is deprecated. The current result is np.dtype(np.str_) which is not strictly correct. Note that np.character is generally deprecated and 'S1' should be used. Converting np.character to a dtype is deprecated. The current result is np.dtype(np.str_) which is not strictly correct. Note that np.character is generally deprecated and 'S1' should be used.
Note that the two outputs are slightly different. With the difference
inversion (reg_mode == 2) being smoother thant the background
constrained (reg_mode == 1). Indeed, in this last one, we haven’t
inverted the differences in resistivity and the difference(percent)
argument has been computed afterwards by ResIPy.
Download python script: nb_3d-time-lapse.py
Download Jupyter notebook: nb_3d-time-lapse.ipynb
View the notebook in the Jupyter nbviewer
Run this example interactively: 