simulator.rockphysics.standardrp

  1"""Descriptive description."""
  2
  3__author__ = 'TM'
  4
  5# standardrp.py
  6import numpy as np
  7import sys
  8import multiprocessing as mp
  9# internal load
 10from misc.system_tools.environ_var import OpenBlasSingleThread  # Single threaded OpenBLAS runs
 11
 12
 13class elasticproperties:
 14    """
 15    Calculate elastic properties from standard
 16    rock-physics models, specifically following Batzle
 17    and Wang, Geophysics, 1992, for fluid properties, and
 18    Report 1 in Abul Fahimuddin's thesis at Universty of
 19    Bergen (2010) for other properties.
 20
 21    Example
 22    -------
 23    >>> porosity = 0.2
 24    ... pressure = 5
 25    ... phases = ["Oil","Water"]
 26    ... saturations = [0.3, 0.5]
 27    ...
 28    ... satrock = Elasticproperties()
 29    ... satrock.calc_props(phases, saturations, pressure, porosity)
 30    """
 31
 32    def __init__(self, input_dict):
 33        self.dens = None
 34        self.bulkmod = None
 35        self.shearmod = None
 36        self.bulkvel = None
 37        self.shearvel = None
 38        self.bulkimp = None
 39        self.shearimp = None
 40        # The overburden for each grid cell must be
 41        # specified as values on an .npz-file whose
 42        # name is given in input_dict.
 43        self.input_dict = input_dict
 44        self._extInfoInputDict()
 45
 46    def _extInfoInputDict(self):
 47        # The key word for the file name in the
 48        # dictionary must read "overburden"
 49        if 'overburden' in self.input_dict:
 50            obfile = self.input_dict['overburden']
 51            npzfile = np.load(obfile)
 52            # The values of overburden must have been
 53            # stored on file using:
 54            # np.savez(<file name>,
 55            # obvalues=<overburden values>)
 56            self.overburden = npzfile['obvalues']
 57            npzfile.close()
 58        if 'baseline' in self.input_dict:
 59            self.baseline = self.input_dict['baseline']  # 4D baseline
 60        if 'parallel' in self.input_dict:
 61            self.parallel = self.input_dict['parallel']
 62
 63    def _filter(self):
 64        bulkmod = self.bulkimp
 65        self.bulkimp = bulkmod.flatten()
 66
 67    def setup_fwd_run(self, state):
 68        """
 69        Setup the input parameters to be used in the PEM simulator. Parameters can be a an ensemble or a single array.
 70        State is set as an attribute of the simulator, and the correct value is determined in self.pem.calc_props()
 71
 72        Parameters
 73        ----------
 74        state : dict
 75            Dictionary of input parameters or states.
 76
 77        Changelog
 78        ---------
 79        - KF 11/12-2018
 80        """
 81        # self.inv_state = {}
 82        # list_pem_param =[el for el in [foo for foo in self.pem['garn'].keys()] + [foo for foo in self.filter.keys()] +
 83        #                 [foo for foo in self.__dict__.keys()]]
 84
 85        # list_tot_param = state.keys()
 86        # for param in list_tot_param:
 87        #    if param in list_pem_param or (param.split('_')[-1] in ['garn', 'rest']):
 88        #        self.inv_state[param] = state[param]
 89
 90        pass
 91
 92    def calc_props(self, phases, saturations, pressure,
 93                   porosity, wait_for_proc=None, ntg=None, Rs=None, press_init=None, ensembleMember=None):
 94        ###
 95        # This doesn't initialize for models with no uncertainty
 96        ###
 97        # # if some PEM properties have uncertainty, set the correct value
 98        # if ensembleMember is not None:
 99        #     for pem_state in self.__dict__.keys(): # loop over all possible pem vaules
