DCIP3D manual - running the codes

DCIP3D manual:
3 Executing the DCIP3D Programs

Contents of this page:
5 program files:	DCIPF3D:	performs forward modelling for DC and IP data over a 3D structure.
	DCINV3D:	inverts DC potentials to recover a 3-D conductivity model using a Gauss-Newton method.
	IPSEN3D:	calculates the sensitivity matrix for the 3-D IP inversion.
	IPINV3D:	inverts IP data to recover a 3-D chargeability model.
	MAKE_WDAT:	makes a weighting file for smoothing near surface zones of the model.

3.1 Introduction

All programs in the package can be run by typing the program name followed by a command line argument. With such a format, they can be executed directly on the command line or in a shell script or batch file. When a program is executed without any arguments, it will print a simple message describing the usage. The command format and the control file format are described below.

Command format:
PROGRAM arg_1

PROGRAM = name of the executable program. If the program is not in the current directory, its path must be included also.
arg_1 = a command line argument which is the name of a file. It is usually one of those described in the preceding section or a control file containing input parameters, as follows:

Control file format:

PARAMETER_1   !COMMENT
 
            :
 
            PARAMETER_n   !COMMENT
 
            FILE_NAME_1   !COMMENT
 
            :
 
            FILE_NAME_n   !COMMENT

PARAMETER_n = a constant number which is passed to the program. If more than one parameter is required on the same line, they must be separated by a space.
FILE_NAME_n = name of a file to be passed to the program. There is always only one file name per line.
!COMMENT = optional comments which describe the parameter(s) on the same line. All comments must be in the last fields of each line and be preceded by "!".

3.2 DCIPF3D

This program performs 3D forward modelling of DC resistivity and IP data. As a command line argument, it requires an input file containing all parameters and files needed to carry out the modelling calculations. For users familiar with DCIP3D Version 1.0 (obtained prior to June 2004), this is a major change in usage. The name of all input files are specified by users but the names of output files are fixed. When topo.dat or topo.idx is supplied, it is used by DCIPF3D to define the topographic surface of the model and the forward modelled DC potential data and IP data simulate the data acquired in the presence of the topography; otherwise, the program assumes a flat earth's surface at z = 0 and carries out the modelling accordingly.

Command line usage:
dcipf3d dcipf3d.inp

Format of the control file dcipf3d.inp:

mesh3d.dat
 
            obs3d.loc
 
            model.con
 
            model.chg
 
            topo.dat | NULL | null
 
            0 | 1
 
            tol
 
            nvectors

Control parameters and input files:

NOTE: file formats are detailed in the "Elements" chapter section 2.2.

`dc \| ip \| ipL`	\|	Enter dc to perform only DC forward modelling, or enter ip to perform both DC and IP modelling. Entering ipL calculates the ip data by multiplying the sensitivity matrix by chargeability. When dc is chosen, the chargeability model line is ignored
`mesh3d.dat`	\|	file containing the finite difference mesh.
`obs3d.loc`	\|	file containing the electrode locations.
`model.con`	\|	file containing the cell values of a conductivity model. If the conductivity model is constant, the value in units of S/m can be entered instead.
`model.chg`	\|	file containing the cell values of a chargeability model. Required only if ip or ipL is selected. If the chargeability model is constant, the value can be entered instead.
`topo.dat`	\|	surface topography file (optional). When the electrode locations file obs3d.dat is stored in the surface data format, the topography should be in topo.dat format. When the locations file is in the general data format, the topography should be in topo.idx format. If there is no topography information this line should contain NULL or null.
`0 \| 1`	\|	if 1 is entered on line 7, the potentials for every cell are output in file potentials_???.txt where ??? is the transmitter number. This file can then be viewed using the UBC-GIF viewing utility program Meshtools3d.
`tol`	\|	a real number (such as 1.0e-5). This value indicates how well the forward modelled system is solved. If this line is missing, a default of 1.0e-5 will be used.
`nvectors`	\|	-1 if there is to be no limit on the amount of memory to be used for calculating the sensitivity matrix G. Enter an integer (such as 10) to specify how solution vectors are to be stored in the computer's memory. See "Sensitivity calculations" above under section 2.1. Used only when ipL is chosen.

