Reference Profiles for Compressible Mantle Convection

pnguyen · November 29, 2019, 11:18pm

Hello!

I am using ASPECT for a project in a geodynamics class I am taking this semester. I am interested in reproducing the reference profiles generated by BurnMan from an example cookbook. In particular, it is section 5.3.6 in the manual titled “2D compressible convection with a reference profile and material properties from Burn-Man”. The section mentions a Python script used to created the profiles called “simple_adiabat.py”. I have searched through the BurnMan repository and have been unable to find this file. Does anyone know where I might be able to find it?

Thank you,

Pham

tjhei · November 30, 2019, 5:13pm

Hi Pham,

it looks like the instructions are out of date. There is an ipynb at https://github.com/geodynamics/burnman/blob/master/contrib/ipython/Create_1D_ASPECT_profile.ipynb (part of BurnMan) that might work.
Would you please try to see if you can use it to reproduce the output and report back? If yes, I will update the ASPECT manual.

bobmyhill · December 1, 2019, 9:36pm

Dear Pham,

Sorry for my delayed reply; I was flying yesterday.

In a similar vein to Timo’s response, there is also “example_geodynamic_adiabat.py” in the examples directory. Just in case you want raw code and not an ipython notebook.

If you want any more help, let us know.

Best wishes,

Bob

pnguyen · December 3, 2019, 5:27pm

Hi all,

I checked out the notebook to see what I could find out. Using default values the profiles are too shallow and are only calculated to a depth of ~700 km. I tried increasing the maximum pressure but I get the following error:

ValueError: The set_state pressure (25064003600.0000) is outside the bounds of this rock (0.0001-25.0000 GPa)

I think this error originates from the default_file which is set to ‘…/…/burnman/data/input_perplex/in23_1.tab’ which looks like a file tabulating rock densities, expansion coefficients, seismic velocities etc. I am not sure how this file was generated. Would I have to generate a new file like this to make profiles for greater depths?

Thanks,

Pham

bobmyhill · December 3, 2019, 6:35pm

Dear Pham,