100        #         if pem_state is not 'inv_state': # do not alter the ensemble
101        #             if type(eval('self.{}'.format(pem_state))) is dict:
102        #                 for el in eval('self.{}'.format(pem_state)).keys():
103        #                     if type(eval('self.{}'.format(pem_state))[el]) is dict:
104        #                         for param in eval('self.{}'.format(pem_state))[el].keys():
105        #                             if param in self.inv_state:
106        #                                 eval('self.{}'.format(pem_state))[el][param]=\
107        #                                     self.inv_state[param][:, ensembleMember]
108        #                             elif param + '_' + el in self.inv_state:
109        #                                 eval('self.{}'.format(pem_state))[el][param] = \
110        #                                     self.inv_state[param+'_' + el][:, ensembleMember]
111        #                     else:
112        #                         if el in self.inv_state:
113        #                             eval('self.{}'.format(pem_state))[el] = self.inv_state[el][:,ensembleMember]
114        #             else:
115        #                 if pem_state in self.inv_state:
116        #                     setattr(self,pem_state, self.inv_state[pem_state][:,ensembleMember])
117
118        # Check if the inputs are given as a list (more
119        # than one phase) or a single input (single
120        # phase). If single phase input, make the input a
121        # list with a single entry (s.t. it can be used
122        # directly with the methods below)
123        #
124        if not isinstance(phases, list):
125            phases = [phases]
126        if not isinstance(saturations, list):
127            saturations = [saturations]
128        if not isinstance(pressure, list) and \
129                type(pressure).__module__ != 'numpy':
130            pressure = [pressure]
131        if not isinstance(porosity, list) and \
132                type(porosity).__module__ != 'numpy':
133            porosity = [porosity]
134        #
135        # Load "overburden" into local variable to
136        # comply with remaining code parts
137        overburden = self.overburden
138        if press_init is None:
139            p_init = self.p_init
140        else:
141            p_init = press_init
142        #
143        # Check that no. of phases is equal to no. of
144        # entries in saturations list
145        #
146        assert (len(saturations) == len(phases))
147        #
148        # Make saturation a Numpy array (so that we
149        # can easily access the values for each
150        # phase at one grid cell)
151        #
152        # Transpose makes it a no. grid cells x phases
153        # array
154        saturations = np.array(saturations).T
155        #
156        # Check if we actually inputted saturation values
157        # for a single grid cell. If yes, we redefine
158        # saturations to get it on the correct form (no.
159        # grid cells x phases array).
160        #
161        if saturations.ndim == 1:
162            saturations = \
163                np.array([[x] for x in saturations]).T
164        #
165        # Loop over all grid cells and calculate the
166        # various saturated properties
167        #
168        self.phases = phases
169
170        self.dens = np.zeros(len(saturations[:, 0]))
171        self.bulkmod = np.zeros(len(saturations[:, 0]))
172        self.shearmod = np.zeros(len(saturations[:, 0]))
173        self.bulkvel = np.zeros(len(saturations[:, 0]))
174        self.shearvel = np.zeros(len(saturations[:, 0]))
175        self.bulkimp = np.zeros(len(saturations[:, 0]))
176        self.shearimp = np.zeros(len(saturations[:, 0]))
177
178        if ntg is None:
179            ntg = [None for _ in range(len(saturations[:, 0]))]
180        if Rs is None:
181            Rs = [None for _ in range(len(saturations[:, 0]))]
182        if p_init is None:
183            p_init = [None for _ in range(len(saturations[:, 0]))]
184
185        for i in range(len(saturations[:, 0])):
186            #
187            # Calculate fluid properties
188            #
189            # set Rs if needed
190            densf, bulkf = \
191                self._fluidprops(self.phases,
192                                 saturations[i, :], pressure[i], Rs[i])
193            #
194            denss, bulks, shears = \
195                self._solidprops(porosity[i], ntg[i], i)
196            #
197            # Calculate dry rock moduli
198            #
199
200            bulkd, sheard = \
201                self._dryrockmoduli(porosity[i],
202                                    overburden[i],
203                                    pressure[i], bulks,
204                                    shears, i, ntg[i], p_init[i], denss, Rs[i], self.phases)
205            # -------------------------------
206            # Calculate saturated properties
207            # -------------------------------
208            #
209            # Density
210            #
211            self.dens[i] = (porosity[i]*densf +
212                            (1-porosity[i])*denss)*0.001
213            #
214            # Moduli
215            #
216            self.bulkmod[i] = \
217                bulkd + (1 - bulkd/bulks)**2 / \
218                (porosity[i]/bulkf +
219                 (1-porosity[i])/bulks -
220                 bulkd/(bulks**2))
221            self.shearmod[i] = sheard
222
223            # Velocities (due to bulk/shear modulus being
224            # in MPa, we multiply by 1000 to get m/s
225            # instead of km/s)
226            #
227            self.bulkvel[i] = \
228                100*np.sqrt((abs(self.bulkmod[i]) +
229                             4*self.shearmod[i]/3)/(self.dens[i]))
230            self.shearvel[i] = \
231                100*np.sqrt(self.shearmod[i] /
232                            (self.dens[i]))
233            #
234            # Impedances (m/s)*(Kg/m3)
235            #
236            self.bulkimp[i] = self.dens[i] * \
237                self.bulkvel[i]
238            self.shearimp[i] = self.dens[i] * \
239                self.shearvel[i]
240
241    def getMatchProp(self, petElProp):
242        if petElProp.lower() == 'density':
243            self.match_prop = self.getDens()
244        elif petElProp.lower() == 'bulk_modulus':
245            self.match_prop = self.getBulkMod()
246        elif petElProp.lower() == 'shear_modulus':
247            self.match_prop = self.getShearMod()
248        elif petElProp.lower() == 'bulk_velocity':
249            self.match_prop = self.getBulkVel()
250        elif petElProp.lower() == 'shear_velocity':
251            self.match_prop = self.getShearVel()
252        elif petElProp.lower() == "bulk_impedance":
253            self.match_prop = self.getBulkImp()
254        elif petElProp.lower() == 'shear_impedance':
255            self.match_prop = self.getShearImp()
256        else:
257            print("\nError in getMatchProp method")
258            print("No model output type selected for "
259                  "data match.")
260            print("Legal model output types are "
261                  "(case insensitive):")
262            print("Density, bulk modulus, shear "
263                  "modulus, bulk velocity,")
264            print("shear velocity, bulk impedance, "
265                  "shear impedance")
266            sys.exit(1)
267        return self.match_prop
268
269    def getDens(self):
270        return self.dens
271
272    def getBulkMod(self):
273        return self.bulkmod
274
275    def getShearMod(self):
276        return self.shearmod
277
278    def getBulkVel(self):
279        return self.bulkvel
280
281    def getShearVel(self):
282        return self.shearvel
283
284    def getBulkImp(self):
285        return self.bulkimp
286
287    def getShearImp(self):
288        return self.shearimp
289
290    #
291    # ===================================================
292    # Fluid properties start
293    # ===================================================
294    #
295    def _fluidprops(self, fphases, fsats, fpress, Rs=None):
296        #
297        # Calculate fluid density and bulk modulus
298        #
299        #
300        # Input
301        #       fphases - fluid phases present; Oil
302        #                 and/or Water and/or Gas
303        #       fsats   - fluid saturation values for
304        #                 fluid phases in "fphases"
305        #       fpress  - fluid pressure value (MPa)
306        #       Rs      - Gas oil ratio. Default value None
307
308        #
309        # Output
310        #       fdens - density of fluid mixture for
311        #               pressure value "fpress" inherited
312        #               from phaseprops)
313        #       fbulk - bulk modulus of fluid mixture for
314        #               pressure value "fpress" (unit
315        #               inherited from phaseprops)
316        #
317        # -----------------------------------------------
318        #
319        fdens = 0.0
320        fbinv = 0.0
321
322        for i in range(len(fphases)):
323            #
324            # Calculate mixture properties by summing
325            # over individual phase properties
326            #
327            pdens, pbulk = self._phaseprops(fphases[i],
328                                            fpress, Rs)
329            fdens = fdens + fsats[i]*abs(pdens)
330            fbinv = fbinv + fsats[i]/abs(pbulk)
331        fbulk = 1.0/fbinv
332        #
333        return fdens, fbulk
334    #
335    # ---------------------------------------------------
336    #
337
338    def _phaseprops(self, fphase, press, Rs=None):
339        #
340        # Calculate properties for a single fluid phase
341        #
342        #
343        # Input
344        #       fphase - fluid phase; Oil, Water or Gas
345        #       press  - fluid pressure value (MPa)
346        #
347        # Output
348        #       pdens - phase density of fluid phase
349        #               "fphase" for pressure value
350        #               "press" (kg/m³)
351        #       pbulk - bulk modulus of fluid phase
352        #               "fphase" for pressure value
353        #               "press" (MPa)
354        #
355        # -----------------------------------------------
356        #
357        if fphase.lower() == "oil":
358            coeffsrho = np.array([0.8, 829.9])
359            coeffsbulk = np.array([10.42, 995.79])
360        elif fphase.lower() == "wat":
361            coeffsrho = np.array([0.3, 1067.3])
362            coeffsbulk = np.array([9.0, 2807.6])
363        elif fphase.lower() == "gas":
364            coeffsrho = np.array([4.7, 13.4])
365            coeffsbulk = np.array([2.75, 0.0])
366        else:
367            print("\nError in phaseprops method")
368            print("Illegal fluid phase name.")
369            print("Legal fluid phase names are (case "
370                  "insensitive): Oil, Wat, and Gas.")
371            sys.exit(1)
372        #
373        # Assume simple linear pressure dependencies.
374        # Coefficients are inferred from
375        # plots in Batzle and Wang, Geophysics, 1992,
376        # (where I set the temperature to be 100 degrees
377        # Celsius, Note also that they give densities in
378        # g/cc). The resulting straight lines do not fit
379        # the data extremely well, but they should
380        # be sufficiently accurate for the purpose of
381        # this project.
382        #
383        pdens = coeffsrho[0]*press + coeffsrho[1]
384        pbulk = coeffsbulk[0]*press + coeffsbulk[1]
385        #
386        return pdens, pbulk
387
388    #
389    # =======================
390    # Fluid properties end
391    # =======================
392    #
393
394    #
395    # =========================
396    # Solid properties start
397    # =========================
398    #
399    def _solidprops(self, poro, ntg=None, ind=None):
400        #
401        # Calculate bulk and shear solid rock (mineral)
402        # moduli by averaging Hashin-Shtrikman bounds
403        #
404        #
405        # Input
406        #       poro -porosity
407        #
408        # Output
409        #       denss  - solid rock density (kg/m³)
410        #       bulks  - solid rock bulk modulus (unit
411        #                inherited from hashinshtr)
412        #       shears - solid rock shear modulus (unit
413        #                inherited from hashinshtr)
414        #
415        # -----------------------------------------------
416        #
417        # From PetroWiki (kg/m³)
418        #
419        densc = 2540
420        densq = 2650
421        #
422        # From "Step 1" of "recipe" in Report 1 in Abul
423        # Fahimuddin's thesis.
424        #
425        vclay = 0.7 - 1.58*poro
426        #
427        # Calculate solid rock (mineral) density. (Note
428        # that this is often termed \rho_dry, and not
429        # \rho_s)
430        #
431        denss = densq + vclay*(densc - densq)
432        #
433        # Calculate lower and upper bulk and shear
434        # Hashin-Shtrikman bounds
435        #
436        bulkl, bulku, shearl, shearu = \
437            self._hashinshtr(vclay)
438        #
439        # Calculate bulk and shear solid rock (mineral)
440        # moduli as arithmetic means of the respective
441        # bounds
442        #
443        bulkb = np.array([bulkl, bulku])
444        shearb = np.array([shearl, shearu])
445        bulks = np.mean(bulkb)
446        shears = np.mean(shearb)
447        #
448        return denss, bulks, shears
449    #
450    # ---------------------------------------------------
451    #
452
453    def _hashinshtr(self, vclay):
454        #
455        # Calculate lower and upper, bulk and shear,
456        # Hashin-Shtrikman bounds, utilizing that they
457        # all have the common mathematical form,
458        #
459        # f = a + b/((1/c) + d*(1/e)).
460        #
461        #
462        # Input
463        #       vclay - "volume of clay"
464        #
465        # Output
466        #       bulkl  - lower bulk Hashin-Shtrikman
467        #                bound (MPa)
468        #       bulku  - upper bulk Hashin-Shtrikman
469        #                bound (MPa)
470        #       shearl - lower shear Hashin-Shtrikman
471        #                bound (MPa)
472        #       shearu - upper shear Hashin-Shtrikman
473        #                bound (MPa)
474        #
475        # -----------------------------------------------
476        #
477        # From table 1 in Report 1 in Abul Fahimuddin's
478        # thesis (he used GPa, I use MPa), ("c" for clay
479        # and "q" for quartz.):
480        #
481        bulkc = 14900
482        bulkq = 37000
483        shearc = 1950
484        shearq = 44000
485        #
486        # Calculate quantities common for both bulk and
487        # shear formulas
488        #
489        lb = 1 - vclay
490        ub = vclay
491        le = bulkc + 4*shearc/3
492        ue = bulkq + 4*shearq/3
493        #
494        # Calculate quantities common only for bulk
495        # formulas
496        #
497        bld = vclay
498        bud = 1 - vclay
499        blc = bulkq - bulkc
500        buc = bulkc - bulkq
501        #
502        # Calculate quantities common only for shear
503        # formulas
504        #
505        sld = 2*vclay*(bulkc + 2*shearc)/(5*shearc)
506        sud = 2*(1 - vclay)*(bulkq + 2*shearq)/(5*shearq)
507        slc = shearq - shearc
508        suc = shearc - shearq
509        #
510        # Calculate bounds utilizing generic formula;
511        # f = a + b/((1/c) + d*(1/e)).
512        #
513        # Lower bulk
514        #
515        bulkl = self._genhashinshtr(bulkc, lb, blc, bld,
516                                    le)
517        #
518        # Upper bulk
519        #
520        bulku = self._genhashinshtr(bulkq, ub, buc, bud,
521                                    ue)
522        #
523        # Lower shear
524        #
525        shearl = self._genhashinshtr(shearc, lb, slc,
526                                     sld, le)
527        #
528        # Upper shear
529        #
530        shearu = self._genhashinshtr(shearq, ub, suc,
531                                     sud, ue)
532        #
533        return bulkl, bulku, shearl, shearu
534
535    #
536    # ---------------------------------------------------
537    #
538    def _genhashinshtr(self, a, b, c, d, e):
539        #
540        # Calculate arbitrary Hashin-Shtrikman bound,
541        # which has the generic form
542        #
543        # f = a + b/((1/c) + d*(1/e))
544        #
545        # both for lower and upper bulk and shear bounds
546        #
547        #
548        # Input
549        #       a - see above formula
550        #       b - see above formula
551        #       c - see above formula
552        #       d - see above formula
553        #       e - see above formula
554        #
555        # Output
556        #       f - Bulk or shear Hashin-Shtrikman bound
557        #           value
558        #
559        # -----------------------------------------------
560        #
561        cinv = 1/c
562        einv = 1/e
563        f = a + b/(cinv + d*einv)
564        #
565        return f
566    #
567    # =======================
568    # Solid properties end
569    # =======================
570    #
571
572    #
573    # ============================
574    # Dry rock properties start
575    # ============================
576    #
577    def _dryrockmoduli(self, poro, poverb, pfluid, bulks, shears, ind=None, ntg=None, p_init=None, denss=None, Rs=None, phases=None):
578        #
579        #
580        # Calculate bulk and shear dry rock moduli,
581        # utilizing that they have the common
582        # mathematical form,
583        #
584        #     --                                -- ^(-1)
585        #     |(poro/poroc)      1 - (poro/poroc)|
586        # f = |------------   +  ----------------|   - z.
587        #     |   a + z              b + z       |
588        #     --.                               --
589        #
590        #
591        # Input
592        #       poro   - porosity
593        #       poverb - overburden pressure (MPa)
594        #       pfluid - fluid pressure (MPa)
595        #       bulks  - bulk solid (mineral) rock bulk
596        #                modulus (MPa)
597        #       shears - solid rock (mineral) shear
598        #                modulus (MPa)
599        #
600        # Output
601        #       bulkd  - dry rock bulk modulus (unit
602        #                inherited from hertzmindlin and
603        #                variable; bulks)
604        #       sheard - dry rock shear modulus (unit
605        #                inherited from hertzmindlin and
606        #                variable; shears)
607        #
608        # -----------------------------------------------
609        #
610        # Calculate Hertz-Mindlin moduli
611        #
612        bulkhm, shearhm = self._hertzmindlin(poverb,
613                                             pfluid)
614        #
615        # From table 1 in Report 1 in Abul Fahimuddin's
616        # thesis (I assume \phi_max in that
617        # table corresponds to \phi_c in his formulas):
618        #
619        poroc = 0.4
620        #
621        # Calculate input common to both bulk and
622        # shear formulas
623        #
624        poratio = poro/poroc
625        #
626        # Calculate input to bulk formula
627        #
628        ba = bulkhm
629        bb = bulks
630        bz = 4*shearhm/3
631        #
632        # Calculate input to shear formula
633        #
634        sa = shearhm
635        sb = shears
636        sz = (shearhm/6)*((9*bulkhm + 8*shearhm) /
637                          (bulkhm + 2*shearhm))
638        #
639        # Calculate moduli
640        #
641        bulkd = self._gendryrock(poratio, ba, bb, bz)
642        sheard = self._gendryrock(poratio, sa, sb, sz)
643        #
644        return bulkd, sheard
645    #
646    # ---------------------------------------------------
647    #
648
649    def _hertzmindlin(self, poverb, pfluid):
650        #
651        # Calculate bulk and shear Hertz-Mindlin moduli
652        # utilizing that they have the common
653        # mathematical form,
654        #
655        # f = <a>max*(peff/pref)^kappa.
656        #
657        #
658        # Input
659        #       poverb - overburden pressure
660        #       pfluid - fluid pressure
661        #
662        # Output
663        #       bulkhm  - Hertz-Mindlin bulk modulus
664        #                 (MPa)
665        #       shearhm - Hertz-Mindlin shear modulus
666        #                 (MPa)
667        #
668        # -----------------------------------------------
669        #
670        # From table 1 in Report 1 in Abul Fahimuddin's
671        # thesis (he used GPa for the moduli, I use MPa
672        # also for them):
673        #
674        bulkmax = 3310
675        shearmax = 2840
676        pref = 8.8
677        kappa = 0.233
678        #
679        # Calculate moduli
680        #
681        peff = poverb - pfluid
682        if peff < 0:
683            print("\nError in _hertzmindlin method")
684            print("Negative effective pressure (" + str(peff) +
685                  "). Setting effective pressure to 0.01")
686            peff = 0.01
687   #         sys.exit(1)
688        common = (peff/pref)**kappa
689        bulkhm = bulkmax*common
690        shearhm = shearmax*common
691        #
692        return bulkhm, shearhm
693    #
694    # ---------------------------------------------------
695    #
696
697    def _gendryrock(self, q, a, b, z):
698        #
699        # Calculate arbitrary dry rock moduli, which has
700        # the generic form
701        #
702        #     --                                -- ^(-1)
703        #     |     q                1 - q       |
704        # f = |------------   +  ----------------|   - z,
705        #     |   a + z              b + z       |
706        #     --.                               --
707        #
708        # both for bulk and shear moduli
709        #
710        #
711        # Input
712        #       q - see above formula
713        #       a - see above formula
714        #       b - see above formula
715        #       z - see above formula
716        #
717        # Output
718        #       f - Bulk or shear dry rock modulus value
719        #
720        # -----------------------------------------------
721        #
722        afrac = q/(a + z)
723        bfrac = (1 - q)/(b + z)
724        f = 1/(afrac + bfrac) - z
725        #
726        return f
727    # ===========================
728    # Dry rock properties end
729    # ===========================
730
731
732if __name__ == '__main__':
733    #
734    # Example input with two phases and three grid cells
735    #
736    porosity = [0.34999999, 0.34999999, 0.34999999]
737#    pressure = [ 29.29150963, 29.14003944, 28.88845444]
738    pressure = [29.3558, 29.2625, 29.3558]
739#    pressure = [ 25.0, 25.0, 25.0]
740    phases = ["Oil", "Wat"]
741#    saturations = [[0.72783828, 0.66568458, 0.58033288],
742#                   [0.27216172, 0.33431542, 0.41966712]]
743    saturations = [[0.6358, 0.5755, 0.6358],
744                   [0.3641, 0.4245, 0.3641]]
745#    saturations = [[0.4, 0.5, 0.6],
746#                   [0.6, 0.5, 0.4]]
747    petElProp = "bulk velocity"
748    input_dict = {}
749    input_dict['overburden'] = 'overb.npz'
750
751    print("\nInput:")
752    print("porosity, pressure:", porosity, pressure)
753    print("phases, saturations:", phases, saturations)
754    print("petElProp:", petElProp)
755    print("input_dict:", input_dict)
756
757    satrock = elasticproperties(input_dict)
758
759    print("overburden:", satrock.overburden)
760
761    satrock.calc_props(phases, saturations, pressure,
762                       porosity)
763
764    print("\nOutput from calc_props:")
765    print("Density:", satrock.getDens())
766    print("Bulk modulus:", satrock.getBulkMod())
767    print("Shear modulus:", satrock.getShearMod())
768    print("Bulk velocity:", satrock.getBulkVel())
769    print("Shear velocity:", satrock.getShearVel())
770    print("Bulk impedance:", satrock.getBulkImp())
771    print("Shear impedance:", satrock.getShearImp())
772
773    satrock.getMatchProp(petElProp)
774
775    print("\nOutput from getMatchProp:")
776    print("Model output selected for data match:",
777          satrock.match_prop)