Output files generated by DCIPF3D are listed next. These filenames are generated by the program and can not be modified. Therefore any new modelling run carried out in a directory that contains previous results will result in overwriting the earlier files.

dc3d.dat = the DC potential data.
ip3d.dat = the IP data.
ip3d_lin.dat = IP data when ipL is chosen.
obs.loc = location of the electrodes specified by three coordinates. Produced only if surface-data format is used in the input.

3.3 DCINV3D

DCINV3D performs the inversion of the DC resistivity data in file obs3d.dat. The program requires a control file as the argument on the command line. The control file contains the control parameters and the names of input files. As input, the program requires the mesh file, potential data file, initial and reference models, the optional topography and special weighting files. In addition, it requires parameters specifying the choice of regularization parameter, and the wavelet and threshold used to compress the sensitivity matrix. It outputs the inverted model to a file named dcinv3d.con and generates a predicted data file dcinv3d.pre and other auxiliary information files.

After reading in the dc data, it is checked to see if some data might have a wrong sign. A file 'check_sign.txt' is created. This is just a warning about the potential of an incorrect sign.

Command line usage:
dcinv3d dcinv3d.inp

Format of the control file dcinv3d.inp:

MODE, PAR
 
            OBS3D.DAT
 
            MESH3D.DAT
 
            TOPO.DAT | TOPO.IDX | NULL
 
            INI_MOD.CON | CONST | NULL
 
            REF_MOD.CON | CONST | NULL
 
            ACTIVE CELL FILE | NULL
 
            Le Ln Lz | a_s, a_x, a_y, a_z | NULL
 
            WVLT | NULL | NONE
 
            ITOL PAR | NULL
 
            W.DAT
 
            IDISK
 
            tol
 
            nvectors

Control parameters and input files:

