{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Julia-Perple_X interface demo\n", "\n", "This notebook runs the StatGeochem.jl package, which implements an interface for interacting with Perple_X from the Julia programming language, including from Jupyter notebooks such as this.\n", "\n", "StatGeochem.jl also includes some of the codes and utilities used in Keller & Schoene 2012, Keller et al. 2015 and Keller & Schoene 2018.\n", "\n", "\"Launch \n", "

If running this notebook as an online Binder notebook and the webpage times out, click the badge at left to relaunch (refreshing will not work). Note that any changes will be lost!

\n", "\n", "Hint: `shift`-`enter` to run a single cell, or from the `Cell` menu select `Run All` to run the whole file. Any code from this notebook can be copied and pasted into the Julia REPL or a `.jl` script.\n", "***\n", "### Load required Julia packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Recompiling stale cache file /Users/cbkeller/.julia/compiled/v1.1/StatGeochem/Ht7Cf.ji for StatGeochem [df4de05a-b714-11e8-3c2a-c30fb13e804c]\n", "└ @ Base loading.jl:1184\n" ] } ], "source": [ "## --- Load (and install if neccesary) the StatGeochem package which has the resampling functions we'll want\n", "\n", "try\n", " using StatGeochem\n", "catch\n", " using Pkg\n", " Pkg.add(PackageSpec(url=\"https://github.com/brenhinkeller/StatGeochem.jl\"))\n", " using StatGeochem\n", "end\n", "\n", "using Statistics, StatsBase, DelimitedFiles\n", "using Plots; gr();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Try to download and install Perple_X" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "## --- Configure\n", "\n", "# Absolute paths to perplex resources\n", "perplexdir = joinpath(resourcepath,\"perplex-stable\")\n", "scratchdir = \"./scratch/\" # Location of directory to store output files\n", "\n", "# Attempt to install perplex, if not already extant\n", "if !isfile(joinpath(perplexdir,\"vertex\"))\n", " # Make sure resourcepath exists\n", " run(`mkdir -p $resourcepath`)\n", "\n", " # Try to compile PerpleX from source; if that fails, try to download linux binaries\n", " try\n", " # Check if there is a fortran compiler\n", " run(`gfortran -v`)\n", "\n", " # Download Perplex v6.8.7 -- known to work with interface used here\n", " file = download(\"https://storage.googleapis.com/statgeochem/perplex-6.8.7-source.zip\", joinpath(resourcepath,\"perplex-stable.zip\"))\n", "\n", " # # For a more updated perplex version, you might also try\n", " # file = download(\"https://petrol.natur.cuni.cz/~ondro/perplex-sources-stable.zip\", joinpath(resourcepath,\"perplex-stable.zip\"))\n", "\n", " run(`unzip -u $file -d $resourcepath`) # Extract\n", " system(\"cd $perplexdir; make\") # Compile\n", " catch\n", " @warn \"Failed to compile from source, trying precompiled linux binaries instead\"\n", " run(`mkdir -p $perplexdir`)\n", " file = download(\"https://petrol.natur.cuni.cz/~ondro/Perple_X_6.8.7_Linux_64_gfortran.tar.gz\",\"perplex-6.8.7-linux.tar.gz\")\n", " run(`tar -xzf $file -C $perplexdir`)\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Configure Perple_X options" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "9-element Array{Float64,1}:\n", " 50.0956\n", " 15.3224\n", " 8.5103\n", " 9.252 \n", " 9.6912\n", " 2.5472\n", " 0.8588\n", " 2.0 \n", " 0.6 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## --- # # # # # # # # # # # # # Initial composition # # # # # # # # # # # # # #\n", "\n", "## McDonough Pyrolite\n", "#elements = [ \"SIO2\", \"TIO2\", \"AL2O3\", \"FEO\", \"MNO\", \"MGO\", \"CAO\", \"NA2O\", \"K2O\", \"H2O\", \"CO2\",]\n", "#composition = [45.1242, 0.2005, 4.4623, 8.0723, 0.1354, 37.9043, 3.5598, 0.3610, 0.0291, 0.1511, 0.0440,]\n", "\n", "## Kelemen (2014) primitive continental basalt. H2O and CO2 are guesses\n", "#elements = [ \"SIO2\", \"TIO2\", \"AL2O3\", \"FEO\", \"MNO\", \"MGO\", \"CAO\", \"NA2O\", \"K2O\", \"H2O\", \"CO2\",]\n", "#composition = [50.0956, 0.9564, 15.3224, 8.5103, 0.1659, 9.2520, 9.6912, 2.5472, 0.8588, 2.0000, 0.6000,]\n", "\n", "# Kelemen (2014) primitive continental basalt excluding Mn and Ti since most melt models can\"t handle them..\n", "elements = [ \"SIO2\", \"AL2O3\", \"FEO\", \"MGO\", \"CAO\", \"NA2O\", \"K2O\", \"H2O\", \"CO2\",]\n", "composition = [50.0956, 15.3224, 8.5103, 9.2520, 9.6912, 2.5472, 0.8588, 2.0000, 0.6000,]\n", "\n", "## Average Archean basalt (EarthChem data)\n", "#elements = [ \"SIO2\", \"TIO2\", \"AL2O3\", \"FEO\", \"MNO\", \"MGO\", \"CAO\", \"NA2O\", \"K2O\", \"H2O\", \"CO2\",]\n", "#composition = [49.2054, 0.8401, 12.0551, 11.4018, 0.2198, 12.3997, 9.3113, 1.6549, 0.4630, 1.8935, 0.5555,]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"\"" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## --- # # # # # # # # # # # Some solution model options # # # # # # # # # # # #\n", "\n", "# Emphasis on phases from Green (2016) -- developed for metabasites, includes what is probably the best (and most expensive) amphibole model. Use with hp11ver.dat\n", "G_solution_phases = \"Augite(G)\\nOpx(JH)\\ncAmph(G)\\noAmph(DP)\\nO(JH)\\nSp(JH)\\nGrt(JH)\\nfeldspar_B\\nMica(W)\\nBio(TCC)\\nChl(W)\\nCtd(W)\\nCrd(W)\\nSa(WP)\\nSt(W)\\nIlm(WPH)\\nAtg(PN)\\nT\\nB\\nF\\nDo(HP)\\nScap\\nChum\\nNeph(FB)\\n\"\n", "G_excludes =\"ged\\nfanth\\ngl\\nilm\\nilm_nol\\n\"\n", "\n", "# Emphasis on phases from White (2014) -- developed for metapelites. Use with hp11ver.dat (though can apparenty run with hp02ver.dat without crashing)\n", "W_solution_phases = \"Omph(HP)\\nOpx(W)\\ncAmph(DP)\\noAmph(DP)\\nO(JH)\\nSp(JH)\\nGt(W)\\nfeldspar_B\\nMica(W)\\nBi(W)\\nChl(W)\\nCtd(W)\\nCrd(W)\\nSa(WP)\\nSt(W) \\nIlm(WPH)\\nAtg(PN)\\nT\\nB\\nF\\nDo(HP)\\nScap\\nChum\\nPu(M)\\n\"\n", "W_excludes = \"andr\\nts\\nparg\\ngl\\nged\\nfanth\\n\"\n", "\n", "# Emphasis on phases from Jennings and Holland (2015) -- developed for mantle melting. Use with hp11ver.dat\n", "JH_solution_phases = \"Cpx(JH)\\nOpx(JH)\\ncAmph(DP)\\noAmph(DP)\\nO(JH)\\nSp(JH)\\nGrt(JH)\\nfeldspar_B\\nMica(W)\\nBio(TCC)\\nChl(W)\\nCtd(W)\\nCrd(W)\\nSa(WP)\\nSt(W)\\nIlm(WPH)\\nAtg(PN)\\nT\\nB\\nF\\nDo(HP)\\nScap\\nChum\\nNeph(FB)\\n\"\n", "JH_excludes = \"ts\\nparg\\ngl\\nged\\nfanth\\n\"\n", "\n", "# Emphasis on phases from Holland and Powell -- all phases can be used with hp02ver.dat.\n", "HP_solution_phases = \"Omph(HP)\\nOpx(HP)\\nGlTrTsPg\\nAnth\\nO(HP)\\nSp(HP)\\nGt(HP)\\nfeldspar_B\\nMica(CF)\\nBio(TCC)\\nChl(HP)\\nCtd(HP)\\nSapp(HP)\\nSt(HP)\\nIlHm(A)\\nDo(HP)\\nT\\nB\\nF\\n\"\n", "HP_excludes = \"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run Perple_X" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 54.792925 seconds (464.08 k allocations: 23.311 MiB, 0.02% gc time)\n", "\n", "----------------------------------------\n", "\n", "Stable phases at:\n", " T(K) = 1723.15 \n", " P(bar) = 10000.0 \n", "\n", "Phase Compositions (molar proportions):\n", " wt % vol % mol % mol H2O CO2 FEO MGO CAO NA2O K2O SIO2 AL2O3\n", " melt(G) 99.35 98.28 97.33 0.590 0.18379 0.00000 0.20071 0.38896 0.29283 0.06964 0.01545 1.41272 0.25463\n", " F 0.65 1.72 2.67 0.162E-01 0.15752 0.84248 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n", "\n", "Phase speciation (molar proportions):\n", "\n", " melt(G) foTL: 0.11411, faTL: 0.05889, abL: 0.16344, silL: 0.02543, kspL: 0.03626, woL: 0.17010, q8L: 0.04254\n", " h2oL: 0.21569, oanL: 0.17354\n", " F H2O: 0.15752, CO2: 0.84248\n", "\n", "Molar Properties and Density:\n", " N(g) G(J) S(J/K) V(J/bar) Cp(J/K) Alpha(1/K) Beta(1/bar) Cp/Cv Density(kg/m3)\n", " melt(G) 166.44 -2819391 506.13 6.3017 256.47 0.87739E-04 0.40948E-05 1.0865 2641.2 \n", " F 39.92 -647615 214.51 4.0329 60.884 0.26245E-03 0.43779E-04 1.2189 989.75 \n", " System 98.88 -1674426 302.18 3.7844 152.35 0.90752E-04 0.47792E-05 1.0796 2612.7 \n", " System - fluid 98.23 -1663946 298.71 3.7192 151.37 0.87739E-04 0.40948E-05 1.0865 2641.2 \n", "\n", "Seismic Properties:\n", " Gruneisen_T Ks(bar) Mu(bar) V0(km/s) Vp(km/s) Vs(km/s) Poisson ratio\n", " melt(G) 0.57200 0.26533E+06 88443. 3.1695 3.8093 1.8299 0.35000 \n", " F 0.48403 27842. 0.0000 1.6772 1.6772 0.0000 0.50000 \n", " System 0.57049 0.24627E+06 43459. 3.0701 3.4123 1.2897 0.41667 \n", " System - fluid 0.57200 0.26533E+06 88443. 3.1695 3.8093 1.8299 0.35000 \n", "\n", "Bulk Composition:\n", "\n", " Complete Assemblage Solid+Melt Only\n", " mol g wt % mol/kg mol g wt % mol/kg\n", " H2O 0.111 2.000 2.023 1.123 0.108 1.954 1.989 1.104\n", " CO2 0.014 0.600 0.607 0.138 0.000 0.000 0.000 0.000\n", " FEO 0.118 8.510 8.607 1.198 0.118 8.510 8.664 1.206\n", " MGO 0.230 9.252 9.357 2.322 0.230 9.252 9.419 2.337\n", " CAO 0.173 9.691 9.801 1.748 0.173 9.691 9.866 1.759\n", " NA2O 0.041 2.547 2.576 0.416 0.041 2.547 2.593 0.418\n", " K2O 0.009 0.859 0.869 0.092 0.009 0.859 0.874 0.093\n", " SIO2 0.834 50.096 50.664 8.432 0.834 50.096 50.997 8.488\n", " AL2O3 0.150 15.322 15.496 1.520 0.150 15.322 15.598 1.530\n", "\n", "Other Bulk Properties:\n", "\n", " Enthalpy (J/kg) = -.116682E+08\n", " Specific Enthalpy (J/m3) = -.304862E+11\n", " Entropy (J/K/kg) = 3056.09 \n", " Specific Entropy (J/K/m3) = 0.798480E+07\n", " Heat Capacity (J/K/kg) = 1540.80 \n", " Specific Heat Capacity (J/K/m3) = 0.402573E+07\n", "\n", "\n", " Solid Enthalpy (J/kg) = -.116992E+08\n", " Solid Secific Enthalpy (J/m3) (2) = -.309002E+11\n", " Solid Entropy (J/K/kg) = 3040.85 \n", " Solid Specific Entropy (J/K/m3) = 0.803158E+07\n", " Solid Heat Capacity (J/K/kg) (1) = 1557.70 \n", " Solid Specific Heat Capacity (J/K/m3) (1) = 0.406988E+07\n", "\n", "\n", "N.B.: Aggregate properties represent the entire stable assemblage.\n", "Solid aggregate properties represent solid and melt properties,\n", "but do not include molecular fluid properties.\n", "\n", "\n", "Chemical Potentials (J/mol):\n", "\n", " H2O CO2 FEO MGO CAO NA2O K2O SIO2 AL2O3\n", " -509890. -673366. -478356. -738298. -864198. -0.102378E+07 -0.111639E+07 -0.106581E+07 -0.194481E+07\n", "\n", "Variance (c-p+2) = 9\n", "\n", "\n", "----------------------------------------\n", "\n" ] } ], "source": [ "## --- # # # # # # # # # # # # # Isobaric example # # # # # # # # # # # # # # # #\n", "\n", " # Input parameters\n", " P = 10000 # Pressure, bar\n", " T_range = [500+273.15, 1500+273.15] # Temperature range, Kelvin\n", " melt_model = \"melt(G)\"\n", "\n", " # Configure (run build and vertex)\n", " @time perplex_configure_isobar(perplexdir, scratchdir, composition, elements,\n", " P, T_range, dataset=\"hp11ver.dat\", npoints=100, excludes=G_excludes,\n", " solution_phases=melt_model*\"\\n\"*G_solution_phases)\n", "\n", "## --- Query all properties at a single temperature -- results returned as text\n", "\n", " T = 1450+273.15\n", " data_isobaric = perplex_query_1d(perplexdir, scratchdir, T) |> print" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Warning: Perplex seems to be reporting mole fractions instead of weight percentages, attempting to correct\n", "└ @ StatGeochem /Users/cbkeller/Documents/julia/packages/StatGeochem.jl/src/utilities/Geochemistry.jl:649\n" ] } ], "source": [ "## --- Query the full isobar -- results returned as dict\n", "\n", "bulk = perplex_query_system(perplexdir, scratchdir) # Get system data for all temperatures. Set include_fluid = \"n\" to get solid+melt only\n", "modes = perplex_query_modes(perplexdir, scratchdir) # || phase modes\n", "melt = perplex_query_phase(perplexdir, scratchdir, melt_model) # || melt data\n", "\n", "# Melt wt.% seems to be slightly inaccurate; use values from modes instead\n", "melt[\"wt_pct\"] = modes[melt_model]\n", "\n", "# Create dictionary to hold solid composition and fill it using what we know from system and melt\n", "solid = Dict()\n", "solid[\"wt_pct\"] = 100 .- melt[\"wt_pct\"]\n", "for e in [\"SIO2\",\"AL2O3\",\"FEO\",\"MGO\",\"CAO\",\"NA2O\",\"K2O\"]\n", " solid[e] = (bulk[e] - (melt[e] .* melt[\"wt_pct\"]/100)) ./ (solid[\"wt_pct\"]/100)\n", "end\n", "renormalize!(solid,[\"SIO2\",\"AL2O3\",\"FEO\",\"MGO\",\"CAO\",\"NA2O\",\"K2O\"],total=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot results" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "25\n", "\n", "\n", "50\n", "\n", "\n", "75\n", "\n", "\n", "100\n", "\n", "\n", "0\n", "\n", "\n", "10\n", "\n", "\n", "20\n", "\n", "\n", "30\n", "\n", "\n", "40\n", "\n", "\n", "50\n", "\n", "\n", "60\n", "\n", "\n", "\n", "\n", "melt(G) + G_solution_phases, 10000 bar\n", "\n", "\n", "Percent melt\n", "\n", "\n", "Wt. % in melt\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "SIO2\n", "\n", "\n", "\n", "AL2O3\n", "\n", "\n", "\n", "FEO\n", "\n", "\n", "\n", "MGO\n", "\n", "\n", "\n", "CAO\n", "\n", "\n", "\n", "NA2O\n", "\n", "\n", "\n", "K2O\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Plot melt composition as a function of melt percent\n", "\n", "h = plot(xlabel=\"Percent melt\", ylabel=\"Wt. % in melt\", title=\"$melt_model + G_solution_phases, $P bar\")\n", "i = 0\n", "for e in [\"SIO2\",\"AL2O3\",\"FEO\",\"MGO\",\"CAO\",\"NA2O\",\"K2O\"]\n", " plot!(h, melt[\"wt_pct\"], melt[e], label=e, color=lines[global i += 1])\n", " plot!(h, melt[\"wt_pct\"], bulk[e], label=\"\", color=lines[i], linestyle=:dot)\n", "end\n", "plot!(h,fg_color_legend=:white, framestyle=:box)\n", "# savefig(h,\"MeltComposition.pdf\")\n", "display(h)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "25\n", "\n", "\n", "50\n", "\n", "\n", "75\n", "\n", "\n", "100\n", "\n", "\n", "0\n", "\n", "\n", "10\n", "\n", "\n", "20\n", "\n", "\n", "30\n", "\n", "\n", "40\n", "\n", "\n", "50\n", "\n", "\n", "\n", "\n", "melt(G) + G_solution_phases, 10000 bar\n", "\n", "\n", "Percent melt\n", "\n", "\n", "Wt. % in solid\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "SIO2\n", "\n", "\n", "\n", "AL2O3\n", "\n", "\n", "\n", "FEO\n", "\n", "\n", "\n", "MGO\n", "\n", "\n", "\n", "CAO\n", "\n", "\n", "\n", "NA2O\n", "\n", "\n", "\n", "K2O\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Plot solid composition as a function of melt percent\n", "\n", "h = plot(xlabel=\"Percent melt\", ylabel=\"Wt. % in solid\", title=\"$melt_model + G_solution_phases, $P bar\")\n", "i = 0\n", "for e in [\"SIO2\",\"AL2O3\",\"FEO\",\"MGO\",\"CAO\",\"NA2O\",\"K2O\"]\n", " plot!(h, melt[\"wt_pct\"], solid[e], label=e, color=lines[global i +=1])\n", "end\n", "plot!(h,fg_color_legend=:white, framestyle=:box, legend=:topleft)\n", "# savefig(h,\"SolidComposition.pdf\")\n", "display(h)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "600\n", "\n", "\n", "800\n", "\n", "\n", "1000\n", "\n", "\n", "1200\n", "\n", "\n", "1400\n", "\n", "\n", "0\n", "\n", "\n", "25\n", "\n", "\n", "50\n", "\n", "\n", "75\n", "\n", "\n", "100\n", "\n", "\n", "\n", "\n", "melt(G) + G_solution_phases, 10000 bar\n", "\n", "\n", "T (C)\n", "\n", "\n", "Weight percent\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Augite(G)\n", "\n", "\n", "\n", "cAmph(G)\n", "\n", "\n", "\n", "Chl(W)\n", "\n", "\n", "\n", "Mica(W)\n", "\n", "\n", "\n", "Bio(TCC)\n", "\n", "\n", "\n", "Do(HP)\n", "\n", "\n", "\n", "zo\n", "\n", "\n", "\n", "ab\n", "\n", "\n", "\n", "q\n", "\n", "\n", "\n", "F\n", "\n", "\n", "\n", "feldspar_B\n", "\n", "\n", "\n", "Opx(JH)\n", "\n", "\n", "\n", "melt(G)\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Plot modes of all phases as a function of temperature\n", "\n", "h = plot(xlabel=\"T (C)\", ylabel=\"Weight percent\", title=\"$melt_model + G_solution_phases, $P bar\")\n", "for m in modes[\"elements\"][3:end]\n", " plot!(h, modes[\"T(K)\"] .- 273.15, modes[m], label=m)\n", "end\n", "plot!(h,fg_color_legend=:white, framestyle=:box)\n", "# savefig(h,\"PhaseModes.pdf\")\n", "display(h)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "ename": "MethodError", "evalue": "MethodError: no method matching perplex_configure_geotherm(::String, ::String, ::Array{Float64,1}, ::Array{String,1}, ::Array{Int64,1}, ::Float64, ::Float64; dataset=\"hp02ver.dat\", excludes=\"\", solution_phases=\"Omph(HP)\\nOpx(HP)\\nGlTrTsPg\\nAnth\\nO(HP)\\nSp(HP)\\nGt(HP)\\nfeldspar_B\\nMica(CF)\\nBio(TCC)\\nChl(HP)\\nCtd(HP)\\nSapp(HP)\\nSt(HP)\\nIlHm(A)\\nDo(HP)\\nT\\nB\\nF\\n\", npoints=200, index=2)\nClosest candidates are:\n perplex_configure_geotherm(::String, ::String, ::Array{#s47,N} where N where #s47<:Number; elements, P_range, T_surf, geotherm, dataset, solution_phases, excludes, index, npoints) at /Users/cbkeller/Documents/julia/packages/StatGeochem.jl/src/utilities/Geochemistry.jl:427", "output_type": "error", "traceback": [ "MethodError: no method matching perplex_configure_geotherm(::String, ::String, ::Array{Float64,1}, ::Array{String,1}, ::Array{Int64,1}, ::Float64, ::Float64; dataset=\"hp02ver.dat\", excludes=\"\", solution_phases=\"Omph(HP)\\nOpx(HP)\\nGlTrTsPg\\nAnth\\nO(HP)\\nSp(HP)\\nGt(HP)\\nfeldspar_B\\nMica(CF)\\nBio(TCC)\\nChl(HP)\\nCtd(HP)\\nSapp(HP)\\nSt(HP)\\nIlHm(A)\\nDo(HP)\\nT\\nB\\nF\\n\", npoints=200, index=2)\nClosest candidates are:\n perplex_configure_geotherm(::String, ::String, ::Array{#s47,N} where N where #s47<:Number; elements, P_range, T_surf, geotherm, dataset, solution_phases, excludes, index, npoints) at /Users/cbkeller/Documents/julia/packages/StatGeochem.jl/src/utilities/Geochemistry.jl:427", "", "Stacktrace:", " [1] top-level scope at util.jl:156", " [2] top-level scope at In[10]:8" ] } ], "source": [ "## --- # # # # # # # # # # # Geothermal gradient example # # # # # # # # # # # #\n", "\n", "# Input parameters\n", "P_range = [280, 28000] # Pressure range to explore, bar (roughly 1-100 km depth)\n", "T_surf = 273.15 # Temperature of surface (K)\n", "geotherm = 0.1 # Geothermal gradient of 0.1 K/bar == about 28.4 K/km\n", "melt_model = \"\"\n", "\n", "# Configure (run build and vertex)\n", "@time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,\n", " P_range, T_surf, geotherm, dataset=\"hp02ver.dat\", excludes=HP_excludes,\n", " solution_phases=HP_solution_phases, npoints=200, index=2)\n", "\n", "# # Alternative configuration, using hpha02ver.dat\n", "# @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,\n", "# P_range, T_surf, geotherm, dataset=\"hpha02ver.dat\", excludes=\"qGL\\n\"*HP_excludes,\n", "# solution_phases=HP_solution_phases, npoints=200, index=2)\n", "\n", "# # Alternative configuration, using hpha02ver.dat and new phases for metapelites\n", "# @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,\n", "# P_range, T_surf, geotherm, dataset=\"hpha02ver.dat\", excludes=\"qGL\\n\"*W_excludes,\n", "# solution_phases=W_solution_phases, npoints=200, index=2)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "500\n", "\n", "\n", "1000\n", "\n", "\n", "1500\n", "\n", "\n", "2000\n", "\n", "\n", "2500\n", "\n", "\n", "0\n", "\n", "\n", "10\n", "\n", "\n", "20\n", "\n", "\n", "30\n", "\n", "\n", "40\n", "\n", "\n", "50\n", "\n", "\n", "\n", "\n", "T (C)\n", "\n", "\n", "Weight percent\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Bio(TCC)\n", "\n", "\n", "\n", "Ctd(HP)\n", "\n", "\n", "\n", "Gt(HP)\n", "\n", "\n", "\n", "GlTrTsPg\n", "\n", "\n", "\n", "sid\n", "\n", "\n", "\n", "acti\n", "\n", "\n", "\n", "q8L\n", "\n", "\n", "\n", "Do(HP)\n", "\n", "\n", "\n", "cz\n", "\n", "\n", "\n", "mic\n", "\n", "\n", "\n", "ky\n", "\n", "\n", "\n", "ab\n", "\n", "\n", "\n", "Mica(CF)\n", "\n", "\n", "\n", "q\n", "\n", "\n", "\n", "Anth\n", "\n", "\n", "\n", "zo\n", "\n", "\n", "\n", "feldspar_B\n", "\n", "\n", "\n", "cc\n", "\n", "\n", "\n", "Chl(HP)\n", "\n", "\n", "\n", "F\n", "\n", "\n", "\n", "Opx(HP)\n", "\n", "\n", "\n", "Omph(HP)\n", "\n", "\n", "\n", "O(HP)\n", "\n", "\n", "\n", "abL\n", "\n", "\n", "\n", "faGL\n", "\n", "\n", "\n", "kspL\n", "\n", "\n", "\n", "anL\n", "\n", "\n", "\n", "diL\n", "\n", "\n", "\n", "enL\n", "\n", "\n", "\n", "anl\n", "\n", "\n", "\n", "qGL\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Plot modes of all phases as a function of temperature\n", "\n", "# Get phase modes\n", "modes = perplex_query_modes(perplexdir, scratchdir, index=2)\n", "\n", "h = plot(xlabel=\"T (C)\", ylabel=\"Weight percent\")\n", "for m in modes[\"elements\"][3:end]\n", " plot!(h, modes[\"T(K)\"] .- 273.15, modes[m], label=m)\n", "end\n", "plot!(h,fg_color_legend=:white, framestyle=:box)\n", "# savefig(h,\"GeothermPhaseModes.pdf\")\n", "display(h)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "5.0×10\n", "\n", "\n", "3\n", "\n", "\n", "1.0×10\n", "\n", "\n", "4\n", "\n", "\n", "1.5×10\n", "\n", "\n", "4\n", "\n", "\n", "2.0×10\n", "\n", "\n", "4\n", "\n", "\n", "2.5×10\n", "\n", "\n", "4\n", "\n", "\n", "1\n", "\n", "\n", "2\n", "\n", "\n", "3\n", "\n", "\n", "4\n", "\n", "\n", "5\n", "\n", "\n", "6\n", "\n", "\n", "7\n", "\n", "\n", "Pressure\n", "\n", "\n", "Property\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "vp,km/s\n", "\n", "\n", "\n", "vp/vs\n", "\n", "\n", "\n", "rho, g/cc\n", "\n", "\n", "\n", "T(K)/1000\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Plot seismic properties\n", "\n", "# Query seismic properties along the whole profile\n", "seismic = perplex_query_seismic(perplexdir, scratchdir, index=2)\n", "seismic[\"vp/vs\"][seismic[\"vp/vs\"] .> 100] .= NaN # Exclude cases where vs drops to zero\n", "\n", "h = plot(xlabel=\"Pressure\", ylabel=\"Property\")\n", "plot!(h,seismic[\"P(bar)\"],seismic[\"vp,km/s\"], label=\"vp,km/s\")\n", "plot!(h,seismic[\"P(bar)\"],seismic[\"vp/vs\"], label=\"vp/vs\")\n", "plot!(h,seismic[\"P(bar)\"],seismic[\"rho,kg/m3\"]/1000, label=\"rho, g/cc\")\n", "plot!(h,seismic[\"P(bar)\"],seismic[\"T(K)\"]/1000, label=\"T(K)/1000\")\n", "# savefig(h,\"GeothermSeismicProperties.pdf\")\n", "display(h)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.5.3", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.3" } }, "nbformat": 4, "nbformat_minor": 2 }