Yes, this example uses output from PerpleX, which deals efficiently with the tricky problem of thermodynamic equilibration (http://www.perplex.ethz.ch/). In order to keep the size of the burnman package to a minimum, we only include a small PerpleX-generated table.

The programs provided by PerpleX are fairly straightforward to use, but there is a bit of a learning curve. A first step will be to download PerpleX and maybe try using build and vertex (example here: http://www.perplex.ethz.ch/perplex_66_seismic_velocity.html). If you want to use the Stixrude dataset, then I have some build files where you can just change the composition and P-T bounds, which will save you some time.

The steps to creating your own tables are as follows:

Make a PerpleX input file (with PerpleX-build, or by modifying an existing file)
Run PerpleX-vertex (this calculates equilibrium phase assemblages as a function of P and T).
Make the burnman-readable tables. burnman/misc/create_burnman_readable_perplex_table.py makes this step nearly-automatic.

Let me know if you want any help.

Best wishes,

Bob

pnguyen · December 3, 2019, 7:43pm

Hi Dr. Myhill,

Thank you for the detailed explanation! Since I am still new to using BurnMan, I think for now I would be interested in the build files you have for the Stixrude data set. I will look into PerpleX in more detail as well.

Pham

bobmyhill · December 3, 2019, 8:09pm

Hi Pham,

Ok, great. You’ll still need to run PerpleX to be able to use the build files, so let me know when you have things up and running and we can go from there.

Bob

pnguyen · December 3, 2019, 9:28pm

Hi Dr. Myhill,

Ok, I think I have PerpleX up and running. What’s the next step?

Pham

bobmyhill · December 3, 2019, 10:01pm

Ok, now you need to create a directory containing all the files you need to run vertex.
You will need a copy of the correct endmember file, solution model file and preferably also an option file.

If you want to use Lars Stixrude’s dataset, stx11ver.dat is the endmember file, stx11_solution_model.dat is the solution model file, and copy the perplex_option.dat file too (you may wish to change some of the options at some point, and it’s always good to store a new copy of this file in each directory).

Then, you will either need to create a new input file using PerpleX-build, or edit an existing file. I’ve attached a file for a pyrolite to this message (if you use this, you will need to rename pyrolite.txt to pyrolite.dat). I strongly recommend that you take the time to understand the contents of this file, so that you know which bits to edit. This file will be used to calculate properties between 1-0.3e6 bar (0.0001-30 GPa) and 1400-2200 K.

After doing this, run PerpleX-vertex to compute all of the properties you need. This will take a few minutes.

Finally, run create_burnman_readable_perplex_table.py in the same directory. This will create an output file that burnman can read and process.

And that should be everything! PerpleX is pretty awesome software.

Let me know if you have any problems. Sometimes PerpleX-werami spits out NaNs in the file created by create_burnman_readable_perplex_table.py. You may need to change these to some sensible number (I usually copy-paste the values from the preceding or following rows; the NaNs are usually not in important parts of P-T space, so it doesn’t matter too much which numbers you use).

pyrolite.txt (2.3 KB)

pnguyen · December 5, 2019, 8:48pm

I have successfully ran PerpleX-vertex and created output files using pyrolite as an input. I am having some issues running the create_burnman_readable_perplex_table.py script. I get the following error:

Traceback (most recent call last):
File “create_burnman_readable_perplex_table.py”, line 38, in
create_perplex_table(args.werami_path, args.project[0], args.outfile[0], args.n_pressures[0], args.n_temperatures[0], args.pressure_range, args.temperature_range)
File “/Users/Pham/anaconda3/lib/python3.7/site-packages/burnman/perplex.py”, line 67, in create_perplex_table
stdout = p.communicate(input=stdin)[0]
File “/Users/Pham/anaconda3/lib/python3.7/subprocess.py”, line 939, in communicate
stdout, stderr = self._communicate(input, endtime, timeout)
File “/Users/Pham/anaconda3/lib/python3.7/subprocess.py”, line 1666, in _communicate
input_view = memoryview(self._input)
TypeError: memoryview: a bytes-like object is required, not ‘str’

My terminal input looks like this:

python create_burnman_readable_perplex_table.py --werami_path /Users/path_to_werami --project pyrolite --outfile results --n_pressures 100 --n_temperatures 100

I have tried changing the strings to a byte-like object but this just raises another error. Could this be an issue with the Python version I am using? My terminal defaults to Python version 3.7.3.

bobmyhill · December 5, 2019, 10:50pm

Oh darn. Yes, the code is written for python 2.

If you replace that line with stdout = p.communicate(input=stdin.encode(‘utf-8’))[0].decode(“utf-8”) then the code should work.

It’s probably a good idea to quickly check the columns in the output file for the right properties in the correct order. This is what the python code assumes:
2 - Density (kg/m3)
4 - Expansivity (1/K, for volume)
5 - Compressibility (1/bar, for volume)
10 - Adiabatic bulk modulus (bar)
11 - Adiabatic shear modulus (bar)
12 - Sound velocity (km/s)
13 - P-wave velocity (Vp, km/s)
14 - S-wave velocity (Vs, km/s)
17 - Entropy (J/K/kg)
18 - Enthalpy (J/kg)
19 - Heat Capacity (J/K/kg)
22 - Molar Volume (J/bar)

Best wishes,

Bob

pnguyen · December 9, 2019, 9:53pm

Thanks for the fix! Unfortunately, another issue came up after I modified the create_burnman_readable_perplex_tab.epy file:

Select the grid resolution (to use an arbitrary grid set sample_on_grid to F):
1 - 40 x 40 nodes
2 - 79 x 79 nodes
3 - 157 x 157 nodes
4 - 313 x 313 nodes [default]

Program received signal SIGILL: Illegal instruction.

Backtrace for this error:
#0 0x1040f5c6c
#1 0x1040f5393
#2 0x7fff6e71cb1c
#3 0x1041cd672

Traceback (most recent call last):
File “create_burnman_readable_perplex_table.py”, line 38, in
create_perplex_table(args.werami_path, args.project[0], args.outfile[0], args.n_pressures[0], args.n_temperatures[0], args.pressure_range, args.temperature_range)
File “/Users/Pham/anaconda3/lib/python3.7/site-packages/burnman/perplex.py”, line 70, in create_perplex_table
out = [s for s in stdout.split(’\n’) if “Output has been written to the” in s][0].split()[-1]
IndexError: list index out of range

Any idea what might be causing this? Thank you!

bobmyhill · December 10, 2019, 6:59am

Hi Pham,

PerpleX 6.8.8 has changed its input structure slightly, and seems to be in a bit of flux at the moment. I’m not having any problem with today’s version (collected from here: https://petrol.natur.cuni.cz/~ondro/perplex). Nevertheless, maybe this is the cause of the problem.

Three suggestions if you are still having problems:

Try changing sample_on_grid to F in perplex_option.dat (as suggested in the output you posted)
if that doesn’t work, try running werami manually, and go through the instructions to see where the input changes. The instructions given to werami by create…py are as follows:

mode 2

a set of properties (given as numbers, and then “n” to individual phase properties)

“y” to changing the default P-T range

Range values for P

Range values for T

Number of nodes

0 to Exit

Try today’s build of PerpleX

I think you’re almost there. Good luck!

Bob

pnguyen · December 18, 2019, 5:53pm

Hello Dr. Myhill,

I changed the sample_on_grid to F in the perplex_option.dat and it worked! I was able to create reference profiles using the Jupyter Notebook as well. Going back to the original question of this thread, can the pyrolite file be extended into the pressure and temperature range to recreate the reference profile found in the ASPECT manual?

Pham

bobmyhill · December 18, 2019, 9:05pm

Dear Pham,

That’s excellent news!

Where are you currently stuck? It sounds like you have everything working. These are the steps (I’m guessing you’ve done most of these already):

Modify the pyrolite PerpleX-build file so that PerpleX-vertex computes over your chosen P-T range.
Run vertex
Create the pyrolite tab file using create_burnman_table_…py
Modify the 1D adiabat ipython or python file to read in this new tab file rather than in23_1.tab.
Change maxP and anything else you want in the 1D adiabat ipython or python file.
Everything should now work.

This should allow you to create consistent profiles for any pressure range, adiabat and bulk composition. You should be able to calculate at any resolution you like.

Best wishes,

Bob

pnguyen · December 19, 2019, 6:37pm

Hello Dr. Myhill,

Everything is in working order now! Here is a sample reference profile from BurnMan after modifying the pyrolite file for a larger P-T range that is more in line with the profile seen in the reference manual (although not quite as deep).

bobmyhill · December 19, 2019, 9:10pm

That’s brilliant, thanks for letting me know!

Best wishes for your continued research,

Bob

vivianotieno · May 11, 2023, 12:48pm

Dear Bob

I have been trying to use the newer version of burnman to create a burnman readable file, and I keep getting the below error: I am using the newer version of perple_x of v.7.0.9. I have tried the various versions of burnman, and non seem to work. Thank you.

Keeping P-T range the same as the original project range.

Perple_X release 7.0.9, March 27, 2023.

Enter the project name (the name assigned in BUILD) [default = my_project]:

Reading Perple_X options from: stx21_perplex_option.dat
Writing Perple_X option summary to: not requested

Perple_X computational option settings for WERAMI:

Keyword:               Value:     Permitted values [default]:

Input/Output options:

aqueous_output          T         [F] T
aqeuous_species         20        [20] 0-150
aq_solvent_composition  y         [y] m: y => mol fraction, m => molality
aq_solute_composition   m         y [m]: y => mol fraction, m => molality
spreadsheet             T         [F] T
logarithmic_p           F         [F] T
logarithmic_X           F         [F] T
bad_number             NaN        [NaN]
composition_constant    F         [F] T
composition_phase       mol       [mol] wt
composition_system      wt        [wt] mol
proportions             vol       [vol] wt mol
absolute                F         [F] T
cumulative              F         [F] T
fancy_cumulative_modes  F         [F] T
interpolation           on        [on] off 
melt_is_fluid           T         [T] F
solution_names          mod       [model] abbreviation full
structural_formulae     T         [T] F
output_species          T         [T] F
output_species_props    F         [F] T
seismic_output          som       [some] none all
pause_on_error          T         [T] F
poisson_test            F         [F] T
interim_results         aut       [auto] off manual
sample_on_grid          T         [T] F

Information file output options:

option_list_files       F         [F] T; echo computational options

Thermodynamic options:

approx_alpha            T         [T] F
Anderson-Gruneisen      F         [F] T
hybrid_EoS_H2O             4      [4] 0-2, 4-7
hybrid_EoS_CO2             4      [4] 0-4, 7
hybrid_EoS_CH4             0      [0] 0-1, 7
fd_expansion_factor     2.0       [2] >0
finite_difference_p     0.1E+05   [1d4] >0; fraction = 0.1E-02   [1d-2]

Seismic wavespeed computational options:

seismic_data_file       T         [F] T
bounds                  VRH       [VRH] HS
vrh/hs_weighting        0.5       [0.5] 0->1
explicit_bulk_modulus   T         [T] F
poisson_ratio           on        [on] all off; Poisson ratio = 0.35
seismic_output          som       [some] none all
poisson_test            F         [F] T
Tisza_test              F         [F] T
fluid_shear_modulus     T         [T] F
phi_d                   0.36      [0.36] 0->1

To change these options see: www.perplex.ethz.ch/perplex_options.html

warning ver099 The P derivative of mu for qtz is < 0, this is improbable.

Select operational mode:
1 - properties at specified conditions
2 - properties on a 2d grid
3 - properties along a 1d path
4 - as in 3, but input from file
0 - EXIT

Select a property:
1 - Specific Enthalpy (J/m3)
2 - Density (kg/m3)
3 - Specific heat capacity (J/K/m3)
4 - Expansivity (1/K, for volume)
5 - Compressibility (1/bar, for volume)
6 - Composition (Mol, Mass, or Wt%) of the system
7 - Mode (Vol, Mol, or Wt proportion) of a phase
8 - Composition (Mol, Mass, or Wt%) of a solution phase
9 - Grueneisen thermal ratio
10 - Adiabatic bulk modulus (bar)
11 - Adiabatic shear modulus (bar)
12 - Sound velocity (km/s)
13 - P-wave velocity (Vp, km/s)
14 - S-wave velocity (Vs, km/s)
15 - Vp/Vs
16 - Specific entropy (J/K/m3)
17 - Entropy (J/K/kg)
18 - Enthalpy (J/kg)
19 - Heat Capacity (J/K/kg)
20 - Specific mass of a phase (kg/m3-system)
21 - Poisson ratio
22 - Molar Volume (J/bar)
23 - Dependent potentials (J/mol, bar, K)
24 - Assemblage Index
25 - Modes of all phases
26 - Sound velocity T derivative (km/s/K)
27 - P-wave velocity T derivative (km/s/K)
28 - S-wave velocity T derivative (km/s/K)
29 - Adiabatic bulk modulus T derivative (bar/K)
30 - Shear modulus T derivative (bar/K)
31 - Sound velocity P derivative (km/s/bar)
32 - P-wave velocity P derivative (km/s/bar)
33 - S-wave velocity P derivative (km/s/bar)
34 - Adiabatic bulk modulus P derivative (unitless)
35 - Shear modulus P derivative (unitless)
36 - All phase &/or system properties
37 - Absolute amount (Vol, Mol, or Wt) of a phase
38 - Multiple property output
39 - Heat capacity ratio (Cp/Cv)
40 - Lagged or back-calculated aqueous solute chemistry

Calculate individual phase properties (y/n)?