`niter`	\|	maximum number of iterations to be performed.
`irest`	\|	restarting control parameter: 0: begin inversion from scratch 1: restart inversion from the previous iteration. This requires that the files dcinv3d.con, dcinv3d.out, and dcinv3d.log be present.
`mode, par`	\|	an integer and a real number that specify one of the three choices to determine the regularization parameter (see the flowchart in Figure 5 in Section 1.4). The parameter par has different meanings in different modes. mode=1: the program chooses a regularization parameter : to produce a user specified target misfit. par is the misfit parameter chifact such that the final Chi² = N x chifact, where N is the number of data. Usually chifact = 1.0 . mode=2: the user inputs the regularization parameter. In this mode, : = par. mode=3: the program calculates the regularization parameter according to the L-curve criterion. par is ignored.
`obs3d.dat`	\|	name of the file containing observed potential data, electrode configuration and errors.
`mesh3d.dat`	\|	name of the file containing finite difference mesh.
`topo`	\|	name of the file containing surface topography. When NULL or null is entered, the surface is assumed flat. It is essential that the information about the topography is specified in the format that is consistent with the information of electrode positions that are supplied in the obs3d.dat file. Therefore: When the observations file obs3d.dat is stored in the surface data format, the topography should be in topo.dat format. When the observations file is in the general data format, the topography should be in topo.idx format.
`ini_mod.con`	\|	name of the file containing initial conductivity model. If the initial model is constant, the filename can be replaced by the constant value. If NULL or null is entered, the initial model will be equal to the reference model. Stored in model.con format.
`ref_mod.con`	\|	name of the file containing reference conductivity model. If the reference model is constant, the filename can be replaced by the constant value. If NULL or null is entered, the reference model will be the best fitting halfspace. Stored in model.con format.
`active cells`	\|	A file indicating which cells are active (1) and which are inactive (0) . The values in this file are ordered the same as the model file. The inactive cells are set to the reference model. If NULL or null is entered, all cells are active.
`Le, Ln, Lz \| a_s, a_x, a_y, a_z`	\|	either length scales or alpha values can be specified here, or NULL. When three values are entered, it is assumed they are length scales: (`Le Ln Lz`). When four values are entered, they are assumed to be alpha values: (`a_s, a_x, a_y, a_z`). The .log file would have both sets of values printed out. Length scales are in easting, northing, and depth directions respectively. These parameters determine the weighting coefficients (`a_s, a_x, a_y, a_z`) in the model objective function. The relationship between length scales and the weighting coefficients is given in Section 1.3. The recommended value of the length scale is two to five cell widths in the corresponding direction. If `NULL` or `null` is entered, the length scale will be equal to two times the maximum cell width in the middle of the mesh. For example, if the cells are 50m x 50m x 25m at the centre of the mesh, then the default values are `Le` `= Ln = Lz = 100m`.
`wvlt`	\|	a five-character string identifying the type of wavelet used to compress the sensitivity matrix. The types of wavelets available are Daubechies wavelet with 1 to 6 vanishing moments (daub1, daub2, and so on) and Symmlets with 4 to 6 vanishing moments (symm4, symm5, symm6). Note that daub1 is the Haar wavelet and daub2 is the Daubechies-4 wavelet. The Daubechies-4 wavelet should be used for most inversions, while the others are provided for users' experimentation. If NULL or null is entered, the default Daubechies-4 wavelet (daub2) is used. If NONE or none is entered, no compression is performed and the program generates a dense matrix in its original form.
`itol, par`	\|	an integer and a real number that specify how the wavelet threshold level is to be determined. itol=1: par is the relative reconstruction error r* of the sensitivity. A reconstruction error of 0.05 is usually adequate, but 0.02 is recommended. The relative threshold value e for compressing each sensitivity is calculated by the program. itol=2: par is the relative threshold e to be used. A value of par=0.001 is recommended. If null is entered on this line, a default relative reconstruction error of 0.02 is used and the relative threshold level is calculated (i.e. itol=1, par=0.02). The detailed explanation of threshold level and reconstruction error can be found in the background section (section 1.6) of this manual.
`w.dat`	\|	name of the file containing weighting matrix. If NULL is entered, default values of unity are used.
`idisk`	\|	=0: store the entire sensitivity matrix in memory. =1: access the sensitivity matrix from disk when needed.
`tol`	\|	a real number (such as 1.0e-5). This value indicates how well the forward modelled system is solved. If this line is missing, a default of 1.0e-5 will be used.
`nvectors`	\|	-1 if there is to be no limit on the amount of memory to be used for calculating the sensitivity matrix G. Enter an integer (such as 10) to specify how many solution vectors are to be stored in the computer's memory. See "Sensitivity calculations" above under section 2.1.

Output files:
DCINV3D Versions 2.0 and later save a model after each iteration. The models are ordered: dcinv3d_01.con, dcinv3d_02.con, dcinv3d_03.con, etc.

dcinv3d.con = conductivity model of the latest iteration. The model is stored in the model.con format. This file is overwritten at the end of each iteration.

dcinv3d.out = convergence information of past iterations. This file lists the values of data misfit, the model objective function, and the regularization parameter as functions of iteration. In Versions 2.0 and later of the program library, a header line has been added to this output file.

dcinv3d.log = log file containing the detailed information about each iteration.

dcinv3d.pre = predicted potential data from the inverted model in the latest iteration. The predicted data is in the format of obs3d.dat with the column corresponding to data error omitted. True data are included as a separate column (in program library Versions 2.0 and later). This file is overwritten at the end of each iteration.

obs.loc = electrode locations with three coordinates specified. This is output only when the obs3d.dat is stored in the surface data format.

sensitivity.txt = the average absolute value of the sensitivity matrix for each cell. This file is stored in the same format as the model files and can be displayed using Meshtools3d (use a log scale).

Example DCINV3D Control File: The following is an example control file. The inversion is started from scratch and a maximum 25 iterations are to be performed. The inversion will use mode=1 to determined the regularization parameter and the chifact is equal to 1.0. The observed potential data are in data.dc. The model discretization is defined by the finite difference mesh in mesh3d.dat and the surface topography in topo.dat. The initial and reference models are equal to 0.005 S/m. The sensitivity matrix is compressed using Daubechies-4 wavelet and a relative threshold of 0.001 is supplied by the user. The default weighting function (unity) is used for the model objective function.

25 0       ! maxit, irest 
 
            1 1.       ! mode, par
 
            data.dc    ! data file
 
            mesh3d.dat ! mesh
 
            topo.dat   ! topography | null
 
            null       ! initial model | null
 
            0.005      ! reference model | null
 
            null       ! active cells
 
            100.0, 100.0, 100.0 ! Le, Ln, Lz | null
 
            daub2      ! wavelet | null |none
 
            2 0.001    ! itol, par | null
 
            null       ! 3D weighting | null
 
            0          ! idisk 0|1
 
            1.0e-5     ! forward modelling toleranc
 
            -1         ! no limit on number of solution vectors stored in memory

NOTE-1: A sample input file can be obtained by executing: dcinv3d -inp.

NOTE-2: In mode=1, DCINV3D will terminate before the specified maximum number of iterations is reached if the expected data misfit is achieved and if the model norm has plateaued. However, if the program exits when the maximum iteration is reached, the file dcinv3d.out should be checked to see if the desired misfit (equal to chifact times the number of data) has been reached and if the model norm is no longer changing. If either of these conditions has not been met then the program should be restarted.

NOTE-3: In mode=1, if the desired misfit level cannot be achieved, DCINV3D will determine a reasonable regularization parameter and carry out the inversion to completion. This gives a valid model from the inversion unless the program is terminated at the maximum number of iterations.

NOTE-4: The restart option is only used to continue the inversion after the program has been terminated abnormally, such as after it is interrupted or when the maximum number of iterations has been reached. If the user wants to re-invert the data by changing the mode or inversion parameters, a different inversion must be performed independently.

3.4 IPSEN3D

This program calculates the sensitivity file for the IP inversion.

Command line usage:
ipsen3d ipsen3d.inp

Format of the control file ipsens3d.inp:

OBS3D_IP.DAT 
 
            MESH3D.DAT
 
            DCINV3D.CON | CONST
 
            TOPO.DAT | TOPO.IDX | NULL
 
            ACTIVE CELL FILE | NULL
 
            WVLT | NULL | NONE
 
            ITOL PAR | NULL
 
            tol
 
            nvectors

Control parameters and input files:

`obs3d_ip.dat`	\|	name of the file containing observed chargeability data, electrode configuration and errors.
`mesh3d.dat`	\|	name of the file containing finite difference mesh.
`dcinv3d.con`	\|	name of the file containing the conductivity model used for sensitivity calculation. If the model is constant, the filename can be replaced by the constant value.
`topo`	\|	name of the file containing surface topography. When NULL or null is entered, the surface is assumed flat. It is essential that the information about the topography is specified in the format that is consistent with the information of electrode positions that are supplied in the obs3d.dat file. Therefore: When the observations file obs3d.dat is stored in the surface data format, the topography should be in topo.dat format. When the observations file is in the general data format, the topography should be in topo.idx format.
`active cells`	\|	A file indicating which cells are active (1) and which are inactive (0) . The values in this file are ordered the same as the model file. The inactive cells will be set to 0 in the inversion. If NULL or null is entered, all cells are active.
`wvlt`	\|	a five-character string identifying the type of wavelet used to compress the sensitivity matrix. The types of wavelets available are Daubechies wavelet with 1 to 6 vanishing moments (daub1, daub2, and so on) and Symmlets with 4 to 6 vanishing moments (symm4, symm5, symm6). Note that daub1 is the Haar wavelet and daub2 is the Daubechies-4 wavelet. The Daubechies-4 (daub2) wavelet should be used for most inversions, while the others are provided for users' experimentation. If NULL or null is entered, the default Daubechies-4 wavelet (daub2) is used. If NONE or none is entered, no compression is performed and the program generates a dense matrix in its original form.
`itol,`	\|	par an integer and a real number that specify how the wavelet threshold level is to be determined. itol=1: par is the relative reconstruction error r* of the sensitivity. A reconstruction error of 0.05 is usually adequate, but 0.02 is recommended. The relative threshold value e for compressing each sensitivity is calculated by the program. itol=2: par is the relative threshold e to be used. A value of par=0.001 is recommended. If NULL or null is entered on this line, a default relative reconstruction error of 0.02 is used and the relative threshold level is calculated (i.e.,itol=1, par=0.02). The detailed explanation of threshold level and reconstruction error can be found in the Section 1.6 of this manual.
`tol`	\|	a real number (such as 1.0e-5). This value indicates how well the forward modelled system is solved. If this line is missing, a default of 1.0e-5 will be used.
`nvectors`	\|	-1 if there is to be no limit on the amount of memory to be used for calculating the sensitivity matrix G. Enter an integer (such as 10) to specify how many solution vectors are to be stored in the computer's memory. See "Sensitivity calculations" above under section 2.1.

Output files:
ipinv3d.mtx = sensitivity matrix file to be used by ipinv3d to perform the IP inversion.

ipsens3d.log = log file containing information about the sensitivity calculation.

obs.loc = electrode locations. This is output only when the obs3d.dat is stored in the surface data format.

obs_ip.dat  ! data file 
 
            mesh        ! mesh file
 
            dcinv3d.con ! conductivity file | constant
 
            topo.dat    ! topography
 
            null        ! active cells
 
            daub2       ! wavelet
 
     ! itol, par | null
 
            1.0e-5      ! forward modelling tolerance
 
            -1          ! no limit on number of solution vectors stored in memory

3.5 IPINV3D

IPINV3D performs the inversion of the IP data in file obs3d_ip.dat. The program requires a control file as the argument on the command line. The control file contains the control parameters and the names of input files. The program requires as input the mesh file, apparent chargeability data file, initial and reference models, and possibly a special weighting file. It outputs the inverted model in a file named ipinv3d.chg. The command format and the control file format are described below.

Command line usage:
ipinv3d ipinv3d.inp

Format of the control file ipinv2d.inp:

IREST 
 
            MODE, PAR
 
            OBS3D_IP.DAT
 
            IPINV3D.MTX
 
            INI_MOD.CHG
 
            REF_MOD.CHG
 
            Le Ln Lz | a_s, a_x, a_y, a_z | NULL
 
            W.DAT
 
            IDISK

Control parameters:

`irest`	\|	restarting flag: =0: start inversion from scratch. =1: restart inversion after it is interrupted. Restart requires two files written out by IPINV3D at the end of each iteration: ipinv3d.aux and ipinv3d.eta.
`mode, par`	\|	an integer and a real number that specify three choices to determine the regularization parameter (see the flowchart in Figure 7 in Section 1.5). The parameter par has different meanings in different modes. mode=1: the program chooses a regularization parameter m to produce a user specified target misfit. par is the misfit parameter chifact such that the final Chi² = N x chifact, where N is the number of data. Usually chifact = 1.0 . mode=2: the user inputs the regularization parameter. In this mode, m = par. mode=3: the program calculates the regularization parameter using the GCV criterion. par is ignored.
`obs3d_ip.dat`	\|	observed IP data and electrode configuration. This should be the same file as that input to ipsen3d. (standard deviations can change.)
`ipinv3d.mtx`	\|	sensitivity matrix computed by IPSEN3D.
`ini_mod.chg`	\|	initial chargeability model. This is either a file name or a real value. If NULL is entered, the default initial model will be 0.05 of the maximum data if IPTYPE=1, or 0.01 if IPTYPE=2. See description of the obs3d.dat file in section 2.2. Since the model in the IP inversion cannot be exactly on the zero bound, the initial model must not be zero. Stored in model.chg format. If the initial chargeability model is constant, the file name can be replaced by the constant value in the parameter file.
`ref_mod.chg`	\|	reference chargeability model. This is either a file name or a real value. If NULL is entered, a default value of 0.0 will be used. Stored in model.chg format. If the reference chargeability model is constant, the file name can be replaced by the constant value in the parameter file.
`Le, Ln, Lz \| a_s, a_x, a_y, a_z`	\|	either length scales or alpha values can be specified here, or NULL. When three values are entered, it is assumed they are length scales: (`Le Ln Lz`). When four values are entered, they are assumed to be alpha values: (`a``_s, a_x, a_y, a_z`). The .log file would have both sets of values printed out. Length scales are in easting, northing, and depth directions respectively. These parameters determine the weighting coefficients (`a``_s, a_x, a_y, a_z`) in the model objective function. The relationship between length scales and the weighting coefficients is given in Section 1.3. The recommended value of the length scale is two to five cell widths in the corresponding direction. If NULL or null is entered, the length scale will be equal to two times the maximum cell width in the middle of the mesh. For example, if the cells are 50m x 50m x 25m at the centre of the mesh, then the default values are Le = Ln = Lz = 100m.
`w.dat`	\|	user supplied weighting functions. If NULL is entered, default values of unity are used.
`idisk`	\|	=0: store the entire sensitivity matrix in memory. =1: access the sensitivity matrix from disk when needed.

Output files:
ipinv3d.chg = chargeability model of latest iteration. This file is overwritten at the end of each iteration.

ipinv3d.log = log file containing detailed information about each iteration.

ipinv3d.pre = predicted apparent chargeability data from the inverted model in the latest iteration. The predicted data are in the format of obs3d_ip.dat with the field corresponding to the data error removed. True data are included as a separate column (in program library Versions 2.0 and later). This file is overwritten at the end of each iteration.

ipinv3d.aux = auxiliary file storing information on the progress of the inversion. This file is used for restarting the inversion after interruption.

ipinv3d.eta = auxiliary file storing the temporary chargeability model at each iteration. This file is used for restarting the inversion after interruption. At the end of the inversion, it is identical to the file ipinv3d.chg, so it can be deleted.

Example of IPINV3D Control File: The following is an example control file for IPINV3D.

0           ! irest 
 
            1, 1.0      ! mode (1, 2, or 3), par
 
            obs_ip.dat  ! data file
 
            ipinv3d.mtx ! mtx file
 
            0.0         ! initial model | null
 
            0.0         ! reference model | null
 
            100 100 100 ! Le, Ln, Lz | null
 
            w.dat       ! 3D weighting | null
 
            0           ! idisk (0 or 1)

NOTE-1: A sample input file can be obtained by executing: ipinv3d -inp.

NOTE-2: The default conductivity of a uniform halfspace for IP inversions should only be used for preliminary examination of the data. When there is little structure in the background conductivity, the inversion using a uniform conductivity can yield a reasonable chargeability model and it is justifiable to fit the data close to the expected misfit value. However, when the background conductivity deviates greatly from a uniform halfspace, reproducing the data to within the assumed errors will certainly result in overfitting the data. If the halfspace conductivity is assumed and inversion is run with mode=1, then it is prudent to assign a value greater than 1.0 for chifact when the background conductivity is structurally complex. The judgement can be made based upon the complexity of the apparent resistivity pseudo-section.

NOTE-3: In general, when any form of approximate conductivity other than the one obtained from full nonlinear inversion of the DC resistivity data is used, it is highly recommended that the IP inversion be run with mode=3 so that an optimal regularization parameter is determined by the GCV analysis.

NOTE-4: The restart option is only used to continue the inversion after the program has been terminated abnormally, for instance, after it is interrupted. If the user wants to re-invert the data by changing the mode or inversion parameters, a different inversion must be performed independently.

3.6 MAKE_WDAT.EXE

This is a utility used to make a 'w.dat' file that has smoothing in the X and Y directions for the first few layers that underlie the topography surface. This suppresses the tendency of the algorithm to make highly variable structure in these top layers of the model.

Command line usage:
make_wdat make_wdat.inp

Format (by example) of make_wdat.inp:

mesh.txt     ! mesh file
 
            topo.dat     ! topography file | null
 
            6            ! # of layers where extra smoothing will be applied
 
            64. 32. 16. 8. 4. 2.   ! weighting values to place in top layers of 'w.dat'

In this example, the top 6 layers below air will have extra smoothing in the X and Y directions. The first layer will have a smoothing value of 64, the second layer down will have weighting value 32, etc. The first "layer" is defined as the first cell below topography, the second "layer" is defined as the next cell down, etc. No weighting is applied to any cells that are in the air above topography.

Weighting values in all other cells of the mesh will remain at 1.0 .